Intro Aux Elements Finis GAIE(1)

Département de génie civil Laboratoire de mécanique des sols

ÉCOLE POLYTECHNIQUE FÉDÉRALE DE LAUSANNE

Cycle postgrade: Géologie Appliquée à l'Ingénierie et à l'Environnement Michel Dysli

B3-5: Introduction aux éléments finis Documents distribués

2ème édition, juin 1997

Liste des Cycle postgrade : Géologie Appliquée à l'Ingénierie et à l'Environnement documents Module B3-5: Introduction aux éléments finis M. Dysli mise à jour: 1997-06-10 Document No 1 2 2a 3 4 4a 5 5a 6 7 8 8a 9 9a 10 11 11a 12 13 14 15 16 Annexe

Titre Buts du cours Bibliographie Scalaire, vecteur, tenseur Equations de contraintes et déformations (cours mécanique des sols B2-2) Exemple d'un tenseur: le tenseur des contraintes Les équations et leur solution Solution équations de contraintes et déformations (cours mécanique des sols B2-2) Différences finies Calcul d’un écoulement souterrain avec un tableur Eléments finis Autres éléments finis Exemples de discrétisation en éléments finis Exemples de discrétisation en éléments finis (suite) Fonctions de transformation et d'interpolation Fonctions de transformation et d'interpolation (suite) Matrices d'élasticité Relations entre les noeuds et l'intérieur de l'élément (1ère partie) Relations entre les noeuds et l'intérieur de l'élément (suite) Détermination de la matrice déformation - déplacement Application numérique Application numérique (suite) Principe de la résolution par la méthode des déplacements FEM, BEM, KEM Rappel du calcul matriciel

Cycle postgrade : Géologie Appliquée à l'Ingénierie et à l'Environnement Module B3-5: Introduction aux éléments finis M. Dysli

Buts du cours

Buts du cours "Introduction aux éléments finis" (8 heures):

• Comprendre les principes de la méthode des éléments finis appliquée au calcul des contraintes et déformations et à celui des écoulements souterrains. • Avoir les bases nécessaires pour s'en servir dans des cas simples. Et non pas : Etre capable d'écrire des programmes sur ordinateur relatifs à cette méthode... Il faudrait pour cela un cours de 100 à 200 heures!

Méthode: approche physique par des cas particuliers


Document No 1 ajout 1997-06-24

Bibliographie 1. Acquisitions conseillées BATHE K.-J. (1996). Finite element procedures . Prentice-Hall. BRITTO A. M., GUNN M. J. (1987). Critical State Soil Mechanics via Finite Elements. Ellis Horwood Limited, Halsted Press: a division of John Wiley & Sons, DHATT G., TOUZOT G. (1984). Une présentation de la méthode des éléments finis. Maloine SA Paris. PRAT M. (1995). La modélisation des ouvrages, Hermès, Paris. ZIENKIEWICZ O. C. (1979). La méthode des éléments finis. 3e éd., McGraw Hill. ZIENKIEWICZ O. C. (1977). The Finite Element Method. 3rd ed., McGraw Hill.

2. Autres CHEN W. F., MIZUNO E. (1990). Nonlinear Analysis in Soil Mechanics. Elsevier DESAI C. S., SIRIWARDANE H. J. (1984). Constitutive laws for engineering materials. PrenticeHall.

Des éditions plus récentes pourraient exister.

Document No 2


Scalaire, vecteur, tenseur En physique, lorsqu'un simple nombre suffit à définir une grandeur (température, densité, module élastique, par ex.), ce nombre est un scalaire. Si la grandeur a non seulement une dimension mais aussi une direction, on l'appellera alors vecteur. Une force, une accélération, un déplacement sont des vecteurs. Dans un système d'axes cartésien à 3 dimensions, il faut 3 composantes au moins pour définir un vecteur, soit par exemple deux angles et une intensité, ou 3 projections sur les axes du système choisi. Une modification du système de référence induit une modification des 3 composantes définissant le vecteur. Il est cependant possible de dériver de ce vecteur une quantité scalaire qui est indépendante de sa direction. Par exemple, si l'on considère un vecteur déplacement de composantes d1, d 2, d3, la grandeur de ce vecteur: d = d12 + d22 + d32 est indépendante de sa direction; c'est l'invariant du vecteur. La notion de tenseur est plus abstraite. On pourrait lui donner la signification physique d'une représentation d'un champ: champ de vitesses, champ de contraintes par exemple. Dans la définition des composantes d'un tel champ, nous devons nous référer deux fois à des directions, par exemple, pour un champ de contrainte: premièrement à l'orientation du solide élémentaire et deuxièmement à l'orientation des contraintes proprement dites sur les faces de ce solide. Il faudrait donc 3 . 3 = 9 composantes pour définir un tel tenseur. Cependant un tenseur cartésien est symétrique et seules 6 composantes suffisent à le définir complètement. Si pour un vecteur il est possible de dériver un invariant, un tenseur donne trois invariants dans un système à 3 dimensions. Type

Scalaire

Vecteur

Tenseur

0

1

2

poids volumique

déplacement

contrainte

Notations

γ

 d1    d2  d   3

Notations indicielles

γ

di

σij

Nombre de composantes dans un système de coordonnées à 3 dimensions

30

31 = 3

32 = 9

Valeurs indépendantes

1

3

9 en général

Invariants

1

Ordre de la matrice Exemple

 σ11 σ12 σ  21 σ 22 σ 31 σ 32

σ13  σ 23  σ 33 

6 avec symétrie: σij = σji 1

3

Document No 02a


Les équations des contraintes et déformations Notation indicielle (d'Einstein)

Nbre Nbre équations inconnues

Notation xyz

∂σ'xx ∂σyz ∂σzx ∂u + + + + Fx = 0 ∂x ∂y ∂z ∂x

1. Equations d'équilibre i = 1,2,3 = direction ∂σ'ij ∂u des normales + + Fj = 0 aux faces ∂xi ∂xi j = 1,2,3 = direction des axes σ = σ' + u u = pression interstitielle Fj = force par unité de volume En général, F x (ou F 1) = F y (ou F 2) = 0 sauf si forces d'inertie (séisme par ex.) Fz (ou F3) = -ρ.g car sol = milieu pesant

qui s'écrit habituellement :

3 z ou 3

∂σ'x ∂τyz ∂τzx ∂u σz + + + + Fx = 0 ∂x ∂y ∂z ∂x τzx τzy ∂τxy ∂σ'y ∂τzy ∂u τyz + + + Fy = 0 τ + σy xz ∂x ∂y ∂z ∂y σx τxy τyx y ou 2 ∂τzx ∂τyz ∂σ'z ∂u + + + Fz = 0 + ∂x ∂y ∂z ∂z x ou 1

2. Relations déformation-déplacement (équations géométriques) ∂ux ∂ux ∂ux ∂u εx = 1 + = = ∂ui ∂u j i = 1,2,3 2 εij = 1 + ∂x ∂x ∂x ∂x j = 1,2,3 2 ∂x j ∂xi ux = u ; uy = v ; uz = w ui,j = vecteur des déplacements ∂w ∂v εz = εy = ∂z Formulation de Lagrange : la position d'un ∂y point est décrite en fonction de ses coordonγ ∂u ∂v ∂u ∂v nées initiales (avant déformation). εxy = xy = 1 + γxy = + ∂y ∂x 2 2 ∂y ∂x ∂v ∂w Par opposition à la formulation d'Euler où la γyz = + ∂z ∂y position d'un point est décrite en fonction de ses coordonnées après déformation. ∂w ∂u + γzx = ∂x ∂z

εx εy ε εij = εz xy εyz εzx

i = direction des σ, τ j = direction des ε

loi linéaire loi non linéaire

sij est la matrice du matériau et son inverse sij-1 la matrice d'élasticité . sij est symétrique, elle a donc, pour le cas le plus général d'un matériau anisotrope : 36 - 15 = 21 termes. Exemples (élastiques linéaires) :

fil

σx

σx

εx = s 11 σx = 1/E σx σz

σz εx = s31 σ 3 = -ν/E σz εy = s32 σ 3 = -ν/E σz εz = s33 σ3 = 1/E σz σz

σy

σx σz

voir § 2.2.3 du cours polycopié

σx σy σ σij = z σxy σyz σzx

sij =

9

9

15

il manque 6 équations

(lois constitutives, équations d'état, équations de déformabilité)

R(εij, σij) = 0 ou εij = f(σij) = sij . σij

3σ 3τ τij = τji

u, v, w εx , εy , εz γxy, γyz, γzx

Total intermédiaire

3. Relations contrainte-déformation

sij = cte => sij = f(σ, ε, t) =>

6

6

6

0

15

15

s11 s12 s13 s14 s15 s16 s21 s31 etc. s41 s51 s61

Dans le cas de la mécanique des sols, les σ des exemples seraient négatifs (compression = +)

Total final

Document No 3


Exemple d'un tenseur: le tenseur des contraintes 1, 2, 3 = axes des contraintes principales τ = 0 sur les faces z

z σz τzx τxz σx

τxy

τzy τyz σy

' z' σ z

τ z'x'

y

τ x'z'

τyx

σ x'

x

z 3 σ 3

τ z'y'

y' ' τ y'zσ y'

τ x'y'

y

σ2

2 y

τ y'x'

x

σ1

x

x' 1

tenseur des contraintes :

σx τ yx τ zx

τ xy σy τ zy

tenseur des contraintes principales : σ1 0 0 σ2 0 0 σ3 0 0

τ xz τ yz σz

τ xy = τ yx avec: τ xz = τ zx τ zy = τ yz

Invariants du tenseur des contraintes (sans démonstration) I1 = σ x + σ y + σ z

I1 = σ1 + σ 2 + σ3

I2 = σ xσ y + σ y σ z + σ zσ x − τ xy2 − τ yz2 − τ zx2

I2 = σ 1σ 2 + σ 2 σ 3 + σ 3 σ1 I3 = σ1 σ 2 σ3

I3 = σ x σ y σ z + 2 τ xy τyz τzx − σ x τyz2 − σ y τzx2 − σ z τxy2

σ' = σ - u σx' τ xy ' ' σ y' τ yx τ zx ' τ zy '

τ xz' σx τ yz' = τ yx σ z' τ zx

τ xy σy τ zy

σ ' = contrainte effective σ = contrainte totale u

τ xz u 0 0 τ yz - 0 u 0 σz 0 0 u

pression interstitielle = tenseur sphérique

= pression interstitielle

σ1' 0 0

0 σ 2' 0

σ1 0 0 = 0 σ 3' 0

0 σ2 0

0 u 0 0 0 - 0 u 0 σ3 0 0 u


Document No 4

Les équations et leur solution Ecoulement souterrain 1 ∂h ∇h = c v ∂t avec

h = charge hydraulique t = temps cv = coefficient de

Contraintes-déformations ∂σ ' ∂u + + Fx = 0 équilibre des forces ∂x ∂x avec

consolidation

σ'= u= Fx = x=

Conduction thermique Champs électriques etc.

contraintes effectives pression interstitielle force volumique axe cartésien quelconque

∂vi ∂vj relation déformation + 2 ∂xj ∂xi déplacement avec vi,j = vecteur des déplacements

εij = 1

εij =

ε = f(σ ')

déformations

relation contraintes - déformations

.. M .v + K .v = R - F + U avec

M= K= v= R= F= U=

équilibre des déplacements (document No 15)

matrice des masses. matrice de rigidité. vecteur des déplacements. vecteur des forces extérieures. vecteur des forces intérieures vecteur des pressions interstitielles.

Cas géométriquement simple : solutions analytiques

Cas plus complexe: différences finies ou éléments finis Excavation avec paroi moulée WS

puits filtrant

couche 1

WS

couche 2

Différences finies : document No 5

Eléments finis : document No 6


Document No 4a

Solution des équations de contraintes et déformations Charges simples pouvant conduire à une solution analytique :

Semi-infini élastique

Solutions

Abaques

analytiques

et

formules

Recordon Ed. (1980). Abaques du cours polycopié de mécanique des sols de l'EPFL. Giroud J.-P. (1975). Tables pour le calcul des fondations, Tomes 1 et 2. Dunod, Paris. Poulos H. G., Davis E. H. (1974). Elastic solutions for soil and rock mechanics. John Wyley & Sons.


Document No 5

Différences finies Ecoulement souterrain

Contraintes-déformations WS

puits filtrant

couche 1

WS

couche 2

y Fi+1,j - Fi-1,j ∂F = ∂x 2∆

j+2

Fi,j+1 - Fi,j-1 ∂F = ∂y 2∆

j+1

∂2 F ∂x2 ∂2 F ∂y2

=

=

i-2

i-1

i

i+1

i+2

x

j

Fi+1,j - 2Fi,j + Fi-1,j 2∆2

j-1

Fi,j+1 - 2Fi,j + Fi,j-1 2∆2

∆

∆

j-2

F = contraintes, déformations, déplacements, températures, pressions, moments, charge hydraulique, etc. Exemple: écoulement souterrain z

palplanche h = 24

24

h0 = h2/4 + h3/2 + h 4/4

2∆2 h = 18

h = 24

h6

à savoir: hi+1,j - 2hi,j + hi-1,j

h = 18

h5

2 2 Equation pour régime permanent: ∂ h + ∂ h = 0 2 ∂x ∂ y2 h = charge hydraulique = z + u / γw avec u = pression interstitielle

hi,j =

+

hi,j+1 - 2hi,j + hi,j-1 2∆ 2

hi-1,j + hi,j-1 + hi+1,j + h i,j+1 4

=0

ou h 0 =

h1 + h 2 + h 3 + h 4 4 h4

Cas particuliers aux limites :

h1

0

h0

imperméable h2

h0 = h1/4 + h3/4 + h4/2 h0 = h1/2 + h2/4 + h4/4

h0 = h 1/4 + h 2/4 + h3/4 + h 5/8 + h 6/8

Excellent moyen = tableur

h3

Calcul d'un écoulement souterrain avec un tableur 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

A

B

C

D

E

F

G

z 24 22 20 18 16 14 12 10 8 6 4 2 0

0 24 24 24 24 24 24 24 24 24 24 24 24 24

2 24 23.92 23.83 23.74 23.66 23.59 23.54 23.50 23.48 23.47 23.46 23.46 23.46

4 24 23.83 23.66 23.48 23.31 23.16 23.05 22.98 22.94 22.92 22.91 22.91 22.91

6 24 23.75 23.49 23.21 22.94 22.70 22.53 22.44 22.39 22.37 22.36 22.36 22.36

8 24 23.67 23.32 22.94 22.53 22.15 21.94 21.85 21.82 21.81 21.81 21.82 21.82

10 24 23.61 23.19 22.70 22.10 21.42 21.23 21.20 21.22 21.24 21.26 21.28 21.29

=(E9+F10+G9+F8)/4

H x 12 24 23.59 23.14 22.57 21.75 20.21 20.37 20.51 20.61 20.68 20.72 20.75 20.76

I

J

K

L

M

N

14

16

18

20

22

24

18 18.36 18.81 19.53 19.86 20.04 20.14 20.20 20.23 20.24

18 18.32 18.68 19.08 19.36 19.54 19.65 19.71 19.74 19.75

18 18.25 18.51 18.76 18.96 19.10 19.20 19.25 19.28 19.29

18 18.17 18.33 18.49 18.62 18.71 18.78 18.82 18.84 18.85

18 18.08 18.16 18.24 18.30 18.35 18.39 18.41 18.42 18.42

18 18 18 18 18 18 18 18 18 18


Document No 6

Eléments finis Ecoulement souterrain

Contraintes-déformations WS

puits filtrant

couche 1

WS

couche 2

WS

uy,Fy hm

hj hn

hi

N qM-N h p ho

hk

O

hi = charge hydraulique qi = débit

j n

m

i

N

p

M

o hl

qO

M qM

qO-N

Eléments rectangulaires

ux,Fx

k

x l

O ui = déplacements Fi = forces

y

Document No 7


ajout 1997-06-24

Autres éléments finis exemples parmi de nombreux autres

barres

seul l'effort de traction / compression est transmis

poutre

une unité

2D

Efforts dans le plan

3D

σ

joint cinématique

corps 1

Elément de cisaillement

τ

τ = c + σ tg Φ Intersection ?

joints

corps 2

Document No 8


Exemples de discrétisation en éléments finis z

paroi moulée

axe de symétrie

couche 1

seulement déplacement vertical

couche 2

2D couche 3

Réseau initial couche 4

y aucun déplacement possible

z

naissance du buton après étape excavation 2 e1

couche 1

couche 2

e = étapes d'excavation

e2 e4 e5

e3

Réseau après plusieurs pas d'excavation

couche 3

couche 4

y

De par la symétrie seul 1/4 de l'excavation est modélisée

3D Elément isoparamétrique à 21 noeuds


Document No 8a

Exemples de discrétisation en éléments finis (suite)


Document No 9 modification 1997-06-24

Fonctions de transformation et d'interpolation (1ère partie) Approche très simple n'utilisant que les forces, les déplacements, les contraintes et les déformations. Les mêmes principes s'appliquent par analogie aux écoulements souterrains et à tous les autres champs. 1) Compatibilité des déplacements entre les côtés de deux éléments voisins

mauvais

2) Fonction de transformation géométrique: Transformation des coordonnées réelles en coordonnées de référence Exemple pour un élément triangulaire

η

uy

y

ux xk

xj xi

x

Elément réel Coordonnées locales

h(ξ) = fonction de transformation géométrique Par exemple, pour l'élément triangulaire :

0, 1

x( ξ , η) = h 1 ( ξ , η)⋅ xi + h 2 (ξ , η) ⋅ x j + h 3 ( ξ , η) ⋅ xk 1, 0

ξ

0, 0 Elément de référence Coordonnées naturelles

 xi    x( ξ , η) = {h }  xj     xk 

De la même manière :  yi    y( ξ , η) = {h }  yj     yk 

Document No 9a


modification 1997-06-24

Fonctions de transformation et d'interpolation (suite) Approche très simple n'utilisant que les forces, les déplacements, les contraintes et les déformations. Les mêmes principes s'appliquent par analogie aux écoulements souterrains et à tous les autres champs. 3) Fonctions d'interpolation: Passage des valeurs aux noeuds (variables nodales) à celles en tout point d'un élément Exemple pour un élément triangulaire ui,j,k = variables nodales = variables attachées aux noeuds de l'élément réel comme des déplacements, des forces ou une pression interstitielle. Passage des valeurs aux noeuds (variables nodales u) à celles en tout point x,y ou ξ,η du triangle :  u   i  ui  η η η η u( , ) = h ( , ) h ( , ) h ( , ) ξ {1ξ }uj  2 ξ 3 ξ     u( x, y) = {h1( x, y) h 2 ( x, y) h 3 ( x, y)} u j    h(ξ,η) = fonction d'interpolation sur uk  u  k l'élément de référence h(x) = fonction d'interpolation sur l'élément réel Pour le triangle (interpolation linéaire) : Pour le triangle (interpolation linéaire) : 1 h 1 (ξ , η ) = 1 − ξ − η h 1 ( x, y) = (y − y j )(x j − x) − (xk − x j )(y j − y) 2A k 1 h 2 ( x, y) = h 2 ( ξ , η) = ξ [( y − yk )( xk − x) − ( xi − xk ) (yk − y)] 2A i 1 h 3 ( x, y) = h 3( ξ , η) = η (y − y )( xi − x) − (x j − xi )( yi − y) 2A j i

[

]

[

]

2 A = (xk − x j)(yi − y j )− (xi − x j )(yk − y j)

Elément: isoparamétrique sub-paramétrique super-paramétrique

h(ξ) = h(ξ) ordre h(ξ) < ordre h(ξ) ordre h(ξ) > ordre h(ξ)

4) Jacobien Toutes les expressions des fonctions d'interpolation qui impliquent des dérivées de u (variable nodale) en x, y sont transformées en dérivées de ξ et η grace à la matrice de transformation dite matrice jacobienne ou Jacobien J : u = quelque chose ∂u   ∂ x ∂ξ  ∂ξ  =  ∂u   ∂ x ∂η  ∂η

∂ y  ∂u  ∂ξ ⋅ ∂ x    ∂ y  ∂u  ∂η  ∂ y 

{∂ξ}= [ J ] ⋅{∂ x}

{∂ x} = [ J] −1⋅ {∂ξ} Problème lors de l'inversion si, par exemple, le déterminant de J est nulle

Document No 10


ajout 1997-06-24

Matrices d'élasticité Forme générale de la relation σ - ε pour un corps élastique: Première solution: décomposition du tenseur des contraintes en en tenseur sphérique σm et en un tenseur déviatorique σo. Les tenseurs correspondant des déformations sont εm et εo. L'indice o de σ et ε est la lettre o et non pas le chiffre 0. élasticité linéaire: σ = 3K . ε m

σ

élasticité non linéaire : σ = 3K . ε

m

m

σo = 2G . εo

avec : K = coefficient de compressibilité [kPa par ex.] G = module de glissement ou de cisaillement [kPa[

m

σo = 2G . εo

ε

Deuxième solution: usage de vecteurs de contraintes et de déformations εx εy εz εxy εyz εzx

σx σy σz = σxy σyz σzx

εij = sij . σij

z

s11 s12 s13 s14 s15 s16 s21 s31 etc. s41 s51 s61 et

σij = sij-1 . εij

σz

τzx τxz σx

τxy

τzy τyz σy τyx

y

x

Rappel: σxy = τxy 2εxy = γxy

forme la plus utilisée

sij-1 est l'inverse de la matrice sij

y

Simplifications ( 2D ) εz ≠ 0

σij =

1 E σij = ν (1-ν 2) 0

=

ν 1

0 0 εij = [ D]⋅ ε ij (1-ν) 0 D = matrice d'élasticité 2

z

Déformations planes (εz = 0) 1-ν ν ν 1-ν 0

0

0 εij = [ D] ⋅ ε ij 0 1-2ν 2

σij =

εz εij = εr εθ γrz

z z

Symétrie de révolution σz σr σθ τrz

0

E σij = (1+ν) (1-2ν)

=

εx εij = εy γxy

y

1-ν E σij = (1+ν) (1-2ν)

ν 0 1-ν 0 0 0

0

ν ν 1-ν 0

0 0 εij 0 1-2ν 2

x

εz

σij =

σx σy τxy

x

z

εx εij = εy γxy

σ

σx σy τxy

0

Contraintes planes (σz = 0)

σz ≠ 0

σz θ

r

τzr τrz

σθ

σr


Document No 11 modification 1997-06-24

Relations entre les noeuds et l'intérieur de l'élément (1ère partie) Approche très simple n'utilisant que les forces, les déplacements, les contraintes et les déformations. 1)

Relation entre les déformations à l'intérieur de l'élément et les  ∂u  déplacements aux noeuds : Matrice des déformations - déplacements  ∂x  v1 v2 v5  ∂v  u2 u  ∂y  ε (x,y)= document No 35 module B22 : 1 u5   1 v  ∂u ∂v  8 Une matrice B différente pour chaque point x,y 5 v6 2  u 8 ∂ y + ∂ x εx,y u6 8 6 (1) {ε (x,y)} = [B (x,y)] uv noeuds v4 v3 noeuds 7 u 4 4 Dimensions (lignes . colonnes): u3 3 v7 u, v ε B u .6 élément triangulaire à 3 noeuds 7 3 3 6 x .

{

y

3 3

38 3.16

8 16

}

élément rectangulaire à 4 noeuds élément rectangulaire à 8 noeuds

2)

Relation entre les contraintes à l'intérieur de l'élément et les forces aux noeuds : Matrice des forces aux noeuds - contraintes σx Fy2 Fy1 Fy5 σ ( x,y) = σy τxy Fx2 Fx1 Fx5 (2) F = [ C(x,y) ] ( x, y) σ { noeuds} { } 1 Fy8 2 5 Fy6 Dimensions (lignes . colonnes): σx,y Fx6 F 8 x8 σ F C 6 F y4 élément triangulaire à 3 noeuds 6. 3 6 3 Fy3 7 Fx4 4 élément rectangulaire à 4 noeuds 8 3 8. 3 3 y Fx3 élément rectangulaire à 8 noeuds 16 3 16.3 Fx7 F Pas utilisée par la suite y7 x

3) Points de Gauss et méthode des résidus pondérés Dans le cas de l'élément isoparamétrique avec h = h > 1 (par exemple interpolation cubique) on peut utiliser, pour passer des déplacements aux noeuds aux déformations à l'intérieur de l'élément, la méthode d'intégration numérique de Gauss qui permet de calculer, avec une très bonne approximation, les déformations et les contraintes en des points déterminés: les points de Gauss, dont les positions permettent d'obtenir la précision maximale. La méthode de Gauss fait partie d'un ensemble de méthodes dénommé: méthode des résidus pondérés. Cette méthode générale d'intégration numérique – qui utilise des fonctions de pondérations – permet de passer d'un système d'équations aux dérivées partielles à une formulation intégrale en utilisant les fonctions de d'interpolation (polynomiale par exemple). La précision de la formulation intégrale dépend de la position du point considéré dans l'élément. Ces formulations intégrales sont nombreuses, par exemple: • méthode de Gauss déjà citée  • méthode de Ritz   = environ Gauss • méthode de Legendre  • méthode de Galerkine (à la mode aujourd'hui!)  • méthode de Newton-Cotes => différente des autres et moins précises.

{

}

Dans la relation {ε (x,y)} = [B (x,y)] uv noeuds , la matrice B est formée de termes provenant des noeuds fonctions de pondérations, des fonctions d'interpolation et des fonctions de transformation géométrique. Un exemple de constitution de cette matrice est donné sur le document no 12.

Document No 11a



Relations entre les noeuds et l'intérieur de l'élément (suite) Approche très simple n'utilisant que les forces, les déplacements, les contraintes et les déformations.

4) Relation forces - déplacements avec usage de la matrice d'élasticité D et de la matrice B

{Fnoeuds} = [K ] {uvnoeuds } noeuds [K ] =

∫S [B ]T [ D] [B] dS

dS = det [ J] dξ dη

(5)

[K ] = ∫ [A] dξ dη (6) S

[K ] =

∑ i,j

avec [K] = matrice de rigidité de l'élément

(3)

(4)

Dimensions de K (lignes . colonnes): 6.6 élément triangulaire à 3 noeuds élément rectangulaire à 4 noeuds 8.8 16.16 élément rectangulaire à 8 noeuds

sans démonstration, intégration sur le surface de l'élément

ξ, η = coordonnées naturelles, J = Jacobien (document 9a), det = déterminant avec [Α] = [Β]Τ [D] [B] det [J]

αi,j [A ] i,j (7) ou [A ] i,j est évaluée à chaque point ξi et ηi

[α ] i,j sont des constantes qui dépendent des valeur de ξi et ηi (intégration de Gauss)

Il reste à assembler tous les éléments. La méthode générale utilisée pour cela fait l'objet du document 15. Nous allons avant dire encore quelques mots sur les méthodes d'intégration numérique, dont celle de Gauss (document 12).

Document No 12



Détermination de la matrice déformation - déplacement Deux dimension, élément à 4 noeuds avec 4 points d'intégration (degré n = 2) y, v

η

noeud 2

noeud 1

0, +1 4

2 4 -0,577

ξ

4 +0,577

+1, 0 déplacements v4 u4

4 -0,577 -1, 0

3 4 -0,577

Fonctions d'interpolation: (sans démonstration) Noeuds : 4 5 6 7 8 -1/2h 8 h1 = 1/4 (1 + ξ)(1 + η) -1/2h 5 h2 = 1/4 (1 - ξ)(1 + η) -1/2h 5 -1/2h 6 h3 = 1/4 (1 - ξ)(1 - η) -1/2h 6 -1/2h 7 h4 = 1/4 (1 + ξ)(1 - η) -1/2h 7 -1/2h 8 h5 = 1/2 (1 - ξ2)(1 + η) h6 = 1/2 (1 - ξ)(1 - η2) h7 = 1/2 (1 - ξ2)(1 - η) h8 = 1/2 (1 + ξ)(1 - η2) h9 = 1/2 (1 - ξ2)(1 - η 2)

Point de Gauss

1

y4

noeud 4 0, -1

x4

noeud 3 ainsi:

x = 1/4 (1 + y = 1/4 (1 + u = 1/4 (1 + v = 1/4 (1 +

ξ)(1 + η)x1 + 1/4 (1 ξ)(1 + η)y1 + 1/4 (1 ξ)(1 + η)u1 + 1/4 (1 ξ)(1 + η)v1 + 1/4 (1 -

x, u

z = h1z1 + h2z2 + h3z3 + ... + hizi avec z = x, y, u ou v

ξ)(1 + η)x2 + 1/4 (1 ξ)(1 + η)y2 + 1/4 (1 ξ)(1 + η)u2 + 1/4 (1 ξ)(1 + η)v2 + 1/4 (1 -

ξ)(1 - η)x3 + 1/4 (1 + ξ)(1 - η)y3 + 1/4 (1 + ξ)(1 - η)u3 + 1/4 (1 + ξ)(1 - η)v3 + 1/4 (1 +

Nous savons déjà que (doc. No 11 ou module B22) : ∂u ∂v εx = ∂ u εy = ∂ v γxy = + ∂x ∂y ∂y ∂x

∂x /∂ξ = 1/4 (1 + η)x1 - 1/4 (1 + ∂x /∂η = 1/4 (1 + ξ)x1 + 1/4 (1 ∂y /∂ξ = 1/4 (1 + η)y1 - 1/4 (1 + à savoir en dérivant (1) : ∂y /∂η = 1/4 (1 + ξ)y1 + 1/4 (1 ∂u /∂ξ = 1/4 (1 + η)u1 - 1/4 (1 + ∂u /∂η = 1/4 (1 + ξ)u1 + 1/4 (1  ∂  ∂v /∂ξ = 1/4 (1 + η)v1 - 1/4 (1 + ∂  ∂ x -1  ∂ξ  ∂v /∂η = 1/4 (1 + ξ)v1 + 1/4 (1    = [J ] ⋅  ∂ ∂     ∂ y  ∂η 

et que (Jacobien, document No 9a) :  ∂   ∂x  ∂ξ   ∂ξ  =  ∂   ∂x  ∂η   ∂η

∂y   ∂  ∂ξ   ∂ x  ⋅  ∂y   ∂  ∂η   ∂ y 

 ∂  ∂  et  ∂ξ   ∂ x  son   = [ J ]⋅    ∂  inverse:  ∂  ∂ y  ∂η 

 ∂ u ∂ x 1 -1 1 + ηj ainsi :   = [ J ] ⋅ ij 1 + ξi  ∂ u 4 ∂ y à η = η   i

-(1 + ηj) 1 - ξi

0 0

0 0

-(1 - ηj) -(1 - ξi)

0 0

ξ)(1 - η)x4 ξ)(1 - η)y4 ξ)(1 - η)u4 ξ)(1 - η)v4

1 + ηj 1 + ξi

0 0

-(1 + ηj) 1 - ξi

0 0

-(1 - ηj) -(1 - ξi)

1 - ηj

0 . u    -(1 + ξi) 0  v 

(3)

1 - ηj . u  -(1 + ξi)  v 

(4)

0 0

(1)

η)x2 - 1/4 (1 - η)x3 + 1/4 (1 - η)x4 ξ)x2 - 1/4 (1 - ξ)x3 - 1/4 (1 + ξ)x4 η)y2 - 1/4 (1 - η)y3 + 1/4 (1 - η)y4 ξ)y2 - 1/4 (1 - ξ)y3 - 1/4 (1 + ξ)y4 (2) η)u2 - 1/4 (1 - η)u3 + 1/4 (1 - η)u4 ξ)u2 - 1/4 (1 - ξ)u3 - 1/4 (1 + ξ)u4 η)v2 - 1/4 (1 - η)v3 + 1/4 (1 - η)v4 ξ)v2 - 1/4 (1 - ξ)v3 - 1/4 (1 + ξ)v4

ξ = ξj

 ∂ v ∂ x 1 -1 0   = [J ] ⋅ ∂ v 4 ij 0    ∂ y  à η = ηi

9 -1/4h 9 -1/4h 9 -1/4h 9 -1/4h 9 -1/4h 9 -1/4h 9 -1/4h 9 -1/4h 9

 u1 v  1  u2 u   v2 u =   = u  v  3 v  u3  4  v4

          

ξ = ξj

En évaluant les relations (3) et (4), on peut établir le matrice de transformation déformation - déplacement, à savoir la matrice B, en un point ξi, η j : εij = Bij . u. Les indices i et j indiquent que la transformation est évaluée au point ξi, ηj. Par exemple, si les axes x et y de l'exemple à 4 noeuds sont confondus avec les axes ξ et η et si l'élément est un carré de dimension 2 de coté, le Jacobien est une matrice unitaire et ainsi:  ∂ u /∂ x   εx  1 + ηj 0 -(1 + ηj) 0 0 -(1 - η j) 0 1 - ηj     1 εij =  εy  =  ∂ v / ∂ y  d'où : B = 0 1 + ξi 0 1 - ξi 0 -(1 + ξi) 0 -(1 + ξi) 4   γ  xy 1 + η 1 ξ -(1 + η ξ η ξ ) ) -(1 ) -(1 ) -(1 + 1 + ξ i j i j i j i 1 - ηj ∂ u /∂ y + ∂ v /∂ x    On peut constater que les valeurs de la matrice B dépendent des coordonnées naturelles ξi et ηj

Document No 13


Application numérique y, v noeud 2

noeud 1

η

2

4

ξ

4m

En appliquant les relations (1) du document No 12: x = 3ξ et y = 2η

1

x, u

3

 ∂x  ∂ξ Jacobien = J=   ∂x  ∂η

∂y  ∂ 3ξ ∂ 2 η ∂ξ   ∂ξ ∂ξ   3 0 = = ∂y  ∂ 3ξ ∂ 2 η  0 2  ∂η   ∂η ∂η 

Dans ce cas particulier J est indépendant des ξ et η noeud 3

noeud 4 6m

 0,333 0  Son inverse, J-1 =   0 0,5 

Avec les relations (3) et (4) du document No 12, on peut maintenant déterminer la matrice B, par exemple pour le point de Gauss No 1 (ξi = -0,5774, ηj = -0,5774) :  ∂ u  ∂ x  1  0,333 0  1 + ηj 0 -(1 + ηj) 0 -(1 - ηj) 0 1 - ηj 0   =   0 0,5  1 + ξi 0 1 - ξi 0 -(1 - ξi) 0 -(1 + ξi) 0  ∂ u 4 ∂ y

u  1  =  v 4

 ∂ v  ∂ x  1  0,333 0  0 1 + ηj 0 -(1 + ηj) 0 -(1 - ηj) 0 1 - ηj  =    ∂ v  4  0 0,5  0 1 + ξi 0 1 - ξi 0 -(1 - ξi) 0 -(1 + ξi) ∂ y

u  1  =  v 4

et ainsi B1 = 1 4

0.14 0 -0.14 0.21 0 0.79

0 0

-0.53 -0.79

0 0

0.53 -0.21

0  u 0  v 

+ 0 0

0.14 0.21

0 0

-0.14 0.79

0 -0.53 0 -0.79

0 0.53 u 0 -0.21  v 

0.14 0 -0.14 0 -0.53 0 0.53 0 0 0.21 0 0.79 0 -0.79 0 -0.21 0.21 0.14 0.79 -0.14 -0.79 -0.53 -0.21 0.53

Calculons enfin les déformations au point de Gauss No 1 (ε1) pour le déplacement des noeuds suivant (cas trivial pour contrôle) : noeud 2

noeud 1 0,5 m

4m 1 noeud 3

noeud 4 6m

 u1   0  v   0.5   1   u2   0   εx  0.14 0 -0.14 0 -0.53 0 0.53 0        ε1 =  εy  = B  uv2  = 14 0 0.21 0 0.79 0 -0.79 0 -0.21  0.5 0  3 γ  0.21 0.14 0.79 -0.14 -0.79 -0.53 -0.21 0.53  v3   0   xy u   0   4    v4   0    εx   0     ε1 =  εy  =  -0.125   γ     xy  0

c.q.f.d. ! Suite sur document No 14

Document No 14


ajout 1997-06-24

Application numérique (suite)

Il reste maintenant à constituer la matrice de rigidité K de l'élément au moyen de la relation (7) du document No 11. Les coefficients α sont donnés sur le tableau ci-dessous (sans démonstration).

∑

[K] =

i,j

αi,j [A ] i,j = α1 [A] 1 + α2 [A] 2 + α3 [A] 3 + α4 [A] 4

Nbre de pts de Gauss

Dans notre cas α1 = α2 = α3 = α4 = 1

2 4 6

No des points de Gauss

Pour calculer A, il faut connaître la matrice d'élasticité D. En admettant E = 20'000 et ν = 0,3 et des contraintes planes (doc. 10) : 21'978 6'593 0 6'593 21'978 0 0 0 7'692

[D] =

292.3 159.5 317.2 188.8 159.5 425.2 246.9 1316.3 317.2 246.9 1957.9 -595.2 188.8 1316.3 -595.2 5184.0 -1091.2 -595.2 -1183.9 -704.9 -595.2 -1587.2 -921.5 -4913.1 481.7 188.8 -1091.2 1111.2 246.9 -154.3 1269.8 -1587.2

[A] 1 = [B]1T [D] [B]1 =

2407.4 595.2 -1797.8 704.9 -1091.2 -1111.2 481.7 -188.8

[A] 2 = etc.

595.2 1165.5 921.5 576.0 -1269.8 -1587.2 -246.9 -154.3

 Fx1  F   y1   Fx2  F   y2   Fx3   Fy3  F   x4   Fy4 

=

k11

k12

8730.5 3571.4 3571.4 12699.0 -2961.3 274.7 -274.7 3784.6 -4364.7 -3571.4 -3571.4 -6348.7 -1404.5 -274.7 274.7 -10134.9

-2961.3 274.7 8730.5 -3571.4 -1404.5 -274.7 -4364.7 3571.4

704.9 576.0 -2221.6 5924.3 921.5 -4913.1 595.2 -1587.2

α

0,0 ±0,5774 ±0,7746 0,0000 ±0,8611 ±0,3400

2.0000 1.0000 0,5555 0,8888 0,3479 0,6521

-1091.2 -595.2 481.7 246.9 -595.2 -1587.2 188.8 -154.3 -1183.9 -921.5 -1091.2 1269.8 -704.9 -4913.1 1111.2 -1587.2 4072.9 2221.6 -1797.8 -921.5 2221.6 5924.3 -704.9 576.0 -1797.8 -704.9 2407.4 -595.2 -921.5 576.0 -595.2 1165.5 -1091.2 -1111.2 481.7 -188.8 -1269.8 -1587.2 -246.9 -154.3 -1183.9 704.9 -1091.2 595.2 921.5 -4913.1 595.2 -1587.2 1957.9 595.2 317.2 -246.9 595.2 5184.0 -188.8 1316.3 317.2 -188.8 292.3 -159.5 -246.9 1316.3 -159.5 425.2

{Fnoeuds} = [K] {uv noeuds } noeuds

[K] = 1.[A] 1 + 1.[A] 2 + 1.[A] 3 + 1.[A] 4 k21

-1797.8 921.5 4072.9 -2221.6 -1183.9 704.9 -1091.2 595.2

8

ξ

[K] =

matrice de rigidité d'un élément

symétrique kij = kji

-274.7 3784.6 -3571.4 12699.0 274.7 -10134.9 3571.4 -6348.7

-4364.7 -3571.4 -1404.5 274.7 -3571.4 -6348.7 -274.7 -10134.9 -1404.5 -274.7 -4364.7 3571.4 274.7 -10134.9 3571.4 -6348.7 8730.5 3571.4 -2961.3 -274.7 3571.4 12699.0 274.7 3784.6 -2961.3 274.7 8730.5 -3571.4 -274.7 3784.6 -3571.4 12699.0

 u1 v  1  u2 v  2  u3  v3 u  4  v4

          

F   x1   Fy1 

u 

u 

u 

u 

i = no noeud où l'on désire connaître les forces

 Fy2 

 v1 

 v2 

 v3 

 v4 

connus

= k11  1  + k12  2  + k13  3  + k14  4   v1   v2   v3   v4  etc. Les sous-matrices k sont ij j = no noeud où les F   u1   u2   u3   u4  aussi elles-mêmes symétriques x2 déplacements sont   = k21   + k22   + k23   + k24   Appuis : noeud 2

c1, c 2, c 3, c 4 = constantes

noeud 1

k14 = k41 =

c1 0 0 0

k24 = k42 =

c2 0 0 0

k34 = k43 =

c3 0 0 0

k44 =

c4 0 0 0

noeud 4 noeud 4 noeud 3

v4 = 0 d'où :

u4 = v4 = 0 d'où : k14 = k41 = k24 = k42 = k34 = k43 = k44 = 0 0 0 0

Document No 15



Principe de la résolution par le méthode des déplacements Point de calcul

ne figure pas dans l'équation

Equation d'équilibre des déplacements .. . v = R - F+ U M.v +

C

,ε

R2

2

1

3

4

K

K

U

M

v

des σ

R1

R k

k31 k

k

k

k11 k 21 k

k

k13 k12 k 22 k

{σ}=

k23 k

k

k

k

k

k

k

k

k

k

k32 k

k

k

k24 k34

k44 k

k

k

k

k 42 k

k

k43 k

k

k

k

k

k

k

k

k33 k

Relation entre le vecteur des forces aux noeuds F de l'élément et celui des déplacements v correspondants :

{F} =

eb

an

k

{ε}

voire élasto-plastique D = loi de comportement

k

k14 k

D

D = matrice d'élasticité

k41

de

k

{v}

rd

e signifie "de l'élément"

rg

eu

symétrique

Ke

la

matrice K : voir document No 11a

Matrice de rigidité de l'ensemble de la structure : arrangement des sous-matrices kij en fonction de la numérotation globale de tous les noeuds du réseau. Dans cette matrice les termes kij sont en fait la somme des termes occupant la même position (provenant des éléments voisins).

k 11 k 12 k13 k 14 e

k 21 k 22 k23 k 24

K = k k k k 31 32 33 34 k 41 k 42 k 43 k 44

avec : M = K= v= R= F=

matrice de rigidité de l'élément qui établit la relation entre le vecteur des forces F et le vecteur des déplacements v aux noeuds 1,2,3,4. F = Ke . v sous-matrice relative à chaque noeud, de dimension k ij = . dl dl, dl = degrés de liberté du noeud

matrice des masses. matrice de rigidité. u  vecteur des déplacements, dénommé u ou   dans les précédents documents.  v vecteur des forces extérieures. vecteur des forces aux noeuds du système; forces équivalentes aux contraintes effectives. U = vecteur de la résultante aux noeuds des pressions interstitielles.


Document No 16

FEM, BEM, KEM

FEM = Finite Element Method = MEF = Méthode des Eléments finis

Un seul élément avec résolution à sa frontière qui sera "tué" pour simuler l'excavation

Simulation creuse tunnel

BE(M) = Boundary Element (Method)

Elément eau avec résolution à sa frontière

Eléments avec résolution à leur frontière

∞

∞ ∞

∞

= (Méthode des) Eléments aux frontières (Méthode des éléments finis avec éléments spéciaux où les déplacements ne sont connus qu'à leurs "frontières")

∞ Effet séisme sur un barrage poids

t = 0,26 sec.

t = 0,51 sec.

t = 0,75 sec.

Tous les éléments sont indépendants avec joints cinématiques entre eux.

KEM = Kinematic Element Method = Méthode des Eléments cinématiques (Méthode des éléments finis avec éléments indépendants et joints cinématiques) Méthode dénommée aussi: DEM = Discret Element Method

Intro Aux Elements Finis GAIE(1)

Recommend Documents