Examen09 Corrige
Examen09 Corrige
Examen09 Corrige
C OURS D ’A NALYSE DES S TRUCTURES M ÉCANIQUES PAR LA M ÉTHODE DES E LEMENTS F INIS
(MEC 568)
contrôle non classant (24 mars 2009, 2 heures)
Documents autorisés : livre du cours ; documents et notes de PC
On se place dans le cadre élastique linéaire HPP isotherme, de sorte que N (x, t) est relié à la déformation
∂u(x, t)/∂x d’extension de la barre par
∂u
N (x, t) = EA (x, t)
∂x
où E est le module de Young (constant) du matériau. On admettra que la formulation faible pour la dynamique
des barres est de la forme
∂u ∂w ∂2u
Z n o
∀w, trouver u tel que EA (x, t) (x) + ρA 2 (x, t)w(x) dx = F (w)
barre ∂x ∂x ∂t
où la solution u(x, t) et le champ virtuel w(x) vérifient des conditions d’admissibilité cinématique en cas de
liaisons cinématiques (non spécifiées pour l’instant), et F (w) est la forme linéaire sur w exprimant la puissance
virtuelle des efforts extérieurs appliqués (non spécifiés ici) dans le champ virtuel w(x).
1. Matrice de rigidité d’un élément de barre à 3 nœuds. On considère l’utilisation d’éléments linéiques à 3
noeuds, l’élément de référence étant le segment −1 ≤ a ≤ 1 et reposant sur les fonctions d’interpolation
1 1
N1 (a) = a(a − 1) , N2 (a) = 1 − a2 , N3 (a) = a(a + 1) (−1 ≤ a ≤ 1)
2 2
Les trois nœuds de l’élément sont supposés régulièrement espacés, deux noeuds consécutifs étant donc
séparés de h/2 (h : longueur de l’élément considéré).
h/2
a=−1 a=0 a=1
1 2 3
h
(i) Donner le jacobien dx/da de l’élément.
(ii) Exprimer la déformation ∂uh /∂x(x, t) sur l’élément en fonction de a, h et des valeurs nodales
u(1) , u(2) , u(3) de u sur l’élément.
(iii) Formuler la matrice de rigidité élémentaire [Ke ] de l’élément, et calculer explicitement [Ke ]. On
pourra pour cela utiliser les identités
Z 1 Z 1 Z 1
7 8
N1′2 (a) da = N3′2 (a) da = N2′2 (a) da =
−1 −1 6 −1 3
Z 1 Z 1 Z 1
4 1
N1′ (a)N2′ (a) da = N2′ (a)N3′ (a) da = − N1′ (a)N3′ (a) da =
−1 −1 3 −1 6
(iv) Pour quelles combinaisons de valeurs nodales {Ure } = {u(1) , u(2) , u(3) } a-t’on [Ke ]{Ue } = {0} ?
Commenter.
(v) Combien faudrait-il utiliser de points de Gauss pour une évaluation exacte de la rigidité élémentaire
[Ke ] par intégration numérique ?
(vi) Il est parfois avantageux (moyennant certaines précautions !) de mettre en œuvre des méthodes de
sous-intégration consistant à utiliser délibérément un nombre a priori insuffisant de points de Gauss
pour le calcul par intégration numérique de [Ke ].
Quelle est, dans le cas présent, la valeur approchée de [Ke ] obtenue au moyen d’un seul point de
Gauss (abscisse ag = 0, poids wg = 2) ? Montrer que [Ke ]{Ue } = {0} pour (a) {Ure } trouvé au
(iv), et (b) une autre combinaison indépendante {Ube } que l’on déterminera. Montrer (au moyen de la
valeur correcte de [Ke ]) que {Ube } a une énergie de déformation non nulle sur l’élément.
(vii) Formuler la matrice de masse élémentaire [Me ] de l’élément générique, et calculer explicitement
[Me ]. On pourra pour cela utiliser les identités
1 1 1
4 16
Z Z Z
N12 (a) da = N32 (a) da = N22 (a) da =
−1 −1 15 −1 15
1 1 1
2 1
Z Z Z
N1 (a)N2 (a) da = N2 (a)N3 (a) da = N1 (a)N3 (a) da = −
−1 −1 15 −1 15
(viii) Proposer une approximation [M̃e ] de [Me ] obtenue par condensation sur la diagonale. Quelle est
l’utilité pratique d’une telle condensation ?
2. Algorithme d’intégration en temps. On considère maintenant une barre de longueur h modélisée par un
unique élément à 3 nœuds, tel que 0 ≤ x ≤ h, et encastrée en x = 0, soit
u(0, t) = 0
(i) Donner en fonction des paramètres h, E, A, ρ du problème, pour ce modèle à un seul élément, la
forme explicite des matrices de raideur et de masse dans l’équation du mouvement
telle que {U} ne contienne que les déplacements nodaux restant inconnus après prise en compte de
l’encastrement (le mode de sollicitation n’étant pas précisé, on ne cherchera pas à spécifier les forces
nodales {F}).
(ii) On introduit une discrétisation temporelle de pas de temps constant ∆t, à laquelle on applique le
schéma de Newmark avec β = 1/6, γ = 1/2. Ce schéma, correspondant à une précision optimale,
est conditionnellement stable (cf chapitre 9), le pas de temps devant vérifier ∆t ≤ ∆tstab . Définir
précisément une procédure permettant d’évaluer la limite de stabilité ∆tstab de ce schéma de New-
mark sur cet exemple. Donner la valeur de ∆tstab en fonction des divers paramètres du problème.
Exercice 2 : échantillon élastoplastique sous sollicitation antiplane
On considère un échantillon de section rectangulaire S = {−A ≤ x ≤ A, 0 ≤ y ≤ B} dans le plan Oxy et de
longueur infinie dans la direction Oz (Oxyz désignant un système de coordonnées cartésiennes orthonormées).
L’unique sollicitation appliquée à cet échantillon consiste en une densité d’effort T D = T (x)ez imposée sur sa
face y = B de façon à créer un état de déformation antiplane, c’est-à-dire correspondant à un déplacement de la
forme
u(x) = u(x, y)ez ,
les seules composantes non nulles de déformation et de contrainte étant alors εxz , εyz et σxz , σyz .
T D= T(x) ez y a2
4 3
A A
a1
B (S)
x 1 2
L’objet de cet exercice est d’examiner la recherche de solutions approchées en régime élastique puis élasto-
plastique, en s’appuyant sur la modélisation de la section S par un unique élément fini quadrangulaire à 4 nœuds
(cadre de travail retenu pour la totalité de cet exercice). Les fonctions de forme pour cet élément sont données
au tableau 2.1 du chapitre 2 du livre ; elles sont utilisées pour l’interpolation de la géométrie et du déplacement
(élément isoparamétrique). On respectera la numérotation des noeuds définie par la figure.
1. Interpolation géométrique. Exprimer un point générique x de l’élément S en fonction des coordonnées pa-
ramétriques a1 , a2 . Donner les expressions de la matrice jacobienne J et du déterminant jacobien J(a1 , a2 )
de l’élément.
2. Sollicitation. Donner la valeur du vecteur élémentaire {Fe } des forces généralisées associées au chargement
décrit plus haut, en fonction de T (x). Expliciter complètement {Fe } pour le cas d’une densité d’effort
T (x) = T uniforme.
3. Régime élastique. Le matériau est dans cette partie supposé purement élastique (propriétés isotropes ca-
ractérisées par les constantes de Lamé λ, µ).
(a) Exprimer les déformations εxz , εyz sur l’élément, pour des valeurs quelconques des déplacements no-
daux, en fonction des coordonnées paramétriques a1 , a2 .
(b) Evaluer la matrice de rigidité élémentaire [Ke ] de l’élément.
(c) On complète la condition précédemment introduite de densité d’effort uniforme T D = T ez imposée
sur sa face y = B par les conditions aux limites suivantes : blocage en déplacement de la face y = 0, et
faces x = ±A libres de contraintes. Donner le système d’équations vérifié par les déplacements nodaux
inconnus. Calculer la solution élastique approchée.
4. Régime élastoplastique. On suppose maintenant que le matériau constitutif est élastique (λ, µ) et plastique
avec écrouissage cinématique linéaire (limite d’élasticité en traction σ0 , module d’écrouissage h, règle de
normalité). La limite d’élasticité est définie par le critère de von Mises, qui compte tenu des hypothèses
faites ici s’écrit
σ eq − R(p) ≤ 0, R(p) = σ0 + hp,
où p désigne la déformation plastique cumulée. Le chargement appliqué et les conditions aux limites sont
définis de la même manière que précédemment.
(a) Déterminer la valeur T0 de la densité uniforme d’effort correspondant à l’apparition de la plasticité
dans l’échantillon.
(b) On considère, à partir du niveau de chargement T = T0 , un pas de chargement en régime élastoplastique
tel que le chargement (uniforme) imposé en fin de pas soit T0 + ∆T (avec ∆T > 0).
– Donner le champ de contraintes sur l’élément prédit par l’algorithme de retour radial pour un
incrément de déplacement dont les valeurs nodales sont {0, 0, ∆u, ∆u}. Etablir l’équation vérifiée
par le scalaire ∆u issue de la formulation faible de l’équilibre.
– Donner le module tangent cohérent associé à la méthode de Newton pour la résolution en ∆u de
l’équation d’équilibre précédemment établie.
– Donner la valeur de ∆u en fonction de ∆T et des paramètres de géométrie et de comportement.
C ORRIG É
(viii) On peut formuler une approximation diagonale [M̃e ] de [Me ] en affectant à chaque coefficient diago-
nal la somme des coefficients de la ligne correspondante, ce qui donne ici :
1 0 0
ρAh
[M̃e ] = 0 4 0
6
0 0 1
Une telle condensation conduit à une approximation diagonale de la matrice de masse globale [M],
ce qui permet d’accélérer les schémas « explicites » tels que celui des différences centrées en évitant
la résolution d’un système d’équations de matrice [M].
2. Algorithme d’intégration en temps.
(i) Le déplacement nodal u(1) (t) étant imposé à zéro, on incorpore cette liaison dans les équations
discrètes, et on ne conserve que celles associées aux champs virtuels de valeur nodale w(1) nulle.
Cela revient à écrire (discrétisation avec un seul élément) l’équation [Ke ]{Ue } + [Me ]{Üe } = {Fe }
en supprimant la première ligne de l’équation et la première colonne des matrices, soit :
EA 16 −8 ρAh 8 1
[K]{U} + [M]{Ü} = {F} avec [K] = , [M] =
3h −8 7 15 1 2
(ii) Pour le schéma proposé (Newmark avec β = 1/6, γ = 1/2), la limite de stabilité ∆tstab est donnée
(chapitre 9) par √
1 2 2 3
∆tstab = min √ =
J ωJ 2γ − 4β maxJ ωJ
où ω1 , ω2 sont les racines de l’équation caractéristique
Det [K] − ω 2 [M] = 0
p √
2
Comptep tenu des expressions de [K], [M] données en (i) et posant ω = ω̄ 5E/(ρh ) = ω̄ 5 c/h
(c = E/ρ étant la célérité des ondes de compression dans une barre élastique), ω̄ est solution de
16 −8 2 8 1
Det − ω̄ = 0 soit 15ω̄ 4 − 104ω̄ 2 + 48 = 0
−8 7 1 2
dont les racines sont √
2 52 ± 8 31
ω̄ =
15
La limite de stabilité ∆tstab est donc donnée, tous calculs faits, par
s √
13 − 2 31
c∆tstab = h ≈ 0, 61h
5
Elle correspond donc au temps que met une onde de compression pour parcourir (environ) 0, 61
longueur d’élément.
Exercice 2 : échantillon élastoplastique sous sollicitation antiplane
1. Interpolation géométrique. Les fonctions de forme de l’élément sont données (chapitre 2) par
1 1
N1 (a) = (1 − a1 )(1 − a2 ), N2 (a) = (1 + a1 )(1 − a2 ),
4 4
1 1
N3 (a) = (1 + a1 )(1 + a2 ), N4 (a) = (1 − a1 )(1 + a2 )
4 4
La position d’un point x de l’élément est donnée par interpolation des positions nodales. En posant l’inter-
polation et effectuant les calculs, on trouve :
x −A A A −A a1 0 A
= N1 (a) + N2 (a) + N3 (a) + N4 (a) =
y 0 0 B B 0 (1 + a2 )/2 B
La matrice jacobienne J et le déterminant jacobien J(a1 , a2 ) de l’élément sont donc constants (cela résulte
de la forme rectangulaire de l’élément), et donnés par :
A 0 1
J= , J(a1 , a2 ) = AB
0 B/2 2
2. Sollicitation. Le vecteur élémentaire {Fe } des forces généralisées est défini ici à travers la relation
Z
{We }T {Fe } = w(x)T (x) dsx avec w(x) = [Ne (a)]{We }
bord supérieur y = B
ce qui donne (en remarquant ici que dsx = A da1 sur le bord supérieur, lui-même défini par a2 = 1), en
exploitant les expressions des fonctions de forme :
N1 (a1 , 1) 0
Z 1 Z 1
N2 (a1 , 1) 0
{Fe } = T (x(a1 , 1))A da1 = A T (x(a1 , 1)) da1
−1
N3 (a1 , 1)
−1
(1 + a1 )/2
N4 (a1 , 1) (1 − a1 )/2
Pour le cas d’une densité d’effort T (x) = T uniforme, on effectue les intégrations et trouve :
{Fe } = {0 0 AT AT }T (1)
3. Régime élastique.
(a) Compte tenu du caractère diagonal et constant de la matrice jacobienne J, on a
∂ ∂ ∂x ∂ ∂ ∂ ∂y B ∂
= =A , = =
∂a1 ∂x ∂a1 ∂x ∂a2 ∂y ∂a2 2 ∂y
Les déformations associées à un champ de déplacement u(x) = [Ne (a)]{Ue } défini par interpolation
de déplacements nodaux sur l’élément sont donc données par :
∂u 1 ∂u 1 h ∂Ne i ∂u 2 ∂u 2 h ∂Ne i
2εxz = = = {Ue }, 2εyz = = = {Ue } (2)
∂x A ∂a1 A ∂a1 ∂y B ∂a2 B ∂a2
les matrices des dérivées des fonctions de forme étant données, après calcul, par
h ∂N i
e
= − 41 ((1 − a2 ) 1 1
− 41 ((1 + a2 )
∂a1 4 (1 − a2 ) 4 (1 + a2 )
h ∂N i
e
= − 41 ((1 − a1 ) − 14 (1 + a1 ) 1
− 41 ((1 − a1 )
∂a2 4 (1 + a1 )
(b) La matrice de rigidité élémentaire [Ke ] de l’élément est définie à travers la relation
Z Z
T
{We } [Ke ]{Ue } = σ[u] : ε[w] dV = 2 σxz [u]εxz [w] + σyz [u]εyz [w] dV
élément élément
Z
= 4µ εxz [u]εxz [w] + εyz [u]εyz [w] dV
élément
En reportant les expressions (2) pour u et w, on trouve par identification (en remarquant que dV =
J da1 da2 = (AB/2) da1 da2 )
(c) Les conditions de blocage conduisent, pour ce modèle à un seul élément, à écrire l’équation [Ke ]{Ue } =
{Fe } dans laquelle les lignes 1 et 2, ainsi que les colonnes 1 et 2 de [Ke ], sont supprimées. Compte
tenu des expressions (1) et (3), l’équilibre de l’échantillon est gouverné par l’équation matricielle d’in-
connues u(3) , u(4)
µB 1 −1 µA 2 1 u(3) T A
+ =
6A −1 1 3B 1 2 u(4) TA
La solution élastique approchée (qui est en fait exacte pour T uniforme) est donnée par résolution du
système ci-dessus (on cherche directement u(3) = u(4) par symétrie) :
TB
u(3) = u(4) = (4)
µ
4. Régime élastoplastique. Il est utile de noter pour la suite que les déformations et contraintes sont, dans le
cadre antiplan considéré ici, toujours égales à leurs parties déviatoriques.
(a) La solution élastique approchée définie par (4) est telle que
T
εxz = 0, εyz = , σxz = 0, σyz = T
2µ
La contrainte équivalente pour cette solution est donnée par
√
σ eq = 3|T |
elas ∆u
σxz = selas
xz = 0,
elas
σyz = selas
yz = T0 + 2µ∆εyz = T0 + µ
B
√ √
La contrainte équivalente associée σ elas,eq valant 3(T0 + µ∆u/B) > 3T0 (pour ∆u > 0), le
prédicteur élastique n’est pas plastiquement admissible et l’incrément de contrainte est élasto-plastique.
L’algorithme de retour radial (chapitre 6) se traduit ici par les étapes suivantes :
– Détermination de l’incrément de déformation plastique cumulée ∆p par résolution de l’équation de
cohérence : √
elas,eq 3µ∆u
σ − 3µ∆p − σ0 + h∆p = 0 soit ∆p =
(3µ + h)B
– Détermination de l’incrément de déformation plastique : on trouve
3∆p elas 3µ∆u
∆εPyz = elas,eq
syz =
2σ 2B(3µ + h)
L’inconnue principale étant le scalaire ∆u, il suffit de construire une équation scalaire à partir de cette
formulation faible. Prenant w(3) = w(4) , remarquant que σyz [u + ∆u] donné par (5) et εyz [w] sont
constants sur S, on obtient l’équation scalaire sur ∆u :
w(4) ∆u hµ
R(∆u; w) = 2AB T0 + − 2Aw(4) (T0 + ∆T ) = 0 ∀w(4)
B B 3µ + h
qui donne après simplification :
hµ
R(∆u; w) = 2Aw(4) ∆u − ∆T = 0
3µ + h
Cette équation est linéaire en ∆u. Son module tangent global Aep est donc constant (et scalaire), il est
défini par
dR (4) 2Ahµ
w = Aep w(4) , soit Aep =
d∆u 3µ + h
La résolution de l’équation d’équilibre donne
3µ + h
Aep ∆u = 2A∆T, soit ∆u = ∆T
hµ