Méthodes numériques en génie civil et géologique – Introduction à la méthode des éléments finis Illustration des concepts de base sur un exemple élémentaire
Introduction générale • Méthode des éléments finis : Méthode numérique de résolution approchée des équation différentielles décrivant les phénomènes physiques de l’ingénierie (pas CAO !!)
• Bref historique : Mise au point chez Boeing (calcul d’ailes d’avion) en 1953 Baptisée ‘Méthode des éléments finis’ en 1960 Dès 1970, application progressive à tous les créneaux de l’ingénierie
• Domaines d’application : Analyse des structures, transferts de chaleur, mécanique des fluides, électromagnétisme, écoulements souterrains, combustion, diffusion des polluants…
Introduction générale • Méthode des éléments finis : Méthode numérique de résolution approchée des équation différentielles décrivant les phénomènes physiques de l’ingénierie (pas CAO !!)
• Bref historique : Mise au point chez Boeing (calcul d’ailes d’avion) en 1953 Baptisée ‘Méthode des éléments finis’ en 1960 Dès 1970, application progressive à tous les créneaux de l’ingénierie
• Domaines d’application : Analyse des structures, transferts de chaleur, mécanique des fluides, électromagnétisme, écoulements souterrains, combustion, diffusion des polluants…
Exemple élémentaire – définition du problème • Exemple : Tassement d’un sol homogène, isotrope, élastique soumis à son poids propre.
x, u h
γ
Exemple élémentaire – solution analytique Inconnue de base : Déplacement vertical: u(x,y,z) = u(x) Relations fondamentales : Déformations: ε x = Contraintes:
du dx
1 −ν σx = E εx = K εx ( 1 − 2ν )( 1 + ν ) σ dA
Equilibre:
γ dx dA x (σ + dσ) dA
dσ +γ = 0 dx
Exemple élémentaire – solution analytique • Equation différentielle du tassement d 2u K 2 +γ = 0 dx
Conditions aux limites u( h ) = 0
σ(0 ) = K
Déplacement nul au niveau du rocher
du ⎞ ⎟ =0 dx ⎠ x =0
Contrainte nulle à la surface libre
• Solution exacte
4
u( x)
u( x ) =
γ
( h2 − x2 )
4
2
0
2K σ ( x ) = −γ x
0
0.5
0
0
s( x)
−2
0
1
1.5
x
0.5
1
2 h
1.5
2
1
2 0
x
h
Résolution approchée – résidus pondérés Équilibre en volume
d 2u K 2 +γ=0 dx
Équilibre en surface
K
Equations à résoudre
Evaluation du résidu pondéré (forme faible)
⎛ d 2u* ⎞ du* ( 0 ) + γ ⎟⎟ dx + w( 0 )K WR = ∫ w( x ) ⎜⎜ K 2 dx 0 ⎝ dx ⎠
Fonction de pondération ayant * la même forme que u ( x )
w( h ) = 0
du( 0 ) =0 dx
h
Solution approchée respectant les conditions limites cinématiques
u* ( h ) = 0
Résolution approchée – résidus pondérés h Evaluation du résidu ⎛ d 2 u* ⎞ du* ( 0 ) + γ ⎟⎟ dx + w( 0 )K pondéré (forme faible) WR = ∫ w( x ) ⎜⎜ K 2 dx dx 0
⎝
⎠
h
h dw du* h ⎡ du ⎤ du* ( 0 ) WR = ⎢ Kw( x ) dx + ∫ w( x )γ dx + Kw( 0 ) ⎥ −K∫ dx ⎦ 0 dx 0 dx dx 0 ⎣ *
h dw du* h du* ( h ) dx + ∫ w( x )γ dx = Kw( h ) −K∫ dx 0 dx dx 0
=
w( h ) = 0 Annulation du résidu pondéré
h dw du* dx + ∫ w( x )γ dx −K∫ 0 dx dx 0 h
h⎡
⎤ dw( x ) du* ( x ) − w( x ) γ ⎥ dx = 0 ∫ ⎢K dx dx 0⎣ ⎦
Résolution approchée – résidus pondérés Choix du déplacement approché : champ linéaire Choix de la fonction de pondération : Même forme que le champ inconnu Annulation du résidu pour évaluer le paramètre a
WR = 0 h
u* ( x ) = a x + b = a ( x − h ) w( x ) = w0 ( x − h )
Æ « Meilleure » solution approchée
u* = h
⇒ ∫ Ka w0 dx − ∫ w0 ( x − h )γ dx = 0 0
0
h2 ⇒ w0 ( Ka h + γ ) = 0 2 hγ ⇒ a=− 2K
4
hγ (h−x) 2K
σ* = −
hγ 2
4
u( x) u1( x)
2
0
0
0.5
0
0
0
1
1.5
x
0.5
1
2 h
1.5
2
s ( x) s1( x)
−2
1
2 0
x
h
Résolution approchée – minimisation de l’énergie Justification physique des particularisations de la méthode des résidus pondérés : Approche énergétique Energie dans un volume infinitésimal : u
Changement d’énergie de déformation
dΠ = d U + d P =
γ dx dA
Energie totale sur l’épaisseur de sol (pour une surface unitaire) :
1 σ ε dx dA − ( γ dx dA ) u( x ) 2 Changement d’énergie potentielle
2
h ⎛ ⎞ K du ( x ) * ⎟⎟ dx − ∫ γ u* ( x ) dx Π ( u ) = ∫ ⎜⎜ 2 0 ⎝ dx ⎠ 0 h
*
Résolution approchée – minimisation de l’énergie Choix du déplacement approché : champ linéaire
u* ( x ) = a ( x − h )
Æ Evaluation de l’énergie approchée 2
h K ⎛ du* ( x ) ⎞ * ⎟⎟ dx − ∫ γ u* ( x ) dx Π ( u ) = ∫ ⎜⎜ 2 0 ⎝ dx ⎠ 0 K 2 h2 * ⇒ Π (u ) = a h+γ a 2 2 h
La meilleure approximation du champ de déplacement est celle qui minimise l’énergie totale dΠ h2 hγ = 0 ⇒ K ah+γ =0 ⇒ a=− da 2 2K A comparer avec la relation obtenue par les résidus pondérés :
γ h2 w0 ( Ka h + )=0 2
Résolution approchée – principe des travaux virtuels Justification physique des particularisations de la méthode des résidus pondérés : principe des travaux virtuels Formulation générale :
δ WI = δ WE ∀δ w
Pour le cas du tassement : h
h
∫ σ δε dx = ∫ γ δ w dx 0
0
avec
δε =
dδw dx
du( x ) d δ w( x ) K∫ dx = ∫ γ δ w dx dx dx 0 0 h
h
Si u(x) est le champ de déplacement exact, vrai quel que soit δw cinématiquement admissible ÅÆ Equilibre
Résolution approchée – principe des travaux virtuels Si le champ inconnu est approché (u*), le PTV approché (PTV*) ne peut plus être vérifié pour tout δw
du* ( x ) d δ w( x ) K∫ dx = ∫ γ δ w dx dx dx 0 0 h
h
?!
Æ Vérifié pour un champ particulier, « homothétique » au champ inconnu approché ⎧ u = a( x − h ) ⎨ ⎩δw = δa ( x − h ) *
PTV *
h
h
0
0
⇒ K ∫ a δa dx − γ ∫ δa ( x − h ) dx = 0
Formellement identique à l’approche résidu pondéré si on remplace w0 par δa Æ La fonction de pondération particulière peut être interprétée comme un champ de déplacement cinématiquement admissible
Amélioration de la solution approchée Formulation du problème approché (3 méthodes équivalentes) •
Particularisation de la méthode des résidus pondérés
•
Minimisation de l’énergie potentielle totale
•
Application du principe des travaux virtuels, avec choix d’un champ cinématiquement admissible de même forme que le champ approché
Amélioration de la solution •
Enrichir le champ inconnu : u* linéaire Æ u* parabolique
•
Pas toujours évident !! (Equations différentielles à solution non polynomiales ou sans solution, problèmes bi- ou tri-dimensionnels…)
Æ Approche « éléments finis » : Solution simple sur des sous-domaines de l’espace étudié
Principe de la méthode des éléments finis U0 h
x, u
I
U1
II
U2
Choix du champ « simple » Variation linéaire du déplacement sur couches I et II 2x 2x ) + U1 h h
I : entre U0 et U1
uI ( x ) = U 0 ( 1 −
II : entre U1 et U2 (=0)
u II ( x ) = 2 U1 ( 1 −
x ) h
⎡ h⎤ x ∈ ⎢0 , ⎥ ⎣ 2⎦ ⎡h ⎤ x ∈ ⎢ , h⎥ ⎣2 ⎦
Intérêt de ce choix: continuité du déplacement assurée à l’interface Epaisseur des couches arbitraires Æ ici, deux couches d’épaisseurs égales
Application des éléments finis à l’exemple élémentaire Approche « travaux virtuels » - la plus courante en mécanique des solides et mécanique des structures Champ approché : uI ( x ) = U 0 ( 1 −
2x 2x ) + U1 h h
u II ( x ) = 2 U1 ( 1 −
x ) h
⎡ h⎤ x ∈ ⎢0 , ⎥ ⎣ 2⎦ ⎡h ⎤ x ∈ ⎢ , h⎥ ⎣2 ⎦
du I ⎧ = K σ I ⎪ dx ⎨ du ⎪σ II = K II dx ⎩
Champ virtuel associé : d δu I ⎧ δε = ⎪ I dx ⎨ d δu II ⎪δε II = dx ⎩
2x 2x ⎧ u ( ) U δ 1 δ δU1 = − + 0 ⎪ I h h ⎨ x ⎪ δu II = 2( 1 − )δU1 h ⎩
PTV* :
h
h2
h
0
h2
δWI = ∫ σ ( x )δε ( x ) dx = ∫ σ I ( x )δε I ( x ) dx + ∫ σ II ( x )δε II ( x ) dx *
0
h
h2
h
0
0
h2
δWE = ∫ γ δu( x ) dx = ∫ γ δu I ( x ) dx + ∫ γ δu II ( x ) dx
Application des éléments finis à l’exemple élémentaire PTV* : δWI = δWE ( 2 δU 0
K K γh γh K K − 2 δU1 )U 0 + ( −2 δU 0 + 4 δU1 )U1 = δU 0 + δU1 h h h h 4 2
PTV* - Forme matricielle K⎤ ⎡ K ⎧γ h ⎫ − 2 ⎥ ⎧U ⎫ ⎢2h ⎪ 4 ⎪ 0 h ∀ [δU 0 δU1 ] = δ δ [δU 0 δU1 ] ⎢ K [ ] U U 0 1 ⎨ γ h⎬ K ⎥ ⎨⎩U1 ⎬⎭ ⎪ ⎪ ⎢− 2 4 ⎥ h h ⎩ 2 ⎭ ⎦ ⎣ K ⎡ 1 − 1⎤ ⎧U 0 ⎫ γ h ⎧1 ⎫ ⇒ 2 ⎢ ⎨ ⎬= ⎨ ⎬ h ⎣− 1 2 ⎥⎦ ⎩U1 ⎭ 4 ⎩2⎭ ⎧U ⎫ γ h h ⎧4⎫ ⇒ ⎨ 0⎬ = ⎨ ⎬ ⎩U1 ⎭ 4 2 K ⎩3⎭ ⎧ γ h2 ⎪⎪ U 0 = 2K ⇒ ⎨ 2 ⎪U = 3 γ h ⎪⎩ 1 8 K
Application des éléments finis à l’exemple élémentaire PTV* : δWI = δWE ( 2 δU 0
K K γh γh K K − 2 δU1 )U 0 + ( −2 δU 0 + 4 δU1 )U1 = δU 0 + δU1 h h h h 4 2
PTV* - Forme matricielle K⎤ ⎡ K ⎧γ h ⎫ − 2 ⎥ ⎧U ⎫ ⎢2h ⎪ 4 ⎪ 0 h ∀ [δU 0 δU1 ] = δ δ [δU 0 δU1 ] ⎢ K [ ] U U 0 1 ⎨ γ h⎬ K ⎥ ⎨⎩U1 ⎬⎭ ⎪ ⎪ ⎢− 2 4 ⎥ h h ⎩ 2 ⎭ ⎦ ⎣ Matrice de rigidité K ⎡ 1 − 1⎤ ⎧U 0 ⎫ γ h ⎧1 ⎫ ⇒ 2 ⎢ ⎨ ⎬= ⎨ ⎬ h ⎣− 1 2 ⎥⎦ ⎩U1 ⎭ 4 ⎩2⎭
Vecteur de forces nodales énergétiquement équivalentes
⎧U ⎫ γ h h ⎧4⎫ ⇒ ⎨ 0⎬ = ⎨ ⎬ ⎩U1 ⎭ 4 2 K ⎩3⎭ ⎧ γ h2 ⎪⎪ U 0 = 2K ⇒ ⎨ 2 ⎪U = 3 γ h ⎪⎩ 1 8 K
4
4
0
u( x)
0
0.5
1
1.5
2
s( x)
u1( x)
s1( x)
2
u2( x)
0
1
s2( x)
0 0
0.5
1 x
1.5
2 h
−2
2 0
x
h
Résumé des concepts de base Déplacement
Contraintes s ( x) 0
s1( x) 1
2
4
0
0
II 0.5
0.5
1
x
x
1
Noeuds
s2( x)
4
−2
u( x)
2
U2
u1( x)
0
γ
I
0
U1
0
x, u
h
u2( x)
U0
Eléments 1.5
1.5
[K].{U} = {F(γ)}
Vecteur des forces nodales équivalentes
2
Vecteur des inconnues nodales
h
h
2
Matrice de rigidité
Généralisation Concept applicable à différents types de « sous-domaines »
Elément « poutre » : linéaire à 2 noeuds Elément « plaque-membrane » : plan à 4 noeuds
Méthodes numériques en génie civil et géologique – Introduction à la méthode des éléments finis Discrétisation, interpolation, matrice de rigidité
Problème à résoudre Equation différentielle
d 2u K 2 +γ=0 dx du( 0 ) =0 dx
Condition limite faible
K
Condition limite forte
u( h ) = 0
A[ u( x )] = 0 B [ u( 0 )] = 0
⎛ d 2 u* ⎞ du* ( 0 ) WR = ∫ w( x ) ⎜⎜ K + γ ⎟⎟ dx + w( 0 )K =0 2 dx dx 0 ⎝ ⎠ h
Résidu pondéré
h
WR = ∫ w( x ) A[ u( x )] dx + w( 0 )B [ u( 0 )] = 0 0
h
WR = ∫ δu( x ) A[ u( x )] dx + δu( 0 )B [ u( 0 )] = 0 0
Problème à résoudre ⎛ d 2u ⎞ du( 0 ) WR = ∫ δu( x ) ⎜⎜ K 2 + γ ⎟⎟ dx + δu( 0 )K =0 dx 0 ⎝ dx ⎠ h
Résidu pondéré
h
WR = ∫ δu( x ) A[ u( x )] dx + δu( 0 )B [ u( 0 )] = 0 0
⎡ 0⎣ h
Intégration par partie WR = ∫ ⎢ K
du ( x ) dδu( x ) ⎤ − δu( x ) γ ⎥ dx = 0 dx dx ⎦
E [ u( x )] F [ δu( x )]
(ici, les opérateurs E et F sont identiques) h
h
0
0
WR = ∫ E [ u( x )] F [ δu( x )] dx − ∫ δu( x ) γ dx = 0
Autre exemple : poutre en flexion C
v( x )
P M d 2v χ= =− 2 EI dx T=
dM dx
p=−
d 4v r( x ) = p − EI 4 dx
x
dT dx
d 2v ⇒ M = − EI 2 dx d 3v ⇒ T = − EI 3 dx d 4v d 4v ⇒ p = EI 4 ⇒ p − EI 4 = 0 dx dx
est le résidu local d’équilibre, positif vers le bas, si v(x) n’est pas la solution exacte
Autre exemple : poutre en flexion C
v( x ) Conditions limites faibles
P
x
Conditions limites fortes
v( 0 ) = 0
dv( 0 ) =0 dx
d 2 v( L ) M ( L ) = C ⇒ − EI −C = 0 2 dx d 3v( L ) +P=0 T ( L ) = P ⇒ EI 3 dx
d 2 v( L ) −C M r ( L ) = − EI 2 dx
est le résidu d’équilibre de rotation en x=L, positif dans le sens horlogique, si v(x) n’est pas la solution exacte
d 3 v( L ) R( L ) = EI +P 3 dx
est le résidu d’équilibre de translation en x=L, positif vers le bas, si v(x) n’est pas la solution exacte
Autre exemple : poutre en flexion L
WR = ∫ r( x ) δv( x ) dx + M r ( L ) 0
L⎛
d 4v WR = ∫ ⎜⎜ p − EI 4 dx 0⎝ L d 4v
dδv( L ) + R( L ) δv( L ) = 0 dx
⎞ ⎛ ⎞ dδv( L ) ⎛ d 3v( L ) ⎞ d 2 v( L ) ⎟δv dx + ⎜ − EI − C ⎟⎟ + ⎜⎜ EI + P ⎟⎟ δv( L ) = 0 2 3 ⎟ ⎜ dx dx ⎠ ⎝ ⎠ dx ⎝ ⎠ L d 2v
2 d 2 δv d 3v d v dδv L L v dx dx [ v ] [ δ = + δ − ]0 = 0 ∫ 4 ∫ 2 0 2 3 2 dx dx dx dx 0 dx 0 dx
δv( 0 ) = 0
dδv( 0 ) =0 dx
dδv( L ) L d 2 v d 2 δv WR = ∫ p δv dx + Pδv( L ) − C − ∫ EI 2 dx = 0 2 dx dx dx 0 0 L
δWE
δWI
Principes de base DISCRETISATION DU PROBLEME 0
Domaine de résolution
h
x
Principes de base DISCRETISATION DU PROBLEME 1
2
3
4
5
0
6
7 h
Définition des nœuds et des éléments Æ Solution sous forme d’un vecteur de valeurs particulières du champ inconnu u(x) aux nœuds.
⎡ U1 ⎤ { q } = ⎢⎢ ⎥⎥ ⎢⎣U N ⎥⎦
x
Principes de base Description du champ inconnu par INTERPOLATION entre les valeurs particulières aux nœuds u( x ) = ∑ hi ( x )U i = h1( x ) i
hi(x) :
⎡ U1 ⎤ hN ( x ) ⎢⎢ ⎥⎥ = N ( x ) { q } ⎢⎣U N ⎥⎦
Fonction d’interpolation associée au nœud i = forme du champ inconnu si on impose Ui = 1, et toutes les autres inconnues nodales = 0
Æ hi(x) = 1 en x = xi et hi(x) = 0 en x = xj (j = 1…N sauf i)
Définition de la fonction de pondération (ou champ virtuel) associé : ⎡ δU1 ⎤ ⎥ = N ( x ) { δq } δu( x ) = ∑ hi ( x )δU i = h1( x ) hN ( x ) ⎢⎢ ⎥ i
⎢⎣δU N ⎥⎦
Fonctions d’interpolation 1°/ Interpolation linéaire 1 h1 h2 h3 h4 h5
Fonctions d’interpolation 2°/ Interpolation parabolique
h1
h2
h3
h4
h5
Fonctions d’interpolation 3°/ Interpolation cubique
h1
h2
h3
h4
h5
Exemple: tassement du sol ⎡ du ( x ) dδu( x ) ⎤ − δu( x ) γ ⎥ dx = 0 WR = ∫ ⎢ K dx dx ⎦ 0⎣ h
h
h
0
0
WR = ∫ E [ u( x )] F [ δu( x )] dx − ∫ δu( x ) γ dx = 0
u( x ) = ∑ hi ( x )U i = h1( x ) i
dh ( x ) dh1 ( x ) du( x ) =∑ i Ui = dx dx dx i
dh ( x ) dh1 ( x ) dδu( x ) =∑ i δU i = dx dx dx i
⎡ U1 ⎤ hN ( x ) ⎢⎢ ⎥⎥ = N ( x ) { q } ⎢⎣U N ⎥⎦
⎡ U1 ⎤ dhN ( x ) ⎢ ⎥ ⎢ ⎥ = B( x ) { q } dx ⎢⎣U N ⎥⎦
E [ u( x )] ⇒ K B( x ) { q }
⎡ δU 1 ⎤ dhN ( x ) ⎢ ⎥ = B( x ) { δq } F [ δu( x )] ⇒ B( x ) { δq } ⎢ ⎥ dx ⎢⎣δU N ⎥⎦
Exemple: tassement du sol ⎡ du ( x ) dδu( x ) ⎤ − δu( x ) γ ⎥ dx = 0 WR = ∫ ⎢ K dx dx ⎦ 0⎣ h
dδu( x ) du( x ) ⎡h ⎤ { } {q} = K dx δ q B ( x ) K B ( x ) dx ∫ ∫ ⎢ ⎥ dx dx 0 ⎣0 ⎦ h
[K ] ∫ [δu( x ) γ ]dx = δq ∫ γ {N ( x )}dx
h
h
0
0
{Q} δq [K ]{q} = δq {Q} ⇒
matrice de rigidité forces nodales énergétiquement équivalentes aux forces appliquées
[K ]{q} = {Q}
Exemple: tassement du sol En pratique, on effectue les calculs élément par élément h
NEL
0
K =1 K
∫ (∗) dx = ∑ ∫ (∗) dx avec NEL = nombre d’éléments
[K ] = ∫ {B( x )} K h
0
B( x ) dx = ∑ ∫ {B( x )} K B( x ) dx NEL
J =1 J
h
NEL
0
J =1 J
{Q} = ∫ γ {N ( x )}dx =
∑ ∫ γ {N ( x )}dx
Exemple: tassement du sol 1 h1
1
2
3
4
h2 h3 h4 h5
⎧− 1⎫ ⎪+ 1⎪ 1 1 ⎪⎪ ⎪⎪ {B( x )} = ⎨ 0 ⎬ ⎪0⎪ ⎪ ⎪ ⎪⎩ 0 ⎪⎭
h=4
⎧0⎫ ⎪− 1⎪ 2 1⎪ ⎪ ⎪⎪ {B( x )} = ⎨+ 1⎬ ⎪0⎪ ⎪ ⎪ ⎪⎩ 0 ⎪⎭
⎧0⎫ ⎪0⎪ 3 1⎪ ⎪ ⎪⎪ {B( x )} = ⎨− 1⎬ ⎪+ 1⎪ ⎪ ⎪ ⎪⎩ 0 ⎪⎭
⎧ dh1 ( x ) ⎫ ⎪ dx ⎪ ⎪ dh ( x ) ⎪ ⎪ ⎪ 2 ⎪ dx ⎪ {B( x )} = ⎪⎨ dh3 ( x ) ⎪⎬ ⎪ dx ⎪ ⎪ dh4 ( x ) ⎪ ⎪ dx ⎪ ⎪ dh5 ( x ) ⎪ ⎪ ⎪ ⎩ dx ⎭
⎧0⎫ ⎪0⎪ 4 1 ⎪⎪ ⎪⎪ {B( x )} = ⎨ 0 ⎬ ⎪− 1⎪ ⎪ ⎪ ⎪⎩+ 1⎪⎭
Exemple: tassement du sol [K ]
J
= ∫ {B( x )} K B( x ) J
J
dx
J
[K ]
[K ]
1
3
⎡+ 1 − 1 0 0 0⎤ ⎢ − 1 + 1 0 0 0⎥ ⎥ K⎢ = ⎢0 0 0 0 0⎥ ⎢ ⎥ 0 0 0 0 0 ⎢ ⎥ ⎢⎣ 0 0 0 0 0⎥⎦
[K ]
⎡0 ⎢0 K⎢ = ⎢0 ⎢ ⎢0 ⎢⎣0
⎡0 ⎢0 4 K⎢ [K ] = ⎢0 ⎢ ⎢0 ⎢⎣0
0⎤ 0 0 0 0⎥⎥ 0 + 1 − 1 0⎥ ⎥ 0 − 1 + 1 0⎥ 0 0 0 0⎥⎦ 0
0
0
2
0 0 0⎤ ⎡0 0 ⎢0 + 1 − 1 0 0 ⎥ ⎥ K⎢ = ⎢0 − 1 + 1 0 0 ⎥ ⎢ ⎥ 0 0 0 0 0 ⎢ ⎥ ⎢⎣0 0 0 0 0⎥⎦
0 0 0 0 0
0 0 0⎤ 0 0 0 ⎥⎥ 0 0 0⎥ ⎥ 0 + 1 − 1⎥ 0 − 1 + 1⎥⎦
Exemple: tassement du sol [K ] = ⎡+1 ⎢−1 K⎢ [K] = ⎢ 0 ⎢ ⎢0 ⎢⎣ 0
4
∑ [K ]
J
J =1
−1 +2 −1 0 0
0 −1 +2 −1 0
0 0 −1 +2 −1
0⎤ 0 ⎥⎥ 0⎥ ⎥ −1⎥ +1⎥⎦
Exemple: tassement du sol h
NEL
0
J =1 J
{Q} = ∫ γ {N ( x )}dx = ⎧1⎫ ⎪1⎪ 1 γ ⎪⎪ ⎪⎪ {Q} = ⎨0⎬ 2 ⎪ ⎪ 0 ⎪ ⎪ ⎪⎩0⎪⎭
{Q}
3
⎧0⎫ ⎪0⎪ γ ⎪⎪ ⎪⎪ = ⎨1⎬ 2 ⎪ ⎪ 1 ⎪ ⎪ ⎪⎩0⎪⎭
{Q}
{Q}
2
4
⎧0⎫ ⎪1⎪ γ ⎪⎪ ⎪⎪ = ⎨1⎬ 2 ⎪ ⎪ 0 ⎪ ⎪ ⎪⎩0⎪⎭ ⎧0⎫ ⎪0⎪ γ ⎪⎪ ⎪⎪ = ⎨0⎬ 2 ⎪ ⎪ 1 ⎪ ⎪ ⎪⎩1⎪⎭
∑ ∫ γ {N ( x )}dx
⎧1 ⎫ ⎪2 ⎪ γ ⎪⎪ ⎪⎪ {Q} = ⎨2⎬ 2 ⎪ ⎪ 2 ⎪ ⎪ ⎪⎩1 ⎪⎭
Exemple: tassement du sol 0 0 ⎤ ⎧U 1 ⎫ ⎧1⎫ ⎡+ 1 − 1 0 ⎪2 ⎪ ⎥ ⎪U ⎪ ⎢− 1 + 2 − 1 0 0 2 ⎥ ⎪⎪ ⎪⎪ γ ⎪⎪ ⎪⎪ K⎢ ⎢ 0 − 1 + 2 − 1 0 ⎥ ⎨U 3 ⎬ = ⎨2⎬ ⎥⎪ ⎪ 2 ⎪ ⎪ ⎢ 0 − 1 + 2 − 1⎥ U 4 2 ⎢0 ⎪ ⎪ ⎪ ⎪ ⎢⎣ 0 ⎪⎩1⎪⎭ 0 0 − 1 + 1⎥⎦ ⎪⎩U 5 ⎪⎭
⎡ +1 −1 0 0 ⎤ ⎧U1 ⎫ ⎧1 ⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ K ⎢ −1 +2 −1 0 ⎥ ⎪U 2 ⎪ γ ⎪2 ⎪ = ⎨ ⎬ ⎨ ⎬ ⎢ 0 −1 +2 −1⎥ ⎪U 3 ⎪ 2 ⎪2 ⎪ ⎢ ⎥⎪ ⎪ ⎪⎩2 ⎪⎭ 0 0 1 2 − + ⎣ ⎦ ⎩U 4 ⎭
avec U 5 = 0 (condition limite forte)
Exemple: tassement du sol
1
⎡ +1 −1 0 0 0 ⎤ ⎧U1 ⎫ ⎧1 ⎫ ⎧ 0 ⎫ ⎢ −1 +2 −1 0 0 ⎥ ⎪U ⎪ ⎪2 ⎪ ⎪ 0 ⎪ ⎥ ⎪⎪ 2 ⎪⎪ γ ⎪⎪ ⎪⎪ ⎪⎪ ⎪⎪ K⎢ ⎢ 0 −1 +2 −1 0 ⎥ ⎨U 3 ⎬ = ⎨2 ⎬ + ⎨ 0 ⎬ ⎢ ⎥⎪ ⎪ 2 ⎪ ⎪ ⎪ ⎪ 2 0 ⎢ 0 0 −1 +2 −1⎥ ⎪U 4 ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎪⎩1 ⎪⎭ ⎩⎪−V ⎭⎪ ⎣ 0 0 0 −1 +1⎦ ⎩⎪U 5 ⎭⎪
2
3
=0
4
−
K
γ U4 = −V 2
γ K ⇒ V = + U4 2
5
V
Fonctions d’interpolation Structure de la matrice de rigidité
Linéaire
⎡X ⎢X ⎢ ⎢0 ⎢ ⎢0 ⎢⎣ 0
X X X 0 0
0 X X X 0
0 0 X X X
0⎤ 0 ⎥⎥ 0⎥ ⎥ X⎥ X ⎥⎦
Parabolique
⎡X ⎢X ⎢ ⎢X ⎢ ⎢0 ⎢⎣ 0
X X X X 0
X X X X X
0 X X X X
0⎤ 0 ⎥⎥ X⎥ ⎥ X⎥ X ⎥⎦
⎡X ⎢X ⎢ ⎢X ⎢ ⎢X ⎢⎣ 0
X X X X X
X X X X X
X X X X X
0⎤ X ⎥⎥ X⎥ ⎥ X⎥ X ⎥⎦
Cubique
Fonctions d’interpolation Æ Equilibre à trouver entre: •
Amélioration de la représentation sur chaque élément fini (augmentation du degré de l’approximation polynomiale)
•
Lourdeur du système d’équations à résoudre
Æ Deux techniques pour améliorer une solution approchée Æ convergence p
•
Augmentation du degré d’interpolation
•
OU Diminution de la taille des éléments, à degré d’interpolation égal Æ convergence h
! Attention aux critères mathématiques de convergence !
Synthèse de la méthode Construction d’un modèle « éléments finis » •
Identification de la forme du domaine de résolution (1D, 2D, 3D – linéaire ou courbe – …)
•
Découpage en éléments, définition des nœuds
•
Choix des inconnues nodales
•
Choix des fonctions d’interpolation
•
Calcul des intégrales élément par élément (matrices de raideur locales)
•
Assemblage des matrices locales : constitution de la matrice de raideur du problème
Exemple d’application Équation à résoudre : d 2u( x ) + α u( x ) = 0 sur [0 , L] 2 dx Avec u(0) = A et u(L) = 0
Solution analytique : Forme générale : u( x ) = A cos( ω x ) + B sin( ω x ) d 2u( x ) + α u( x ) = − A ω 2 cos( ω x ) − B ω 2 sin( ω x ) + α A cos( ω x ) + α B sin( ω x ) 2 dx = A [ α − ω 2 ] cos( ω x ) + B [ α − ω 2 ] sin( ω x ) = 0 ⇒ω = α
Conditions limites :
⎧u( 0 ) = A cos( 0 ) + B sin( 0 ) = C ⎨ ⎩u( L ) = A cos( ω L ) + B sin( ω L ) = 0 ⎧A = C ⎪ ⇒⎨ cos( ω L ) = − B C ⎪ sin( ω L ) ⎩
Solution : u( x ) = C [cos( x α ) −
cos( α L ) sin( x α )] sin( α L )
Exemple d’application Équation à résoudre : d 2 u( x ) + α u( x ) = 0 2 dx
sur [0 , L] Avec u(0) = A et u(L) = 0
⎛ d 2u( x ) ⎞ ⎟⎟ δu( x ) dx + WR = ∫ ⎜⎜ α u ( x ) 2 dx ⎠ 0⎝ L
d 2u du dδu ⎡ du ⎤ = ∫ 2 δu dx + α ∫ u δu dx = − ∫ dx + ⎢ δu ⎥ + α ∫ u δu dx = 0 dx dx dx ⎣ dx ⎦ 0 0 0 0 0 L
L
L
Problème à résoudre par l’approche EF: du dδu ∫0 dx dx dx − α ∫0 u δu dx = 0 L
L
L
L
Exemple d’application Modèle à 4 EF de même longueur L/4 : U1
U2
U3
Fonctions d’interpolation linéaires : h1 h2 h3 h4 h5
U4
U5
Exemple d’application Coordonnées locales : x varie de 0 à L/4 sur chaque élément h1 h2 h3 h4 h5 1
dhi ( x ) 4x 4 → =− Deux fonctions à définir : hi ( x ) = 1 − L dx L 2 dhi ( x ) 4 4x 2 hi ( x ) = → = L dx L 1
Exemple d’application Calcul de la matrice K : Boucle sur les éléments h1 h2 h3 h4 h5 dh1( x ) dh1( x ) k11_1 = ∫ dx − α ∫ h1( x ) h1( x ) dx dx dx 0 0 L
Element 1:
L
dh 2 ( x ) dh 2 ( x ) k 22 _1 = ∫ dx − α ∫ h 2 ( x ) h 2 ( x ) dx dx dx 0 0 L
L
dh1( x ) dh 2 ( x ) k12 _1 = ∫ dx − α ∫ h1( x ) h 2 ( x ) dx dx dx 0 0 L
L
Exemple d’application Calcul de la matrice K : Boucle sur les éléments h1 h2 h3 h4 h5 k 22 _ 2
dh1( x ) dh1( x ) dx − α ∫ h1( x ) h1( x ) dx =∫ dx dx 0 0
k33 _ 2
dh 2 ( x ) dh 2 ( x ) dx − α ∫ h 2 ( x ) h 2 ( x ) dx =∫ dx dx 0 0
k 23 _ 2
dh1( x ) dh 2 ( x ) =∫ dx − α ∫ h1( x ) h 2 ( x ) dx dx dx 0 0
L
Element 2:
L
L
L
L
L
Exemple d’application Même principe pour éléments 3 et 4 puis assemblage des matrices élémentaires ⎡ k11_1 ⎢k ⎢ 12 _1 [K] = ⎢ 0 ⎢ ⎢ 0 ⎢ 0 ⎣
k12 _1 k 22 _1 + k 22 _ 2
⎡ 4 αL ⎢ L − 12 ⎢ 4 αL ⎢− − ⎢ L 24 [K] = ⎢ 0 ⎢ ⎢ 0 ⎢ ⎢ ⎢ 0 ⎢⎣
k 23 _ 2 0
0
0
k 23 _ 2 k33 _ 2 + k33 _ 3
0
k34 _ 3
k34 _ 3 k 44 _ 3 + k 44 _ 4
0
k 45 _ 4
0 4 αL − L 24 8 αL − L 6 4 αL − − L 24 −
0 0
0 4 αL − L 24 8 αL − L 6 4 αL − − L 24
0
−
0
0 4 αL − L 24 8 αL − L 6 4 αL − − L 24 −
0 ⎤ 0 ⎥⎥ 0 ⎥ ⎥ k 45 _ 4 ⎥ k55 _ 4 ⎥⎦
⎤ ⎥ ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ 4 α L⎥ − − ⎥ L 24 ⎥ 4 αL ⎥ − L 12 ⎥⎦ 0
Exemple d’application Résolution du système d’équations Prise en compte des conditions [ K ] {q} = {0} aux limites
⎧U1 = C ⎫ ⎪ U ⎪ 2 ⎪⎪ ⎪⎪ ⇒ [ K ] ⎨ U 3 ⎬ = {0} ⎪ U ⎪ 4 ⎪ ⎪ ⎪⎩U 5 = 0 ⎪⎭
⎧U 2 ⎫ ⎪U ⎪ Réorganisation de * ⎪⎪ 3 ⎪⎪ ⎡ K la matrice K [ K * ] ⎨U 4 ⎬ = {0} ⇒ ⎢ 11* ⎣ K 21 ⎪C ⎪ ⎪ ⎪ ⎪⎩ 0 ⎪⎭
Condensation statique
* ⎤ ⎧U I ⎫ ⎧ 0 ⎫ K12 ⎬=⎨ ⎬ * ⎥⎨ K 21 ⎦ ⎩U II ⎭ ⎩…⎭
(
* * * * [ K11 ] { U I } + [ K12 ] { U II } = { 0 } ⇒ { U I } = [ K11 ] −1 − [ K12 ] { U II }
)
Exemple d’application Exemples de solution – avec 4 et 8 éléments 10
C =10 L=1 α = 0.01
10
Vi Vbib
5
u( x)
0
0
0.1
0.2
0.3
0.4
0
10
C =10 L=1 α=2
0.5 i⋅
L 4
, ib⋅
0.6 L 8
0.7
0.8
0.9 L
,x
10
Vi Vbib
5
u( x)
0
0 0
0.1
0.2
0.3
0.4
0.5 i⋅
L 4
, ib⋅
0.6 L 8
,x
0.7
0.8
0.9 L
Exemple d’application Exemples de solution – avec 4 et 8 éléments 15.67
C =10 L=1 α=6
20
15 Vi Vbib
10
u( x) 5
0
0
0.1
0.2
0.3
0.4
0
70.862
C =10 L=1 α=9
0.5 i⋅
L 4
, ib⋅
0.6 L 8
0.7
0.8
0.9 L
,x
80
60 Vi Vbib
40
u( x) 20
0
0 0
0.1
0.2
0.3
0.4
0.5 i⋅
L 4
, ib⋅
0.6 L 8
,x
0.7
0.8
0.9 L
Exemple d’application Exemples de solution – avec 4 et 8 éléments 13.587
C =10 L=4 α=4
20
10 Vi Vbib
0
0.5
1
1.5
2
2.5
3
3.5
4
u( x) 10
− 11.413 20 0
23.07
C =10 L=4 α=3
i⋅
L 4
, ib⋅
L 8
L
,x
40
20 Vi Vbib
0
0.5
1
1.5
2
2.5
3
3.5
4
u( x) 20
− 21.93 40 0
i⋅
L 4
, ib⋅
L 8
,x
L
Exemple d’application Exemples de solution – avec 4 et 8 éléments 10
C =10 L=2 α = -2
10
Vi Vbib
5
u( x)
0
0
0.2
0.4
0.6
0.8
0
10
C =10 L=2 α = -8
1 i⋅
L 4
, ib⋅
1.2 L 8
1.4
1.6
1.8
2 L
,x
10
Vi Vbib
5
u( x)
0
0 0
0.2
0.4
0.6
0.8
1 i⋅
L 4
, ib⋅
1.2 L 8
,x
1.4
1.6
1.8
2 L
Méthodes numériques en génie civil et géologique – Introduction à la méthode des éléments finis Convergence, transformation isoparamétrique, intégration numérique, mise en charge, conditions limites
Notion de convergence Tassement sol
⎡ du ( x ) dδu( x ) ⎤ − δu( x ) γ ⎥ dx = 0 WR = ∫ ⎢ K dx dx ⎦ 0⎣ h
d 2u K 2 +γ=0 dx
Opérateurs d’ordre 1
Opérateur d’ordre 2
Poutre
d 4v − EI 4 + p = 0 dx
Opérateur d’ordre 4 Opérateurs d’ordre 2
dδv( L ) L d 2 v d 2 δv WR = ∫ p δv dx + Pδv( L ) − C − ∫ EI 2 dx = 0 2 dx dx dx 0 0 L
Notion de convergence Opérateur d’ordre 2m dans l’équation à résoudre
A[ u( x )] = 0
Æ opérateurs différentiels d’ordre m dans la forme faible L
∫ E [ u( x )] F [ δu( x )] dx
0
Raffinement progressif du modèle EF Æ Convergence vers la solution exacte Nécessité de respecter des critères pour assurer la convergence du modèle Deux types de convergence : •
Convergence h : raffinement du maillage sans modifier l’interpolation
•
Convergence p : enrichissement de l’interpolation sans changer le maillage
Critères de convergence Critère 1 : critère de continuité L’interpolation du champ inconnu doit être: 1.
De classe Cm dans l’élément Æ Forme approchée du champ inconnu et ses dérivées jusqu’à l’ordre m continues dans l’élément
2.
De classe Cm-1 aux frontières de l’élément Æ Forme approchée du champ inconnu et ses dérivées jusqu’à l’ordre m-1 continues aux frontières des éléments (nœuds du maillage dans le cas 1D)
1.
Garantit qu’on peut calculer les matrices de rigidité élémentaires
2.
Garantit qu’on peut effectuer l’assemblage des matrices élémentaires
Critères de convergence Conséquence Interpolations classique de u(x) :
Condition 1 : OK Condition 2 : OK uniquement si m = 1 !!
Pour des équations différentielles d’ordre 2m, il est nécessaire de discrétiser les dérivées du champ inconnu jusqu’à l’ordre m-1 Æ Elément conforme
Critères de convergence Exemple : Poutre en flexion: d 4v p = dx 4 EI 1
2
3
2m = 4 Æ discrétiser les dérivées d’ordre 1 4
5
6
0
L V4 et
⎡ V1 ⎤ ⎢θ ⎥ ⎢ 1⎥ {q}= ⎢ # ⎥ ⎢ ⎥ ⎢VN ⎥ ⎢⎣θ N ⎥⎦
7
dV ⎞ ⎟ = θ4 dx ⎠ 4 h1
Fonctions d’interpolation cubiques (4 fonctions par EF)
h2
h3
h4
x
Critères de convergence Exemple : Poutre en flexion: x3 x2 x3 x2 x x3 x2 x3 x 2 v( x ) = v1 ( 2 3 − 3 2 + 1 ) + θ1A( 3 − 2 2 + ) + v 2 ( −2 3 + 3 2 ) + θ 2 A( 3 − 2 ) A A A A A A A A A
x3 x2 h1 ( x ) = 2 3 − 3 2 + 1 A A x3 x2 x h2 ( x ) = A( 3 − 2 2 + ) A A A
x3 x2 h3 ( x ) = −2 3 + 3 2 A A x3 x 2 h4 ( x ) = A( 3 − 2 ) A A
h1
h2
h3
h4
Critères de convergence Critère 2 : critère de complétude Le champ approché et ses dérivées jusqu’à l’ordre m doivent pouvoir prendre des valeurs constantes arbitraires dans l’élément Æ Le polynôme doit être complet au moins jusqu’au degré m Cas particuliers •
Modes rigides
Si U1 = U2 = A,
•
Modes homogènes
Alors u(x) = A pour x ∈ [0,ΔL]
Possibilité de représenter les modes à « déformation » (=dérivée d’ordre m) constante, pour permettre le passage à la limite quand ΔLÆ 0
Critères de convergence Assouplissement des critères 1.
Vérification du critère uniquement si ΔLÆ 0 : convergence non monotone
2.
Éléments non conformes : éléments ne vérifiant pas les critères de continuité Éléments surcompatibles : éléments avec continuité supérieure à (m-1) aux frontières Æ Patch test : découpage quelconque du domaine et vérification globale sur les cas de déformations homogènes.
Convergence h Convergence h : raffinement du maillage sans changer le degré d’interpolation
Augmentation du nombre de nœuds, mais conservation de la simplicité des calculs.
Convergence p Convergence p : enrichissement de l’interpolation sans changer le maillage Interpolation linéaire 1
h1
2 U2
U1
u( x ) = U1 h1( x ) + U 2 h2 ( x ) = U1 ( 1 −
x x ) +U2 L L
h2
Interpolation parabolique : nécessité d’introduire un paramètre supplémentaire UA (par exemple valeur du champ inconnu à la moitié de l’élément) 1 U1
2 UA
U2
Convergence p 1
2
U1
UA
U2
u( x ) = U1 h1( x ) + U 2 h2 ( x ) + U A hA ( x )
h1
DDL connectés h2
hA
Mode « bulle » (DDL interne)
Convergence p Complication des calculs + augmentation du nombre de DDL Mais possibilité de condensation statique des DDL internes
[K ]{q} = {R} EL
⎡ k11 ⇒ ⎢⎢ k 21 ⎢⎣k A1
k12 k 22 k A2
k1 A ⎤ ⎧U1 ⎫ ⎧ R1 ⎫ ⎪ ⎪ ⎪ ⎪ k 2 A ⎥⎥ ⎨U 2 ⎬ = ⎨ R2 ⎬ k AA ⎥⎦ ⎪⎩U A ⎪⎭ ⎪⎩ 0 ⎪⎭
⎡ K Ext K' ⎤ ⎧U Ext ⎫ ⎧ R ⎫ ⇒ ⎢ T ⎥ ⎨ U ⎬ = ⎨0⎬ K ' k ⎩ ⎭ AA ⎦ ⎩ A ⎭ ⎣ ⎧[ K Ext ] U Ext + [ K' ] U A = { R } ⇒⎨ T [ K ' ] U Ext + k AA U A = 0 ⎩
⇒ UA = −
1 k AA
[ K' ] T { U ext }
⎛ [ K' ] [ K' ] T ⎞ ⎟ { U ext } = { R } ⇒ ⎜⎜ [ K Ext ] − ⎟ k AA ⎝ ⎠ ⇒ [ K Ext ] * { U ext } = { R }
Puis assemblage des matrices de rigidité élémentaires modifiées
Transformation isoparamétrique - Coordonnées naturelles Cas d’un domaine d’intégration curviligne Æ
Changement de variable
Abscisse curviligne s y
ξ
sj si
-1 x
Interpolation entre si et sj : linéaire, parabolique… Æ approximation de la forme du domaine Linéaire :
⎧x j ⎫ ⎧ xi ⎫ ⎧x⎫ = h1( ξ ) ⎨ ⎬ + h2 ( ξ ) ⎨ ⎬ ⎨ ⎬ y ⎩ ⎭ sur ij ⎩ yi ⎭ ⎩yj ⎭ 1 avec h1( ξ ) = ( 1 − ξ ) 2 1 et h2 ( ξ ) = ( 1 + ξ ) 2
y
1
Transformation isoparamétrique - Coordonnées naturelles Cas d’un domaine d’intégration curviligne Æ
Changement de variable
Abscisse curviligne s y
A
ξ
sj
si
-1 x
Interpolation entre si et sj : linéaire, parabolique… Æ approximation de la forme du domaine ⎧x j ⎫ ⎧ xi ⎫ ⎧ xA ⎫ ⎧x⎫ = ξ ξ h ( ξ ) h ( ) h ( ) + + ⎨ ⎬ ⎨ ⎬ 2 ⎨ ⎬ ⎨ ⎬ A 1 Parabolique : y y y ⎩ ⎭ sur ij ⎩ yA ⎭ ⎩ i⎭ ⎩ j⎭ Ajout d’un point intermédiaire
1 1 avec h1( ξ ) = ξ ( ξ − 1 ); h2 ( ξ ) = ξ ( ξ + 1 ) 2 2 et hA ( ξ ) = ( 1 + ξ )( 1 − ξ )
1
Transformation isoparamétrique - Coordonnées naturelles Transformation isoparamétrique = Utilisation des mêmes fonctions d’interpolation pour le champ inconnu que pour la forme du domaine de définition du problème s2
ξ
sA -1
0
1
s1
h1(ξ) h2(ξ) hA(ξ)
⎧s( ξ ) = h1( ξ ) s1 + h2 ( ξ ) s2 + hA ( ξ ) s A ⎨ ⎩u( ξ ) = h1( ξ ) u1 + h2 ( ξ ) u2 + hA ( ξ ) u A ⎧ x( ξ ) = h1 ( ξ ) x1 + h2 ( ξ ) x 2 + h A ( ξ ) x A ⎪ y( ξ ) = h ( ξ ) y + h ( ξ ) y + h ( ξ ) y ⎪ 1 1 2 2 A A ⎨ ⎪u( ξ ) = h1 ( ξ ) u1 + h2 ( ξ ) u 2 + h A ( ξ ) u A ⎪⎩v( ξ ) = h1 ( ξ ) v1 + h2 ( ξ ) v 2 + h A ( ξ ) v A
Transformation isoparamétrique - Coordonnées naturelles Calcul de la matrice de rigidité
Tassement sol
L
∫ E [ u( s )] F [ δu( s )] ds
h
0
dhi ( ξ ) d n hi ( ξ ) E [ u ] = E0 ∑ hi ( ξ )U i + E1 ∑ U i + " + En ∑ Ui n dx dx i i i ⎛ dN ( ξ ) d n N( ξ ) = ⎜⎜ E0 N ( ξ ) + E1 + " + En dx dx n ⎝ = B( ξ ) { q }
⎞ ⎟{ q } ⎟ ⎠
⎛ dN ( ξ ) d n N( ξ ) + " + Fn F [ δu ] = ⎜⎜ F0 N ( ξ ) + F1 dx dx n ⎝ = C( ξ ) { δq } = δq { C( ξ )}
⎞ ⎟ { δq } ⎟ ⎠
∫K
0
du( x ) dδu( x ) dx dx dx
E 0 = 0 ; E1 = K ; E 2 = ..... = 0
F0 = 0 ; F1 = 1; F2 = ..... = 0
Poutre
d 2 v d 2 δv dx ∫ EI 2 2 dx dx 0
L
L
E0 = E1 = 0 ; E2 = EI ; E3 = ..... = 0
0
F0 = F1 = 0 ; F2 = 1; F3 = ..... = 0
δq [ ∫ { C( ξ )} B( ξ ) ds ] { q } δq [ K ] { q }
Transformation isoparamétrique - Coordonnées naturelles L
L
[ K ] = ∫ { C( ξ )} B( ξ ) ds ] 0
K ij = ∫ Ci ( ξ ) B j ( ξ )ds 0
Jacobien de la transformation isoparamétrique :
ds( ξ ) J(ξ ) = dξ
J ( ξ ) > 0 ∀ξ !!!!
L
K ij = ∫ Ci ( ξ ) B j ( ξ )ds 0
+1
= ∫ Ci ( ξ ) B j ( ξ ) J ( ξ ) dξ −1
⎛ dhi ( ξ ) ⎞ = ∫ Ci ( ξ ) B j ( ξ ) ⎜⎜ ∑ si ⎟⎟ dξ −1 ⎝ N dξ ⎠ +1
Intégration numérique ⎛ dhi ( ξ ) ⎞ K ij = ∫ Ci ( ξ ) B j ( ξ ) ⎜⎜ ∑ si ⎟⎟ dξ −1 ⎝ N dξ ⎠ +1
+1
Nip
−1
IP =1
∫ f ( ξ ) dξ = ∑ wIP f ( ξ IP ) f ( ξ IP ) = Ci ( ξ IP ) B j ( ξ IP
⎛ dhi ( ξ IP ) ⎞ ) ⎜⎜ ∑ si ⎟⎟ dξ ⎠ ⎝N
Méthodes possibles • Rectangle • Trapèze • Simpson (parabole) Æ Nbre d’IP important Plus efficaces: Gauss : intègre exactement un polynôme de degré 2n-1 avec n points Lobatto : intègre exactement un polynôme de degré 2n-3 avec n points
Intégration numérique Choix du nombre de points d’intégration 1.
Évaluation du degré du polynôme à intégrer et choix de NbIP pour obtenir l’intégrale exacte
2.
Si NbIP devient très grand, on peut se limiter.
3.
•
Méthode des EF = méthode approchée !
•
Constatation: les erreurs de l’approximation EF et de l’intégration approchée se compensent souvent.
NbIP minimum: être capable de calculer le volume de l’EF (ici, sa longueur) Æ en fonction du degré d’interpolation 1
L = ∫ ds = ∫ J ( ξ ) dξ −1
Mise en charge Tassement sol
⎡ du ( x ) dδu( x ) ⎤ − δu( x ) γ ⎥ dx = 0 WR = ∫ ⎢ K dx dx ⎦ 0⎣ h
d 2u K 2 +γ=0 dx
Mise en charge
d 4v − EI 4 + p = 0 dx
Poutre
dδv( L ) L d 2 v d 2 δv WR = ∫ p δv dx + Pδv( L ) − C − ∫ EI 2 dx = 0 2 dx dx dx 0 0 L
Mise en charge
Mise en charge Ajout d’un terme indépendant du champ inconnu dans le problème de base L
L
L
WR = ∫ A [ u( x )] w( x ) dx = ∫ A [ u( x )] δu( x ) dx = ∫ [A[ u( x )] − f ( x )]δu( x ) dx *
*
0
0
0
L
L
0
0
⇒ WR = ∫ E [ u( x )] F [ δu( x )] dx − ∫ f ( x ) δu( x ) = 0
Fonction de pondération
⎡ δU1 ⎤ δu( x ) = ∑ hi ( x )δU i = h1( x ) " hN ( x ) ⎢⎢ # ⎥⎥ = N ( x ) { δq } = δq {N ( x )} i ⎢⎣δU N ⎥⎦ L
L
0
0
WR = ∫ F [ δu( x )] E [ u( x )] dx − ∫ δu( x ) f ( x ) dx = 0 L
→
∫ δq { C( x )} 0
→
L
B( x ) { q } dx − ∫ δq { N ( x )} f ( x ) dx = 0 0
L ⎡ L ⎤ δq ⎢[ ∫ { C( x )} B( x ) dx ] { q } − [ ∫ { N ( x )} f ( x ) dx ] ⎥ = 0 ∀ δq 0 ⎣ 0 ⎦
Mise en charge L
L
0
0
[ ∫ { C( x )} B( x ) dx ] { q } = ∫ f ( x ){ N ( x )} dx [ K ]{ q }= { F }
Matrice de rigidité
Vecteur des force nodales équivalentes Vecteur des inconnues nodales L
Fi = ∫ f ( x ) hi ( x ) dx = 0
∑ ∫ f ( x ) h ( x ) dx i
NbEF EF
Fi : force généralisée équivalente = Grandeur nodale représentative des contributions des forces généralisées agissant sur les éléments aboutissant au nœud considéré
Conditions aux limites Conditions naturelles : •
Imposition de forces généralisées
•
Soit on impose Fi, soit on impose f(x) et on en déduit les Fi
Conditions essentielles : •
Valeurs imposées des inconnues nodales
•
Pour un problème d’ordre 2m (Æ forme faible d’ordre m), au minimum m conditions essentielles, sinon [K] est singulière (mécanisme)
•
m conditions = système « isostatique »
•
+ de m conditions = système « hyperstatique »
Conditions aux limites Quelques remarques •
On ne peut imposer une condition essentielle que sur un degré de liberté discrétisé
•
On ne peut pas imposer une condition essentielle et une condition naturelle au même DDL;
•
Imposition de conditions essentielles : soit par élimination, soit par pénalisation * ⎡ k11 ⎢k ⎢ 21 ⎢k31 ⎢ ⎣k 41
k12 k 22 k32 k 42
k13 k 23 k33 k 43
k14 ⎤ ⎧U1* ⎫ ⎧ R1 ⎫ ⎪ ⎪ k 24 ⎥⎥ ⎪U 2 ⎪ ⎪⎪ F2 ⎪⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
k12 k 22
k13 k 23
k32
k33
k 42
k 43
k14 ⎤ ⎧U1 ⎫ ⎧ R1 ⎫ ⎪ ⎪ k 24 ⎥⎥ ⎪U 2 ⎪ ⎪⎪ F2 ⎪⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
⎡k11 + k11P ⎢ ⎢ k 21 ⎢ k31 ⎢ ⎢⎣ k 41
k12
k13
k 22 k32
k 23 k33
k 42
k 43
⎡ k11 ⎢k ⎢ 21 ⎢ k31 ⎢ ⎣k 41
k14 ⎤ ⎧U1* ⎫ ⎧ R1 + R1P ⎫ ⎥⎪ ⎪ ⎪ ⎪ k 24 ⎥ ⎪U 2 ⎪ ⎪ F2 ⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎥⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
Conditions aux limites Quelques remarques •
On ne peut imposer une condition essentielle que sur un degré de liberté discrétisé
•
On ne peut pas imposer une condition essentielle et une condition naturelle au même DDL;
•
Imposition de conditions essentielles : soit par élimination, soit par pénalisation * ⎡ k11 ⎢k ⎢ 21 ⎢k31 ⎢ ⎣k 41
k12 k 22 k32 k 42
k13 k 23 k33 k 43
k14 ⎤ ⎧U1* ⎫ ⎧ R1 ⎫ ⎪ ⎪ k 24 ⎥⎥ ⎪U 2 ⎪ ⎪⎪ F2 ⎪⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
k12 k 22
k13 k 23
k32
k33
k 42
k 43
k14 ⎤ ⎧U1 ⎫ ⎧ R1 ⎫ ⎪ ⎪ k 24 ⎥⎥ ⎪U 2 ⎪ ⎪⎪ F2 ⎪⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
⎡k11 + k11P ⎢ ⎢ k 21 ⎢ k31 ⎢ ⎢⎣ k 41
k12
k13
k 22 k32
k 23 k33
k 42
k 43
⎡ k11 ⎢k ⎢ 21 ⎢ k31 ⎢ ⎣k 41
k14 ⎤ ⎧U1* ⎫ ⎧ R1 + R1P ⎫ ⎥⎪ ⎪ ⎪ ⎪ k 24 ⎥ ⎪U 2 ⎪ ⎪ F2 ⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪U 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎥⎦ ⎪⎩U 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
Conditions aux limites Cas particulier : éléments infinis Pseudo-transformation isoparamétrique ξ -1
x1
1
x( ξ ) = x1 + α
x
1+ ξ 1− ξ
→ J(ξ ) =
dx 2α = dξ ( 1 − ξ )2
1− ξ 1+ ξ U1 + U2 2 2 2 du du dξ ⎛ 1 ⎞ (1 − ξ ) → = = ⎜ ( U 2 − U1 ) ⎟ dx dξ dx ⎝ 2 ⎠ 2α → dans WR
u( x ) =
x2 → ∞
Exemple 1 Écoulement en milieu poreux (Loi de Darcy) h1 h2
Potentiel h(x) – débit q(x) – Flux extrait Q(x) dx (q+dq) dA
q dA
dq +Q = 0 dx
q = −α Q dV
dh dx
d 2 h( x ) α = Q( x ) 2 dx
Exemple 1 Écoulement en milieu poreux (Loi de Darcy) h1
Hi h2
d 2 h( x ) α = Q( x ) dx 2
⎧ H1 ⎫ ⎪ ⎪ h( x ) → ⎨ # ⎬ ⎪H ⎪ ⎩ N⎭
Discrétisation du potentiel 2m = 2 Æ m = 1 Æ Interpolation linéaire [ K ]{ H } = { F }
Avec Fi = ∫ Q( x ) hi ( x ) dx L
Fi : résultante pondérée de Q(x) au voisinage du nœud i
Exemple 1 Écoulement en milieu poreux (Loi de Darcy) h1
Hi h2
Conditions naturelles Imposer Q(x) Æ en déduire les Fi associées OU Imposer directement des Fi Conditions essentielles •
m = 1 Æ au moins une condition, sinon [K] singulière (mode rigide : le modèle ne peut pas fixer le niveau de référence du potentiel)
•
Imposer le potentiel par exemple aux extrémités
Exemple 1 Conditions essentielles - élimination K' ⎤ ⎧ H1 ⎫ ⎧ F1 ⎫ ⎧ k1 H1 + K' h = F1 ⎡ k1 ⎢ K' T K * ⎥ ⎨ h ⎬ = ⎨ F ⎬ → ⎨ K' T H + K * h = F ⎦⎩ ⎭ ⎩ ⎭ 1 ⎣ ⎩ → [ K * ] { h } = { F } − H1 { K' T } Cas particulier : H1 = 0 → F1 = k1 H1 + { K' } { h } ≈ Q1
Conditions essentielles - pénalisation HA
H1
H2
H3
H4
H5
EF supplémentaire, sans modification de la fonction de pondération
H6
H7
δh( x ) = h1( x )δH1 d δh( x ) dh1( x ) dx
=
dx
δH1
Exemple 1 Modification à la matrice de rigidité dh( x ) H1 − H A = dx L
∫α
EF
α dh( x ) dδh( x ) dh ( x ) δH1 dx dx = ∫ ( H1 − H A ) 1 dx dx L EF dx
⎡α dh ( x ) ⎤ ⎡ α dh ( x ) ⎤ =⎢ ∫ 1 dx ⎥ H1 δH1 − ⎢ H A ∫ 1 dx ⎥ δH1 = k11_ S H1 δH1 − F1_ S δH1 L dx L dx EF ⎣ ⎦ ⎣ EF ⎦
⎡k11 + k11_ S ⎢ k 21 ⎢ ⎢ k31 ⎢ ⎣ k 41 H1 ≈ H A q = −α ⇒
α L
k12 k 22 k32
k13 k 23 k33
k 42
k 43
k14 ⎤ ⎧ H1 ⎫ ⎧ F1 + F1_ S ⎫ k 24 ⎥⎥ ⎪⎪ H 2 ⎪⎪ ⎪⎪ F2 ⎪⎪ ⎨ ⎬=⎨ ⎬ k34 ⎥ ⎪ H 3 ⎪ ⎪ F3 ⎪ ⎥ k 44 ⎦ ⎪⎩ H 4 ⎪⎭ ⎪⎩ F4 ⎪⎭
∀q
dh( x ) α = ( H1 − H A ) dx L >> ⇒ k11_ S >>
Limites pratiques à k11_S
Exemple 2 Poutre en flexion
d 4v p = dx 4 EI
2m = 4 Æ m = 2 Æ Discrétisation des dérivées
Conditions d’appui : •
On peut imposer des conditions sur v (Æ déplacements) et sur v’ (Æ rotations)
•
Minimum 2 conditions essentielles pour éviter les modes rigides
Exemple 2 Poutre en flexion Charges équivalentes Par EF, 4 DDLs Æ 4 efforts équivalents Fi ( x ) =
∫ p( x ) h ( x ) dx i
EF
