Méthode Des Éléments Finis M1 GC KADA Abdelhak [PDF]

Zitiervorschau

République Algérienne Démocratique et Populaire Ministère de l’Enseignement Supérieur

‫وزارة التعليــــم العــــالي و البحث العلمــــي‬

et de la Recherche Scientifique

Université Hassiba Benbouali de Chlef

FACULTE DE GENIE CIVIL ET D’ARCHITECTURE DEPARTEMENT DE GENIE CIVIL

POLYCOPIE Par : KADA Abdelhak

METHODES Octobre 2017 DES ELEMENTS FINIS

Master 1 LMD Génie Civil 𝐾 𝑈 = 𝐹

Par : KADA Abdelhak

"Les lois fondamentales de la physique, qui sont la base de l’analyse statique et dynamique, sont vieilles de plusieurs centaines d’années. De ce fait, toute personne qui croit avoir découvert un principe fondamental de calcul des structures nouveau est tout simplement victime de sa propre ignorance." [ Pr E. L. Wilson 1998 – Berkeley University]

A la mémoire de mon Père

AVANT- PROPOS Le contenu de ce polycopié s’adresse en premier lieu aux étudiants de première année de master de génie civil, leur permettant de s’imprégner des différents aspects théoriques et approches de la méthode des éléments finis (MEF). Il est question de présenter la MEF selon le programme de Master académique, unité d’enseignement UEM 1.2 avec le but de faire comprendre à l’étudiant le déroulement de la méthode pour pouvoir assimiler les procédures et les étapes souvent dissimulées dans les logiciels. Le polycopié est organisé en quatre chapitres, d’abord une introduction avec des énoncés de base, l’approximation par des éléments finis à une, deux et trois dimensions et enfin une partie réservée à un certain nombre d’applications de la méthode. Ces derniers sont dédiées à la

résolution de

problèmes de barres, de poutres et des problèmes d’élasticité plane dans le but de faciliter la compréhension des étapes de la MEF pour le calcul des structures. Ce polycopié représente un outil de travail pour les étudiants de master 1 et les chapitres correspondent à des cours dispensés par l’auteur. Je remercie, Pr Lamri Belkacem, Dr Benarous Abdellah (enseignants à l’Université Hassiba Benbouali de Chlef) et Pr Kerdal Djamel (enseignant à l’Université des Sciences et de la Technologie d’Oran Mohamed Boudiaf), pour la relecture du polycopié, ainsi que pour les nombreuses remarques constructives qui m’ont été prodiguées.

L’auteur : KADA Abdelhak Université Hassiba Benbouali de Chlef Octobre 2017 .

TABLE DES MATIERES Page

Chapitre1

Introduction

1.1. Introduction

1

1.2. Objectifs du cours de la MEF

1

1.3. Idée de base de la Méthode des Eléments Finis

1

1.4. Pourquoi la MEF ?

3

1.5. Domaines d’application de la MEF

3

1.6. Problèmes de conditions aux limites

4

A- Equations gouvernantes B- Autre problème par comparaison 1.7. Un bref historique de la MEF

6

1.8. Concepts de base et procédures de la MEF

7

A- Concepts B- Procédures C- Support informatique et logiciels incorporant la MEF D- Avantages de la MEF 1.9. Notion de modélisation et comparaison de la MEF à d’autres méthodes

13

A- Notion de modélisation - system continu et system discret B- La MEF et les autres méthodes d’analyse C- Illustration du principe de la modélisation par un exemple

Chapitre2

Éléments finis à une dimension

2.1. Introduction

18

2.2. Méthode Directe pour les structures à éléments discrets

18

A- Elément ressort linéaire B- Elément fini barre (Treillis plan de barres) C- .Formulations pour les éléments finis 2.3. Application à des barres prismatiques uniformes A. Hypothèses et position du problème

25

B. Formation du système d’équations d’équilibre C. Transformation d’une base locale vers une base globale D. Matrice de rigidité dans un espace bidimensionnel –2D E. Contraintes et forces internes élémentaires 2.4 Application aux ossatures planes de poutres

31

A. Position du problème B. Matrice de rigidité élémentaire C. Matrice de rigidité d’un élément poutre à orientation arbitraire 2.5 Généralisation

35

A. Elément de poutre bidimensionnelle- 2D B- Elément de poutre tridimensionnelle -3D 2.6 Influence de la numérotation des nœuds sur la largeur de la bande de K

37

A- Matrice de connectivité C- Semi – bande 2.7 Méthode Formelle -Fonction d’Interpolation & Fonction de Forme

42

A. Introduction B. Position du problème et définition C. Approximation nodale D. Choix de la fonction de déplacement E. Principe de l’approximation par éléments finis E- Types d’Eléments Finis à une dimension 2.8 Matrices Caractéristiques par minimisation de l’énergie potentielle

51

A. Energie potentielle totale minimum B- Formulation de l’énergie potentielle stationnaire C. Energie potentielle totale 2.9 Formulation de la matrice de rigidité et du vecteur charge équivalent

54

A. Déformation et contrainte moyenne dans un E.F unidimensionnel 1-D B. Matrice de rigidité élémentaire C. Vecteur de charges équivalentes aux nœuds 2.10 Construction des Fonctions de Forme A. Principe de la méthode B. La méthode des déterminants

56

Chapitre 3 Éléments finis à deux et trois dimensions 3.1 Idée de base

59

3.2 Types d’Eléments Finis 2D et 3D

59

A- Eléments bidimensionnels - 2D B- Triangle de Pascal 3.3 Construction des fonctions de forme pour L’EF ‘T3’

65

3.4 Transformation géométrique

66

A- Relation entre les variables généralisées et les variables nodales B- Règles et propriétés de la matrice de transformation Te 3.5 Système de coordonnées de référence

72

A- Système de coordonnées locales normalisé B- Système de coordonnées naturel 3.6 Cas de problèmes d’éléments finis bidimensionnels de classe C0

75

A- Rappel et formulation de la théorie de base B- Contraintes planes et déformations planes C- Relation contrainte- déformation D- Formulation de la matrice de rigidité E- Equations d’équilibre et conditions aux limites 3.7 Intégration numérique – Méthode de Gauss

84

A- Intégration dans un espace unidimensionnel 1-D B- Intégration dans un espace à deux dimensions (2-D) C- Intégration dans un espace à trois dimensions (3-D) D- Transformation du domaine d’intégration (base globalebase locale) 3.8 Intégration de la matrice de rigidité

92

Chapitre4- Applications des éléments finis de base, aux problèmes de classe C1 et C0 4.1 Applications au problème ‘1D’ de barre

94

4.2 Applications au problème de classe C1 élément fini type Hermitien

99

4.3 Applications au problème ‘2D’ de Classe C0 E.F. Triangle linéaire à 3 nœuds (T3)

103

4.4 Applications au problème ‘2D’ élément fini rectangle linéaire à 4 nœuds – Q4

118

Chapitre 1 Introduction

Stade d’Oran (en réalisation) : Toiture réalisée par des barres tubulaires en acier tubulaires

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

Chapitre1

Introduction

APERÇU Ce chapitre, présente l’historique et la nécessité d’étudier la méthode des éléments finis.Il énumère brièvement les différents types de problèmes pour lesquels elle est appliquée ainsi que ses concepts de base. Il s’agit en particulier de discerner la différence entre les autres méthodes et la MEF et d’examiner la notion de l’approximation par les éléments finis. 1.1. Introduction : Des problèmes, qui dans un passé récent ont été considérés comme insolvables par les méthodes analytiques classiques, sont maintenant aisément résolus par les méthodes numériques dont la plus utilisée est la Méthode des Eléments Finis ou ‘MEF’. De ce fait, la complexité des calculs n’est plus d’actualité scientifique, surtout par l’avènement de l’ordinateur qui a amené les sciences de l’ingénieur au summum jamais atteint auparavant. 1.2. Objectifs du cours de la MEF : 

Comprendre les idées fondamentales de la MEF.



Connaître le comportement et l’utilité et de chaque type d’élément.



Comprendre le comportement physique du problème.



Être capable de préparer un modèle EF convenable pour un problème donné.



Connaître les limites du modèle par la MEF (outil numérique !).



Avoir un bagage théorique, pour être à même de comprendre

comment

fonctionnent les logiciels (Boites Noires), et de pouvoir réagir efficacement vis à vis des messages d’erreurs sur ordinateur. 

Pouvoir consulter sans difficulté les livres et les recueils sur le ‘Net’, concernant la MEF.

_____________________________________________________________________ ! Un logiciel commercial dont lequel la MEF est complètement dissimulée, est un système CAE (Computer Aided Engineering) qui ne peut être utilisé comme un outil de dessin autoCAD, donc une connaissance en la matière est nécessair.! Extrait de propos de BATH –MIT University

!.... bien que les logiciels de MEF peuvent rendre un bon ingénieur excellent, ils peuvent cependant rendre un mauvais ingénieur dangereux ! Commentaire de Cook & Ass. dans son livre ″ Concepts and Applications of Finite Element Analysis"

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-1-

Chapitre1

Introduction

1.3. Idée de base de la Méthode des Eléments Finis: On apprend comme ingénieur débutant à calculer les surfaces et volumes de corps de forme quelconque en les décomposant en un ensemble de corps élémentaires de formes connues pour ensuite l’appliquer au calcul de moments d’inertie ou de centres de gravité. Ce mode de pensé à conduit à la Méthode des Eléments Finis (MEF), ou analyse par élément finis (AEF), qui est basée sur l’idée de construire un objet compliqué avec des blocks simples, c’est à dire diviser l’objet compliqué en un petit nombre de pièces facilement manipulables. On peut rencontrer l’application de cette idée simple aussi bien dans la vie de tous les jours qu’en technologie et pour tous les problèmes de l’ingénieur. Exemples : -

Lego (jeux d’enfants de construction)

-

Bâtiments

-

L’approximation de la circonférence et de la surface d’un cercle en est un simple exemple : S(sup) S

S(inf)

Elément Si

)

i R

(a)

(b)

Figure1.1 (a) Limite supérieure et inférieure de la circonférence d’un cercle (b) Approximation de l’aire d’un cercle par des triangles La figure1.1(a) est une illustration du concept même des éléments finis par les mathématiques anciennes, pour trouver la circonférence d’un cercle par approximation de polygones. En utilisant l’appellation moderne on peut qualifier chaque côté du polygone d’élément finis. On remarquera que lorsque le nombre de polygones augmente la valeur d’approximation converge vers la valeur exacte.

La figure1 (b) illustre aussi l’idée bien ancienne d’un calcul approché bien avant que le nom ‘’ élément finis’’ ne soit courant.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-2-

Chapitre1

Introduction

La surface d’un seul triangle est :

S  12 R sin 2

i

i

2 2  1 R N sin(2 ) R lorsque N  S N  S i 2 N i 1 N

La surface du cercle :

N - nombre total des triangles considérés dans l’approximation (éléments).

On remarque que des objets compliqués peuvent être représentés par des parties simples (éléments) 1.4. Pourquoi choisir la MEF ? : Les procédés de conception de calcul et d’analyse reposent sur le calcul à la main, l’expérience ou le calcul automatique et la simulation par ordinateur.  La MEF c’est la méthode de simulation par ordinateur la plus utilisée par les ingénieurs. C’est une technique essentiellement numérique à partir de laquelle les équations gouvernantes (systèmes d’équations différentielles), sont représentées sous une forme matricielle très adaptée à une solution automatique par ordinateur.  Elle est intégrée dans toutes les applications de logiciels commerciaux de calcul des structures à interface graphique (GUI). 1.5. Domaines d’application de la MEF : Les principaux domaines d’application de la MEF sont au nombre de trois :  Problèmes d’équilibre et statique : dans lequel le comportement du système ne varie pas avec le temps,  Problèmes de dynamique et de stabilité (valeurs propres): ce sont des extensions des problèmes d’équilibre pour lesquelles des valeurs spécifiques ou critiques de certains paramètres sont déterminés,  Problèmes de propagation : ils concernent les problèmes où les phénomènes dont le comportement est dépendant du facteur temps.

Le tableau 1 suivant résume ces types d’applications : Spécialité Génie Civil structure

Problèmes d’équilibre

Problèmes Problèmes de valeurs propres de propagation Analyse statique de Fréquences et modes Réponse des structures : treillis, portiques, propres et stabilité structures à des plaque, coques, voiles, ponts, des structures. charges accidentelles béton précontraint (séisme, incendie)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-3-

Chapitre1 Géotechnique

Hydraulique

Génie Mécanique

Biomédical

Introduction Analyse des excavations, stabilité des talus, murs de soutènement, interaction sol-structure , analyse des contraintes dans les sols

Fréquence et modes propres des ouvrages enterrés et semienterrés , et problèmes d’interaction solstructure Analyse d’écoulements Périodes et modes potentiels, d’écoulements à propres de bassins surface libre, écoulement superficiels, digues, visqueux. mouvements des Analyse de structures liquides dans des bacs hydrauliques et barrages, etc. (conteneurs) rigides ou flexibles Problèmes de concentration Fréquences propres de contraintes. de vibrations et Analyse de contrainte de stabilité des machines pistons, de matériaux composites, etc. Analyse de contraintes dans les os , les dents. Capacité portante des systèmes des implants et -----_________-----prothèses. Mécanique des valves du cœur artificiel

Problèmes de solstructure dépendant du temps. Propagation de contraintes dans les sols et les roches Analyse de problèmes d’écoulements turbulents et propagation d’ondes. Ecoulements hydrodynamiques Problèmes de fissures et de fractures sous charges dynamiques Analyse d’impact Sur le crane. Dynamiques de structures anatomiques

1.6. Problèmes de conditions aux limites : A- Equations gouvernantes Les différents problèmes d’ingénieurs présentent tous une écriture en équations différentielles, dont la solution dépend des conditions aux limites (problème à valeur initiale), qui pour un ingénieur de génie civil, sont spécifiques aux appuis d’une structure quelconque. Ces équations peuvent prendre plusieurs formes : *

l’exemple le plus simple est peut être celui de la barre tendue d’un système à treillis figure 1.2:

du N dx d  du   A.E.   0 dx  dx  u (a )  u (b)  0 AE

N (a) Figure 1.2 : Barre tendue

(b)

x

a xb

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(1.1)

-4-

Chapitre1

Introduction

Le plus souvent les équations différentielles de ces problèmes sont identifiées, comme étant de type équations de Poisson ou de type équations de Laplace à une ou deux dimensions. *

Un autre exemple étant celui de l’équation de mouvement harmonique simple : d 2u   2 . u 2 dx u (a)  

a xb

(1.2)

u b    /

L’équation (1.2) définie un problème à valeurs initiales, unidimensionnel linéaire du second ordre. Le terme unidimensionnel fait référence au domaine qui est défini par un segment a  x  b . Le terme "second ordre" tient de l’équation (1.2) dont laquelle l’ordre le plus élevé de la dérivée est 2 (second ordre). Le terme linéaire provient du faite que le degré de chaque terme faisant intervenir la variable u est égal à 1. Bien que la MEF possède des applications très larges, l’objet de notre cours traitera dans la plus part des cas des problèmes à valeurs initiales linéaires du 2eme et 4eme ordre, qui sont importants en génie civil, pour le calcul des structures. *

Un exemple d’un problème à valeurs initiales de 4 eme ordre est celui de la flexion d’une poutre, figure 1.3. Q L

Figure 1.3 : Poutre encastrée aux extrémités *

 d 2u   EI . 2   F x  dx   du u (0)  (0)  0 dx du u (l )  (l )  0 dx d2 dx 2

0 xl (1.3)

Un autre exemple de problème bidimensionnel linéaire est celui de torsion d’une poutre (type équation de Poisson)  2u  2u   2 x 2 y 2 u0

*

 x, y   D x, y  sur  )

(domaine)

(1.4)

(limite du domaine D)

Un dernier exemple est celui de l’analyse de la flexion de plaque (dalle) simplement appuyée sur ses côtés  sous sollicitation F(x,y) :

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-5-

Chapitre1

Introduction   2u  2u   2   2u  2u   2  2   2  2  2   F x, y   x, y   D y  y  x y   x u0 et ( 2 u x 2   2 u y 2 )  0 sur  2 x 2

(1.5)

B- Autre problème par comparaison Toutes les formules de calcul des structures possèdent leurs vis-à-vis dans d’autres domaines proches de celui du génie civil, tel que le transfert thermique (conduction) largement utilisé dans le bâtiment que ce soit pour le calcul incendie ou l’isolation thermique des parois. Le tableau 2, ci-dessous dresse une comparaison des équations de base dans les deux domaines : Tableau 2 : variables de calcul et équations de base pour une analyse par les éléments finis Paramètres Variables indépendantes

structure Coordonnées x, y, z

Conduction thermique Coordonnées x, y, z

Variable dépendante(s)

Déplacements u, v, w

Température T

Champ de gradient

Déformations  x ,  y ,  xy

T=[ T,x T,y T,z]

Matrice de compatibilité

Constantes élastiques D

Conductivité thermique k

Champ induit

Contraintes = D.

Flux de chaleur f=-k. T

Charge de surface

Forces reparties en surface t

Flux normal fn aux limites

Effort interne

Forces internes Fx, Fy, Fz

Chaleur interne diffusée Q

    F   0

Fx,x+fy,y+fz,z – Q = 0

Equation d’équilibre

T

On notera au passage, que la conduction de chaleur est un problème scalaire car le paramètre champ, qui est la température T, n’est associé à aucune direction. Par opposition au champ de déplacement en structure qui est un champ de vecteur ayant des composantes orientées selon un système de coordonnées. 1.7. Un bref historique de la MEF : Quoi que le label Méthode des Eléments Finis fut utilisé pour la première fois en 1960(4), lorsqu’il a été utilisé par Ray Clough dans une publication sur les problèmes d’élasticité plane, l’idée de l’analyse par éléments finis date des années d’avant. En effet, aux questions qui est à l’origine de le MEF ? et en quelle période ?, il existe trois réponses selon que l’on s’adresse à un physicien( 1), mathématicien(2),ou un ingénieur ( 3). Le Tableau 2 montre l’évolution chronologique de la méthode. (2)

Une approche similaire à la méthode des éléments finis, utilisant un développement de

fonctions continues (une approximation à un certain ordre) défini pour des régions A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-6-

Chapitre1

Introduction

triangulaires a été pour la première fois proposé par Courant en 1943 dans la littérature des mathématiques appliquées pour la solution au problème de torsion de St Venant. (3)

La méthode des éléments finis tel que reconnue aujourd’hui a été présentée par Turner,

Clough, Martin et Topp en 1956 dans ‘Journal of the aeronautical sciences), sous le titre "Stiffness and Deflection Analysis of Complex Structures". Tableau 3 : Publications et périodes ayant marquées l’évolution de la MEF 1941 (1)

Hrenikoff : Division d’un problème d’élasticité au domaine continue, en un certain nombres d’éléments.

1943 (2)

Courant : Méthodes variationelles

1956 (3)

Turner,Clough,Martin, Topp : Rigidité-Methode Directe Clough : Finite Element, problèmes plans (le terme Elément Fini utilisé pour la 1er fois)

1960 (4)

1970 70)

(années

1980 (années 80) 1990 (années 90)

Applications sur gros ordinateurs Applications sur micro-ordinateurs Possibilité d’analyse de gros systèmes de structures

La publication est une représentation systématique de la méthode des déplacements. C’est une contribution clé pour la solution aux problèmes de contraintes planes tel quelle se présente aujourd’hui en utilisant des éléments finis dit "triangulaires" dont les propriétés sont déterminé à partir des équations de la théorie de l’élasticité . L’ordinateur a offert un moyen très rapide pour tous les calculs multiples et divers (essentiellement numériques) qu’exige la MEF, ce qui a fait que la méthode est devenue très intéressante dans sa pratique. En même temps que le développement d’ordinateurs de plus en plus rapides, l’application de la MEF a aussi progressée à grands pas. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-7-

Chapitre1

Introduction

1.8. Concepts de base, procédures de la MEF: A- Concepts

C’est l’intuition physique qui pour la première fois, a mis en évidence les concepts de la MEF, pour les ingénieurs. Pour un ingénieur de structure, le problème d’un treillis par exemple figure 5 (a), sera un ensemble de barres dont il combine les caractéristiques individuelles selon les lois d’équilibre pour ensuite résoudre le système d’équations pour tout le système. Cette procédure marche bien pour un certain nombre de points de connections (nœuds), mais qu’en est il cependant d’une structure continue comme la plaque figure 5 (b) qui possède un nombre infini de points de connexion (nœuds) ? Le problème aurait été plus difficile sans l’intuition de Hrenikoff ((1) Tableau 2) qui proposa de diviser la structure continue en un certain nombre d’éléments de sections de structure (barres ou poutres) interconnectés à un nombre fini de points nodaux, figure 1.4 .

Charge

Charge

(a)

(b)

Figure 1.4 : Exemple (a) d’un treillis et (b) d’une plaque de forme similaire supportant la même charge  Concept 1 Les problèmes d’ingénieurs sont le plus souvent exprimés en terme :  

d’équations gouvernantes (différentielles essentiellement) et des conditions aux limites

Exemple : Barre en traction (figure 1.5) du A.E  N  Cst .a  x  b  (1.6) dx u (a)  u (b)  0 Plus généralement :

(a )

1

(b) N 2 x (u)

Figure 1.5 : Barre tendue

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-8-

Chapitre1

Introduction

L(u) + f = 0 eqt. diff.

Problèmes d’Elasticité Problèmes de transfert de chaleur Ecoulement hydraulique Etc.…



(1.7) B(u)+ g =0 Condition aux limites

L et B sont des opérateurs : d2( . )/dx2 ; . [d2( . )/dx2+. d2( . )/dy2] ; d/dx  Concept 2 On connaît toutes les équations mais on ne peut pas résoudre manuellement. L(u) + f = 0 Eqt. Diff.

MEF 

(Système d’équations algébriques) K.U=F (1.8)

B(u)+ g =0 Condition aux limites -Approximation(Géométrie du domaine compliquée….)

K – Matrice de rigidité du système U – Vecteur de déplacement F – Vecteur force

 Concept 3

K.U=F Propriété

Comportement

 Résolution U = K-1. F

(1.9)

Sollicitation

La nature des différents paramètres pour certaines applications est décrite dans le tableau 3.

Tableau 4 : Matrices et vecteurs caractéristiques pour certaines applications Propriété de K

Comportement U

Sollicitation F

Elasticité

Rigidité

Déplacement

Force

Transfert de chaleur

Conductivité thermique

Température

Source de chaleur

Fluide

Viscosité

Vitesse

Force interne

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

-9-

Chapitre1

Introduction

 Concept 4 Il est difficile d’obtenir un système d’équations pour la totalité du domaine. 

Diviser le domaine en un nombre d’éléments petits et simples,



Choisir un paramètre de champ (déplacement, température, vitesse) représentant les degrés de liberté (d.d.ls.) via une interpolation (approximation) polynomiale sur chaque élément,



Les éléments adjacents partagent les mêmes degrés de liberté aux points d’intersection (nœuds).

Un exemple de plaque est présenté dans la figure 1.9  Concept 5 Obtenir les équations algébriques pour chaque élément (Facile !) 

Mettre tous les éléments ensembles (processus d’assemblage de tous les éléments,

K

e

.U e   F e



K .U  F

(1.10)

 Concept 6 Résoudre le system d’équations pour obtenir les variables inconnues (déplacements) aux nœuds.

K.U=F



Résolution

U = K-1. F

(1.11)

B- Procédures 

Diviser la structure en un certain nombre de pièces (formant des éléments avec des nœuds)



Décrire le comportement des quantités physiques pour chaque élément (champs de déplacement - choix d’une forme d’approximation)



Relier (assembler) les éléments aux nœuds pour former un système d’équations approximatif pour toute la structure.



Résoudre le system d’équations comportant les inconnues aux nœuds (en génie civil les inconnues sont en général des déplacements).



Calculer les quantités désirées (en génie civil, il s’agit de contraintes et déformations) pour des éléments choisis.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 10 -

Chapitre1

Introduction

Exemples :

Eléments Barres (raidisseurs) Elément de panneaux rectangulaires Elément de plaque triangulaire (tôle de couvertures) Figure 1.6 : Décomposition de la structure d’une aile d’avion (premier exemple de discrétisation par Eléments finis de la publication (3)

P2 P1

(a)

(b)

Figure 1.7: Poutre à caisson (a) structure originale (b) Modèle éléménts finis Elément typique

(a)

Nœud typique

(b)

Figure 1.8 Plaque rectangulaire avec ouverture (a) Model de base (b) Discrétisation par EF de tailles différentes A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 11 -

Chapitre1

Introduction

C- Support informatique et logiciels incorporant la MEF: Exemples de logiciels La possibilité d’une application presque généralisée de la MEF, lui a permis d’être un outil puissant et versatile pour résoudre une large variété de problèmes. De ce fait, un grand nombre de logiciels professionnel ont été développés pour résoudre un grand nombre de problèmes d’ingénieur et de génie civil en particulier. On peut citer un certain nombre de logiciels pouvant être utilisés sur ordinateurs personnels:  SAP, ETAPS, ROBOT - Analyse statique/dynamique, linéaire/non linéaire, effet de température, sols et structures etc. ayant acquis une popularité d’utilisation dans les bureaux d’études de génie civil. 

ANSYS -Utilisation pour une analyse plastique, fracture mécanique, hautes températures et retrait, coques et systèmes tubulaires,  ABAQUS -Utilisation pour des analyses non linéaires et dynamiques,  PAFEC -Utilisation pour l’analyse des structures avec différents types de propriétés .Analyses de contraintes, déplacements, stabilité et fréquences, fatigue, optimisation et transfert de chaleur,  NAG library : Bibliothèque de sous-programmes d’Eléments Finis en langage Fortran pouvant êtres assembler et appelé par le biais d’un programme principale pour résoudre un problème précis. (Bon pour l’apprentissage didactique de la mise en œuvre de la MEF), -MATLAB

(Utilisation pour tous les problèmes d’ingénieurs, ayant acquis une

popularité parmi les ingénieurs de mécanique, électrotechniques et d’électronique) -MATHEMATICA (Ensemble de programmes et sous programmes pour l’analyse matricielle des structures et par MEF) On notera enfin que pour une utilisation rationnelle et économique des matériaux dans un processus de conception et de calcul, on est amené plus souvent, à utiliser les propriétés non linéaires (i.e. des équations non linéaires). Ceci nécessite des utilisateurs d’une certaine expérience qui retrouvent dans la plus part de ces logiciels la possibilité d’analyse de contraintes non linéaires pour les problèmes de grandes déformations, plasticité, retrais et stabilité. Support informatique de mise en œuvre de la MEF : Pour l’analyse par la MEF, la mémoire, l’espace disque, la vitesse du processeur, et la résolution graphique sont très importantes . A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 12 -

Chapitre1

Introduction

On peut cependant dire que l’analyse par la MEF demande beaucoup plus de puissance matérielle que de programmes. La MEF dans les packages des logiciels Multiphysic repose sur trois phases: 

Construire le modèle éléments finis, charges et conditions aux limites (prétraitement)



Assembler et résoudre le système d’équations (résolution)



Sélectionner et afficher les résultats (post-traitement)

D- Avantages de la MEF: La MEF étant versatile, elle présente des avantages par rapport aux autres méthodes numériques :  Elle est applicable à tout type de problème dit de champs, d’analyses de contraintes, de transfert thermique, etc,  Elle n’impose aucune restriction géométrique, le corps ou le domaine à modéliser peut avoir une forme quelconque,  Elle n’impose aucune restriction sur les conditions aux limites et le type de chargement,  Elle n’impose aucune restriction sur les propriétés du matériau. Ces propriétés ne sont donc pas réduites à l’isotropie, et peuvent changer d’un élément à un autre,  Les éléments possédant différents comportements (modèles mathématiques) peuvent être combinées (ex : voile-portique dans une structure 3D), donc un seul modèle d’élément fini (EF) peut contenir des barres, des poutres, des plaques, etc.  La structure modélisée par les éléments finis représente le possible à la structure réelle ou au domaine à analyser.  L’approximation peut être facilement améliorée en développant la taille du maillage en augmentant le nombre d’éléments. 1.9. Notion de modélisation et comparaison de la MEF à d’autres méthodes: A- Notion de modélisation - System continu et system discret: La figure 1.9 résume la notion fondamentale de modélisation dans un domaine prédéfini au préalable pour pouvoir faire une approximation du champ de déplacement avec une source d’erreur appréciable. Alors que le system continue est le domaine (structure) réel qui possède une infinité de degré de libertés (d.d.ls.), le system discret possède un nombre fini de d.d.ls. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 13 -

Chapitre1

Introduction Sans Intérêt !

Model idéal Modèle mathématique (Model physique)

Discrétisation

Solution discrète

Model Numérique

(Modèle discret) (Système d’équations algébriques)

U (Champ de déplacement) Approchée (Interpolation)

Erreur de la simulation=Modélisation + erreur de la solution

Figure 1.9: Schéma fondamental de modélisation B- La MEF et les autres méthodes d’analyse : On a donc le choix entre la recherche d’une solution précise pour un problème trop simplifié ou d’opter pour une solution approximative pour un problème pris généralement dans sa forme la plus complexe. Dans les deux cas on a créé une source d’erreur, cependant dans le cas de la MEF on a la possibilité de construire une bonne précision et une bonne exactitude du résultat, avec plus de travail sur le modèle pour son amélioration. Les méthodes d’analyses existantes pour les problèmes de champs à équations différentielles tel que les problèmes d’élasticité, des écoulements de fluides et de transfert de chaleur peuvent être classés comme suit, figure 1.10:

Méthodes d’analyses Méthodes Analytiques Méthodes exactes : Les méthodes des variables séparables , et de la transformation de Laplace

Méthodes Numériques Méthodes approximatives : Les méthodes de Rayleigh-Ritz , et de Galerkin- Bubnov

Solution numérique pour les équations

différentielles

Méthode des Eléments Finis

MEF Intégration numérique

Différences Finis

Figure1.10: Méthodes de résolution A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 14 -

Chapitre1

Introduction MEF (méthode numérique)

Méthode Analytique

Problème (domaine, structure) Faire des hypothèses de simplification

Laisser toutes les complications (au moins les termes importants même si le problème reste difficile)

Processus de Modélisation

!

Modèle simplifié

Modèle complexe Processus de résolution

Résoudre exactement (Solution Analytique)

Résoudre ! approximativement (Solution numérique) Solution approximative (pour un model ± "Exacte")

Solution Exacte (pour un modele approximatif)

!

Niveaux d’introduction d’approximations (source d’erreur)

Figure 1.11: Schéma comparatif des processus de modélisation C- Illustration du principe de la modélisation par un exemple

Structure

Modèle

P

Modèle discret

P Ah

 t . At  P

Ah

P

1

A1

2

A2

3 Poteau Sol

Ab

A3

Ab Support rigide

Représentation physique

Représentation par Elements Finis

Figure 1.12: Etapes pour la modélisation et l’analyse par la MEF d’un poteau à section variable A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 15 -

Chapitre1

Introduction

Pour l’illustration, l’application qui suit considère une section variable égale à "7A" à l’encastrement et "A" à l’extrémité libre, figure 1.13. u P

6A

x

1 L x/

LT

u3

4A

2A

2

A 7A

u2

P 4

3 L

L

x/

x/

Figure 1.13 : Exemple de discrétisation à base de 3 éléments de sections différents L’essence de la méthode est l’approximation du champ de déplacement "u"par un ensemble de déplacements "ui" pris dans des intervalles, formés par des éléments dont les extrémités sont des nœuds. C’est un problème isostatique, de contrainte uniaxiale. Chaque élément est considéré uniforme, linéaire élastique de longueur L, figure 1.12 et le déplacement "u" pour l’élément du milieu par exemple peut écrit en fonction des déplacements nodaux u2 et u3 selon l’équation suivante :

 x/  x/ u   1   u 2  . u3 L L 

(1.12)

x/ est une coordonnée axiale prise le long de l’élément L’équation (1.12) exprime une variation linéaire de u en fonction de x/ et u= u2 pour x/=0 et u= u3 pour x/=L. Il y va de même pour l’élément le plus à droite, avec des déplacements nodaux u3 et u4 et de manière similaire pour l’élément le plus à gauche, avec le déplacement u1=0 car l’extrémité gauche de la structure est fixée. A partir de la relation contrainte-déformation et la définition de la déformation comme une variation de la longueur par rapport à la longueur initiale, on obtient les expressions des contraintes axiales comme suit :

 12  E.

u2 L

;

 2 3  E.

u3  u 2 ; L

 3 4  E .

u 4  u3 L

(1.13)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 16 -

Chapitre1

Introduction

Le problème bien sûr est élémentaire et ne nécessite pas des écritures matricielles, (comme on le verra au fil des chapitres qui vont suivrent, et les déplacements aux nœuds peuvent être déduites comme suit :

u1  0 u3  u 2 

u2  PL 4A E

PL 6A E u 4  u3 

PL 2A E

(1.14)

Les déplacements peuvent être aussi exprimés en fonction de la longueur totale LT. les contraintes de l’élément obtenues à partir des équations (1.12) peuvent être vérifiées en divisant la charge "P" par les sections droites. Les résultats sont présentés dans la figure 1.14 et donnent une idée sur l’évolution de la solution par la MEF, par rapport à la solution exacte.

u



Déplacement axial Solution EF Solution exacte

Contrainte axiale Solution EF Solution exacte

x

x LT

LT

Figure 1.14 : Barre à section variable, discrétisée par éléments uniformes à 2 noeuds Pour cet exemple, la représentation graphique de la contrainte axiale en "escalier" montre que les contraintes sont exactes au milieu des éléments. Ailleurs, les contraintes sont représentées avec moins d’exactitude que les déplacements, car la plus part des éléments finis se basent sur une forme d’interpolation des champs de déplacement, et les contraintes sont généralement calculées (déduction) à partir des variations des déplacements.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 17 -

Chapitre 2

Éléments finis à une dimension

Projet du stade de Baraki-Alger : Toiture réalisée par des barres tubulaires en acier tubulaires

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

Chapitre2

Éléments finis à une dimension

APERÇU On examine dans ce chapitre, l’approche des éléments finis par la méthode directe en considérant des structures composés d’éléments discrets et ensuite par la méthode formelle en se basant sur le principe de l’énergie potentielle minimum. On présente des éléments simples dans une structure à une dimension pour étudier les différentes étapes pour la résolution par la MEF. Il s’agit d’établir d’une manière simple, la matrice de rigidité élémentaire, l’assemblage, l’introduction des conditions d’appuis et de charge et enfin la résolution. 2.1 Méthode directe pour les structures à éléments discrets: A. Introduction La méthode directe est une approche pour des systèmes discrets, basée sur la méthode des rigidités. Elle est de loin le procédé le plus simple pour introduire les concepts de base de MEF et présente les avantages suivants : 

Application de concept physique (équilibre des forces, conservation d’énergie, conservation de masse,…) directement à des éléments discrets.



Facile dans son interprétation physique.



Ne demande pas de concept ou de manipulation mathématique sophistiquée.

Son application est limitée à un certain nombre de problèmes pour lesquelles les lois d’équilibre et de conservation peuvent être facilement exprimées en termes des quantités physiques que l’on désire obtenir (déplacements). En général, les systèmes à treillis et les portiques sont constitués d’éléments discrets d’euxmêmes selon le sens physique et sont de parfaits exemples pour illustrer la méthode. B. Elément fini barre (Treillis plan de barres) : Le système à treillis est formé d’un ensemble membrures appelé barres, qui sont sollicités par des efforts agissant le long de leur axe moyen. Le schéma statique de la barre impose la présence de rotules aux deux extrémités (nœuds). Ce sont les premiers éléments présentés par la MEF suivi des éléments poutres assemblés en ossature.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 18 -

Chapitre2

Éléments finis à une dimension

Elément barre de référence

Eléments barre

(b) Noeud Support s (a) Figure 2.1 : (a) Structure à treillis représentant une ferme de toiture. (b) Modèle d’élément barre

C. Elément ressort linéaire : L’élément barre possède des caractéristiques similaires à celles d’un ressort élastique.

k

Figure 2.2 : Analogie barre- ressort On considère que chaque élément de structure se comporte comme un ressort élastique c'està-dire que la relation charge déplacement est linéaire.

i

k

j

F

u Figure2.3 : Déformation d’une ressort élastique On appelle k la raideur (rigidité) et correspond à la pente du graphe charge –déplacement. Charge F

Pente k Déplacement u de l’extrémité j Figure2.4 : Relation charge –déplacement d’un ressort élastique A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 19 -

Chapitre2

Éléments finis à une dimension

Connaissant la valeur de la rigidité et de la charge appliquée on a les relations :

F=k.u

1 u  .F k

et

(2.1)

C. Formulations en éléments finis : D.1. Matrice de rigidité élémentaire

x u1, u2: Déplacements aux nœuds 1 et 2

k

1

2

f1

f2 u1

f1, f2 : Forces (internes) aux nœuds 1 et 2

u2

Figure 2.5 : Ressort équivalent d’une barre à deux rotules Convention de signes : Cette même convention est adoptée pour les charges et les déplacements.

- ( f, u )

+ ( f, u ) Figure 2.6 : Convention de signes

Équilibre : Nœud 1

f1= -k. ( u2 – u1 )

Nœud 2

f2= k. ( u2 – u1 )

(2.2)

Pour l’équilibre des forces f1=- f2 Ecriture matricielle pour un élément :

k  k 

 k  u1   f1      k  u2   f 2 



K .U  F e

e

(2.3) e

Ke- Matrice de rigidité élémentaire Fe- Vecteur résultant de forces nodales A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 20 -

Chapitre2

Éléments finis à une dimension

Ue- Vecteur de déplacements aux nœuds Remarque1: Une colonne de Ke représente le vecteur charge qui doit être appliqué aux nœuds de l’élément pour obtenir un état de déformation où le degré de liberté nodal est égal à 1 alors que les autres sont nuls. Exemple : u1=0 et u2=1

F   

k  k

 k  0  1    k .   k  1  1

(2.4)

Donc le produit Ke . Ue = Fe représente la 2eme colonne de K. D.2. Matrice de rigidité d’un ensemble d’éléments (Assemblage) : Il est question dans cette étape de décider comment les matrices de rigidités élémentaires sont combinées ensembles pour former la matrice d’une structure composée de plusieurs éléments. Pour la simplicité de l’analyse étudions dans un premier temps le système de deux ressorts colinéaires [1] et [2].

f2(1)

f1(2)

K1

f1(1) u1

K2

[1]

F1

f2(2)

u2 [ 2 ]

u3

F2

F3

Figure 2.7 : Système à deux ressorts co-linéaires Elément 1:  k1  k  1

 k1   u1   f11      k1  u2   f 21 

 k2  k  2

 k2   u1   f12      k2  u2   f 2 2 

(2.5a)

Elément 2:

(2.5b)

Avec fie- Forces internes agissant au nœud i de l’élément e (i=1,2) En utilisant l’équilibre des charges aux nœuds 1, 2 et 3 respectivement, on aboutit à l’assemblage de la matrice de rigidité pour tout le system.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 21 -

Chapitre2

Éléments finis à une dimension F1 = k1 . u1 – k1 . u2 

F2 = - k1 . u1 + (k1 + k2 ) u2– k2 . u3

(2.6)

F3 = -k2 . u2 + k2 . u3 Soit : F1 = f11 F2 = f21+f12 F3= f22 Ecriture matricielle :

Elément 1

Soit

 k1 k1  k 2

 k1  k  1  0

 k2

Elément 2

 u1   F1   . u    F    2  2 k 2  u 3   F3 

0  k2

(2.7)

K.U=F

avec

 k1 K   k1  0

 k1 k1  k2  k2

0   k2  k2 

(2.8)

K – Matrice de rigidité globale du système. Remarque2 : Les matrices K et Ke sont symétriques et, moyennant une numérotation appropriée des nœuds K prennent la forme d’une "matrice bande". On verra par la suite, que ces deux caractéristiques en plus, de la propriété d’une matrice à être définie positive, vont contribuer efficacement à une optimisation dans la résolution du système d’équations. D.3. Technique d’assemblage par superposition (addition) : Bien que l’assemblage de la matrice de rigidité K ne soit pas difficile pour ce cas particulier, il le sera si la structure comporte un nombre important de ressorts. On se demande s’il est possible d’obtenir la matrice de rigidité globale K à partir des matrices de rigidité Ke de chaque élément ?

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 22 -

Chapitre2

Éléments finis à une dimension

Les équations (2.5) et (2.8) possèdent certaines similarités, cependant la procédure n’est pas immédiate. Les équations (2.5a) et (2.5b) sont du même ordre, mais ne peuvent être sommées directement car elles font référence à des déplacements (degré de libertés) différents. La technique est de procédé par insertion de zéros en termes de lignes et de colonnes, tel que chaque matrice soit étendue pour englober tous les degrés de liberté u1, u2, u3. Soit :

 f11   1  f2   0  

 k1  k  1  0

 k1 k1 0

0  u1    0  u2  0  u3 

0  2  f1    f 2  2

et

0 0  0

0 k2  k2

 u1   u  (2.9)  2 0  u3 

0  k2

On utilise ensuite la règle de l’addition de matrices. C’est une opération identique à l’utilisation de la superposition pour obtenir l’équation (2.7), sauf que les déplacements de la structure sont considérés élément par élément au lieu de nœud par nœud.

 k1  k  1  0

 k1 k1  k 2  k2

0  u1   f11       k 2  . u 2    f 21  f12  k 2  u 3   f 22 

(2.10)

C’est la même écriture matricielle obtenue par le concept de l’équilibre des forces. D.4. Procédure de résolution : La matrice de rigidité (2.8) est singulière, sont déterminant est nul et donc mathématiquement elle ne possède pas de matrice inverse. Le system d’équations ne peut être résolu ! L’explication physique est que la structure est dépourvue d’attaches ; l’application de n’importe quel chargement résulte en un "mouvement de corps rigide" du système. Pour y remédier, on doit imposer au système des conditions aux limites suffisantes.  Conditions aux limites : Condition géométrique – Condition limite essentielle Déplacement : u1=0 spécifiée

u2, u3 inconnues

Condition de charge – Condition aux limite naturelle Charges : F2= F3= F spécifiées

F1 réaction inconnue

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 23 -

Chapitre2

Éléments finis à une dimension

Remarque3 : Pour chaque nœud, on n’admet qu’un seul des deux conditions aux limites, déplacement ou force soit spécifié. Il n’est pas de nature de spécifier les deux à la fois. Si aucun des deux n’est connu, alors le problème est mal posé, ou tout simplement on n’a pas de problème à résoudre. On notera que si on n’a pas de condition géométrique, il n’y a pas de solution unique et dans notre cas de notre système, le ressort peut se déplacer selon l’axe des x sans élongation. Dans ce cas la matrice de rigidité globale devient singulière avec : ∑ lignes de la matrice =0



éléments kij linéairement dépendants

 Résolution : L’équation (2.9) peut être écrite sous la forme partitionnée suivante :  k1 k1  k 2

 k1  k  1  0

 k2

  0   F1   . u    P    2   k 2  u3   P 

0  k2

(2.10a)

A ce stade nous utiliserons la méthode dite "brutale" en enlevant les lignes et les colonnes de K, correspondant aux déplacements nuls. K se réduit alors à :

 k1  k 2  k 2 

k 2  u 2   p  .     k 2  u3   p 

et

F1= - k1. u2

(2.11)

u  Les inconnues sont : u  2  et la force de réaction F1 (si l’on veut). u 3  Remarque 4 : Une fois les déplacements calculés, les forces internes dans chaque élément peuvent être déterminées à l’aide des relations forces déplacements des ressorts. Si f1 et f2 représentent les forces internes dans les ressorts [1] et [2] respectivement, alors :

f1= k1 . (u2-u1)

et

f2= k2. (u3 – u2)

(2.12)

Ce qui complète le processus de résolution.  Résumé des Opérations pour une analyse par la méthode directe: Il est important de présenter les différentes étapes de l’analyse (Tableau 2.1), car elles restent inchangées quel que soit le problème et le type d’éléments. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 24 -

Chapitre2

Éléments finis à une dimension

Tableau 2.1 : Etapes d’analyse par la méthode directe Etape

Opération

1

-Former la matrice de rigidité Ke pour chaque élément.

2

-Assembler la matrice de rigidité globale K du système à partir de Ke pour chaque élément,

3

-Appliquer les conditions aux limites,

4

-Résoudre afin d’obtenir les déplacements, et les réactions si nécessaire, -Utiliser les relations charge –déplacement de l’élément pour obtenir

5

les forces internes. 2. 2. Application à des barres prismatiques uniformes : A. Hypothèses et position du problème On a présenté la méthode d’analyse d’une série de ressorts colinéaires par analogie à un système de barres à rotules. Nous voulons maintenant attribuer une valeur à la quantité k pour le cas d’une barre prismatique, figure 2.8, à partir de la relation contrainte – déformation par une analyse statique linéaire basée sur les suppositions suivantes : 

petites déformations



matériau élastique



charge statique

L’étude statique en élasticité linéaire consiste, connaissant les caractéristiques des éléments – barres (matériau et section) et la géométrie de la structure, à calculer : o les déplacements des nœuds non liés au support ; o les efforts de liaison exercés par le support sur les nouds concernés ; o les contraintes dans les barres. Cette étude peut apporter la plupart des informations sur le comportement de la structure, et elle peut être une bonne approximation pour beaucoup d’analyses. La barre prismatique uniforme possède les caractéristiques suivantes, longueur (L), section (A), module d’élasticité (E).

ui

uj

fi i

j

x

fj

L Figure 2.8 : Barre prismatique uniforme à deux rotules A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 25 -

Chapitre2

Éléments finis à une dimension

B. Formation du système d’équations d’équilibre u=u(x)

 =  (x)

déplacement ;

 (x) contrainte

déformation ;

Relation déformation – déplacement

  du

(2.13)

dx

Relation contrainte –déplacement (2.14)

  E La rigidité de la barre peut être obtenue à partir de (2.14) :

fi 

A.E .u i L



k

A.E L

Pour la barre de section uniforme l’équation (2.3) s’écrit sous la forme :

 fi  A.E  1  1  ui     .    f j  L  1 1 u j 

(2.15)

La relation (2.15) représente la relation entre les forces agissant aux extrémités de la barre et les déplacements au niveau des nœuds dans un système de coordonnées propre à l’élément dit "local ". Cependant, puisque les éléments d’une structure plane se présentent sous des angles différents les uns par rapport aux autres, il est nécessaire d’en tenir compte. Dans l’élaboration de la matrice de rigidité globale K du system, il est nécessaire d’écrire la matrice de rigidité pour chaque élément de la structure non pas en fonction du système de coordonnées local (x,y), mais en fonction du système de repère globale (X,Y), référence pour l’ensemble de la structure. C. Transformation d’une base locale vers une base globale Cette transformation est très importante avant l’assemblage des matrices et vecteurs caractéristiques pour tenir compte de l’aspect vectoriel du champ de déplacement. Toutes les caractéristiques élémentaires sont aisément programmées dans un système de coordonnées locales, pour minimiser l’effort de programmation.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 26 -

Chapitre2

Éléments finis à une dimension

En définitif si un système de coordonnées local est utilisé, il est nécessaire avant d’assembler les équations élémentaires de les transformer de façon à ce qu’elles fassent référence à un système de coordonnée global. La figure 2.9 montre, une barre inclinée d’un angle  par rapport au système d’axe global et les conventions de signes adoptés avec  pris positif dans le sens contraire des aiguilles d’une montre.

j

x

y

Local

X Global

+

ui’ vi’ v 'i’ i

i 

ui

+

+

Y

Figure 2.9 : Systèmes de coordonnées de transformation Le tableau 2.2 présente les degrés de libertés à considérer dans les deux cas de repères, local et global. Tableau 2.2 : Degrés de libertés selon les deux différents axes de coordonnées local et global. Base locale à l’élément x, y

Base Globale à la structure X,Y

u’i ,( v’i)

ui ,vi

1 ddl par nœud

2 ddls par nœud

A noter que le déplacement latéral vi’ ne contribue pas à l’allongement de la barre, conformément à la théorie linéaire.

 Relations de transformation : ui’ = ui cos + vi sin  =< c

ui s> vi

vi’ = -ui sin  + vi cos  =

ui

(2.16)

vi A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 27 -

Chapitre2 où 

Éléments finis à une dimension c= cos

s= sin 

;

Forme matricielle :

u i’

c

s

ui

-s

c

vi

=

vi’

(2.17a)

soit , ui’ = T ui

(2.17b)

La matrice de transformation étant :

T =

c

s

-s

c

(2.18)

Elle est orthogonale, soit T-1 = TT Pour les deux nœuds de l’ EF barre, on a :

u i’ vi’ = u j’ vj’

Soit

u’ = T u

c -s 0 0

avec

s c 0 0

0 0 c -s

0 0 s c

ui vi uj vj

T

0

0

T

T =

(2.19)

(2.20)

Les charges aux nœuds sont transformées de la même façon. f’ = T . f

(2.21)

D. Matrice de rigidité dans un espace bidimensionnel –2D Un élément est caractérisé par une matrice de rigidité reliant les efforts Ni et Nj exercés par les nœuds (i,j), aux déplacements de ces nœuds.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 28 -

Chapitre2

Éléments finis à une dimension

Nj

y E,A L

x Ni

Figure 2.10 : Elément barre sous sollicitation axiale Dans un système de base locale, la relation (2.15) s’écrit :

EA L

1 -1

-1 1

u i’

fi’ =

u j’

Ni =

fj’

Nj

(2.22)

En étalant cette écriture à l’ensemble des degrés de liberté on obtient :

EA L

1 0 -1 0 0 0 -1 0 1 0 0 0

0 0 0 0

u i’ vi’ u j’ vj’

=

Ni 0 Nj 0

Soit,

k’ . u’ = f’

(2.23a)

(2.23b)

En utilisant les transformations (2.20), on obtient :

k’. T . u = T . f En multipliant les deux côtés de l’égalité par T T et sachant que TT T = I , on obtient : TT. k’. T . u = f On obtient alors l’écriture de la matrice de rigidité Ke dans le système de coordonnée globale, sous la forme suivante : Ke = TT. k’ .T

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(2.24)

- 29 -

Chapitre2

Éléments finis à une dimension

qui représente une matrice 4x4 symétrique La forme explicite de Ke devient:

EA L

Ke =

ui

vi

uj

vj

c2

c.s

-c2

-c.s

c.s

s2

-c2 -c.s

-c.s

-s2

-c.s

c2

c.s

-s2

c.s

s2

(2.25)

Calcul des cosinus directeurs c et s : c= cos =

X j Xi , L

s =sin  =

Y j Yi L

La matrice de rigidité de la structure est assemblée en utilisant les matrices de rigidité élémentaires de la même manière que dans le cas des barres colinéaires. E. Contraintes et forces internes élémentaires :

ui’

e =E e =E B

u j’

=E L

c s

0

0

0

c

s

0

ui vi uj vj

Soit :  Contraintes  e=

E  -c -s L

c

s>

ui vi uj vj

(2.26)

 Forces internes dans chaque élément [e] aux nœuds i,j

f

e

 Ni j

 A.E   . c  L 

s

i j

u i  u j  .  vi  v j 

e

(2.27)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 30 -

Chapitre2

Éléments finis à une dimension

2.3 Application aux ossatures planes de poutres : A. Position du problème Pour une poutre droite d’une ossature plane Figure2.11, qui se différencie de l’élément précédent par le fait que les nœuds peuvent transmettre des moments aux extrémités, la matrice de rigidité fait aussi intervenir le moment quadratique I.

y D h2 H

x

A

H

B

E, I2

A-A

A

(a)

E, I3

E, I1

h1

B

L b

A

A

(b)

C

Figure 2.11 : (a) Portique plan à éléments poutres AB, BC, CD. (b) Cantilever fixé en A et chargé en B B. Matrice de rigidité élémentaire : En négligeant l’influence de l’effort tranchant, la matrice de rigidité peut être obtenue en utilisant les formules de "Bresse".

Tj z j Ti

Mj

j z i

vi

Mi i

y x

vj j

y (a)

x

L i

(b)

Figure 2.12 : Elément poutre dans le plan xy, ses charges nodales (a) et ses degrés de liberté (b) . Conventions de signes : Les rotations nodales contiennent l’indice z pour marquer leurs représentations vectorielles selon l’axe z, normal au plan xy. Les charges nodales sont prises chacune positive s’ils agissent dans la même direction que leurs degrés de libertés respectifs ; Figure 2.12. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 31 -

Chapitre2

Éléments finis à une dimension

 6.E.I   4.E.I   6.E.I   2.E.I  M i   2  . vi  . z i   . v j   . z j  L   L   L   L  et

 6.E.I   2.E.I   6.E.I   4.E.I  M j   2  . vi  . z i   . v j   . z j  L   L   L   L  À partir de l’équilibre statique :

 M i  M j   12.E.I     3  . vi  6.E2.I . z i   12.E3 .I . v j   6.E2.I . z j Ti   T j   L  L   L   L     L  Ces équations peuvent être écrites sous la forme matricielle suivante :  Ti   M    i  E.I   .   3  Tj  L   M j   Résumé sous la forme suivante :

12 6.L  12 6.L 6.L 4.L2  6.L 2.L2  12  6.L 12  6.L 6.L 2.L2  6.L 4.L2

F  K .U e

e

  vi     .  z i    vj      z j 

(2.28)

e

e

Ceci défini la matrice de rigidité K pour un élément poutre horizontal. Remarque 5: Si l’extrémité gauche de l’élément poutre est fixée tel que vi=0 et zi=0, on obtient une structure avec vj=0 et zj=0 comme des degrés de libertés actives et la matrice de rigidité de l’élément poutre cantilever figure2.11b, est la sous matrice inférieure droite (2.29) de dimension 2 par 2

  E.I e K = 3 . L   

12 6.L  12 6.L 2 6.L 4.L  6.L 2.L2  12  6.L 12  6.L 6.L

2.L2

 6.L

 vi  . zi  vj  4.L2   z j

(2.29)

C. Matrice de rigidité d’un élément poutre à orientation arbitraire. Pour permettre à un élément poutre de s’allonger en même temps que de fléchir, on ajoute les déplacements ui et uj au vecteur degré de liberté Ue et étendre Ke à une dimension de 6x6 en introduisant la rigidité axiale AE/L à partir de (2.15). Après combinaison avec la matrice de rigidité axiale (élément barre), on obtient Ke : A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 32 -

Chapitre2

Éléments finis à une dimension

ui

vi

i

uj

vj

j

EA L

0

0

- EA

0

0

0

12EI

6EI

0

- 12EI 3

6EI

0

L 6EI

0

L - 6EI

L 2EI

L

L

3

2

L 4EI

2

2

2

L

L

- EA

0

0

EA L

0

0

- 12EI 3

- 6EI2

0

12EI

0

L 6EI

L 2EI

0

L - 6EI

L 4EI

L

L

e

K =

L

L

2

L

(2.30)

0 3

- 6EI2

2

L

Pour un élément de poutre orienté arbitrairement tel que l’élément BC de la figure 2.11a, il est nécessaire d’étendre la matrice K

e

pour permettre la transformation des degrés de liberté

selon les deux déplacements u et v dans le système de base globale. Les moments ne sont pas concernés. Les relations de transformation utilisées pour le cas de l’élément barre (2.24) restent valable pour ce cas d’élément poutre sauf que l’équation (2.30) doit être considérée comme k’ et la matrice de transformation T est :

T=

c

s

0

0 0

-s

c

0

0 0

0

0

0

1

0 0

0

0

0 0

c s

0

0

0 0

-s

c

0

0

0

0 0

1

0

0

soit

t 0 T   0 t  (2.31)

avec

 c s 0 t   s c 0   0 0 1 

En utilisant (2.30) et (2.31) et en effectuant K =TT.Ke.T, on obtient l’écriture matricielle suivante à utiliser dans le cas d’éléments de portique

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 33 -

Chapitre2

Éléments finis à une dimension

  A.c 2  12I .S 2 ( A  12I ).c.s  L2 L2   12I  A.s 2  2 .c 2  L   E K L      SYMETRIQUE    

6I 12I  .s  ( A.c 2  2 ).s 2 L L 6I .c L

 (A 

12I .c.s L2 )

6I .s L

4I

A.c 2 

12I 2 .s L2

6I  .s L    6I  .c 12I  ( A.s 2  2 .c 2 ) L   L  2I  6I  .c  L  6I  .s 12I ( A  2 ).c.s L  L   6I   .c 12I A.s 2  2 .c 2 L  L   4I  12I  ( A  2 )c.s L



Remarque 6 : Cas particulier : Matrice de rigidité d’une poutre arbitrairement orienté dans le plan :

v’j

Y y v’i ’z i

’z j

x 

j

X

i

La relation entre les degrés de libertés de la base locale et ceux de la base globale s’écrit :

v    v  

' i ' zi ' j ' zj

  s    0   0   0 

c 0 0 1

0 0

0 0

0 0 s c 0 0 0 0

 ui    0  vi  0  zi    0  u j   1  v j     zj 

  E.I ' K = 3 . L   

 s 0 T 0  0

avec

c 0 0 0 0 0 1 0 0 0 0 0  s c 0  0 0 0 0 1

12 6.L  12 6.L 2 6.L 4.L  6.L 2.L2  12  6.L 12  6.L 6.L

2.L2

 6.L

  .  2  4.L 

En utilisant la relation (2.24) on obtient la matrice de rigidité de cette élément poutre: A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 34 -

Chapitre2

Éléments finis à une dimension

 12 s 2  12 s.c  6 L.s  12 s 2 12 s.c  6 L.s    2 6 L.c 12 s.c  12 c 2 6 L.c   12 s.c 12 c 6 L.c 4 L2 6 L.s  6 L.c 2 L2  E.I   6 L.s e K  3 .  L   12 s 2 12 s.c 6 L.s 12 s 2  12 s.c 6 L.s   12 s.c  12 c 2  6l.c  12 s.c 12 c 2  6 L.c    6 L.c 2 L2 6 L.s  6 L.c 4 L2    6 L.s

2.4 Généralisation : A. Elément de poutre bidimensionnelle- 2D On modifie les termes de la rigidité de flexion pour tenir compte des déformations transversales de cisaillement, produisant ainsi un élément poutre de Timoshenko. Pour une déformation dans le plan xy la matrice K s’écrit sous la forme :  X  0   0 K   X  0   0

0 Y1 Y2

0 Y2 Y3

X 0 0

0  Y1  Y2

0  Y1 Y2

0  Y2 Y4

X 0 0

0 Y1  Y2

 ui v  i   zi  0 uj  Y2  v j  Y3   zj 0 Y2 Y4

(2.32)

avec X 

A.E L

Y1 

12E.I z (1   y ).L3

Y2 

6 E.I z (1   y ).L2

(2.33) Y3 

(4   y ).EI z (1   y ) L

Y4 

(2   y ) EI z

y 

(1   y ) L

12EI z .k y A.G.L2

G  Module de cisaillement ; et A/ky est la surface effective de cisaillement pour la déformation de cisaillement transversal dans la direction y. Des exemples de valeurs usuelles adoptées sont :

k y  1.2

Pour une section rectangulaire

k y  2.0

Pour une section tubulaire

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 35 -

Chapitre2

Éléments finis à une dimension

Remarque 7 : Lorsqu’un élément devient de plus en plus élancé, y tend vers zéro, et les coefficients de flexion Yi se réduisent aux coefficients de la relation (2.29), ou la déformation de cisaillement transversale a été négligée. Dans la relation (2.29) le degré de liberté de rotation défini la valeur nodale, de la pente dv/dx et la rotation de la section droite ensembles. Si la déformation transversale de cisaillement est présente, la pente et la rotation de la section droite ne sont pas les mêmes, alors zi et zj doivent être considérées comme des rotations des sections droites de la poutre aux nœuds. B- Elément de poutre tridimensionnelle -3D On considère six degrés de liberté par nœud : trois déplacements et trois rotations, comme le montre la figure 2.13. y L xi

yi vi

ui

yj j

i wi z

zi

wj

vj xj uj

x

zj

Figure 2.13 : Elément fini poutre selon l’axe x d’un système de coordonnées orthonormées, avec des degrés de libertés définissant, un déplacement, une rotation et une déformation latérale selon y et z. Les degrés de liberté w et y tiennent compte des déformations latérales dans le plan zx. Le degré de liberté x tient compte de la torsion autour de l’axe x, pour lequel le coefficient de rigidité est G.K/L, avec K une propriété tenant compte de la forme et de la dimension de la section droite Remarque 8 : K peut être égal à J, le moment d’inertie polaire d’une section droite par rapport à son centre, seulement pour les sections circulaires ou tubulaires. Pour les sections à parois minces, tel que tel que les sections standards de profile I ou U, K représente une petite fraction de J.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 36 -

Chapitre2

Éléments finis à une dimension

Matrice de rigidité K, généralisée à un système 3D :

X         K         

0

X

0

0

0

0

0 0 0  Z2

Y2 0

0 0

 Y1 0

0  Z1

0 0

0  Z2

S

0

0

0

0

0

S

0

Z3

0 Y3

0 0

0  Y2

Z2 0

0 0

Z4 0

X

0 Y1

0 0 Z1

0 0 0 S

0 0 Z2 0 Z3

0

0

0

Y1

0 Z1

0

Symétrique

0  ui Y2  vi 0  wi  0   xi 0   yi  Y4   zi 0  uj   Y2  v j 0  wj  0   xj  0   yj Y3   zj

(2.34)

où S=G.K/L , les termes X et Yi, sont définis selon (2.33) et les termes Zi sont toujours définies selon (2.33) mais avec une permutation des indices. Exemple:

Z1 

12E.I y (1  z ).L

3

et

z 

12EI y .k z A.G.L2

Remarque 9 : Cette forme généralisée d’écriture de K pour un élément fini de type poutre est adoptée par la plupart des logiciels de calcul des structures. La relation (2.34) doit être considérée comme k’ car prise dans sa base locale. La relation de transformation T pour une base globale utilisée pour calculer ke est :

  T    

   , avec   0 0 0 t  t 0 0 0 0 t 0 0 0 0 t 0

 c t   s  0

s c 0

0 0  1 

2.5 Influence de la numérotation des nœuds sur la largeur de la bande de K : En remarque à partir des exemples que la matrice de rigidité globale possède des termes dans la diagonale principale ii et jj et des termes dans des positions en dehors de la diagonale ij et ji. Pour garder une 'largeur de bande' de la matrice K la plus petite possible, les nœuds doivent être numérotés tel que la différence maximale dans la numérotation des nœuds soit la plus petite possible. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 37 -

Chapitre2

Éléments finis à une dimension

On est à la recherche d’une "topologie" qui permet de réduire la capacité de stockage de la matrice de rigidité du système dans la mémoire de l’ordinateur en vue d’une exécution rapide lors de la résolution du système d’équations. Ceci est d’autant plus important que le nombre de degrés de liberté devient de plus en plus élevé. Toutes ces considérations sont discutées par une illustration d’exemple de structure, figure 2.14, formée d’élément finis barres à deux nœuds avec deux degrés de libertés u et v par nœud. P

9 7

6

8 5 Numéro de l’élément

3

2

4

1 9

5

3

9

3

8

7

6 5

4 3

2 1

6

7

6

(b)

(a)

(c)

2

5

4

1

3

2 2

1

1

4 8 5 4 6

Figure 2.14 : (a) structure à élément de barres sous charge horizontale P supportée par deux bases à appuis doubles (b) topologie 1 : numérotation horizontale (c) topologie 2 : numérotation verticale A- Matrice de connectivité : Elle donne une information sur la topologie de la structure, car elle contient le numéro de l’élément et les numéros des nœuds qui lui correspondent. Elle est importante pour l’assemblage des éléments finis du modèle et elle est présentée sous forme de tableau pour ces deux cas.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 38 -

Chapitre2

Éléments finis à une dimension

Tableau2.3 : matrice topologique 1 Elément 1 2 3 4 5 6 7 8 9

Tableau2.4 : matrice topologique 2

Nœuds 1 2 1 2 1 3 2 3 2 4 3 4 3 5 4 5 4 6 5 6

Elément 1 2 3 4 5 6 7 8 9

Nœuds 1 2 1 6 1 2 2 6 5 6 2 5 2 3 3 5 4 5 3 4

Pour une raison de commodité de présentation on représente les différents termes des matrices de rigidités élémentaires par le numéro de l’élément (e). Un exemple de matrice rigidité élémentaire de chaque topologie est considéré pour l’illustration (Elément (1)):

Topologie 2 : Elément (1)

Topologie 1 : Elément (1)

K

(1)

U1 V1  ( 1) (1)  (1) (1)   (1) (1)   (1) (1)

U2 (1)

(1) (1) (1)

V2 (1)  ( 1 )  (1)   (1) 

U1 U1 V1 U2 V2

  (1)  K   

V1

U6

V6

( 1) (1) (1) (1) (1) (1)

(1) (1) (1)

(1)

(1)

(1)  ( 1 )  (1)   (1) 

(1)

U1 V1 U6 V6

En généralisant cette forme de représentation à l’ensemble des éléments pour chaque type de topologie et en procédant à l’assemblage de la matrice globale pour chaque système, on obtient les matrices de rigidité pour la topologie1 et, la toplogie2 respectivement (figures, 2.15 et 2.16). B- Semi – bande : En se référant aux figures 2.15 et 2.16, il y a un total de 12 degrés de libertés (6x 2ddls) et si on n’exploite pas l’avantage de la matrice symétrique bande, on aura besoin de réserver une dimension de matrice pour 122 termes. En tenant compte de la semi – bande, on aura besoin de 6x12 termes seulement pour le modele de topologie1 dont la numérotation des nœuds a été bien choisie, et 12x12 pour le model topologie2 dont la numérotation est mal adapté.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 39 -

Chapitre2

Éléments finis à une dimension

La semi – bande est calculée comme suit : B= f. (d +1)

(2.35)

f- nombre de degré de liberté par nœud d- différence maximale entre numéros de nœud (en inspectant tous les éléments de l’assemblage)

Topologie 1

Topologie 2

B1= 2.( 2+1) =6

B2= 2.( 5+1) =12

Profile de la matrice de rigidité de la structure : On définit une limite dite "ligne de ciel" en englobant les termes non nuls de chaque colonne (figure 2.15). A cause de la symétrie, on n’a besoin de stocker que la matrice triangulaire supérieure pour avoir toute l’information sur la matrice globale. Aussi pour la résolution des équations on ne traite que les termes au-dessous de la "ligne de ciel". On peut donc stocker tous les coefficients nécessaires dans un seul vecteur, en prenant les colonnes dans l’ordre tous en incluant les coefficients de la ligne de ciel jusqu’au bas de la diagonale. Les zéros entre la ligne de ciel et la diagonale sont retenues car les méthodes de résolutions directes (non itératives) transforment ces zéros situés au-dessous de la ligne de ciel en valeurs non nulles. Le nombre de termes stocker sous ce format est appelé "profil" de la matrice de rigidité de la structure. Les logiciels on besoin de savoir quels sont les termes dans ce vecteur de stockage, qui représentent les termes de la diagonale de la matrice carrée (nxn). Cette information est donnée sous forme d’un vecteur auxiliaire à n nombres qui donnent la position de chaque colonne. Comparaison des profils des deux matrices de rigidité respectives: Topologie 1 [ 1 2 3 4 5 6 5 6 5 6 5 6 ]  ∑ (nombres) = 54 Topologie 2 [ 1 2 3 4 3 4 3 4 7 8 11 12 ] ∑ (nombres) = 62 La matrice à toplogie2 procède un profil plus grand, ce qui reflète la mauvaise numérotation adoptée pour ce cas.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 40 -

Chapitre2

Éléments finis à une dimension

Matrice de rigidité : Ligne de ciel U1

V1

U2

V2

U3

V3

U4

V4

U5

V5

U6

V6

(1) (2)

(1) (2)

(1)

( 1)

(2)

(1) (2)

(1) (2)

( 1)

( 1)

( 1 )

(1)

(3)

(4)

(4)

(1)

(1) (3) (4) (1) (3) (4)

(3)

(1)

(1) (3) (4) (1) (3) (4)

( 3 )

( 3 )

( 4 )

( 4 )

( 2 )

( 2 )

( 3 )

( 3 )

( 5 )

( 6 )

( 6 )

U3

( 2 )

( 3 )

( 3 )

( 5 )

( 5 )

( 6 )

( 6 )

V3

( 4 )

( 4 )

(2) (3) (5) (6) (2) (3) (5) (6) ( 5 )

( 5 )

( 2 )

(2) (3) (5) (6) (2) (3) (5) (6) ( 5 )

( 7 )

( 8 )

( 8 )

( 4 )

( 5 )

( 5 )

( 7 )

( 7 )

( 8 )

( 8 )

( 6 )

( 6 )

( 7 )

(4) (5) (7) (8) (4) ( 5) (7) (8) ( 7 )

( 7 )

( 4 )

(4) 5) (7) (8) (4) (5) (7) (8)

( 9 )

( 6 )

( 7 )

( 7 )

(6) (7) (9) (6) (7) (9)

( 9 )

( 6 )

(6) (7) (9) (6) (7) (9)

( 9 )

( 8 )

( 8 )

(9)

( 9 )

(8) (9)

( 9 ) V5 (8) (9) U6

( 8 )

( 8 )

( 9 )

( 9 )

(8) (9)

(8) (9)

(2)

(2)

U1

(2)

V1 U2 V2

U4 V4 U5

V6

Diagonale

Figure 2.15 : Topologie 1 U1

V1

(1) (2)

(1) (2)

(2)

U2

( 2)

(1)

(1)

(1) (2)

(1) (2)

( 2)

( 2)

(1)

(1)

( 2 )

(2)

(2)

(2)

(2) (5) (2) (5)

(3) (6) (3) (6)

V2

(2) (5) (2) (5)

(3) (6) (3) (6)

( 6 )

( 6 )

( 6 )

( 6 )

( 5 )

( 5 )

U3

V3

U4

V4

U5

V6

(6)

( 5 )

( 5 )

(3)

(3)

( 6 )

( 6 )

( 5 )

( 5 )

(3)

(3)

(6) (7) (9) (6) (7) (9) ( 9 )

(6) (7) (9) (6) (7) (9) ( 9 )

( 9 )

( 9 )

( 7 )

( 7 )

( 9 )

( 9 )

( 7 )

( 7 )

( 9 )

( 9 )

( 7 )

( 7 )

(8) (9) (8) (9)

(8 ) (9) (8 ) (9)

( 8 )

( 8 )

( 8 )

( 8 )

( 8 )

( 8 )

(4) (5) (7) (8) (4) (5) (7) (8) ( 4 )

(1) (4) (3)

( 4 )

(1) (4) (3)

( 1 )

( 1 )

( 3 )

( 3 )

( 1 )

( 1 )

( 3 )

( 3 )

( 4 )

( 5 )

U6

(6)

(4) (5) (7) (8) (4) (5) (7) (8) (4)

( 5 )

V5

( 7 )

( 7 )

( 8 )

( 8 )

U1 V1 U2 V2 U3 V3 U4 V4

(4)

(4)

U5 (4)

(4)

(1) (4) (3) U6 (1) (4) (3) V6

Figure 2.16 : Topologie 2 A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

V5

- 41 -

Chapitre2

Éléments finis à une dimension

2.6 Méthode formelle -fonction d’interpolation & fonction de forme A. Introduction Les formulations de matrices de rigidité des éléments discrets étudiés au chapitre 2 par la méthode directe, peuvent être déterminées par la construction d’un champ de déplacement par le biais de "fonction d’interpolation". L’application de la méthode qui représente l’essence même de la MEF, sera mise en décrite dans un premier temps, pour les problèmes d’élasticité unidimensionnels, afin de confirmer les écritures des matrices de rigidité de l’élément barre et de l’élément poutre, adoptées lors de l’utilisation de la méthode directe. Son champ d’application sera ensuite étendu, à des éléments finis plans pour une approximation des domaines continus. On procède dans un premier temps, à la discrétisation du domaine ou structure en éléments finis de forme géométriques simples figure 2.15. (Domaines Ωe)

(1D) x=a

Ω Domaine

Elément

x x=b

(a)

x=a

Noeuds

x x=b

Figure 2.15 : Représentation de domaines pour un problème unidimensionnel -1D Chaque élément est décrit par un champ de déplacement représenté par des fonctions approchées dites "d’interpolation". B. Position du problème et définition Une fonction d’interpolation défini la variation du déplacement sur l’élément fini et doit constituer une approximation raisonnable de la réalité (comportement structural) Les questions à poser sont les suivantes:  Quelle fonction choisir ?  Comment construire cette fonction ? La réponse par la MEF à la première question est de définir une fonction de déplacement approchée u pour représenter le champ de déplacement réel v, qui n’est pas difficile à trouver dans la plupart des problèmes complexes. Cette fonction est dite d’interpolation car v, est égale à u au niveau des nœuds des éléments finis seulement ce qui représente pour la MEF une forme d’approximation dite nodale.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 42 -

Chapitre2

Éléments finis à une dimension

Il existe un grand nombre de fonctions possible toutefois on préfère les polynômes pour les raisons suivantes : 

Ils sont très maniables pour la programmation.



Ils donnent la possibilité d’augmenter la précision en augmentant l’ordre du polynôme.

On considère pour une raison de simplicité de l’illustration, une fonction approchée unidimensionnelle u(x) qui peut être définie comme suit : u(x)= a1+a2.x+a3.x2+……. +an.xn-1

Polynôme simple :

(2.36)

C. Approximation nodale La fonction u(x) (2.36), construite sur la base de fonctions polynomiales peut-être écrite sous forme de vecteurs:

u(x)=   . .   a n 

(2.37)

Les coefficients ai, sont les paramètres de l’approximation et à ce stade, ils n’ont aucune signification physique, sauf si on coïncide u(x) avec la solution exacte v(x) au niveau des nœuds. On aura alors, le système d’équation suivant en tenant compte de la position (coordonnées) des nœuds:

u( x1 )  a1  a2 x1  a3 x12  .......  an x1n1  V ( x1 )  V1 u( x2 )  a1  a2 x2  a3 x22  .......  an x2n1  V ( x2 )  V2

…

…

…

…

…

…

(2.38)

…

u( xn )  a1  a2 xn  a3 xn2  .......  an xnn1  V ( xn )  Vn Ecriture matricielle : 1  1 .  . .  . 1 

x1 x2 . . . . xn

x12 x 22 . . . . x n2

. . . . . . .

. . . . . . .

. x1n 1   a1  V1      . x 2n 1  a 2  V 2  . .  .   .      . . . .    .  Soit . .  .   .      . .  .   .      . x nn 1  a n  V n 

A.a=V

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(2.39)

- 43 -

Chapitre2

Éléments finis à une dimension a = A-1.V

Si A n’est pas singulière

(2.40)

En remplaçant a dans la relation (2.35) on obtient l’écriture de la fonction sous la forme suivante : u(x)=A-1.V

Ou encore sous l’écriture suivante : u(x) = < N1(x) N2(x) …………………… Nn(x) > V

(2.41)

Ni(x) pour i=1,....,n sont appelées fonctions de forme Soit avec

u(x) = < N > V

=. A-1

x

Cette forme d’approximation est dite approximation nodale, avec comme variables nodales les déplacements Vi, et comme fonctions d’interpolations nodales les fonctions de formes Ni(x). Remarque 1 : Sachant que u(xi) = Vi on peut déduire une des propriétés les plus importantes de la fonction de forme à partir de la relation (3.3) qui est la suivante :

0 N j ( xi )   1

si

ij i j

(2.43)

C’est-à-dire que la fonction de forme possède une valeur égale à 1 au nœud considéré et elle est nulle sur les autres nœuds. La fonction de forme possède des propriétés intrinsèques, qui sont utilisées lors de la simplification des calculs intégrales, dans l’étape de construction de la matrice de rigidité (voir chapitre 4). D. Choix de la fonction de déplacement: La réponse à la deuxième (P42) question représente une étape importante de la MEF, car elle détermine la performance de l’élément fini dans l’analyse des structures. Le choix est fait avec précaution pour la sélection de cette fonction qui : 

doit avoir le nombre de coefficients inconnues (ai) égal au nombre total de degrés de liberté de l’élément,



ne doit privilégier aucun sens ou direction,



doit permettre un mouvement de corps rigide (sans déformation interne).



doit être capable de représenter des états de contraintes planes et de déformations planes,

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 44 -

Chapitre2 

Éléments finis à une dimension

doit satisfaire la compatibilité des déplacements (continuité) le long des limites d’interconnections entre les éléments finis.

Remarque 2: Le polynôme d’approximation doit comporter autant de termes (a i) de la série que de nombre de nœuds par élément. E. Principe de l’approximation par éléments finis: La construction d’une fonction de déplacement pour la totalité du domaine continu, conduit à un nombre très élevé de points xi et la solution qui en découle est dite "instable". Pour contourner ce problème, on construit la fonction u(x) par morceaux. On procède à la division du domaine Ω en un nombre fini de sous domaines Ω e, figure 2.15, sur lesquels la construction de la fonction u est plus simple. Le procédé est donc le suivant : 

Discrétisation du domaine Ω en sous domaines Ω e (éléments).



Définition d’une fonction de déplacement ue par la méthode d’approximation nodale, qui peut être différente sur chaque élément.

La fonction u pour l’ensemble du domaine Ω sera : u(x)=∑ue(x)

(2.44)

Les points où ue(x) =V (solution exacte) sont les nœuds d’interpolation. Les coordonnées xi représentent les coordonnées nodales. Les valeurs Vi = u(xi) représentent les valeurs nodales. Remarque 3: La fonction approchée ue(x) sur chaque élément, doit être continue sur Ω e et doit satisfaire la condition de continuité entre les différents éléments finis. Approximation par élément fini unidimensionnel linéaire à 2 nœuds: Pour illustration, on considère une approximation linéaire, figure 2.16, applicable aux problèmes de barres, solution exacte v(x)-inconnue V1 V2

x=x1

[e]

x=x2

x

Figure 2.16 : Approximation d’un champ de déplacement inconnu par un polynôme du 1er ordre. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 45 -

Chapitre2

Éléments finis à une dimension

En considérant que la solution exacte est approchée par un polynôme de premier ordre, on a une représentation sur l’élément par la fonction suivante : u(x) = a1 + a2.x

(2.45)

Les coefficients a1 et a2 définissent le segment de droite et si celui-ci passe par (x1, u1) et (x2, u2) alors : u1= a1 + a2.x1 u2= a1 + a2.x2

(2.46)

Après résolution du système d’équations (2.46) et en substituant dans (2.45) on obtient : u ( x) 

u1 x 2  u 2 x1 u  u1 ( 2 ). x x 2  x1 x 2  x1

(2.47)

Cette écriture est la plus simple, cependant en éléments finis on complique un peu la chose en mettant en évidence les valeurs nodales (déplacements nodaux) u 1 et u2, de la fonction u(x), qui sont plus importants dans l’analyse des structures. On réarrange (2.47) sous une forme équivalente : u ( x)  u1 . (

x2  x x  x1 )  u2 .( ) x 2  x1 x 2  x1

(2.48)

Les fonctions de formes et leurs propriétés : Du point de vu mathématique, le choix du couple de paramètres a 1, a2 ou u1, u2, dans l’approximation linéaire, importe peu. On remarquera qu’au lieu de multiplier a1, a2 par les termes des polynômes ‘1’ et ‘x’ respectivement pour obtenir (2.45), on doit multiplier u1 et u2 par les polynômes x2  x x 2  x1

et

x  x1 pour enfin définir deux polynômes linéaires suivants : x 2  x1

N 1 ( x)  et

x2  x x2  x1

x  x1 N 2 ( x)  x 2  x1

(2.49)

N1(x) et N2(x) sont appelées fonctions de forme et possèdent des propriétés, dont la plus importante a été citée en remarque 1.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 46 -

Chapitre2

Éléments finis à une dimension

Propriété évidente :

1 Ni   0

au noeud i aux autres noeuds (ailleur)

(2.50)

Autres propriétés moins évidentes :

 x1  x  x 2

N1+N2=1 et

(2.51)

x1.N1+x2.N2=x

Remarque 4: Les identités (2.51) confirment par leurs existences, la présence des deux termes, le terme constant " 1 " et, le terme linéaire " x ", de l’expression polynomiale (2.45). Interprétation géométrique: La relation entre u(x), N1(x) et N2(x) est explicité montrée dans la figure 2.17 u(x)=a0+a1.x

u2 u1 u2.N2(x)

u1.N1(x)

1 N1(x)

N2(x) x1

x2

x

Figure 2.17 : Relation entre u(x), N1(x) et N2(x) On déduit de cette relation, une interprétation géométrique, à savoir qu’une ligne droite u(x)=a0+a1.x, peut être considérée comme une combinaison linéaire de lignes standards N1(x) et N2(x) avec des valeurs nodales comme facteurs poids. Généralisation : Une fonction de déplacement v (exacte) qui est définie sur un domaine Ω :[x1,xn] et représenté dans la figure 2.18, est déterminée par éléments finis, par le biais d’une fonction approchée u.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 47 -

Chapitre2

Éléments finis à une dimension

Ve ue(x) V2

V3

Vn

V1

Vn-1 [1]

1 x1

[n-1] x n-1 n xn-1 xn

[2] 2 x2

4 x4

3 x3

Figure 2.18 : Approximation par éléments finis unidimensionnels linéaires On définit d’abord la géométrie des éléments et leurs topologies comme suit : Nœuds :

1, 2, ………………, n-1, n

Coordonnées des nœuds :

x1, x2, ……………., xn-1, xn

Domaine complet :

x1  x  x n

Topologie des éléments :

 1e :

 e2 :

  

 en1 : Déplacements nodaux

V1, V2, ………………, Vn-1, Vn

En prenant comme modèle l’élément fini barre à deux nœuds, les fonctions approchées u e(x) doivent être linéaires en x et leurs fonctions de forme doivent satisfaire les propriétés cités auparavant. L’ensemble est présenté dans le tableau ci-après. Tableau2.5 : Construction de fonction de forme par morceau dans un domaine 1D Elément Fonction approchée ue(x) 1

u1(x)=N1(x).V1+N2(x).V2

Fonctions de forme Propriété évidente

N 1 ( x) 

x2  x x 2  x1

x  x1 N 2 ( x)  x 2  x1

N1(x1)=1 ; N2(x1)=0 N1(x2)=0 ; N2(x2)=1

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 48 -

Chapitre2 2

Éléments finis à une dimension

u2(x)=N1(x).V2+N2(x).V3

…..

……………

n-1

un-1(x)=N1(x).Vn-1+N2(x).Vn

N 1 ( x) 

x3  x x3  x 2

N 2 ( x) 

x  x2 x3  x 2

N1(x2)=1 ; N2(x2)=0 N1(x3)=0 ; N2(x3)=1

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

…………

N 1 ( x) 

xn  x x n  x n 1

N 2 ( x) 

x  x n 1 x n  x n 1

N1(xn-1)=1 ; N2(xn-1)=0 N1(xn)=0 ;

N2(xn)=1

La MEF n’est pas restreinte à l’utilisation d’éléments finis linéaires, et la majorité des logiciels commerciaux disponibles, permettent de choisir entre les éléments finis à fonctions d’interpolation linéaires, quadratiques ou cubiques correspondant à des domaines à une, deux ou trois dimensions. Les éléments finis quadratiques ou cubiques, peuvent aussi représenter des frontières curvilignes ; il suffit que le nombre de nœuds géométriques sur chaque frontière soit compatible avec la forme de la courbe correspondant à la frontière du domaine. Il est donc important d’être capable de choisir le type d’élément qui est le plus approprié au problème étudié et de déterminer les fonctions de forme pour ce type d’élément choisi. Les fonctions de formes du tableau modélise un élément fini barre à deux nœuds et sont un cas particulier des fonctions d’interpolation dite de Lagrange à n points (nœuds) tel que la fonction de déplacement :

0 pour j  i avec N i ( x j )   1 lorsque j  i

n

u ( x)   N i ( x).ui i 1

n

et

N i ( x)   j 1 j i

(x  x j ) ( xi  x j )

(2.52)

(2.53)

On notera que les fonctions N i (x) sont des polynômes d’ordre (n-1). n

A la question : est-ce que

 N ( x)  1 ? , la réponse est oui pour toutes les fonctions de i 1

i

forme de classe C0 avec : N1 ( x1 )  1

N 2 ( x1 )  0

N n ( x1 )  0

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 49 -

Chapitre2

Éléments finis à une dimension N1 ( x 2 )  0

 N1 ( x n )  0

N 2 ( x2 )  1

 N n ( x2 )  0

N n ( x2 )  0

 N n ( xn )  1

Le développement de l’équation (3.11a) donne la relation suivante: Nk 

( x  x1 )( x  x 2 ) [ x  x k ] ( x  x n ) ( x k  x1 )( x k  x 2 ) [ x k  x k ] ( x k  x n )

(2.54)

On ne tient pas compte des termes qui se trouvent à l’intérieur de l’intervalle [….] de l’équation (2.54) pour obtenir la k-ieme fonction de forme, ce qui nous permet de retrouver tous les cas particuliers suivant : 

Pour l’élément fini linéaire à deux nœuds (figure 2.9), les N(s) et les x(s) ayant un indice supérieur à 2 ne doivent pas apparaître et on déduit :

N1 

1

et

1

N2 

( x  x2 ) x1  x 2

(2.55)

( x  x1 ) x 2  x1

ce qui confirme les équations (3.9) u2

u1

u  u ( x)  N . 1  u 2 

avec

x2 2

x1 1 x

Figure 2.19 : Interpolation linéaire et fonctions de forme 

Pour l’élément fini quadratique à trois nœuds figure 2.19(i-b), les N(s) et les x(s) ayant un indice supérieur à 3 ne doivent pas apparaître. N1 

( x  x 2 )( x  x3 ) ( x1  x 2 )( x1  x3 ) N3 

N2 

( x  x1 )( x  x3 ) ( x 2  x1 )( x 2  x3 )

(2.56)

( x  x1 )( x  x 2 ) ( x3  x1 )( x3  x 2 )

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 50 -

Chapitre2

Éléments finis à une dimension

E- Types d’éléments finis à une dimension : Les éléments unidimensionnels - 1D les plus courants, figure 2.20:

a-Linéaire

b-Quadratique -

c-Cubique

Figure 2.20 : Eléments finis 1D 2.8 Matrices caractéristiques par minimisation de l’énergie potentielle A. Energie potentielle totale minimale: Une fois la fonction de déplacement définie, il est possible d’obtenir, toutes les déformations et contraintes dans l’élément fini et de formuler sa matrice de rigidité et celui du vecteur de charges équivalentes (concentrées aux nœuds). Dans le chapitre précédent, la matrice de rigidité élémentaire a été introduite par le biais d’un argument physique direct (méthode directe) et peut être aussi élaborée par la puissante méthode, plus familière aux ingénieurs de génie civil, qu’est le principe des travaux virtuels, qui cependant ne peut être un cadre pour produire des approximations par éléments finis plus générales. La formulation de ces matrices dites "caractéristiques des éléments finis", représente l’étape la plus importante de l’analyse par la MEF, et repose sur le principe fondamental de la mécanique des structures qui est le principe de l’énergie potentielle minimum. C’est Courant (cité au Chapitre 1) qui a développée en 1943 la première application basée sur ce principe pour la solution du problème de torsion de Saint-Venant. On introduit dans un premier temps, les concepts de l’énergie de déformation et de l’énergie potentielle totale et on déduit les équations gouvernantes pour les problèmes d’élasticité. B- Formulation de l’énergie potentielle stationnaire : Les charges externes appliquées à un élément de structure provoque sa déformation et résultent en un travail qui est stocké dans le matériau sous forme d’énergie élastique, appelée énergie de déformation. Pour sa simplicité, prenons l’exemple d’une barre de section constante A, soumise à une charge axiale N, représentée par un ressort linéaire figure 2.21.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 51 -

Chapitre2

Éléments finis à une dimension L

L

y P

(a)

dy' y'

y

PI Longueur non Déformée (Pas de charge)

I

(c)

dy

P

k (b)

Y

x PF

dz

II

dx z

Y

Figure 2.21 : Comportement élastique d’une barre sous sollicitation horizontale (a) structure barre ; (b) élément ressort ; (c) état de contrainte sur un élément de volume

La figure 2.21 (c ), montre aussi un élément de volume sous l’effet de contrainte normale agissant sur ses facettes. L’énergie potentielle totale, comporte deux parties, (a) l’énergie de déformation et (b) le potentiel des charges appliquées. Parce que les forces internes et les forces externes sont les deux conservatives il en est de même pour l’élément de structure. L’énoncé du principe est comme suit: Parmi toutes les configurations possibles d’un système conservatif, celle qui satisfait les équations d’équilibre donne une énergie

potentielle

stationnaire par rapport à des petites variations de déplacements. Si la condition de l’état stationnaire est minimum, l’équilibre est stable. Si le ressort ne dissipe pas d’énergie alors le travail des forces interne (énergie de déformation dans le ressort), dépend seulement de l’allongement L et non pas de son passage par le chemin I ou le chemin II (figure 2.21(b)). De la même manière si la charge externe N possède une valeur et une orientation constante, son travail est égal à N  L indépendamment du chemin choisi pour aller de PI(configuration initiale) à PF (configuration finale) (figure 2.21(b)). L’effort N s’écrit comme suit :

 AE  ' P .L  k y  L 

(2.57)

Et l’énergie emmagasinée d i dans le matériau pour une déformation infinitésimale y / : y/ y/ 1 2 1 (2.58) d i   P.dy /   ky / dy /  ky /  ( ky / ) y / 0 0 2 2 On peut écrire l’équation (2.53) en fonction des contraintes et des déformations normales :

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 52 -

Chapitre2

Éléments finis à une dimension d i 

1 1 1 ( k y / ) y /  ( y .dxdz)   . .dv 2 2 2

Avec

k . y /   y .dxdz : Force élastique

et

y /   .dy

(2.59)

 d i

 Figure 2.22 : Représentation de l’énergie de déformation comme un volume de contrainte D’où l’écriture de l’énergie de déformation di à partir de la figure 2.22, pour l’élément de structure sous sollicitation axiale, est la suivante :  i( e )   d i   V

 . 2

.dV  

avec

V : volume de l’élément

et

E : Module de Young

V

E. 2 .dV 2

(2.60)

C. Energie potentielle totale : L’énergie potentielle totale s’écrit comme suit :  P   i  We

(2.61)

L’application pour notre exemple figure 4.1 donne :

i 

2 1 k. y / 2

et

We   P.y /

La charge en se déplaçant sur une distance y/ produit un travail et en conséquence perd un potentiel de même valeur, justifiant ainsi l’existence du signe négative dans l’expression :

We   P.y / L’énergie potentielle totale

P 

2 1 k . y /  P. y / 2

(2.62)

Peut être considérée comme le travail interne et externe effectué pour un changement de configuration de l’état de référence y /  0 à l’état de déplacement y /  0 .

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 53 -

Chapitre2

Éléments finis à une dimension

Elle est représentée graphiquement pour notre exemple dans la figure 2.23. La position d’équilibre y eq/ obtenue à partir de la valeur stationnaire de P : d P  (kyeq/  P)dy /  0

d’où y eq/  P / k

(2.63)

Energie

We 

2 1 k. y / 2

 P   i  We

Figure 2.23: Représentation graphique des relations d’énergie

y/ Valeur stationnaire

yeq/ 

N k

We   N .y /

2.9 Formulation de la matrice de rigidité et du vecteur charge équivalent A. Déformation et contrainte moyenne dans un E.F unidimensionnel 1-D On utilise les fonctions de formes étudiées au chapitre 3 pour l’écriture du déplacement d’un élément fini aux nœuds i et j soit : u ( e )   N i .u i  N i .ui  N j .u j  N .a

Avec comme Fonctions de Forme : N i 1 

y L

et

Nj 

(2.64)

y L

(2.65)

y est la référence de coordonnée locale de l’élément ayant pour origine le nœud i. N  Ni

Nj

T

y  1 L

y L

T

et

ui  a  u j 

(2.66)

La déformation dans chaque élément peut être calculée selon la relation :

  

du d  ( N .a)  L( N ).a  B.a dy dy

Avec

L : opérateur dit de Laplace dans le cas général

et

B  L(N ) =

1 L

1 L

(2.67)

T

: dérivées des fonctions de formes

(2.68)

Enfin la contrainte dite moyenne dans chaque élément s’écrit comme suit :

     E  E B.a A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(2.69) - 54 -

Chapitre2

Éléments finis à une dimension

B. Matrice de rigidité élémentaire L’énergie de déformation pour un élément fini barre quelconque (e) peut être déduite à partir de l’équation (2.60) :

 i( e) 

1 1  T  .dV   (a T B T E ) B.a.dV  2 V 2 V

(2.70)

L’énergie potentielle s’écrit comme suit :

P 

1 T 1 T T T T a ( B E B).a.dV - a F  a .K .a  a .F  V 2 2

(2.71)

F : le vecteur de charges appliquées aux nœuds

Où et on pose :

  B .E.B.dV T

 K

V

;

la matrice de rigidité élémentaire

(2.72)

La minimisation de l’énergie de déformation par rapport à a s’écrit :

 P 0  a

(e)

K .a

(e)

F

(e)

(2.73)

On retrouve l’écriture familière de la matrice de rigidité élémentaire d’une barre de section uniforme A, sous sollicitation uni-axiale :

K

(e)

0

=

Où



L

1  L

  1    1 L AE 1 .E. L . A.dy = 2 .  1 1   dy 1 L L 10   L

AE  1  1  k = L  1 1   k

k

(2.74)

 k k 

(2.75)

A.E L

C. Vecteur de charges équivalentes aux nœuds En minimisant le travail produit par les charges externes pour un élément fini quelconque (e), e le second terme de l’équation (2.66) résulte en un vecteur charge F , dit de charges ponctuelles (concentrées aux nœuds i et j) pour un problème 1D:

 Pi   (e) (e) ( F .a)  p    a Pj 

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(2.76)

- 55 -

Chapitre2

Éléments finis à une dimension

D’une manière générale, pour des cas de charges de types quelconques, on distingue en plus des forces concentrées ponctuelles p , des forces de volume b (forces de gravité, poids propre) et des forces de traction de surface t . Le travail de ces forces externes s’écrit comme suit : We    u .b.dV   u .t.dA   u . p i t

t

V

t

(2.77)

i

A

D’où l’expression générale du vecteur force équivalente sur chaque élément fini en tenant compte de (2.64) et en minimisant We (2.77) par rapport à a:





   F   N .b.dV  N .t.dA  t

V

t

 .N . p t

i

A

(2.78)

i

On constate que la distribution des charges aux nœuds sur un élément fini quelconque pour l’obtention d’un vecteur de charges équivalentes se fait par le biais des fonctions de forme N. 2.10 Construction des Fonctions de Forme : A. Principe de la méthode Une fois les termes de la fonction de déplacement choisis d’une façon précise à partir du Triangle de Pascal, il convient de construire ensuite les fonctions de forme. Il existe deux procédés possibles : i)

Procédé direct en utilisant les propriétés évidentes des Fonctions de Forme.

ii) Méthode des déterminants La méthode (i), est beaucoup plus rapide pour un calcul manuel, pour des éléments finis simples et prend en compte les propriétés déjà citées auparavant à savoir : 

Pour chaque élément, la fonction de forme est égale à 1 au nœud considéré et 0 ailleurs.



Chaque fonction de forme est un polynôme du même degré que celui de la fonction d’interpolation.

Ce procédé est mieux expliqué par un l’exemple d’élément fini unidimensionnel (1D) quadratique figure 2.24 :

1 Figure 2.24 : Elément fini 1D quadratique

2 L /2 x

3 L /2

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 56 -

Chapitre2

Éléments finis à une dimension

On cherche trois fonctions de forme Ni(x), i=1, 2, 3 et on exprime chacune d’elles comme le produit de deux fonctions : Ni(x)= Fi.Gi

; i=1, 2, 3

(2.79)

Fi est fonction qui nulle aux nœuds autre que le nœud considéré i. Gi est choisie tel que Ni est un polynôme quadratique. Les positions des nœuds : 1: 0 ; Nœud 1 : Et

2: L/2 ;

3: L

F1= (x -L/2).(x -L) est quadratique donc G1=C1 (constante) N1(0)=1 au nœud 1 C1=2/L2 N1(x)= (2/L2)(x -L/2).(x -L)=

Soit

Nœud 2: Et

;car N1(x)=0 aux nœuds 2 et 3.

N1(x) = (x -L/2).(x -L).G1

N2(x) = x.(x -L).G2

2 2 3 x  x 1 L L2

(2.80)

;car N2(x)=0 aux nœuds 1, 3

F2= x.(x -L) est quadratique donc G2=C2 (constante) N2(L/2)=1 au nœud 2 C2=-4/L2 N2(x)= (-4/L2).x.(x -L)=

Soit

Nœud 3: Et

4 2 4 x  x L L2

(2.81)

N3(x) = x.(x –L/2).G3 ;car N3(x)=0 aux nœuds 1, 2

F3= x.(x –L/2) est quadratique donc G3=C3 (constante) N3(L)=1 au nœud 3 C3=2/L2

Enfin

N3(x)= (2/L2).x.(x –L/2)=

2 2 x x  L L2

(2.82)

Les expressions des fonctions de forme (2.80), (2.81) et (2.82), sont les mêmes que celles des équations (2.56) en substituant x1=0, x2=L/2 et x3 =L. On peut rapidement vérifier à chaque fois la propriété évidente ∑N i=1 ; i=1, 2, 3 N1(x) +N2(x) +N3(x)=1 B. La méthode des déterminants Elle est beaucoup plus adaptée au calcul automatique et elle de ce fait la plus évidente pour la programmation sur ordinateur s’agissant des éléments finis d’ordre supérieurs. Cette méthode est plus générale et prend en compte aussi, les propriétés des fonctions de forme. Soit une fonction de déplacement d’ordre n dans un espace 1D : n

u(x)= a1+a2x+a3x2+……+anxn-1 =  N i u i

(2.83)

i 1

Ni- Fonctions de Forme x=xi i=1, 2, ….., n position des nœuds dans l’élément fini A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 57 -

Chapitre2

Éléments finis à une dimension ui=u(x=xi)

i=1, 2, ….., n

On défini un vecteur ligne G tel que : G = < 1, x, x2, ….., xn-1>

(2.84)

On définie en plus n vecteurs lignes Hi tel que : Hi = G (x=xi), i=1, 2, …, n

(2.85)

Soit D le déterminant de la matrice dont les lignes sont H1, H2, ……, Hn prises dans cet ordre respectivement.

1 x1

x12 . . x1n 1

1 x2 D . .

x22 . . x2n 1 . . . .

. . 1 xn

. . . . xn2 . . xnn 1

(2.86)

Enfin soit DiG le déterminant obtenu à partir de D en remplaçant la i eme ligne Hi par le vecteur ligne G.

ieme ligne

1 x1

x12 . . x1n 1

1 x2 . . DiG  1 x

x22 . . x2n 1 . . . . 2 n 1 x . . x

. . 1 xn

. . . . 2 n 1 xn . . xn

(2.87)

Les fonctions de forme se calculent comme suit : N i ( x) 

DiG D

i  1,2,, n

(2.88)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 58 -

Chapitre 3 Éléments finis à deux et trois dimensions

Éléments Finis en Une Dimension

Projet du viaduc ferroviaire, Oued-Tlelat /Sidi-Belabbes / Tlemcen, de 1780 m et de hauteur de 114m

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

Chapitre 3

Éléments finis à deux et trois dimensions

APERÇU On considère dans ce chapitre, les éléments finis à deux et trois dimensions. Une attention particulière est donnée à la dérivation des fonctions de formes des éléments finis triangulaires et quadrilatères en associant la méthode du triangle de Pascal et la méthode des déterminants. On présente l’’importance des formulations par intégrales numériques (méthode de Gauss) dans un système de coordonnées locales normalisées ainsi que la transformation entre le repère local et le repère global. Une application aux problèmes d’élasticité plane de classe C0 est aussi présentée. 3.1 Idée de base Dans ce chapitre on introduit une fonction élément fini pour l’approximation des surfaces à deux dimensions, sur lesquelles on applique des fonctions approchées qui sont généralement définies par des équations différentielles. Pour les problèmes à deux dimensions ‘2D’, bien qu’il existe plusieurs possibilités pour l’approximation (diviser) du domaine d’étude’, en pratique deux formes seulement sont utilisées à savoir le triangle et le quadrilatère figure 3.1. y

(2D)

Domaine (Ω) z

Noeuds x

Eléments (Domaines Ωe)



(b) Figure 3.1 : représentation de domaines pour un problème bidimensionnel- 2D

3.2 Types d’éléments finis 2D et 3D A- Eléments bidimensionnels - 2D 

Eléments triangulaires :

a-Linéaire

b-Quadratique

c-Cubique

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 59 -

Chapitre 3



Éléments finis à deux et trois dimensions

Eléments quadrilatéraux:

a-Linéaire

b-Quadratique

c-Cubique

 Eléments Tridimensionnels - 3D : A-Tétraèdre

b-(Quadratique) c-(Cubique)

B-Hexaèdre

b-(Quadratique) c-(Cubique)

C-Prisme droit (solid element)

b-(Quadratique) c-(Cubique)

Figure 3.2 : Quelques éléments finis classiques les plus utilisés en modélisation numérique Remarque 5 : Le type d’élément dépend : 

de la forme (1-D, 2-D, 3-D, triangulaire, quadrilatère, tétraèdre, etc.…)



du nombre de nœuds (3-nœuds, 4-nœuds, etc.)



du type des variables nodales (degrés de liberté)



du type de fonctions d’interpolation.



On notera aussi que le degré de liberté pour chaque nœud peut varier selon le problème.

B- Triangle de Pascal : Dans le cas d’un choix d’éléments finis d’ordre supérieur où l’ordre de la fonction de déplacement est élevé il y a lieu de prendre garde de l’isotropie géométrique (invariance géométrique). Cela veut dire que l’élément fini modèle doit être indépendant de l’orientation de la base de coordonnées locale, ou globale. Le choix des termes qui composent la fonction d’approximation du champ de déplacement, est plus facilement obtenue par référence au Triangle de Pascal figure 3.3.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 60 -

Chapitre 3

Éléments finis à deux et trois dimensions

a-Élément fini triangulaire : Nombre total des termes Au dessus du niveau

1 x x x x

4

2

3

y y2

xy 2

x y 3

xy 2

x y

x y

2

2

y xy

niveau 1

3

linéaire

niveau 2

6

quadratique

niveau 3 niveau 4

10 cubique 15 quartique

3

3

y4

Figure 3.3 Triangle de Pascal Donc pour un élément fini linéaire, tous les termes, sur et au dessus du niveau 1, sont choisis, alors que pour un élément fini cubique, tous les termes, sur et au dessus du niveau 3 sont nécessaires. On remarque que pour la construction des éléments finis, la distribution des nœuds dans le triangle est dictée par la forme géométrique du Triangle de Pascal. Par exemple, un élément quadratique possède tous les nœuds sur les 3 frontières figure 3.4 (a) ; un élément cubique possède un nœud intérieur figure3.4 (b) ; un élément quartique trois nœuds intérieurs et ainsi de suite. La raison de garder le même nombre de nœuds sur les frontières de l’élément, revient au fait que les éléments seront généralement reliés ensembles le long de leurs côtés et une condition de continuité sur la fonction de déplacement u(x, y) doit être observée. L’illustration de cette approche est faite par l’exemple de l’élément fini triangulaire à trois nœuds, dont la fonction de déplacement u(x,y)= a 1+a2.x+a3.y induit une variation linéaire sur ses frontières, se qui est compatible avec l’existence de deux nœuds par frontière (une seule ligne droite peut être tracée entre deux points). N’importe quel autre élément qui sera attaché à cet élément [1] figure 3.5, aura exactement la même variation le long de la frontière commune, donc la fonction de déplacement sera continue sur la frontière de l’élément.

7

5 8 6 1

4

2 (a)

9

6 10 5

1 3

2 (b)

3

4

Figure 3.4 : Elément fini triangulaire : (a) quadratique à 6 nœuds (b) cubique à nœud intérieur A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 61 -

Chapitre 3

Éléments finis à deux et trois dimensions

3 4 1

[2]

[1]

2

frontière commune

Figure 3.5 : deux éléments finis triangulaires à 3 nœuds ayant une frontière commune L’écriture de la fonction de déplacement en se rapprochant de la frontière commune du côté de l’élément [1] est :

u 1  u1 N1[1]  u 2 N 2[1]  u3 N 3[1]

(3.1)

u 1  u 2 N 2[1]  u3 N 3[1] Car N1[1] le long de la limite commune. Puisque les fonctions de forme sont linéaires, alors u[1] varie linéairement le long de la limite commune ayant comme valeurs, u2 et u3, au nœud 2 et 3 respectivement. Alors que si la frontière est approchée à travers l’élément [2] :

u 2  u 2 N1[ 2]  u 4 N 2[ 2]  u3 N 3[ 2]

(3.2)

u 1  u 2 N 2[ 2]  u3 N 3[ 2] car N2[2]=0 le long de la limite commune. Aussi, la variation des fonctions de forme est linéaire le long de la limite commune, entre les valeurs u2 et u3, au nœud 2 et 3 respectivement. On déduit que u[1] =u[2] le long de la frontière commune, confirmant que la fonction de déplacement u est continue à travers les frontières de l’élément. Remarque 6: Pour les problèmes 2D, un polynôme représentant la fonction de déplacement est de degré n s’il contient un terme de la forme xlym, où l et m sont des entiers non négative et l+m=n. le polynôme est dit "complet" s’il contient toutes les combinaisons de l et m pour lesquelles l+m≤ n. Par exemple un quadratique complet, n=2 ,figure 3.4, possède la forme suivante : u=a1+a2x+a3y+a4x2+a5xy+a6y2

(3.3)

Un polynôme complet de degré n contient (n+1)(n+2)/2 termes (figure 3.7).

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 62 -

Chapitre 3

Éléments finis à deux et trois dimensions

En généralisant au domaine tridimensionnel (3D), les mêmes remarques sont applicables, de telles façons qu’un quadratique complet contient 10 termes, qui comportent un terme constant, des termes linéaires x, y, z, et zx. Un polynôme complet de degré n en 3D contient (n+1)(n+2)(n+3)/6 termes, figure 3.6

z3

2

z z 1

yz2 yz

y

2

zx

y

x

xy

z2x

y2z

2

Y3

zx

xyz

xy2 xy

2

2

x

x3 Figure 3.6 : Pyramide montrant le nombre de termes pour des polynômes complets selon les 3 variables indépendantes x, y, z b- Élément fini quadrilatère: Le Triangle de Pascal est aussi d’une utilité considérable dans la construction des fonctions de forme pour les éléments finis rectangulaires qui sont de deux principaux types : Des exemples d’éléments finis Lagrangiens et les éléments finis "Serendipes" sont représentés dans les figures 3.7 et 3.8 respectivement.

4- noeuds

9- noeuds

16- noeuds

Figure 3.7 : Eléments finis Lagrangiens

4- noeuds

8- noeuds

12- noeuds

Figure 3.8 : Eléments finis ‘Serendipe’ A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 63 -

Chapitre 3

Éléments finis à deux et trois dimensions

Les termes dans l’approximation de la fonction de déplacement u(x,y) pour chacun de ces éléments sont obtenus à partir du Triangle de Pascal comme le montre la figure 3.9 et figure 3.10.

Eléments finis Lagrangiens

1 x x x x

xy

y

xy

xy

x3 y 2 x4 y2

2

2

2 2

xy

x5 y

y

xy

x4 y

4- noeuds

xy 2

3

x5 x6

2

3

4

y

3

y4

x2 y3 x3 y3

9- noeuds 16- noeuds

3

xy 4

y5

x2 y4

xy 5

y6

Figure 3.9 : Triangle de Pascal-Eléments finis Lagrangiens

Eléments finis Serendipes

1 x x2 x x x x6

3

4

5

xy 3

x5 y

xy xy

3 2

y xy

xy x3 y 3

8 noeuds 3

3

2 3

xy x4 y 2

2

2 2

xy xy

y2

xy 2

4

4- noeuds

y

y xy

x2 y 4

12- noeuds 4

4

y5 xy5

y6

Figure 3.10: Triangle de Pascal - Eléments finis ‘Serendipes’ Remarque 7: En général, l’utilisation des Eléments Finis "Serendipes" est plus pratique, du fait que ces derniers présentent des niveaux de précision similaires à ceux des éléments finis Lagrangiens avec moins de nœuds par élément. Les deux types d’éléments produisent une fonction de déplacement u(x,y) qui est continue le long des frontières.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 64 -

Chapitre 3

Éléments finis à deux et trois dimensions

3.3 Construction des fonctions de forme pour L’élément finis ‘T3’: L’élément fini triangulaire à 3 nœuds figure 3.11 est désigné comme T3 et on procède de manière similaire, que pour le cas 1D, par la méthode des déterminants pour déterminer ses fonctions de forme.

(x3, y3) 3

y

1 (x1, y1) 2 (x2, y2)

x

Figure 3.11 : Elément fini 2-D, triangulaire à 3 noeuds

u ( x, y) On pose u ( x, y)    v( x, y) 

et

u ( x, y )  a1  a 2 x  a3 y v( x, y )  b1  b2 x  b3 y

(3.4)

Le champ de déplacement s’écrit, en tenant compte des Fonctions de Forme, comme suit : 3

u ( x , y )  u 1 . N 1 ( x , y )  u 2 . N 2 ( x , y )  u 3 . N 3 ( x, y )   u i

(3.5)

i 1

Avec

u  u1   1  v1 

u  u2   2  v 2 

u  u3   3  v3 

Les Fonctions de Forme Ni s’écrivent comme suit selon la méthode des déterminants : 1 x 1 N 1 ( x, y )  1 x2 D 1 x3

y 1 ( x2 y3  x3 y 2 )  x .( y 2  y3 )  y .( x3  x2 ) y2  2A y3

1 x1 1 N 2 ( x, y )  1 x D 1 x3

y1 1 ( x3 y1  x1 y3 )  x .( y3  y1 )  y .( x1  x3 ) y  2A y3

(3.6)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 65 -

Chapitre 3

Éléments finis à deux et trois dimensions

1 x1 1 N 3 ( x, y )  1 x2 D 1 x

y1 1 ( x1 y 2  x2 y1 )  x .( y1  y 2 )  y .( x2  x1 ) y2  2A y 1 x1 D  1 x2

Avec

y1 y2  2 A

1 x3

Où

(3.7)

y3

A- Aire du Triangle

3.4 Transformation géométrique: A- Relation entre les variables généralisées et les variables nodales: Un élément de référence simple peut être transformé par une transformation T e, en un élément de forme plus complexe figure 3.12.



(0,1)



y (xk) k (1-   =0)

T

 1 2

e

i

 1

ΩE

(xi)



(xj)j

R

Ω

(0,0) (1,0)

 1





1 2

x

Figure 3.12 : Représentation d’un exemple de transformation d’un élément fini de référence vers l’élément parent On prend la notation :    

et

x x

y

La transformation Te définie les coordonnées xe de chaque point de l’élément réel à partir des coordonnées  du point correspondant de l’élément de référence. Elle dépend de la forme et de la position de l’élément. Soit

Te :



    Dans un espace tridimensionnel :       



xe= xe(  ) x   x y  z  

(3.8)

On notera aussi qu’un même élément de référence, un triangle à 3 nœuds, peut se transformer en plusieurs éléments réels de même type par des transformations différentes, figure 3.13.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 66 -

Chapitre 3

Éléments finis à deux et trois dimensions

y



4

2 T1

(0,1)

[1]

T

2

T3

5 1

(0,0)

(1,0)

[3]

[2]

3



x

Figure 3.13: exemple d’un élément fini de référence, représentant plusieurs éléments de même type B- Règles et propriétés de Te Te doit générer des éléments réels et doit être choisie ayant les propriétés suivantes : Doit être bijective en tout point de  situé sur l’élément de référence ou sur sa



frontière. Donc à tout point de l’élément de référence correspond un point et un seul de l’élément réel. 

Les nœuds géométriques de l’élément de référence, définie par les nœuds géométriques de cette frontière, correspondent à la portion de frontière définie par les nœuds correspondant.

Te :





xe= xe(  , x1, x2, … , xn )

Où

x1, x2, …, xn

Et

n- le nombre de nœuds

(3.9)

sont les coordonnées des nœuds de l’élément réel

On utilisera des transformations T linéaires par rapport aux coordonnées xn de l’élément réel. T:



x(  )= < T (  )> xn

(3.10)

De plus les fonctions de transformation sont choisies identiques pour les trois coordonnées. x(  )= < T (  )> x

y(  )= < T (  )> y

z(  )= < T (  )>z (3.11)

avec :

x=< x1 , x2, …, xn>T

y =< y1 , y2, …, yn>T

z=< z1 , z2, …, zn>T

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 67 -

Chapitre 3

Éléments finis à deux et trois dimensions

Applications : 1- exemple de l’élément fini triangulaire à trois nœuds figure 3.11 Position des nœuds : i :< xi , yi > Soit

j : < xj , yj >

k : < xk , yk >

< T (  )>=

(3.12)

On a alors les relations suivantes :

x(  , )= T1(  , )xi + T2(  , )xj + T3(  ,  )xk  xi    x( , )  1         x j  x   k

(3.13)

 yi    y ( , )  1         y j  y   k

Vérification des propriétés :  Les nœuds géométriques dans le system de coordonnées local

1: (0,0)

2: (1,0)

3: (0,1)

se transforment dans la base globale :

xi,

xj,

xk

L’application pour le nœud i comme exemple donne :

x(  =0,  =0) = < 1

0

 xi    0 >  x j   xi x   k

(3.14)

Il en sera ainsi pour les autres nœuds j et k. La correspondance des frontières (limites) : Chaque frontière Rde l’élément de référence ΩR se transforme en une frontière E correspondante de ΩE. On choisi pour l’exemple la frontière qui passe par les nœuds 2: (1,0) et 3: (0,1), et dont l’équation est 1      0 . Elle doit se transformer en la frontière qui passe par les nœuds j: xj et k: xk. En substituant  1   dans les relations (3.8) on obtient : A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 68 -

Chapitre 3

Éléments finis à deux et trois dimensions  xi    x  0  1     x j    x j  (1   ) x k x   k

(3.15)  yi    y  0  1     y j    y i  (1   ) y k y   k

Les relations (3.41) sont donc linéaires et ne dépendent que de xj et xk qui sont situés sur la frontière  E . La transformation est bijective : ( , )   TBijective   ( x, y ) dx.dy  J . d .d Avec J - Jacobien  0

Expression d’une matrice J non singulière.  x   J  x  

y     x j  xi  y   x k  xi  

y j  yi  y k  y i 

(3.16)

det(J)= ( xj-xi )(yk-yi)-(xk-xi)(yj-yi)

k det(J)= 2xA i

A j Figure 3.14: Aire de l’élément fini T3

Ce déterminant est égal à 2 fois l’aire du triangle donc, ne s’annule que si les sommets sont alignés. Les angles du triangle doivent toujours rester inférieur à 180° pour que les det(J)0.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 69 -

Chapitre 3

Éléments finis à deux et trois dimensions 2- exemple d’utilisation des fonctions de forme dans la transformation figure 3.15.

On peut monter que les relations (3.13), transforment l’élément fini triangulaire aux sommets (xi,yi), i=1,2,3, en un élément standard dans le plan  .

1 1 1 x  x1 ( (   ))  x2 ( (1   ))  x3 ( (1   )) 2 2 2 (3.17)

1 1 1 y  y1 ( (   ))  y2 ( (1   ))  y3 ( (1   )) 2 2 2

 3

y (x2,y2)

+1 Transformation (3.17)

-1

+1



(x3,y3) x

1

-1

2

(x1,y1) Figure 3.15: Transformation d’un élément fini triangulaire à 3 nœuds en un élément standard La relation entre (x,y) et , η comme définie (3.17) est linéaire tel que des lignes droites dans le plan , η , restes droites dans le plan xy. Considérant par exemple, le coté

  1,  1    1 du triangle standard. On obtient, en utilisant (3.11) 1 1 x  x1 ( (  1)  x3 ( (1   )) 2 2 ( x  x1 ) x3  x1  3  2 2

(3.18a)

d’une manière similaire y 

y 3  y1 y 3  y1  2 2

(3.18b)

on peut bien constater que ce sont des équations paramétriques de lignes droites passant par les points (x1,y1) lorsque   1 et (x3,y3) lorsque   1 . Il y va de même pour le côté du triangle standard   1 A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 70 -

Chapitre 3

Éléments finis à deux et trois dimensions 1 1 x  x1 ( (  1)  x 2 ( (1   )) 2 2 ( x  x1 ) x 2  x1  2  2 2

.

y 

y 2  y1 y 2  y1  2 2

(3.19a)

(3.19b)

Ce sont aussi des équations paramétriques de lignes droites passant par les points (x 1,y1) lorsque   1 et (x2,y2) lorsque   1 . On remarquera que pour cette transformation on a fait appelle aux fonctions de forme de l’élément fini triangulaire à 3 nœuds, à savoir :

1 N1 ( , )   (   ) 2 N 2 ( , ) 

1 (1   ) 2

N 2 ( , ) 

1 (1   ) 2

(3.20)

et on peut écrire (3.45) sous la forme suivante : 3

x   xi N i ( , ) et i 1

3

y  yi N i ( , )

(3.21)

i 1

3- Elément fini quadrilatère à 4 nœuds. De même que les éléments finis triangulaires à côtés droits on peut utiliser des éléments quadrilatères pour le maillage d’un domaine. On peut facilement montrer, de la même manière que précédemment que la transformation d’un élément fini standard (carré) figure 3.16, aux nœuds suivants : 1 :(-1,-1) ;

2 :(1,-1) ;

3 :(1,1) ;

4 :(-1,1)

en un quadrilatère aux sommets pris à des positions (xi, yi), i=1, 2, 3, 4, s’écrit comme suit selon (3.46): 4

x   xi N i ( , ) i 1

4

et

y   y i N i ( , )

(3.22)

i 1

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 71 -

Chapitre 3

Éléments finis à deux et trois dimensions

 4

y 3

+1

3

4 Transformation (3.11)

+1

-1

2



x -1 1

2

1

Figure 3.16 : Transformation d’un élément fini standard vers un quadrilatère Les fonctions de forme N i ( , )

N1 ( , ) 

N 3 ( , ) 

pour i  1,2,3,4 sont :

1 (  1)(  1) 4

1 N 2 ( , )   (  1)(  1) 4

1 (  1)(  1) 4

1 N 4 ( , )   (  1)(  1) 4

(3.23)

3.5 Système de coordonnées de référence : A- Système de coordonnées locales normalisé En modélisation par la MEF, on a besoin le plus souvent d’utiliser différentes bases de référence. On a besoin d’un système de référence global pour représenter la position des nœuds, orientation de chaque élément fini, et pour l’application des conditions aux limites et charges. Plus important, la solution, tel que les déplacements nodaux sont représentés par rapport au système d’axes de coordonnées globales comme on l’a déjà vu en chapitre 2. D’autre part, on a besoin d’un système de coordonnées dit local à cause de la position de son origine à l’intérieur de l’élément fini. Il présente l’avantage de faciliter la construction de la géométrie ainsi que l’évaluation des calculs des intégrales. Cet avantage devient plus apparent lorsque ces intégrales contiennent des produits de fonctions de forme tel que :



(e)

N1 ( x, y ).N 2 ( x).dxdy ;

dN1 ( x) dN 2 ( x) . .dx dx dx (e)



(3.24)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 72 -

Chapitre 3

Éléments finis à deux et trois dimensions

Pour des éléments finis linéaires les termes de ces intégrales sont de simples polynômes et les intégrales sont faciles à évaluer, ce qui n’est pas le cas pour des éléments fini d’ordre supérieur. Aussi est-il plus facile de construire les fonctions de forme dans un système local dont l’origine sera au centre de l’élément et de normaliser les coordonnées leurs permettant une variation entre -1 et +1 en vue d’une intégration facile selon la méthode de Gauss. On obtient alors des éléments finis plus simples tel que le triangle, le carré, le tétraèdre et le parallélépipède. On appellera ces éléments, "éléments de référence" figure 3.17 pris dans cette base cartésienne et orthonormée dite "normalisée".  Elément Triangulaire 2-D











Linéaire



Quadratique

Cubique

 Elément Rectangulaire 2-D

 (-1, 1)



 (1, 1)

(-1, 1)



(1, -1)

(-1, -1) Linéaire

(-1, 1)

(1, 1)

(1, 1)





(-1, -1) Quadratique

(1, -1)

(1, -1)

(-1, -1) Cubique

Figure 3.17: Eléments finis de Référence (Parents) de type ‘Serendipe’

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 73 -

Chapitre 3

Éléments finis à deux et trois dimensions

B- Système de coordonnées naturel Il représente une base locale dont les propriétés intrinsèques, qui tiennent compte des dimensions de l’élément et une écriture de paramètres invariants, sont proches de celles des fonctions de forme étudiées auparavant.  Elément Fini Triangulaire- 2D : Un system de coordonnées naturel pour un élément fini triangulaire (système barycentrique) est défini sur la figure 3.18 : (x3, y3) 3 A2

A1 L2

y

1 (x1, y1)

L1

L3 A3

2 (x2, y2)

x

Figure 3.18 : Représentation du system de coordonnées barycentriques Li= Ai/A avec

(3.25)

A=A1+ A2 + A3

Li- i=1, 2, 3

représente les coordonnées barycentriques qui définissent une ligne de

coordonnées selon la figure 3.19.

3 L1=0.25 L1=0.75

1

2 L1=0.50

L1=0.00

L1=1.00

Figure 3.19 : Représentation de coordonnée naturelle L1Lignes parallèles au coté 2-3 La relation entre les deux systèmes de coordonnées, cartésien et naturelle : ( L1 , L2 , L3 )

 ( x, y ) ;

 x  Li . xi   y  Li . yi

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(3.25)

- 74 -

Chapitre 3

Éléments finis à deux et trois dimensions

La relation sous forme matricielle s’écrit comme suit :  x  L1 .x1  L2 .x 2  L3 .x3   y  L1 . y1  L2 . y 2  L3 . y 3  1 L  L  L 1 2 3 



 x1 y  1  1

x3   L1   x      y 3 . L2    y  1   L3  1 

x2 y2 1

(3.26)

La résolution selon les Li donne :  L1  ( y 2  y 3 ) ( x 3  x 2 ) ( x 2 y 3  x 3 y 2 )   x    1    ( y 3  y1 ) ( x1  x3 ) ( x3 y1  x1 y 3 ) . y   L2     L  2 A  ( y  y ) ( x  x ) ( x y  x y )  1  2 2 1 1 2 2 1    3  1

(3.27)

A étant l’aire du triangle. Au regard des équations (3.27) on constate que les coordonnées Li possèdent la propriété C0 des Fonctions de Forme Ni en plus on les relations suivantes :

N1  L1

N 2  L2

N 3  L3

(3.28)

Remarque : Il sera facile d’intégrer des polynômes en terme de coordonnées naturelles (barycentriques) on peut utiliser les relations suivantes :

 L .L l





1

2

.dl 

 L .L 



1

2

! ! .l (     1)!

.L3 . dA 

A

!  ! ! .2 A (      2)!

(3.29)

(3.30)

L’équation (3.29) est utilisée pour évaluer une intégrale, qui est fonction de la longueur l, le long de la frontière de l’élément. La distance entre deux nœuds est noté ‘l’. La relation (3.30) est utilisée pour évaluer des intégrales de surface. 3.6 Cas de problèmes d’éléments finis bidimensionnels de classe C0 : A- Rappel et formulation de la théorie de base Pour le cas de problèmes 2D de classe C0, il est important de faire appel à la théorie de base de l’élasticité, qui prend en compte les quantités physiques suivantes :

 Déplacement en un point, représenté par un vecteur colonne. u   u v wT



(3.31)

Matrice de contraintes dans l’élément (matrice carrée et symétrique)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 75 -

Chapitre 3

Éléments finis à deux et trois dimensions  X  XY    YX  Y  ZX  ZY



 XZ   YZ   Z 

(3.32)

Matrice de déformations dans l’élément (matrice carrée et symétrique) X    YX  ZX

 XY Y  ZY

 XZ   YZ   Z 

(3.33)

Dans cette section particulière on fixera notre attention sur le problème d’élasticité plane. Il existe deux méthodes de base par le biais desquelles une réduction dans la théorie tridimensionnelle de l’élasticité linéaire peut être envisagée. En général, les contraintes et les déformations dans une structure se composent de six composants, figure 3.20.



 x,

et

 y,

z ,  xy ,  yz ,  xz

pour les contraintes

 x ,  y ,  z ,  xy ,  yz ,  xz

pour les déformations y yx xy

yz

Y

x

zx z

zx

xz

X Z Figure 3.20 : Volume élémentaire de contraintes pour l’équilibre interne Les déformations en un point sont fonction des déplacements u, v, et w tel que :

X 

u x

 XY 

v u  x y

Y 

v y

 YZ 

w v  y z

(3.34)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 76 -

Chapitre 3

Éléments finis à deux et trois dimensions

Z 

w z

 ZX 

u w  z x

Pour un cas de structure linéairement élastique et dont le matériau est homogène et isotrope, les relations contraintes déformation sont données par la loi de Hooke :  X   X    X 0          Y   Y   Y0    Z   Z    Z 0     C.   0  C.       XY   XY   XY 0   YZ   YZ   YZ 0         XZ   XZ   XZ 0 

(3.35)

Avec C la matrice caractéristique des propriétés élastique du matériau définie par : 1    1   C  . E  0  0   0

 1  0 0 0

  1 0 0 0

0 0 0  0 0 0  0 0 0   2(1   ) 0 0  0 2(1   ) 0   0 0 2(1   )

(3.36)

ε0- vecteur de déformations initiales E–

module de Young

 - coefficient de poisson. Et on note que

G

E 2(1)

module de cisaillement

Dans les cas d’effet de température pour un matériaux isotrope, le vecteur de déformations initiales est donné par :

1  1    1   0   .T   0 0   0

(3.37)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 77 -

Chapitre 3

Éléments finis à deux et trois dimensions

Où

α est le coefficient de dilatation thermique

et

T est la variation de température

Parfois on a besoin des expressions de contraintes en fonction des déformations, ceci peut être obtenu en inversant les équations (3.58). On obtient alors en ajoutant les déformations dues à la température les relations suivantes :

 X   X  1     1  Y  Y     Z    Z  E .T 1 .     D.(    0 )  D.    1  2  XY XY     0  YZ   YZ  0       0  XZ   XZ 

(3.38)

Où la matrice D est donnée par :

  1     1       1   E 0 0 D . 0 (1   )(1  2 )   0 0 0   0 0  0

0

0

0 0 1  2 2

0 0 0

0

1  2 2

0

0

 0  0   0   0   1  2  2  0

(3.39)

Sous des conditions spécifiques, l’état des contraintes et de déformations peuvent être simplifiées. B- Contraintes planes et déformations planes De ce fait une analyse de structure tridimensionnelle 3-D peut être réduite à une analyse bidimensionnelle 2-D avec deux types de distribution de contraintes possibles, contrainte plane et déformation plane.  Contrainte plane Ce cas est applicable pour des structures dont la dimension est petite selon une des directions des axes de coordonnées.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 78 -

Chapitre 3

Éléments finis à deux et trois dimensions

Donc l’analyse de plaque mince chargée dans le plan de la plaque peut être considérée selon le modèle de contraintes planes.

 z=

yz

=

 zx = 0

(  x = 0)

(3.40)

Cas de structure plane mince avec épaisseur constante et le chargement pris dans le plan de la structure (plan x-y) figure 3.21 .

y

y

y

z

z x

Figure 3.21 : Exemple de structures planes à épaisseur constante On peut aussi inclure dans cette catégorie les exemples de problèmes de plaques de poutres à caisson et les plaques formant la partie de l’âme d’un profilé PRS par exemple – figure 3.22.

(a)

(b) Figure 3.22 : (a) Plaque perforée (b) Plaque formant l’âme d’un PRS

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 79 -

Chapitre 3

Éléments finis à deux et trois dimensions

 Déformations planes Ce cas s’applique pour des structures longues dont la géométrie et le chargement ne varie pas d’une façon signifiante selon la direction longitudinale. On considère que les composants de déplacement < u , v , w > T on la forme particulière suivante : u=f(x,y)

v=g(x,y)

w=o

ceci implique que les seules composantes de déformations non nulles sont 

(3.41) x

,

xy ,



yy

tel que les déformations apparaissent dans le plan z=constante ; d’où la terminologie ‘déformation plane’.

z=

yz

=

zx

= 0 ( z = 0 )

(3.42)

L’analyse de barrages, murs de soutènement par exemple figures 3.23 (a), (b), rentrent dans ce cas précis, il faut toutefois prévoir des sections d’études loin des extrémités de la structure pour avoir une bonne approximation.

Section sous Déformation plane

(a) Barrage

(b) Mur de soutènement

Figure 3.23 :(a),(b) Exemples d’ouvrages sous déformations planes

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 80 -

Chapitre 3

Éléments finis à deux et trois dimensions

C- Relation contrainte- déformation :  Etat de contrainte plane

 X   X       1  0       E. .T E   Y  . 1 0 . Y     2 1    1   0 0 1           XY  2   XY          

Soit

 =D  + 

Avec

 0=-D. 0

1      1      0    

(3.43)

0

 Etat de déformation plane plane: Pour l’état de déformation plane on prends pour D l’expréssion suivante :

  1   E D .  1  (1   )(1  2 )  0  0 

 0  0  1  2   2 

(3.44)

D- Formulation de la matrice de rigidité  Relations déformation –déplacement Pour de faibles déformations on a les relations suivantes à partir des équations (3.34):     X   x     Y    0     XY     y

soit

  L.. u

 0    u  .  y  v    x 

(3.45)

(3.46)

A partir de cette relation on reconnaît que les déformations (donc les contraintes aussi) sont d’un degré inférieur par rapport aux déplacements, si ces derniers sont représentés par une interpolation polynomiale. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 81 -

Chapitre 3

Éléments finis à deux et trois dimensions

 Matrice de rigidité Les déplacements (u,v) dans un élément plan sont interpolés à partir des déplacements nodaux (ui , vi ) en utilisant les fonctions de formes Ni comme suit :

u   N u  1 v   0

0 N1

0  Nn N2  0

N2 0

 u1  v   1 u  0  2 .v  N n   2     u n  v   n

(3.47)

U = N . a. N – matrice des fonctions de forme u – vecteur de déplacement a - vecteur de déplacement nodal à partir de la relation déformation-déplacement (3.46) on a le vecteur de déformation :

  L.u  L.N . a soit

  B. a

avec

B L. N

(3.48) matrice déformation-déplacement

On obtient alors et de la même manière que pour le cas du 1-D (2.72), la nouvelle écriture de la matrice de rigidité pour le cas général soit :

 K 

 B .D. B. dv T

(3.49)

v

Remarque2 : similairement à la remarque1, tenant compte du fait qu’on travaille toujours par rapport à un élément fini de référence dans sa base locale (ξ, η) dans un espace 2-D ou (ξ, η, ) dans un espace 3-D, on procède par le Jacobien J à une transformation pour un assemblage dans une base globale. Si l’élément fini possède une épaisseur variable ti définie aux nœuds i, l’épaisseur t en un point arbitraire de l’élément s’écrit en utilisant les fonctions de forme, soit : t   N i .t i

(3.50)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 82 -

Chapitre 3

Éléments finis à deux et trois dimensions

La matrice de rigidité s’écrit alors sous la forme :

   K   B .D. B.t dA   T

1 1

1



1

B .D. B. t.J .d .d T

(3.51)

A

J est donc fonction de ξ, η,  selon le problème et de ce fait B inclura ces coordonnées de base locale et l’intégrale (3.74) ne sera mieux évaluer que numériquement par la méthode de Gauss quadrature. E- Equations d’équilibre et conditions aux limites Selon la théorie d’élasticité, les contraintes dans la structure doivent satisfaire les équations d’équilibre.

 x  xy   fx  0 x y (3.52)

 xy x



 y y

 fy  0

Où fx et fy sont les forces de gravité (forces internes) par unité de volume tel que le poids propre. Il est à la fin important d’appliquer les conditions aux limites essentielles imposée sur la variable u (u ,v) ,par exemple :

u= u ,

v= v

(3.53)

u et v sont des valeurs fixes sur une partie de la frontière 1 du domaine D et sur les forces de traction (contrainte aux limites ) sur une partie de la frontière 1 (figure 3.24):

tx = t x ,

ty = t y

(3.54)

On notera que par la MEF tous les types de charges (uniformément repartie, poids propre, charge concentrée), sont transformé en charges concentrées localisés aux nœuds.

Γ1 Γ2

Figure 3.24: Conditions aux limites A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 83 -

Chapitre 3

Éléments finis à deux et trois dimensions

3.7 Intégration numérique – Méthode de Gauss: La méthode se base sur l’évaluation d’une fonction à un certain nombre de points, en multipliant ensuite le résultat par un coefficient de pondération appelle aussi facteur poids et en sommant enfin l’ensemble des résultats. En application à la matrice de rigidité chaque coefficient k ij

de K peut être considéré

comme une fonction f susceptible d’être intégrée sur la longueur de l’élément, sa surface ou son volume. Cependant afin de satisfaire les conditions de la méthode de Gauss on prendra une autre fonction g, tel que g  g ( ) dans le 1-D, g  g ( , ) dans le 2-D, et

g  g ( ,,  ) pour les problèmes 3-D. A. Intégration dans un espace unidimensionnel 1-D Soit f(x) la fonction figure 3.25 et I son intégrale définie à évaluer pour x1  x  x2 : y f(x)

x1

x2

x

Figure 3.25: Fonction à intégrer dans un espace 1-D tel que x2

I   f ( x ). dx x1

(3.55)

qui doit être transformée sous la forme suivante : 1

I   g ( ). d

(3.56)

1

C’est sous cette nouvelle forme que l’intégrale numérique de Gauss peut être appliquée. La transformation f ( x)  g ( ) , inclus le Jacobien J 

dx qui pour un élément fini 1-D à d

deux nœuds s’écrit comme : J

x 2  x1 2

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

(3.57)

- 84 -

Chapitre 3

Éléments finis à deux et trois dimensions

La manière la plus simple et la plus directe le calcul de I (3.79) est d’évaluer g au milieu de l’intervalle figure 3.26 et de le multiplier par la longueur de celui-ci pour obtenir le résultat :

I  2 g1

(3.58)

Ce résultat sera vrai si seulement la courbe est une ligne droite. g 1

g1 ξ

-1

+1 (a)

I g1+g2

a

1 3 g2

g1 a

b  0 .6 g2

b

b 2

0 (b)

g3

g1

a

1

-1

I5/9.g1+8/9.g2+5/9.g3

ξ +1

1

2

ξ

3

-1

+1 (c)

Figure 3.26 : Intégration de la fonction g=g(ξ) dans un espace 1-D selon 1, 2 et 3 points de Gauss (a),(b),et (c) respectivement En généralisant on obtient l’écriture suivante : 1

n

I   g. d  Wi .g ( i )  W1 .g1  W2 .g 2    Wn .g n 1

(3.59)

i 1

Où les Wi représentent les coefficients de pondération. Les ξi représentant les points d’intégration (points de Gauss). Wi et ξi sont choisi d’une façon optimum à 2n conditions et g(ξ) une fonction polynomiale d’ordre m  (2n  1) . La figure 3.26 (b) et (c) présente des exemples pour n=2 et n=3 On retient aussi les propriétés suivantes des points d’intégration de Gauss : 

 1   i  1 les abscisses  i sont symétriques par rapport à   0 ,



Wi  0,



Pour deux points symétriques les valeurs de Wi sont les mêmes,



Le tableau 3.1 donne les valeurs  i et Wi pour chaque valeur de n.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 85 -

Chapitre 3

Éléments finis à deux et trois dimensions

Tableau 3.1 : Valeurs  i et Wi -Intégration de Gauss Nombre de points Valeurs de de Gauss-n ξi 2 ±0.577350=  1 3 4 5

3

Facteur de pondération Wi 1.0

±0.774597=  0.6 0.0

0.555556 = 5 9 0.888889 = 8 9

±0.861136 ±0.339981

0.347855 0.652145

±0.906180 0.0 ±0.538469

0.236927 0.568889 0.478629

±0.932469 ±0.661209 ±0.238619

0.171324 0.360762 0.467914

6

Exemple : On considère un polynôme de 3ieme degré g  a1  a2  a3 2  a4 3 Où les ai sont des constantes. L’intégrale exacte entre -1 et +1 est : 1 2 I   g.d  2a1  a3 1 3

(3.60)

La règle d’intégration d’ordre 1 donne une première approximation I 1  2a 0 . En prenant 2 points d’intégration avec 1   1 3 et  2   1 3 on obtient :

I 2  1.0(a1  a 2

2

3

2

 1   1   1   1   1   a3    a 4    1.0(a1  a 2    a3    a3   3  3  3  3  3  3

1

2  2a1  a3 3

3

(3.61)

On peut aussi vérifier que 2 est le nombre de points nécessaire pour une intégration exacte à par tir de la relation m  3  2n  1  n  2 B. Intégration dans un espace à deux dimensions (2-D) Les règles d’intégration de Gauss dans un espace multidimensionnel, appelé règles de produit de Gauss, sont obtenues par application successive des règles unidimensionnelles (1-D)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 86 -

Chapitre 3

Éléments finis à deux et trois dimensions

Elément fini quadrilatère : Dans un espace 2-D l’élément fini rectangulaire représente un cas particulier de l’élément quadrilatéral, on considère la fonction g  ,  et on choisit à intégrer par rapport à ξ en premier et ensuite par rapport à η soit : 1  n  g (  ,  ). d  . d   Wi .g ( , ).d 1  1 1 i 1 

I 

1

1

(3.62)

n n  n n  W j Wi .g ( i , j )   Wi .W j .g ( i , j ) j 1  i 1  i 1 j 1



Pour un point d’intégration, g  g1 est déterminée au point ξ = η = 0 tel que I  4g1 dans un espace 2-D.



Pour quatre points d’intégration tel que représenté dans la figure 3.27 (a), WiWj=1 à chaque point de Gauss et l’intégrale (3.62) donne : I  g1  g 2  g 3  g 4

(3.63)

Où gi représentent les valeurs numériques de g au i-eme point de Gauss. 

Pour neuf points de Gauss, figure 3.27 (b) l’équation (3.62) donne :

I

ξ= 

1 3

25 40 64 ( g1  g 3  g 7  g 9 )  ( g 2  g 4  g 6  g 8 )  g9 81 81 81

η ξ=

η

1 3

ξ= -

η= 4

9

3

(a)

η= 

1 3

η=

0.6

6

3

ξ

0.6

0.6

1 3

2

1

ξ=

(3.64)

2

5

1

4

ξ

8

7

η= -

0.6

(b)

Figure 3.26 : Position des points d’intégration d’une fonction g=g(ξ,η) (a) 2 points de Gauss ; (b) 4 points de Gauss Remarque : Les points d’intégration nécessaire pour une intégration exacte d’un polynôme dans un espace à deux dimensions s’illustre bien par une représentation sur le triangle de Pascal figure 3.28. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 87 -

Chapitre 3

Éléments finis à deux et trois dimensions

A considérer qu’il s’agit ici de la fonction à intégrer g et non pas d’une fonction de déplacement. En prenant comme degré du polynôme complet la somme l+m la somme, on constate par exemple, que tous les termes quartique exception faite pour  4 et  4 peuvent être intégrés exactement par 2x2 points de Gauss.

Points de Gauss pour I exacte

Ordre ξl.ηm Constant

(l=m=0)

Linéaire

(l+m=1)

1 ξ ξ

Quadratique (l+m=2) Cubique

(l+m=3)

Quartique

(l+m=4)

I 

1 1



1 1

ξ ξ4

η

ξ η

3

η

2

ξη

2

ξ3 η

ξ 2 η2

ξ3 η2

 . .d .d l

1 point

ξη

2

m

η3

2

ξ η3 η4

2x2

ξ2 η3 ξ3 η3

Figure 3.28 : Précision de la méthode de Gauss en 2-D Tous les termes du polynôme dans le Triangle de Pascal figure 3.28 jusqu’aux limites de forme V peuvent être exactement intégrés par la méthode de Gauss d’ordre 1 et 2 respectivement. Elément Triangulaire: Les intégrales faisant intervenir des Eléments Finis Triangulaires sont mieux exprimées en terme de coordonnées barycentriques (3.25) ou de surfaces Li, introduites dans ce chapitre. La formulation suivante donne une écriture équivalente type intégrale de Gauss : n

I   f ( L1 , L2 , L3 ). dA   Wi . f ( L1( i ) , L(2i ) , L(3i ) ) A



pour n=1 (élément triangulaire linéaire):

W1  1; 

(3.65)

i 1

L1(1)  L(21)  L(31) 

1 3

(3.66)

pour n=2 (élément triangulaire quadratique) :

1 W1  ; 3

L1(1)  L(21) 

1 , L(31)  0 2

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 88 -

Chapitre 3

Éléments finis à deux et trois dimensions

1 W2  ; 3

L1( 2)  0, L(22)  L(32) 

1 W3  ; 3

L1(3)  L(33) 

1 2

(3.67)

1 , L(23)  0 2

Les positions des points d’intégration de Gauss pour ces deux exemples sont montrées sur la figure 3.29. 1 2 1 3 (a) n=1

(b) n=3

Figure 3.29: Position des points d’intégration dans l’élément fini triangulaire (a) Linéaire (b) quadratique C. Intégration dans un espace à trois dimensions (3-D) L’intégrale de Gauss prend la forme suivante : 1 1 1



1 1 1

g ( , ,  ).d d d  Wi .W j .Wk g ( i , j ,  k ) i

j

(3.68)

k

D. Transformation du domaine d’intégration (base globale-base locale)  Dans un espace 1-D : On cherche à évaluer l’intégrale I suivant les bornes d’intégration a et b, figure 3.30 b

I   f ( x). dx ; a

ba

ba ba x .     2   2 

(3.69)

ba dx   . d  2  1  b  a)   b  a)  n I   f x( ).  . d     Wi . f x( i ) 1  2   2  i 1

Pour chaque valeur de  i Substituer

(3.70)

 x(  i )

x( i ) 

f ( x)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 89 -

Chapitre 3

Éléments finis à deux et trois dimensions

g ( )  f x( )

η

f(x)

y Transformation (3.69)

-1

ξ

+1

x

b

a Figure 3.30 : Transformation dans une base 1-D normalisée  Dans un espace 2-D

 f ( x, y).dxdy    A

b

h2 ( x )

a

h1 ( x )



( b  a et d  c )

f ( x, y ).dy .dx

y

(3.71)

y  h2 ( x) A

y  h2 ( x)

x

a

b

Figure 3.31 : Domaine d’intégration 2-D 

Cas particulier : A est un rectangle et l’intégrale possède des bornes d’intégration constantes, figure 3.32. I   f ( x, y ).dxdy  

b

a





d

f ( x, y ).dy .dx

c

(3.72)

A

η

y

+1

d

A +1 -1

ξ

c

-1 Transformation

a

b

x

Figure 3.32: Transformation d’un domaine rectangulaire en un carré standard

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 90 -

Chapitre 3

Éléments finis à deux et trois dimensions

 b  a   d  c  1 1 I  .    f x( ), y ( ).d d  2   2  1 1

(3.73)

 ba ab  x   2 .   2         d c cd  y .     2   2  

(3.74)

L’intégrale de Gauss s’écrit comme suit :  b  a  d  c  n m I   . Wi .W j f  x( i ), y ( j )  2  2  i 1 j 1



(3.75)

Cas général : A est un quadrilatère et l’intégrale I possède des bornes d’intégration variables, figure 3.32.

y

A

x Figure 3.32 : Domaine aux bornes d’intégration variables L’intégration est effectuée selon le changement de variable selon la transformation :

x  x( , )   ( x, y) y  y( , )

TRANSFORMATION

 ( , )

Soit non singulière : det( J )  0



1

 f x( , ), y( , ). det( J )

n

m

f ( x, y ).dxdy  

1 1

(3.76)

1

A



 

  Wi .W j f x( i , j ), y ( i , j ) . det J ( i , j



(3.77)

i 1 j 1

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 91 -

Chapitre 3

Éléments finis à deux et trois dimensions  x   J   x  

tel que

y     y   

(3.78)

le Jacobien : det( J )0

Avec

a/ Dans un espace 3-D :

I 



(3.79)

f ( x, y, z ).dxdydz

x y z



1 1 1

1

n

I  i 1

m

r

 W W j 1 k 1

i

j

  1

1

f ( x, y, z ) det( J ).d d d



(3.80)

 

Wk . f x(  i , j ,  k ), y ( i , j ,  k ), z ( i , j ,  k . det J ( i , j ,  k )

 x    x J     x   

avec

y  y  y 

z    z    z    



(3.81)

(3.82)

3.8 Intégration de la matrice de rigidité : On considère l’intégration numérique qui permet de générer les éléments de la matrice de rigidité élémentaire K

(e )

d’un élément fini quadrilatéral plan.

Le processus étant le même pour un élément fini possédant plus de quatre nœuds, seulement les dimensions des matrices de déformation B et de rigidité élémentaire K

(e )

sont plus larges.

Aussi quel que soit l’ordre de la matrice, chaque cœfficient k i j est traité comme la fonction g dans l’équation (3.85). Pour n’importe quel élément fini plan la matrice Jacobéenne J est d’un ordre 2x2.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 92 -

Chapitre 3

Éléments finis à deux et trois dimensions

L’algorithme pour les séquences d’intégration numérique, peut être établi comme suit : (e )

Initialiser la matrice K , choisir un nom de variable indicée, KE par exemple. Faire une boucle (DO-Loop) sur les points d’intégration, sens de  (pour i=1 à ni ) Introduire (Bloc Data) les points de Gauss    i et leurs cœfficient de pondération Wi Faire une boucle sur les points d’intégration, sens de  (pour j=1 à nj ) Introduire (bloc data) les points de Gauss    j et leurs cœfficient de pondération Wj Faire appel (Call) au sous-programme (fonction) qui calcul la matrice B, l’épaisseur t, et le Jacobien det( J ) , au niveau de chaque points ( i , j ) T

Calculer le produit de matrices B . D. B .t. det( J ).Wi W j et l’ajouter à la matrice Ke Fin de boucle sur l’indice j Fin de boucle sur l’indice i

Etapes pour l’écriture d’un sous-programme de fonctions de forme pour établir: B , t, det ( J ) aux points de Gauss Pour chaque valeur de point de gauss de l’algorithme précédent : 

Calculer les fonctions de forme et leurs dérivées par rapport à ξ et η



Calculer l’épaisseur t 



Calculer la matrice Jacobienne J, son déterminant et son inverse.



Calculer la matrice de déformation B (pour l’élément fini plan selon les formulations données

N t

i i

à partir des valeurs nodales des épaisseurs

dans ce chapitre)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 93 -

Chapitre 4 Application à des problèmes éléments finis ‘ 2D’ & ‘3D’

Projet de la Grande Mosquée d’Alger (en cours de réalisation): Minaret de hauteur de 270 m composé de 37 étages & Coupole en hémisphère de diamètre de 50 mètres, culminant à une hauteur de 70 mètres

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 4.1 Applications au problème ‘1D’ de Barre 1- On reprend l’exemple de la barre figure 4.1 sous l’effet de charge axiale P à son extrémité, en plus elle est soumise à une température uniforme T :. On cherche le déplacement yD produit par P et T en utilisant l’énergie potentielle totale minimum du système. L

y

P

Figure 4.1 :Barre uniforme sous une charge axiale, et plus un effet de température.

T(°C)

avec

 y  E. y   y 0

ou

où

 y  yD L

 y 0   E T

et

 y  E.( y   0 )

(4.1)

L’équation de l’énergie potentielle s’écrit : P  

L

0

 1  yD 2 yD  E. A. y D2 ( E T ). A.dy  y D .P   y D .E. A. T  y D .P    E L 2L  2  L  

(4.2)

On obtient le déplacement à l’extrémité y D à partir de l’équation d P dyD  0 : yD 

PL   .T .L AE

(4.3)

La contrainte  y évaluer à partir de (4.1) donne : P  P    .T    E T   A  E. A 

 y  E

(4.4)

Qui est le résultat escompté.

2- On reprend l’exemple de la barre mais cette fois, ci on applique la force axiale P à la position L/3, figure 4.2 et on cherche le vecteur de charges équivalent ? L 2L/3

L/3 1

P

2

y

T(°C) Figure 4.2 : Barre soumise à charge axiale située au 1/3 de la travée A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 94 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Le vecteur charges équivalentes F  N ( y  L 3) .P  f 0 (effet de T ) T

(4.5)

L

f 0   B .( E .T ) A.dy T

avec

(4.6)

0

et N T et B T sont pris à partir de (2.63) et (2.65) respectivement. L’effet de la températureT peut être considéré comme une variation, telle que :

2 P 3  1 F    E. A. T   P 3  1

D’où

(4.7)

3- Méthode des déterminants de l’élément fini 1-D quadratique figure 2.27 Soit la fonction de déplacement pour un élément fini unidimensionnel quadratique: u(x)= a1+a2x+a3x2

(4.8)

G = < 1, x, x2>

(4.9)

Hi = G (x=xi i=1, 2, 3) ; H2=

H1= ;

(4.10)

H3=

1 0 0 2 L3 L L D 1  2 4 4 2 1 L L

D2G

1 0  1 x

x x2 2 L3 3L2 L L L D1G  1   x  x2 2 4 4 4 2 1 L L2 1

;

0 x 2  L2 x  Lx 2

D3G

2

1 L L

1 0 0 2 L2 L L  1 L   x  x2 2 4 4 2 1 x x2

Les fonctions de forme s’expriment comme :

N 1 ( x) 

D1G 3 2  1 x  2 x2 ; D L L

N 2 ( x) 

D2 G 4 4  x  2 x2 D L L

(4.11) N 3 ( x) 

D3G 1 2   x  2 x2 D L L

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 95 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 On peut aussi constater que les résultats sont similaires à ceux obtenus, par le premier procédé. 4. On reprend l’exemple l’élément fini barre quadratique à 3 nœuds. On cherche à déterminer la nouvelle forme d’écriture des fonctions de forme par rapport aux coordonnées naturelles par la relation de transformation pour un cas plus général figure 4.3. x

2

1

3

1

x1

2

3



x2

1

1

x3 L (b) (a) Figure 4.3: Elément fini barre quadratique à 3 nœuds (a) Elément fini modèle 1-D (b) Elément fini normalisé Utilisons la méthode des déterminants :



 2

G 1 

H1  G(  1)  1  1 1 ;

H 2  G(  0)  1 0 0

H 3  G (  1)  1 1 1

1 1

1

D 1

0

0 2 ;

1

1

1

1  2 D1G  1 0 1 1

0     ; 1 2

D3G

D2 G

1 1 1 0 1



1 1 1  1   2  2(1   2 ) ; 1 1 1

1 0   2 

2

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 96 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0



1

1

N1 ( ) 

D1G 1  (   2 ) D 2

N 2 ( ) 

D2 G  1 2 D

N 3 ( ) 

D3G 1  (   2 ) D 2





1

(4.12)

Figure 4.4: représentation graphique des fonctions de forme La détermination de la position (coordonnée) et du déplacement en n’importe quel point de l’élément dans la base globale, se déduit de:  x1    x  N  x2  x   3

et

Nous avons trois coordonnées nodales xi

u1    u  N u 2  u   3

(4.13)

, i  1,2,3 et trois degrés de liberté représentant des

déplacements axiaux. Les mêmes fonctions de forme dans les deux expressions (4.13) sont utilisées ce qui a permis de donner à cet élément l’appellation d’élément isoparamétrique.

dx qui d

Le passage de la base locale vers la base globale se fait par le biais de l’expression

représente le Jacobien J ; un facteur multiplicateur qui décrit la correspondance de la longueur réelle par rapport à celle de référence dans la base locale :  x1  dx d 1   J  ( N ) x 2   (1  2 )  2 d  2 x   3

 x1  1   (1  2 )  x 2  2 x   3

(4.14)

J est utilisé pour la détermination de la déformation axiale.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 97 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 5 - Une barre de section droite uniforme A modélisée par Le même élément fini quadratique de l’exemple 3 ci-dessus (figure 2.27), qui présente un nœud interne supplémentaire, et on cherche une expression pour sa matrice de rigidité élémentaire Ke ? Les fonctions de forme N dans la base locale étant connues, la déformation axiale dans l’élément est :

u1  du  d    x    N .u 2   B.a dx  dx    u 3 

d d d  . dx dx d

avec

(4.15)

Soit J le Jacobien donné par :

 x1  dx d N . x2   1 (1  2 )  2 J  d d 2 x   3

B

(4.16)

 du  1

 x   .  B.a  d  J

D’où on peut écrire

Avec

 x1  1   (1  2 ) . x 2  2 x   3

1 d 1 1 . ( N )  . (1  2 )  2 J d J 2

1 (1  2 ) 2

(4.17)

La matrice de rigidité élémentaire s’écrit comme suit :

K

(e)

l

1

  B .E.B. Adx   B . AE.B.J .d 0

T

T

1

(4.18)

Remarque1 : le Jacobien (2.84) devient une constante J= L /2 si et seulement si le nœud 2 se trouve au milieu de l’élément sachant que x2-x1=L/2 et x3-x2=L/2. Dans le cas général J est fonction de  et de ce fait B inclura cette variable au numérateur ainsi qu’au dénominateur au niveau de chaque terme, et l’intégrale (4.35) ne peut être évalué sous sa forme compacte. On fait appelle alors à la forme d’intégration la plus pratique pour la MEF à savoir l’intégration numérique par la méthode de Gauss, présentée au niveau du chapitre 3.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 98 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 4.2 Applications au problème de classe C1 élément fini type Hermitien. 1- Fonctions de forme pour un le cas de la flexion d’une poutre Les éléments finis que nous avons considérés, se caractérisent tous par le fait qu’à chaque nœud correspond un seul degré de liberté et ne nécessite qu’une seule condition au limite essentielle, d’où l’appellation de classe C0. Cependant pour les problèmes à une dimension d’ordre quatre , tel que la flexion d’une poutre de longueur L, de moment d’inertie I et de module de Young E, figure 4.5, dont l’équation gouvernante découle de la théorie d’Euler –Bernoulli, on peut avoir jusqu’à deux conditions au limite essentielles imposées à chaque nœud i et j, d’où l’appellation de classe C1. P1 Qx

P2

d2/dx2[EI(x)d2/dx2]=Q(x) i

j A.N

E,I

L

0 ≤x≤l

x

Figure 4.5 : Cas de poutre en flexion- équation différentielle d’ordre quatre Modélisation par Elément fini 1D, Hermitien, à deux nœuds i et j Hermite Des éléments finis à un degré de liberté sont inutiles pour résoudre ce type de problème. On considère alors comme fonction de déplacement, un polynôme cubique à quatre paramètres, comme un ordre le plus réduit possible. u( x)  a1  a2 x  a3 x 2  a4 x 3

(4.19)

Les fonctions de forme qui seront utilisées dans l’écriture de la fonction de déplacement sont appelées, fonctions de forme Hermite et l’élément fini correspondant est dit Hermitien. La méthode des déterminants peut être appliquée pour la détermination de ces fonctions de forme si on définit les vecteurs suivants : G 1 x

H i  G ( x  xi )

x2

x3

dH i  0 1 2 x 3x 2 dx

i  1,2

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil"

- 99 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

On définit aussi :

 D  det  H 1 

dH1 dx

H2

dH 2  dx 

(4.20)

avec DiG comme défini au chapitre 2, par l’équation (2.87) , et les fonctions de forme sont exprimées par : Ni 

DiG D

i  1,2,3,4

1 0 0 1 D 1 L

0 0 L2

0 0  L4 3 L

0 1 2 L 3L2

D1G

D2 G

1 x x2 x3 0 1 0 0   L4  3L2 x 2  2 Lx 3 ; 2 3 1 L L L 0 1 2 L 3L2

1 0 0 x  1 L

0 x2 L2

0 x3  L4 x  2 L3 x 2  L2 x 3 3 L

0 1 2 L 3L2

D3G

1 0 0 1  1 x

0 0 x2

0 0  3L2 x 2  2 Lx 3 3 x

0 1 2 L 3L2

D4G

1 0 0 0 1 0  1 L L2 0

x

x2

0 0  L2 x 3  L3 x 2 3 L x3

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 100 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Les fonctions de forme sont :

N 1( x) 

D1G 3x 2 2 x 3  1 2  3 D L L

D2G 2x 2 x 3 N 2 ( x)   x  2 D L L (4.21)

D 3x 2 2 x 3 N 3 ( x)  3G  2  3 D L L N 4 ( x) 

D4G x 2 x3   2 D L L Au nœud x=0

Fonctions de forme

N1(x)

N2(x)

x=0

x=L

L

1

1

N3(x)

1

1

N4(x)

Ni

dNi/dx

Au nœud x=L Ni

dNi/dx

1

0

0

0

0

1

0

0

0

0

1

0

0

0

0

1

Figure 4.6: Représentation graphique de fonctions de forme à approximation cubique et vérification des propriétés évidentes en x=0 et en x=L.

On peut aussi vérifier que : N1+N3 =1

(4.22)

N2+N3.L +N4 = x Ceci veut dire que u(x) est susceptible capable de représenter un mouvement de corps rigide qui est une caractéristique du principe de l’état complet que doit respecter toute fonction d’interpolation. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 101 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 2- Matrice de rigidité élémentaire : En procédant de manière similaire, on va confirmer les résultats obtenus par la méthode directe pour le cas d’une poutre uniforme (sans déformation de cisaillement transversale) figure 2.12. Dans la formulation d’un élément poutre on fait intervenir le moment M et la courbure k au lieu de la contrainte et de la déformation.

M  EI Z k

k

d 2v dx 2

(4.23a)

v  N .a

k  B.a

(4.23b)

Les fonctions de forme N sont données par les relations (4.21). a  v1  Z 1

Avec

v2  Z 2

T

degrés de liberté nodaux

On déduit : B

d2 6 12 x (N )   2  3 2 dx L L

4 6x  L L2



6 12 x  3 L2 L



2 6x  L L2

(4.24)

avec E et IZ constants, la matrice de rigidité élémentaire est obtenue par :

K

(e)

6 L  12 6 L   12  2 L T  6 L 2 L2  EI Z  6 L 4 L   B .EI Z .B.dx  3 0 L  12  6 L 12  6 L    2  6 L 4 L2   6L 2L

(4.25)

3- Répartition de charge équivalente aux nœuds : On cherche à déterminer les charges équivalentes aux nœuds, produites par une charge uniformément répartie Q, figure 4.6.

Figure 4.7 : Elément poutre uniforme avec une

charge Q uniformément répartie

x

Q

1 Z1

L

V1

2 Z2 V2

On utilise à cet effet la deuxième intégrale de l’équation (2.75), avec t   Q et dA  dx

et

F   N . Q .dx L

0

T

(4.26)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 102 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

 x2 x3    QL  1  3  2  2 3   2  L L            QL2  2x 2 x3   2   x   L L  L  12   F   Q  .dx    0 3   x2    QL x  3 2 2 3    L   L  2       x3 x2   QL2        L2 L   12 

(4.27)

L Q

QL/2

1

2 QL/2 QL2/12

QL2/12

Figure 4.8 : Charges équivalentes nodales pour une charge Q uniformément répartie appliquée entres les nœuds Le signe négatif indique que la charge est dans le sens contraire aux conventions positives des degrés de liberté, figure 4.7. 4.3 Applications au problème ‘2D’ de Classe C0 E.F. triangle linéaire à 3 nœuds ( T3) C’est l’élément le plus simple figure 4.9, pour les problèmes plans (2-D), reconnu aussi sous le nom d’élément triangle linéaire (T3). Il est aussi appelé, élément fini à déformation constante.

v 33

u3

v 2

u 2 Y

v

2

6 degrés de liberté

1

1

u 1

X Figure 4.9: Elément triangulaire à 3 Nœuds A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 103 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 La matrice de rigidité et de masse est obtenue en appliquant l’équation déduite précédemment à partir de l’énergie potentielle minimum. 1 -Fonctions de déplacement et fonctions de forme Les déplacements u et v sont considérés des fonctions linéaires dans l’élément :

 u  a1  a 2 x  a3 y   v  a 4  a5 x  a 6 y ai ( i =1,2,….,6)

(4.28)

sont des constantes et les déformations correspondantes sont :

 y a6

 x  a2

 xy a3 a5

(4.29)

Donc constantes à travers l’élément ; d’où le nom de triangle à déformation constante. Les déplacements doivent satisfaire aux six (6) équations suivantes :

 u1  a1  a 2 x1  a3 y1 u  a a x a y 1 2 2 3 2  2  u 3  a1  a 2 x3  a3 y 3   v1  a 4  a5 x1  a 6 y1  v 2  a 4  a5 x 2  a6 y 2   v3  a 4  a5 x3  a 6 y 3

(4.30)

En résolvant ce système d’équations, on peut trouver les coefficients a1, a2 , …,a6 en fonction des déplacements nodaux et des coordonnées. En remplaçant ces coefficients dans (4.30) et en réorganisant les termes, on obtient :

u   N 1    v   0

0

N2

0

N3

N1

0

N2

0

 u1  v   1 0   u 2    N 3   v 2   u3     v3 

(4.31)

où les fonctions de forme linéaires en x et y sont : N1 

1 ( x2 y3  x3 y 2 )  ( y2  y3 ) x  ( x3  x2 ) y 2A

N2 

1 ( x3 y1  x1 y3 )  ( y3  y1 ) x  ( x1  x3 ) y 2A

(4.32)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 104 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 N3 

1 ( x1 y2  x2 y1 )  ( y1  y 2 ) x  ( x2  x1 ) y 2A

1 x1 1 A  det 1 x 2 2 1 x3

et

y1 y2

(4.33)

y3

A- représentant l’aire du triangle

2. Matrice de déformations B : En utilisant la relation (4.30), les résultats (4.35) et (4.36) on obtient :

 X   y 23 1      Y   B. a   0 2 A    x32  XY 

 y 23 B   0  x32

où

0 x32

y 31 0

0 x13

y12 0

y 23

x13

y 31

x 21

0 x32

y 31 0

0 x13

y12 0

y 23

x13

y 31

x 21

 xi j  xi  x j   yi j  yi  y j

et

 u1  v  0  1  u  x 21 . 2  v y12   2   u3     v3 

0 x 21 . y12 

(4.34)

(4.35)

(i, j =1, 2, 3)

(4.36)

On peut remarquer que les déformations sont constantes dans l’élément d’où l’appellation élément Triangle à Déformations Constantes (TDC). Pour la dérivation de la matrice de déformation, on considère un matériau isotopique. Les cas de contraintes planes et de déformation plane peuvent être combinés pour donner la matrice d’élasticité, qui exprime la relation contraintes déformations.

 C1 D  C1 C2  0

C1 C2 C1 0

0  0  C12 

(4.37)

avec C1  E

(1   2 )

et

C2 

contrainte plane

et

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 105 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 C1 

E (1  ) (1  )(1 2 )

 (1  )

C2 

et

déformation plane

et pour les deux cas, C12  C1 (1 C2 ) / 2

E- désigne le module d’élasticité de Young 3- Matrice de rigidité élémentaire Ke La matrice de rigidité est formulée comme suit :

K

e

T   B . D. B dv = eP . A ( B D B )

T

(4.38)

v

où

e P – épaisseur de l’élément

à noter que K est de dimension (6x6) et que sont calcul relève d’un produit de matrices seulement. K – peut être alors écrite à partir de simples produits de matrices comme suit :

1 D. B  A

 C1 C C  1 2  0

C1 C2 C1 0

0   y23 0   0 C12   x32

0 x23

y31 0

0 x13

y12 0

y23

x13

y31

x21

0 x21  y12 

(4.39)

Soit la matrice de rigidité sous la forme explicite suivante (4.40): C1 . y 23 2  C12 .x32

C1 .C 2 . y 23 .x32

2 C1 .x32

 C12 . y 23 .x32

2  C12 . y 23

C1 . y 23 . y 31

C1 .C 2 . y 31 .x32

C1 .x132

 C12 . y 23 .x13

 C12 .x132

e  C12 .x32 .x13 K= P 4. A C1 .C 2 . y 23 .x13  C12 . y 31 .x32 C1 . y 23 . y12  C12 .x32 .x 21 C1 .C 2 . y 23 .x 21  C12 . y12 .x32

SYMMETRIQUE

(4.40) 2 13

C1 .x32 .x13

C1 .C 2 . y 31 .x13

C1 .x

 C1 . y 23 . y 31

 C12 . y 31 .x13

2  C12 . y 31

C1 .C 2 . y12 .x32

C1 . y 31 . y12

C1 .C 2 . y12 .x13

C1 . y122

 C12 . y 23 .x 21

 C12 .x13 .x 21

 C12 . y 31 .x 21

2  C12 .x 21

C1 .x32 .x 21

C1 .C 2 . y 31 .x 21

C1 .x13 .x 21

C1 .C 2 . y12 .x 21

C1 .221

 C12 . y 23 . y12

 C12 . y12 .x13

 C12 . y 31 . y12

 C12 . y12 .x 21

 C12 . y12

Remarque 1: l’écriture de K pour le T3 est simple ne faisant pas appel à l’intégrale numérique, et son implémentation est donc facile.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 106 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

4- Formulation dans une base locale (, ) Les expressions des fonctions de forme dans les relations (4.32) et dans leurs dérivées sont très lourdes et donnent moins d’informations sur le comportement de l’élément. η=0

ξ=0 3

η= b

ξ= a η=1 (a,b) 2

ξ=1

1 Figure 4.10: Représentation de l’élément T3 en coordonnées naturelles En introduisant les coordonnées naturelles (  ,  ) sur le triangle figure 4.10 , les fonctions de forme se présents comme:

N1(,)

N 2(,)

tel que

N1 + N2 + N3 = 1

et

N i   

N3 (,)1 

(4.41)

1au noeud i 0 autres noeuds

A noter que les fonctions de forme varient linéairement à l’intérieur de l’élément. La figure 4.11 schématise le graphe de la fonction de forme N1 , N2 et N3 présentant des formes similaires.

 0

 1

Figure 4.11 : Fonction de forme N1 pour EDC (T3)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 107 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Nous avons deux systèmes de coordonnées pour cet élément, le système de coordonnées global (x, y) et le système de coordonnées local (ξ, η), liés via les relations suivantes :

 x  N1 x1  N 2 x2  N 3 x3   y  N1 y1  N 2 y 2  N 3 y3 soit

 x  x13  x23  x3   y  y13  y 23  y3

avec

xij = xi – xj

et

yij = yi – xj

(4.42)

comme défini auparavant

Les déplacements u , v dans l’élément peuvent être considérés comme fonction de (x,y) ou (ξ, η) comme:

    

  x y       y   x    

u  u 

    

u x u y

   J  

 u   x   u   y   

(4.43)

J - matrice Jacobienne

x J   13  x23

J

1



 y13  y 23 

1  y 23 2 A  x23

 y13  x13 

avec

dét ( J ) = 2 x13 y23 – x23 y13 = 2A

et

A- Aire du triangle

(4.44)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 108 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 à partir de (4.41),(4.42, (4.43), (4.44) on en déduit :

 u      y13    .  x13    u     

 u   x    1  y 23     u  2 A  x 23    y 



(4.45)

 y13  u1  u3  .  x13  u 2  u3 

1  y 23 . 2 A  x23

(4.46)

De la même manière on a :

 v   x    1  y 23     v  2 A  x 23    y 

v  v   y13   1 3  .  x13   v 2  v3 

(4.47)

Par l’utilisation de (4.41) et (4.42) et les relations suivantes :

  L. u  L. u. a  B. a On obtient

 y 23 1  B . 0 2A   x32

0 x32

y 31 0

0 x13

y12 0

y 23

x13

y 31

x 21

0  x 21  y 21 

(4.48)

Ce qui représente la même expression que celle déduite en (4.35) 5- Relation de transformation - matrice de transformation locale –globale T Dans la formulation qui vient d’être établie, la matrice de rigidité est traitée dans une base de coordonnée générale. Cependant pour compléter l’analyse, il est nécessaire de faire appel à la méthode de transformation de base introduite au chapitre 2. Si les axes dans la base locale ne sont pas parallèles aux axes globaux de la structure toute entière figure 4.12.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 109 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

y

f '3y

y'

v

3

f '3x

u 2 f '2x

x'



1

f '2y

f '1x

f '1y

x Figure 4.12: Orientation de l’élément fini triangulaire à 3 noeuds Pour une correspondance des différentes matrices et vecteurs (déplacements, force, et , rigidité) d’une base locale à une autre globale, on utilise les relations suivantes :

u i  T .u i ;

f T . f ;

/

avec

C  S   0 T   0  0   0

/

S C 0 0 0 0

0 0 C S 0 0

0 0 S C 0 0

0 0 0 0 C S

0 0  0  0 S  C 

K

(e)

 T .K .T T

/

(4.49)

(4.50) C=cos et S=sin

Après résolution et obtention des différents déplacements aux nœuds, on peut déduire les déformations et les contraintes dans l’élément fini comme:

  B. a

et

  D .B . a

(4.51)

6- Traitement du problème de compatibilité inter-éléments On se propose de vérifier par un exemple, la continuité de la fonction de déplacement le long de la frontière commune 2-3 sur un maillage à éléments finis triangulaires à 3 nœuds, figure 4.13.

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 110 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

y 4

(0,-1) 3

(+1,+1)

[2] [1] 1 (0,0)

2 (+1,0)

x

Figure 4.13 Discrétisation d’un domaine rectangulaire, par 2 éléments finis triangles à 3 nœuds ayant une frontière commune On détermine dans un premier temps, les fonctions de forme appropriées pour chaque élément fini, et on écrit ensuite le développement de la fonction de déplacement en tenant compte de celles-ci. Coordonnées des nœuds :

1 : (0,0)

2 : (1,0)

3 : (0,1)

4 : (1,1)

Connectivité : Elément

Nœuds(1,2,3)

[1]

1, 2,3

[2]

2, 4,3

Fonction de déplacement :

u(x,y)= a1+a2 x+a3 y

Méthode des déterminants : G=< 1, x, y> ;

1 x1 D  1 x2 1 x3

y1 1 0 0 y2  1 1 0  1 y3

1 0 0

Elément [1] :

1 x y N 1 ( x, y )  1 1 0  1  x  y ; 1 0 1

1 0 0 N 2 ( x, y )  1 x y  x ; 1 0 1

1 0 0 N 3 ( x, y )  1 1 0  y 1 x

y

La fonction de déplacement dans l’élément fini [1] s’écrit comme suit : u [1] ( x, y)  u1 (1  x  y)  u 2 x  u3 y

(4.52)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 111 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Elément [2] :

1 x

y

1 1 0 N 2 ( x, y )  1 x y  1  x  y ; 1 0 1

N 1 ( x, y )  1 1 1  1  y ; 1 0 1

1 1 0 N 3 ( x, y )  1 1 1  1  x 1 x

y

La fonction de déplacement dans l’élément fini [1] s’écrit comme: u [ 2] ( x, y)  u1 (1  y)  u 4 (1  x  y)  u3 (1  x)

L’équation de la frontière 2-3, commune aux deux éléments finis est 1  x  y  0 et les déplacements le long de celle-ci s’écrivent comme suit : u [1] ( x, y)  u 2 x  u3 (1  x) u [ 2] ( x, y)  u1 (1  (1  x))  u3 (1  x)  u [1] ( x, y)

(4.53)

Ceci prouve que u(x,y) est continue le long de cette frontière. 7- Détermination de la matrice de rigidité élémentaire et des contraintes dans la structure formée d’un élément fini T3, figure 4.14 On prend pour ce cas des conditions de contraintes planes. E=2.1x106daN/cm2, =0.25, eP=1cm (épaisseur). Les déplacements aux nœuds ont été déterminés soit :

u1  0.0 cm ; v1  0.006 cm u 2  0.0030cm ; v2  0.0 cm u 3  0.00 cm ; v3  0.006 cm

y 3

(0,1)

2 (2,0)

x

1 (0,-1) Figure 4.14 : Structure formée d’un élément fini triangulaire à 3 noeuds

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 112 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 On calcule en premier les xij et yij (4.39) : y 23  y 2  y 3  0  1  1

x32  x3  x 2  0  2  2

y 31  y 3  y1  1  (1)  2

x13  x1  x3  0  0  0

y12  y1  y 2  1  0  1

x21  x2  x1  2  0  2

D’où la matrice B : (4.34), (4.47)

 1 0 2 0 1 0  1  B . 0  2 0 0 0 2  ; 2(2)  2  1 0 2 2  1

avec A à partir de (4.33) :

1 0 1 1 A 1 2 0 2 2 1 0 1

La matrice d’élasticité D pour un cas de contraintes planes (4.37) s’écrit:

0.25 0   1 2.1x106  D . 0.25 1 0  2  1  (0.25)  0 0 0.375 Les expressions de B et de D sont substituées dans (4.38) de la matrice de rigidité Ke :   1 0  2  0  2  1   1 0.25 0  1 0 2 0 1 0  6   2 0 1 (2)2.1x10 1  e  0 . 0 2 0 0 0 2  K  4(0.9375) . 2 0 2 . 0.25 1  2(2)   0  2  1 0 2 2  1 0 0.375   1 0 2   2  1   0

En effectuant le produit de ces matrices on obtient : 1.25  2  1.5  0.5 0.25   2.5  1.25 4.375  1  0.75  0.25  3.625      2  1 4 0  2 1 K  0.28x106   daN cm 2  1 . 5  0 . 75 0 1 . 5 1 . 5  0 . 75    0.5  0.25  2 1.5 2.5  1.25     0.25  3.625 1  0.75  1.25 4.375 

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 113 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

Les contraintes planes sont fonctions des déplacements et sont calculées selon (4.51) :  0.0  0.006   x  1 0 . 25 0  1 0 2 0  1 0      0 . 003 1   2.1x106     0 2 0 0 0 .0.25 1 0 . 2 .  y     0.0  2(2)   0.9375  0   2  1 0 2 2  1 0 0.375   xy   0.0    0.006

Les contraintes sont ainsi égales à:  x   3360       y    840  daN cm 2    2520   xy  

On peut aussi calculer, les contraintes principales 1, 2 et l’orientation P :

1 

2 

 x  y 2

 x  y 2

1 2

 x  y   2 

    xy2   max 

 x  y   2 

    xy2   min 

2

2

(4.54)

 2 xy

      y   x

 P  . tan 1 

Soit : 3360  840  3360  840  2 1      (2520)  4917.44 daN 2 cm 2 2   2

3360  840  3360  840  2     (2520)  717.44 daN 2 cm 2 2   2

2 

1 2

 2(2520)    31.7  3360  840

 P  . tan 1 

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 114 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 8- Répartition des forces de volume (poids propre) Soit l’élément fini T3, figure 4.15 dont on cherche à répartir la charge sur les 3 noeuds. On fait appel à l’écriture générale du vecteur force : f   N .b. dV   N t. dA T

T

V



(4.55)

A

Le poids propre W s’écrit: -masse volumique du matériau

W= -g () b3y 3

b3x

b2x 2

y b1x

1 W

b1y

x

b 2y

Figure 4.15: Répartition du poids propre sur l’élément fini -T3

bx   0  b   b y    .g   N1 0  N T f b   N .b. dV  e P .  2 0 V A N3   0

(4.56)

0 0      N1   (  g ) N1    0  0  0 .   e p .  .dA N 2    . g  (  g ) N 2  A   0 0    (  g ) N 3  N 3 

(4.57)

L’intégration se fera simple si on utilise les notions de coordonnées barycentriques, introduite au chapitre 3. Les termes dans l’intégrale

f y i  (  g )  N i dA

i  1, 2, 3

(4.58)

A

sont déterminés selon la relation (3.53) comme: pour

=1

et =0, =0

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 115 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 et par une application numérique de la relation (3.30) on obtient : f y i  (  g )

1! 0! 0! A .2 A  (  g ) (1  0  0  2) ! 3

i=1, 2, 3

La distribution de la force au nœud i due au poids propre W, s’écrit comme:

eP . A  0    3   g 

f bi 

(4.59)

Et le vecteur global de répartition de charge sera alors :

 f b1x   0  f    g   b1 y     f b 2 x  e P . A  0  fb     3   g   f b2 y   f b3 x   0        g   f 3 y 

(4.60)

9- Répartition d’une charge uniformément répartie On considère l’exemple d’un domaine discrétisé par des éléments type T3 (TDC) figure 4.16 et on s’intéresse à la répartition d’une charge uniformément répartie q(daN/m2) agissant entre les nœuds 2 et 3 le long de la frontière de l’élément [1]. y 3 q

b [1]

1

2

a

x

Figure 4.16 : Elément T3 sous effet d’une charge répartie q On prend le second terme du vecteur force généralisé (4.55) soit :

f t   N .t .dA T

A

Le vecteur charge s’écrit sous la forme suivante :

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 116 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

t x  q  t    t y  0

(4.61)

La matrice des fonctions de forme NT est évaluée pour x=a :

 N1  0  N T N  2  0  N3   0

     N2  0   N 3  0 N1 0

(4.62)

Le vecteur général des charges équivalentes s’écrit sous la forme intégrale suivante, à évaluer au point x=a.:

ft 

eP

0



b

0

 N1 0  N 2  0 N3   0

0 N 1  0  q  . .dy dz N 2  0  0  N 3 

(4.63)

Pour une épaisseur constante eP, sachant que q est appliquée le long de la frontière 2-3 de l’élément [1], le vecteur charge f t s’écrit :

 f t1 x  f   t1 y   f t 2 x  ft    eP  ft2y   f t3x     f t 3 y 



b

0

 0   0     N 2 .q   . dy  0   N 3 .q     0 

(4.64)

Encore une fois on utilise les propriétés des coordonnées naturelles et on applique la relation (3.52) pour l’évaluation des intégrales: b

f ti x  ( eP .q).  N i .dy 0

avec On obtient

=1

i  2, 3

et =0 f ti x  (e P .q).

q. b. e P 1! .b  (1  0  1) ! 2

i=2, 3

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 117 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Les forces réparties aux nœuds sont montrées sur la figure 4.17. y

a

ft3 x

3

b

[1] 1

2

x

ft2 x

Figure 4.17 : Charge équivalentes aux nœuds 2 et 3 4.4 Applications au problème ‘2D’ élément fini rectangle linéaire à 4 nœuds – Q4: Dans un programme d’élément finis on doit établir les matrices de déformation et de rigidité pour chaque élément. Cette tâche est bien simple, pour l’élément fini triangulaire, comme on vient de le voir dans ce même chapitre. Étant donné que la matrice de déformation B est constante, les procédures d’intégration deviennent simples et directes. D’où on introduit le second élément fini à contrainte plane figure 4.18, pour montrer comment sont faites ces intégrations. y 4

d.d.l internes u  u  v 

3

b x b 1

2 a

a

d.d.l nodaux u  ri   i  vi 

Figure 4.18 : Elément fini rectangulaire à 4 nœuds 1- Matrice de déformation On définit les coordonnées locales (x,y).Puisqu’il y a 4 nœuds et 2 degrés de liberté par nœud, la fonction de déplacement appropriée/ Fonction de Forme sera :

u  a1  a 2 x  a3 y  a 4 xy  v  a5  a6 x  a7 y  a8 xy

(4.65)

d’où l’appellation d’élément fini bilinéaire

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 118 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

u  N. a

Soit

1 x y N  0 0 0

où

xy 0 0 0 0 1 x y

0  xy 

(4.66)

a  a1 , a2 , a3 , a4 , a5 , a6 , a7 , a8 T

Soit les déplacements aux nœuds r tel que :

r  u1 , v1 , T

u 2 , v2 ,

u3 , v3 ,

(4.67)

u 4 , v4

A noter que la numérotation à l’intérieur de l’élément, est toujours selon le sens adopté, figure 4.18. On établit, r i pour i=1, 2, 3, 4 en utilisant (4.66), puis assembler pour donner rT en fonction de a.

u1   1  a  b ab v   0 0 0 0  1  u 2   1 a  b  ab    0 0 v 2   0 0   b ab u 3   1 a  v3   0 0 0 0    u 4   1  a b  ab v   0 0 0 0  4 

0 1 0 1 0 1 0 1

soit

r  A. a

d’où

a  A .r

et

0   a1  a  b ab   a 2 0 0 0   a3  a  b  ab  a 4 . 0 0 0   a5  a b ab   a 6  0 0 0   a7   a b  ab  a8 0

0

            

(4.68)

(4.69)

1

1

u  N . a  N . A . r  N e. r

(4.70)

(4.71)

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 119 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 La matrice inverse A-1 est donnée explicitement:

0 14 0 14 0 14 0   14   1 4a 0 1 4a 0 1 4a 0  1 4a 0     1 4b 0  1 4b 0 1 4b 0 1 4b 0    1 4ab 0  1 4ab 0 1 4ab 0  1 4ab 0  1  A   0 14 0 14 0 14 0 14     1 4a 0 1 4a 0 1 4a 0  1 4a   0  0  1 4b 0  1 4b 0 1 4b 0 1 4b    1 4ab 0  1 4ab 0 1 4ab 0  1 4ab  0

(4.72)

La matrice de déformation s’écrit comme suit :

  L ( u )  L ( N ) A1 r  B. A1 r

(4.73)

avec

0 1 0 y 0 0 0 0 B   0 0 0 0 0 0 1 x   0 0 1 x 0 1 0 y 

u x v pour y u v pour  y x

pour

(4.74)

On remarque alors que les déformations sont linéaires en x et en y . 2 Conformité de l’élément fini Q4 Puisque N est linéaire en x et y, on ne peut avoir que des mouvements linéaires au niveau des frontières. Ceci implique qu’à partir de la compatibilité nodale, la continuité des déplacements inter élément (C0 continuité), est satisfaite.  La compatibilité interne des déformations est satisfaite. 

Le mouvement de corps rigide pour la translation on prend tous les ai=0 à l’exception de a1=a50

a  à partir de u  N .a ,on déduit u   1  qui représente une translation de corps rigide a5  Considérons la relation de déformation   B. a

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 120 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0  a1   0     0  0 1 0 y 0 0 0 0   0     0 0 0 0 0 0 1 x .   0 Soit a  0 0 1 x 0 1 0 y   5   0     0   0     Pour la rotation on prend tous les ai= 0 à l’exception de a3= -a6 0.

A partir de

u  N. a

On a

 a y u 3    a3 x 

ce qui représente une rotation de corps rigide

a  0 0 a3 T

On vérifie bien que

0 0  a3

0 0

  B. a  0 comme prévu.

 Déformation constante. Prendre tous les ai=0 à l’exception de a20 A partir de (5.36)

u  a 2 x   v0

(4.75)

Ce qui représente un allongement dans le sens de x. Selon que l’on considère une différentiation directe de (4.74) ou la relation   B. a on a

  a 2 , 0, 0

T

(4.76)

ce qui représente une déformation constante dans le sens de x de la même manière  y doit être constant.  Contrainte de cisaillement – rigidité de cisaillement, Il existe une particularité importante à préciser concernant les éléments finis rectangles linéaires (Q4) à savoir qu’ils sont excessivement rigides lorsqu’ils sont appelés à représenter le mode de flexion d’une poutre qui est très important en génie civil pour le calcul des structures. 

On considère un élément en flexion pure selon la figure 4.19.

Les déplacements aux nœuds donnent un mouvement linéaire selon y pour x= a, permettant une déformation linéaire à travers l’élément, d’où une flexion pure. A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 121 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 Cependant, parce que les déplacements aux frontières de l’élément ne peuvent être que linéaires, il ne doit y avoir que des frontières en lignes droites entre les nœuds 1 et 2 et entre les nœuds 3 et 4 comme montré sur la figure ci-dessus. y

y  4

3

 x

x 

1

2



(a)

(b)

(c) Figure 4.19 : (a) Sens des déplacements en vue d’une flexion (b) Quantification de la flexion sur l’élément rectangulaire (c) Comportement réel en flexion Dans un comportement réel en flexion pure, la poutre se déforme et fléchie, en ayant une contrainte de cisaillement nulle. Pour avoir l’état (b) à partir de (c), on doit ajouter un cisaillement qui n’affecte pas le cisaillement, et inversement de (b) à (a). On doit donc corriger notre élément pour permettre l’état (c). 

Vérifions les valeurs trouvées à partir des matrices, des fonctions de forme et de déformation, de l’élément.

Pour les déplacements appliqués on a :

r    1, 0, 1, 0,  1, 0, 1, 0 T

En utilisant l’inverse explicite A

1

(4.77)

T

(4.72) on déduit a à partir de :

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 122 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 1

a  T

Soit

a  A .r

(4.78)

 . 0, 0, 0, 1, 0, 0, 0, 0 ab

(4.79)

On obtient alors, à partir de la relation u  N . a , une écriture du vecteur déplacement

  u   .x y u comme :  ab  v 0

expression prévue !

(4.80)

Et à partir de la relation   B. a on déduit :

D’où

T  

 . y, 0, x ab

x  

 .y ab

Ce qui représente une variation linéaire selon y

Fibre supérieure



x  

Fibre inférieure



x 

 a

 a

y  0

Aussi avec Seulement

(4.81)

!  xy  

. x a .b

(constante par rapport à y)

(4.82) (4.83) (4.84)

Ce qui donne plus de rigidité en flexion pure, et donc absorbe plus d’énergie, alors que la déformation de cisaillement doit être nulle. On conclut que cet élément est trop rigide en cisaillement. Les contraintes de cisaillement absorbent l’énergie de déformation qui ne peut pas être absorbée dans la cas d’une flexion pure réelle (d’où l’appellation de "cisaillement parasite"), ce qui rend l’élément trop rigide d’où l’existence de deux problèmes : 

Elément rigide à cause du cisaillement



Pas de déplacement vertical.

On peut remédier facilement au premier problème en annulant tous simplement les termes linéaires de  x y dans la 3eme ligne de la matrice B .

Soit :

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 123 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0

 0 1 0 y 0 0 0 0 B   0 0 0 0 0 0 1 x   0 0 1 0 0 1 0 0

(4.85)

Cette procédure n’affectera pas les conditions de conformité de l’élément et peut cependant affecter les matériaux anisotropes car le cisaillement est maintenant devenu constant dans l’élément. La deuxième difficulté peut être traitée par l’utilisation de modes de déplacements dits incompatibles. Remarque : Au mieux, il est recommandé de choisir l’élément "Serendipe" Q8, et qui peut représenter un état de flexion pure à condition que l’élément soit rectangle (non quadrilatère) ou l’élément Q9 (Lagrangien) qui peut-être un quadrilatère pour autant que ses cotés soient rectilignes. 3 Matrice de rigidité La formulation des expressions de la matrice de rigidité de l’élément est déduite comme auparavant à partir de l’énergie potentielle et s’écrit comme:

K .r  f (e)

K

(e)

avec

 eP .

a

a



b

b

(4.86) T

B e . D. B e . dy .dx

B e  B. A

(4.87)

1

(4.88)

A ce stade on fait appel à la quadrature de Gauss pour intégrer l’expression : C  B e .D. B e T

(4.89)

avec les éléments de la matrice qui s’expriment comme :

ci j  f (1, x, y, x 2 , y 2 , xy)

(4.90)

On remarque alors, que l’ordre maximum du polynôme dans C est p=2, d’où la nécessité de 2x2 points d’intégration selon le triangle de Pascal, figure 3.27. Bien sûr l’intégration de la matrice, comme on l’a déjà précisé auparavant, est la même que l’intégration de ses termes.

K ( e )  eP .  C. dA

ou bien

A

k i j  e p i   Ci j . dA

(4.91)

A

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 124 -

Chapitre4- Application des éléments finis de base, aux problèmes de classe C1 et C0 L’intégration selon la règle des sommes s’écrit : k i j  e P .

a

a



b

b

2

 e P .

C i j .dy dx  e P .

1 1

1

2

 W W .C m

m 1 n 1

n

ij



1

C i j ( , ). det( J ).d . d

(4.92)

( m ,  n ). det( J m n )

En se référant au système de coordonnées locales de la figure 4.18, on obtient les relations de transformation suivantes :

 x  a.   y  b.

(4.93)

 a 0 J   ; det( J )  a. b 0 b

d’où

 (  , )

(4.94)

On doit toujours calculer les valeurs de x et de y aux points d’intégration qui doivent être ensuite substituées dans la fonction C i j . Donc C i j ne sera pas réécrit en terme de  et  d’où l’expression suivante : 2

2

k i j  (a. b. eP ).  WmWn . Ci j ( xm , y n )

(4.95)

m 1 n 1

C sera est calculée à partir de (4.89) pour une intégration numérique prise avec un certain nombre de points de Gauss (  m , n ) d’où l’expression de la matrice de rigidité K

(e)

2

2



 ( a. b. eP )   Wm .Wn . B e ( x m , y n ). D. B e ( x m , y n ) m 1 n 1

T



(4.97)

L’algorithme suivant donne un aperçu sur la procédure automatique de calcul : Initialiser K  0 Faire une boucle (DO –Loop) pour m=1, 2 xm := a. ξm Faire une boucle (DO –Loop) pour n=1, 2 yn := b. ηn Former Be(xm, yn) Etablir D Multiplier

S  B e .D. B e T

Pondérer

S  S . ( a b e P .Wm .Wn ) Ajouter

K K  S Fin de boucle sur n Fin de boucle sur m

A.KADA – "Introduction à la Méthode des Eléments Finis – Master de Génie Civil" - 125 -

REFERENCES BIBLIOGRAPHIQUES

Dhatt G. & Touzot G. Une présentation de la méthode des éléments finis. Deuxième Edition, Maloine S. A., Paris 1984 Rao S.S. The Finite Element in Engineering. Second Edition Oxford: Pergamon Press Edition, 1989. Lewist P.E. & Ward J. P. The Finite Element Method: Principles and Applications. Addison-Wesley Publishers Ltd., 1991. Rockey K. C., Evans H. R., Griffiths D. W., Nethercot D.W. The Finite Element Method: A Basic Introduction, William Collin Sons & Publishers, UK 1983. Hinton E. & Owen D.R.J. Finite Element Programming. Academic Press, Londre 1977 Khennane A. Méthode des Eléments Finis : Enoncé des principes de base. Office des publications universitaires, OPU, Alger 1997 Ern A. Aide-mémoire Eléments Finis. Dunod, Paris 2005 Herbert Baaser. Development and Application of the Finite Element Method based on MATLAB. Springer-Verlag Berlin Heidelberg 2010