Élaboration d'Éléments Coques Volumiques
Élaboration d'Éléments Coques Volumiques
THÈSE
pour obtenir le grade de
Docteur
de
l’École Nationale Supérieure d'Arts et Métiers
Spécialité “Mécanique ”
présentée et soutenue publiquement
par
Vuong-Dieu TRINH
le 20 avril 2009
Jury :
Alain COMBESCURE, Professeur, LaMCoS, INSA de Lyon ................................. Examinateur
Bruno COCHELIN, Professeur, LAM, Ecole Centrale de Marseille......................... Rapporteur
Pierre VILLON, Professeur, Laboratoire Roberval, UTC.......................................... Rapporteur
Jean-Philippe PONTHOT, Professeur, LTAS, Université de Liège Examinateur
Farid ABED-MERAIM, Maître de Conférences, LPMM, ENSAM de Metz............ Examinateur
Patrick MASSIN, Directeur du LaMSID, UMR EDF-CNRS 2832 ........................... Examinateur
Jean-François BILLAUD, Ingénieur, CETIM ........................................................... Examinateur
X Xavier DESROCHES, Docteur, LaMSID, UMR EDF-CNRS 2832.......................... Invité
Arts et Métiers ParisTech (Ecole Nationale Supérieure d’Arts et Métiers) est un Grand Etablissement
dépendant du Ministère de l’Enseignement Supérieur et de la Recherche, composé de huit centres :
AIX-EN-PROVENCE ANGERS BORDEAUX CHÂLONS-EN-CHAMPAGNE CLUNY LILLE METZ PARIS
Table des matières
TABLE DES MATIÈRES .......................................................................................................1
TABLE DES FIGURES ...........................................................................................................3
LISTE DES TABLEAUX ........................................................................................................4
INTRODUCTION ....................................................................................................................6
CONTEXTE ...............................................................................................................................7
MOTIVATION DE LA THÈSE ......................................................................................................8
OBJECTIFS ET CONTENU DE LA THÈSE ......................................................................................9
CHAPITRE I ÉTUDE BIBLIOGRAPHIQUE ...................................................................11
1. ÉTUDES DE LA MMC ET DE LA MEF ................................................................................12
1.1. Cinématique des milieux continus ........................................................................12
1.2. Contraintes, équations d’équilibre et déformations .............................................13
1.3. Relation de comportement ....................................................................................15
1.4. Principe des travaux virtuels................................................................................16
1.5. Principe variationnel............................................................................................17
1.6. Discrétisation par éléments finis ..........................................................................18
2. MODÉLISATIONS DE COQUES EXISTANTES .....................................................................21
2.1. Modélisation de coque tridimensionnelle dégénérée ...........................................21
2.2. Modélisations de coque en formulation mixte......................................................23
2.3. Modélisations solide–coques................................................................................27
CHAPITRE II ÉLÉMENTS COQUES VOLUMIQUES SHB EN LINÉAIRE.............30
1. ÉLÉMENT COQUE VOLUMIQUE SHB6.............................................................................31
1.1. Modélisation SHB6 en linéaire.............................................................................31
1.1.1. Cinématique et interpolation ............................................................................32
1.1.2. Opérateur gradient discrétisé ............................................................................32
1.1.3. Formulation variationnelle utilisée pour l’élément SHB6................................37
1.1.4. Analyse des modes de « hourglass » pour l’élément SHB6 ............................40
1.1.5. Projection par “Assumed local strain method”.................................................43
1.1.6. Matrice de rigidité géométrique K σ .................................................................45
1.1.7. Forces suiveuses et matrice de pression K P .....................................................46
1.2. Validation de l’élément coque volumique SHB6 en linéaire................................48
1.2.1. Poutre en flexion simple ...................................................................................48
1.2.2. Plaque en flexion et cisaillement dans son plan ...............................................50
1.2.3. Poutre vrillée soumise à un effort tranchant.....................................................53
1.2.4. Coque sphérique pincée....................................................................................55
1.2.5. Coque sphérique pincée avec mélange d’éléments ..........................................56
1.2.6. Coque cylindrique pincée avec diaphragmes ...................................................59
1.2.7. Plaque circulaire soumise à une force ponctuelle.............................................61
1.2.8. Étude fréquentielle d’une poutre libre encastrée ..............................................63
1.2.9. Flambement d’un cylindre libre sous pression externe ....................................65
2. ÉLÉMENTS COQUES VOLUMIQUES SHB15 ET SHB20 ....................................................68
2.1. Modélisation SHB15 en linéaire...........................................................................68
2.1.1. Cinématique et interpolation ............................................................................69
2.1.2. Opérateur gradient discrétisé ............................................................................69
2.2. Modélisation SHB20 en linéaire...........................................................................80
2.2.1. Cinématique et interpolation ............................................................................80
2.2.2. Opérateur gradient discrétisé ............................................................................81
2.2.3. Formulation variationnelle utilisée pour les éléments SHB15 et SHB20 ........92
2.2.4. Matrice de rigidité géométrique K σ .................................................................95
2.2.5. Forces suiveuses et matrice de pression K P .....................................................97
2.3. Validation des éléments coques volumiques SHB15 et SHB20 en linéaire........100
2.3.1. Poutre en flexion simple .................................................................................100
2.3.2. Plaque en flexion et cisaillement dans son plan .............................................101
2.3.3. Poutre vrillée soumise à un effort tranchant...................................................102
2.3.4. Coque sphérique pincée..................................................................................103
2.3.5. Coque cylindrique pincée avec diaphragmes .................................................104
2.3.6. Plaque circulaire soumise à une force ponctuelle...........................................106
2.3.7. Étude fréquentielle d’une poutre libre encastrée ............................................107
2.3.8. Flambement d’un cylindre libre sous pression externe ..................................107
2.3.9. Poutre en flexion avec divers élancements.....................................................109
2.3.10. Flambage d’une coque cylindrique avec raidisseur........................................109
2
Table des figures
Figure 1. Schéma de fonctionnement d’une centrale nucléaire..................................................7
Figure 2. Exemple de maillage 3D et centrales nucléaires.........................................................8
Figure 3. Cinématique des milieux continus ............................................................................12
Figure 4. Facette de normale ....................................................................................................14
Figure 5. Modèles coques 3D d’Ahmad...................................................................................21
Figure 6. Modèles coques 3D de Kim Y.H. .............................................................................22
Figure 7. Hypothèse cinématique de Reissner–Mindlin...........................................................23
Figure 8. Élément de référence de coque triangulaire à 3 nœuds.............................................24
Figure 9. Élément de référence SHB6 et ses points d’intégration............................................31
Figure 10. Géométrie, chargement et conditions aux limites pour le test de la poutre en flexion
simple ; un exemple de maillage (12x2x1)x2...........................................................................49
Figure 11. Géométrie, chargement et déformée de la plaque en flexion et cisaillement dans
son plan ; un exemple de maillage (12x4x1)x2........................................................................50
Figure 12. Poutre vrillée soumise à un effort tranchant ; un exemple de maillage (12x4x1)x2
..................................................................................................................................................54
Figure 13. Hémisphère pincé ; un exemple de maillage 3x(5x5x1)x2.....................................56
Figure 14. Hémisphère pincé (SHB6 à un seul sommet) ; un exemple de maillage 4x(7x7x1)
..................................................................................................................................................57
Figure 15. Maillage mixte de l’hémisphère pincé (SHB6 aux 3 sommets)..............................58
Figure 16. Géométrie, chargement et déformée de la coque cylindrique pincée avec
diaphragmes..............................................................................................................................60
Figure 17. Géométrie, chargement et déformée de la plaque circulaire soumise à une force
ponctuelle ; un exemple de maillage 3x(4x4x1)x2...................................................................61
Figure 18. Géométrie, chargement et modes propres de la poutre libre encastrée ; un exemple
de maillage (30x3x1)x2 ............................................................................................................65
Figure 19. Géométrie, chargement et conditions aux limites du cylindre sous pression
externe ; un exemple de maillage (20x30x1)x2........................................................................66
Figure 20. Courbe de convergence et modes de flambement du cylindre sous pression
externe ; un exemple de maillage (20x30x1)x2........................................................................67
Figure 21. Géométrie de l’élément de référence SHB15 et ses points d’intégration ...............68
Figure 22. Géométrie de l’élément de référence SHB20 et ses points d’intégration ...............80
Figure 23. Géométrie, chargement et conditions aux limites pour le test de la poutre en flexion
avec divers élancements ; un exemple de maillage (10x1x1) SHB20....................................110
Figure 24. Géométrie, conditions aux limites et 1er mode de flambage d’un quart de coque
cylindrique avec raidisseur ; un exemple de maillage mixte de 260 éléments SHB8PS pour le
raidisseur et 360 éléments SHB6 pour la coque principale ....................................................112
Figure 25. Géométrie initiale et 1er mode de flambage de la coque cylindrique avec
raidisseur ; un exemple de maillage de : a) 300 éléments SHB15 ; b) 100 éléments SHB20 113
Figure 26. Géométrie, chargement et déplacement de la poutre console en non-linéaire
géométrique ............................................................................................................................121
Figure 27. Géométrie, chargement, conditions aux limites et déformée du panneau cylindrique
épais sous force ponctuelle ; un exemple de maillage SHB6 de (30x30x2)x2.......................122
Figure 28. Géométrie, chargement, conditions aux limites et déformée du panneau cylindrique
mince sous force ponctuelle ; un exemple de maillage SHB6 de (25x25x2)x2 .....................128
Figure 29. Géométrie, chargement, conditions aux limites et déformée de la coque cylindrique
sous force ponctuelle ; un exemple de maillage SHB20 de (20x20x1)..................................131
3
Liste des tableaux
Tableau 1. Données géométriques, matériau et de chargement pour le test de la poutre en
flexion simple ...........................................................................................................................48
Tableau 2. Déplacement du point A suivant Oz de la poutre en flexion simple ......................49
Tableau 3. Données de géométrie, matériau et de chargement du cas test de la plaque en
flexion et cisaillement dans son plan........................................................................................50
Tableau 4. Résultats du déplacement du point A suivant Oy de la plaque en flexion et
cisaillement dans son plan ........................................................................................................53
Tableau 5. Données de géométrie, chargement et de matériau de la poutre vrillée .................54
Tableau 6. Déplacement du point A suivant Oz de la poutre vrillée........................................54
Tableau 7. Données de géométrie, matériau et de chargement de l’hémisphère pincé............55
Tableau 8. Déplacement du nœud A suivant Ox de l’hémisphère pincé..................................56
Tableau 9. Déplacement du point A suivant Ox, maillage mixte SHB8PS et SHB6 en un seul
sommet .....................................................................................................................................57
Tableau 10. Déplacement du point A suivant Ox maillage mixte SHB8PS et SHB6 aux trois
sommets de l’hémisphère pincé................................................................................................58
Tableau 11. Données de géométrie, chargement et de matériau du test de la coque cylindrique
pincée avec diaphragmes ..........................................................................................................59
Tableau 12. Déplacement du point A suivant Oz de la coque cylindrique pincée avec
diaphragmes..............................................................................................................................60
Tableau 13. Données de géométrie, chargement et de matériau de la plaque circulaire
soumise à un effort ponctuel.....................................................................................................62
Tableau 14. Déplacement du point A suivant Oz de la plaque circulaire soumise à un effort
ponctuel ....................................................................................................................................62
Tableau 15. Données géométriques et matériau de la poutre libre encastrée...........................63
Tableau 16. Fréquences propres de la poutre libre encastrée ...................................................64
Tableau 17. Données géométriques et matériau du cas test de flambement d’un cylindre libre
sous pression externe ................................................................................................................65
Tableau 18. Pressions critiques de flambage du cylindre sous pression externe .....................66
Tableau 19. Déplacement du point A suivant Oz de la poutre en flexion simple ..................100
Tableau 20. Déplacement du point A suivant Oy de la plaque en flexion et cisaillement dans
son plan...................................................................................................................................101
Tableau 21. Déplacement suivant Oz du point A de la poutre vrillée....................................102
Tableau 22. Déplacement du nœud A suivant Ox de l’hémisphère pincé..............................103
Tableau 23. Déplacement du point A suivant Oz du cylindre pincé avec diaphragmes ........105
Tableau 24. Déplacement du point A suivant Oz de la plaque circulaire soumise à un effort
ponctuel ..................................................................................................................................106
Tableau 25. Pression critique de flambage du cylindre libre sous pression externe ..............108
Tableau 26. Données géométriques, matériau et de chargement pour le test de la poutre en
flexion simple avec divers élancements .................................................................................109
Tableau 27. Déplacement du point A suivant Oz de la poutre en flexion simple avec divers
élancements ............................................................................................................................110
Tableau 28. Résultats du test de flambement de la coque cylindrique avec raidisseur..........112
Tableau 29. Paramètres géométriques et matériau de la coque cylindrique avec raidisseur..112
Tableau 30. Données géométriques et matériau de la poutre console....................................119
Tableau 31. Solution de référence pour le cas test de la poutre console en non-linéaire
géométrique ............................................................................................................................120
4
Tableau 32. Résultats obtenus pour le cas test de la poutre console en non-linéaire
géométrique ............................................................................................................................120
Tableau 33. Données de géométrie et de matériau du panneau cylindrique épais sous force
ponctuelle ...............................................................................................................................122
Tableau 34. Résultats obtenus du déplacement du point A suivant Oz du panneau cylindrique
épais soumis à une force ponctuelle .......................................................................................123
Tableau 35. Données de matériau élasto-plastique du panneau cylindrique sous force
ponctuelle ...............................................................................................................................124
Tableau 36. Données de matériau élasto-plastique du panneau cylindrique épais sous force
ponctuelle ...............................................................................................................................124
Tableau 37. Résultats obtenus du déplacement du point A suivant Oz du panneau cylindrique
épais soumis à une force ponctuelle .......................................................................................126
Tableau 38. Données de géométrie et de matériau du panneau cylindrique mince sous force
ponctuelle ...............................................................................................................................127
Tableau 39. Résultats de référence .........................................................................................129
Tableau 40. Données de géométrie et de matériau de la coque cylindrique sous force
ponctuelle ...............................................................................................................................130
Tableau 41. Résultats obtenus des déplacements des points A, B et C de la coque cylindrique
soumise à une force ponctuelle...............................................................................................132
5
Remerciements
Ce travail a été réalisé dans le cadre d’une collaboration entre le Laboratoire de
Physique et Mécanique des Matériaux (LPMM) de l’ENSAM CER de Metz, le Laboratoire de
Mécanique des Contacts et des Solides (LaMCoS) de l’INSA de Lyon et le Laboratoire de
Mécanique des Structures Industrielles Durables (LaMSID), dans le cadre d’un partenariat
avec EDF R&D et le CETIM. Je remercie Monsieur le Professeur El Mostafa DAYA,
directeur du LPMM et Monsieur Patrick MASSIN, directeur du LaMSID, pour leurs accueils
chaleureux et pour avoir accepté que j’effectue mes travaux de recherche au sein de leurs
laboratoires dans de bonnes conditions.
Messieurs les Professeurs Bruno COCHELIN et Pierre VILLON m’ont fait le plaisir et
l’honneur d’être rapporteurs de ma thèse, je les remercie vivement pour l’intérêt qu’ils ont
porté à mon travail et pour avoir accepté la lourde tâche que comporte le travail de rapporteur.
Jean-Michel PROIX et Jean-François BILLAUD ont été des partenaires précieux, voire
indispensables, tant d’un point de vue technique qu’humain, je les remercie sincèrement car
ils ont grandement participé à la réalisation de ce travail.
Je tiens à remercier tous mes collègues du LPMM et du LaMSID qui m’ont aidé en
répondant à mes questions. Que ma famille et tous mes amis soient assurés de mon immense
gratitude et de ma sincère reconnaissance pour leur soutien permanent.
Il me faudrait des siècles entiers pour remercier mon directeur de thèse Monsieur le
Professeur Alain COMBESCURE, car celui qui m’apprend une lettre, je serai pour lui un
esclave. Quant à lui, il était pour moi le bon directeur qui m’a imprégné la mécanique des
structures, et le vrai Professeur reconnu par sa modestie, et ses idées claires, pertinentes et
encourageantes qui ont fait de moi un jeune ingénieur chercheur ambitieux prêt à conquérir le
monde.
Enfin, il est sûr que j’oublie certaines personnes… qu’ils m’en excusent. Trois années
m’ont permis de rencontrer beaucoup de personnes, qui ont toutes eu un rôle dans ma vie et
par conséquent dans la construction de ce travail. Ils se reconnaîtront.
Introduction
Contexte
Le groupe EDF est un leader européen de l’énergie, présent sur tous les métiers de
l’électricité, de la production au négoce et de plus en plus actif sur la chaîne du gaz en
Europe. Acteur principal du marché français de l’électricité nucléaire, il est solidement
implanté en Grande-Bretagne, en Allemagne et en Italie.
7
Introduction
Motivation de la thèse
Toutefois, l’élément hexaédrique SHB8PS ne permet pas de mailler des géométries de
formes complexes quelconques. Le développement d’un élément similaire mais de géométrie
prismatique est donc nécessaire : SHB6.
Plusieurs travaux de recherche sur l’élément de coque volumique SHB6 à six nœuds ont
été réalisés et implantés dans le code de calcul INCA. L’élément SHB6 est également sous-
intégré. Outre la réduction considérable des temps de calculs, la méthode de sous-intégration
permet de réduire nombre des différents verrouillages rencontrés dans la mise en œuvre
numérique des éléments finis. Cependant, cette sous-intégration n’a pas que des avantages :
elle introduit malheureusement des modes parasites associés à une énergie nulle. En statique,
ceci peut conduire à une singularité de la matrice de raideur globale pour certaines conditions
aux limites. En dynamique transitoire, en revanche, cela conduit à des modes en sablier
« hourglass » qui vont déformer le maillage de façon irréaliste et qui finissent par faire
8
Introduction
Les premières études menées sur le SHB6 (Abed-Meraim et al.) ont prouvé que cet
élément ne présentait pas de modes de hourglass, mais après implantation, elles ont aussi
montré que celui-ci présentait un sévère blocage numérique, notamment dans les problèmes
dominés par la flexion où un cisaillement parasite induit un verrouillage en cisaillement. La
méthode « Assumed strain » a ensuite été utilisée pour éliminer certains verrouillages de cet
élément SHB6. Le principe de cette méthode, largement utilisée dans la littérature, consiste à
projeter l’opérateur gradient discrétisé B sur un sous-espace approprié afin d’éviter les
différents problèmes liés au verrouillage. Les auteurs ont réalisé différentes projections pour
trouver celle qui élimine le maximum de verrouillages. Leurs travaux ont apporté quelques
améliorations à l’élément SHB6 qui a montré une assez bonne convergence dans plusieurs
cas-tests adaptés aux éléments de coques. Cependant, certains verrouillages sévères de type
cisaillement ou membrane persistaient encore dans certaines situations montrant que la
formulation de l’élément SHB6 pouvait encore être améliorée du point de vue du verrouillage.
Motivés par ces premiers résultats encourageants, nous poursuivons l’effort mené sur
l’élément SHB6 pour aboutir à une version qui souffre le moins de verrouillage.
Le premier chapitre fait état d’une revue bibliographique qui apporte des éléments
théoriques sur les différentes méthodes de développement des éléments finis coques
volumiques réalisées par des recherches durant ces dernières décennies.
9
Introduction
Par ailleurs, l’élément SHB8PS actuel a été couplé à seulement certaines lois de
comportement telles que élastique ou élasto-plastique avec écrouissage isotrope de type von
Mises. Le second objectif de cette étude est d’élargir le champ d’application de l’élément
SHB8PS ainsi que les autres éléments finis solide–coques SHB6, SHB15 et SHB20 à d’autres
lois de comportement du code ASTER. Le troisième chapitre présente le principe théorique de
ce couplage ainsi que des cas tests pour valider cette approche.
10
Chapitre I
Étude bibliographique
JG Mt
MO U
Ωt
JJJG ΩO
XO
JJG
Xt
z
y
x O
12
Etudes bibliographiques
⎡ x⎤ ⎡ xo ⎤ ⎡ u ⎤
x = ⎢ y ⎥ = xo + u = ⎢⎢ yo ⎥⎥ + ⎢⎢ v ⎥⎥
⎢ ⎥ (2)
⎣⎢ z ⎦⎥ ⎣⎢ zo ⎦⎥ ⎣⎢ w⎦⎥
Nous avons :
⎡ ∂x ∂x ∂x ⎤
⎢ ⎥
⎢ ∂x0 ∂y0 ∂z0 ⎥
⎡ dx0 ⎤ ⎡1 + u, x0 u, y0 u, z0 ⎤ ⎡ dx0 ⎤
⎡ ∂x ⎤ ⎢ ∂y ∂y ∂y ⎥ ⎢ ⎥ ⎢ ⎥
dx = ⎢ ⎥ ⋅ dx o = ⎢ ⎥ ⋅ dy0 = ⎢ v, x0 1 + v, y0 v, z0 ⎥ ⋅ ⎢⎢ dy0 ⎥⎥ = F ⋅ dxo (3)
⎣ ∂xo ⎦ ⎢ ∂x0 ∂y0 ∂z0 ⎥ ⎢ ⎥ ⎢ ⎥
⎢⎣ dz0 ⎥⎦ ⎢ w, x w, y0 1 + w, z0 ⎦⎥ ⎢⎣ dz0 ⎥⎦
⎢ ∂z ∂z ∂z ⎥ ⎣ 0
⎢ ⎥
⎢⎣ ∂x0 ∂y0 ∂z0 ⎥⎦
F = I + L0 ; L 0 = D0 + W0
⎡
⎡ u, x0
⎢ u, x0
1
2
(
u, y0 + v, x0 ) 1
(
u, z + w, x0
2 0
)⎤⎥
u, y0 u, z0 ⎤ ⎢ ⎥
⎢ ⎥ 1 1
L0 = ⎢ v, x0 v, y0 v, z0 ⎥ ; D0 = ⎡⎣ L0 + L 0 ⎤⎦ = ⎢
2
T
⎢
v, y0 (
v, z + w, y0
2 0
) ⎥;
⎥
⎢ ⎥ ⎢ sym ⎥
⎢⎣ w, x0 w, y0 w, z0 ⎥⎦ w, z0
⎢ ⎥
⎢⎣ ⎥⎦
⎡ 0 −θ z θ y ⎤
1 ⎢ ⎥
W0 = ⎡⎣L 0 − L0 ⎤⎦ = ⎢ θ z
T
0 −θ x ⎥
2
⎢ 0 ⎥⎦
⎣ −θ y θ x
1 1 1
( ) ( )
avec θ x = w, y0 − v, z0 ;θ y = u, z0 − w, x0 ;θ z = v, x0 − u, y0
2 2 2
( )
où θ x , θ y ,θ z représentent des rotations infinitésimales autour des axes x, y, z , ne produisant
aucune déformation. D0 est le tenseur des déformations linéarisées, ou petites déformations.
13
Etudes bibliographiques
∆S n
Ωt M
∆f
z
y
x O
⎛ ∆f ⎞
σ ( M , n ) = lim ⎜ ⎟ (4)
∆S →0 ∆S
⎝ ⎠
σ ( i ) = σ xx i + σ yx j + σ zx k
σ ( j) = σ xy i + σ yy j + σ zy k (5)
σ ( k ) = σ xz i + σ yz j + σ zz k
⎡σ xx σ xy σ xz ⎤
⎢ ⎥
σ = ⎢σ yx σ yy σ yz ⎥ (6)
⎢σ zx σ zy σ zz ⎥⎦
⎣
L’équilibre des moments autour des axes passant par M, en l’absence de couples répartis
à l’intérieur et à la surface du solide conduit à :
σ xy = σ yx ; σ xz = σ zx ; σ yz = σ zy (7)
σ = ⎡⎣σ xx σ yy σ zz σ xy σ xz σ zy ⎤⎦ (8)
Le solide dans la configuration actuelle Ωt est soumis à des sollicitations comme des
forces surfaciques f s appliquées sur une partie de la frontière ∂1Ωt , des déplacements
14
Etudes bibliographiques
imposés u d appliqués sur une partie de la frontière ∂ 2 Ωt , et des forces volumiques fv (qui
peuvent contenir des termes d’inertie). La somme des parties de la frontière ∂1Ωt et ∂ 2 Ωt
représente le frontière totale fermée ∂Ωt de Ωt . L’équilibre du système s’écrit de la façon
suivante :
⎧ div ( σ ) + fv = 0 ∀M ∈ Ω t
⎪
⎨u ( M ) = u d ∀M ∈ ∂ 2 Ω t
⎪σ (M ) ⋅n = f ∀M ∈ ∂1Ω t (9)
⎩ s
⎡∂ ∂ ∂⎤
où div = ⎢
⎣ ∂x ∂y ∂z ⎥⎦
Dans notre étude des structures minces, subissant des transformations élastiques
caractérisées par de grands déplacements et de petites déformations, on utilise la mesure des
déformations de Green–Lagrange linéarisé :
⎡ 1 1 ⎤
⎢ u, x 2
( u, y + v, x ) 2
( u, z + w, x ) ⎥
⎢ ⎥ ⎡ ε xx ε xy ε xz ⎤
1
ε=⎢ v, y ( v, z + w, y ) ⎥ = ⎢⎢ ε yy ε yz ⎥⎥ (10)
⎢ 2 ⎥
⎢ sym ⎥ ⎢⎣ sym ε zz ⎥⎦
⎢ w, z ⎥
⎢⎣ ⎥⎦
ou encore :
⎡ ε xx ⎤ ⎡ε xx ⎤ ⎡ u, x ⎤
⎢ ε ⎥ ⎢ε ⎥ ⎢ v ⎥
⎢ yy ⎥ ⎢ yy ⎥ ⎢ , y ⎥
⎢ ε zz ⎥ ⎢ ε zz ⎥ ⎢ w, z ⎥
ε=⎢ ⎥=⎢ ⎥=⎢ ⎥ (11)
⎢ 2ε xy ⎥ ⎢γ xy ⎥ ⎢ u, y + v, x ⎥
⎢ 2ε xz ⎥ ⎢γ xz ⎥ ⎢u, z + w, x ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢⎣ 2ε yz ⎥⎦ ⎢⎣γ yz ⎥⎦ ⎢⎣ v, z + w, y ⎥⎦
σ = C ⋅ ε + σ0 (12)
15
Etudes bibliographiques
où C est un tenseur de comportement d’ordre 4 dont les composantes font intervenir les
caractéristiques physiques du matériau, σ 0 est le tenseur de contrainte à l’état initial (pour
simplifier l’écriture du problème, nous supposons que σ 0 = 0 dans la suite).
σ = C⋅ε (14)
⎡1 −ν ν ν 0 ⎤ 0 0
⎢ 1 −ν ν 0 ⎥⎥ 0 0
⎢
⎢ 1 −ν 0 0 0 ⎥
⎢ ⎥
⎢ 1 − 2ν
C=
E 0 0 ⎥
(1 +ν )(1 − 2ν ) ⎢⎢ 2 ⎥
⎥
(15)
1 − 2ν
⎢ 0 ⎥
⎢ 2 ⎥
⎢ 1 − 2ν ⎥
⎢⎢ sym ⎥
⎣ 2 ⎥⎦
où E est le module d'Young et ν est le coefficient de Poisson
Le principe des travaux virtuels consiste à satisfaire l’équation d’équilibre local (9) sous
forme intégrale, on dit aussi sous forme « faible » :
∫u
*
⋅ div ( σ ) dV = ∫uσ
*
i ij , j dV = ∫ uσ
*
i ij n j dS − ∫ ui*, jσ ij dV (17)
Ωt Ωt ∂Ω t Ωt
16
Etudes bibliographiques
∫ uσ n j dS = ∫ ui*σ ij n j dS + ∫ ui*σ ij n j dS = ∫ u* ⋅ f s dS
*
i ij (18)
∂Ω t ∂1Ω t ∂ 2Ωt ∂1Ω t
⎪⎩ Ωt Ωt ∂1Ωt
T
ε* = ⎡⎣ε xx* ε *yy ε zz* 2ε xy
*
2ε xz* 2ε *yz ⎤⎦ =
(20)
* T
= ⎡⎣u,*x v*
,y w *
,z u +v
*
,y
*
,x u +w *
,z
*
,x v + w ⎤⎦
*
,z ,y
En introduisant les équations (14) et (20) dans (19), on obtient l’expression du principe
des travaux virtuels sous forme intégrale en fonction de u* et u :
⎪ ( )
⎧W u* , u = ε*T ⋅ C ⋅ ε dV − u* ⋅ fv dV − u* ⋅ f s dS = 0
∫ ∫ ∫ ∀u*
⎨ Ωt Ωt ∂1Ωt (21)
⎪ u* = 0 sur ∂ Ω et u = u sur ∂ Ω
⎩ 2 t d 2 t
Nous pouvons définir un modèle général où toutes les relations du problème d’élasticité
(9), (11) et (14) sont représentées sous forme variationnelle comme la forme intégrale
suivante :
⎛ ⎛ 1 ⎞ ⎞
⎝ 2 ⎠
(
W = − ∫ ⎜ u*T ⋅ ( div ( σ ) + fv ) + σ*T ⋅ ⎜ ε − ∇u + ∇uT ⎟ + ε*T ⋅ ( σ − C ⋅ ε ) ⎟ dV )
Ωt ⎝ ⎠
(22)
u ⋅ ( f s − σ ⋅ n ) dS − ⎡⎣σ* ⋅ n ⎤⎦ ⋅ ( u − u d ) dS = 0
T
− ∫ ∫ ∀u , σ , ε
*T * * *
∂1Ωt ∂ 2 Ωt
17
Etudes bibliographiques
⎛1 ⎛ 1 ⎞ ⎞
( )
+ ∇u*T ⎤⎦ ⋅ σ − σ*T ⋅ ⎜ ε − ∇u + ∇uT ⎟ − ε*T ⋅ ( σ − C ⋅ ε ) − u*T ⋅ fv ⎟ dV
T
W= ∫ ⎜⎝ 2 ⎡⎣∇u
*
Ωt ⎝ 2 ⎠ ⎠
(23)
− ∫
∂1Ωt
u ⋅ ( f s − σ ⋅ n ) dS −
*T
∫
∂ 2 Ωt
( T
)
⎡⎣σ* ⋅ n ⎤⎦ ⋅ ( u − u d ) + u*T ⋅ ( σ ⋅ n ) dS = 0 ∀u , σ , ε
* * *
L’expression (23) est un principe variationnel de type mixte (faisant intervenir les
variables mixtes : déplacements, déformations, contraintes).
La méthode des éléments finis est une technique particulière d’approximation des
fonctions solutions par sous-domaines (éléments). Les inconnues notées U sont des valeurs
de ces fonctions en certains points ou nœuds de chaque élément. La forme variationnelle
définie sur le milieu continu est ainsi représentée par une forme variationnelle discrétisée qui
fait intervenir les inconnues nodales U . Nous pouvons résumer les démarches de la MEF
comme suit :
V = ∑V e
x (ξ ,η , ζ ) = ∑ N I (ξ ,η , ζ ) xI = N I (ξ ,η , ζ ) xI
I
où x est la position d’un point quelconque ; xI sont des coordonnées des nœuds I
définissant Ve ; ( ξ ,η , ζ ) sont des coordonnées paramétriques ; N I sont des
fonctions d’interpolation (fonctions de forme) dépendant des variables
paramétriques.
18
Etudes bibliographiques
u (ξ ,η , ζ ) = N I* (ξ ,η , ζ ) uI ; u * (ξ ,η , ζ ) = N I* (ξ ,η , ζ ) uI*
We = U*T (
e ⋅ K e ⋅ U e − fe )
où K e est la matrice de rigidité élémentaire ; fe est le vecteur élémentaire des
sollicitations.
( )
W = ∑ We = ∑ U*eT ⋅ K e ⋅ U e − fe = U*T ⋅ ( K ⋅ U − F ) = 0 ∀ U*
e e
soit K ⋅U = F
• Résoudre les relations K ⋅ U = F en tenant compte des conditions aux limites. Pour
un problème linéaire : U = K −1 ⋅ F .
• Évaluer des quantités relatives à chaque élément :
- extraire U e de U
- calculer des déformations, des contraintes …
- On cherche U t +1 = U t + ∆U vérifiant :
19
Etudes bibliographiques
⎡⎣ K ( U t +1 ) ⎤⎦ ⋅ U t +1 − F = R ( U t +1 ) = 0
⎡ ∂R ⎤
R ( U t +1 ) = R ( U t + ∆U ) = R ( U t ) + ⎢ ⎥ ⋅ ∆U + ο
2
⎣⎢ ∂U ⎥
Ut ⎦
≈ R ( U t ) − ⎡⎣K T ( U t ) ⎤⎦ ⋅ ∆U
20
Etudes bibliographiques
Parmi les premiers à avoir tenté de fournir une réponse au problème de coque, on peut
citer Ahmad [3] dans les années 70. Ces modèles, nommés éléments finis tridimensionnels
dégénérés, se basent sur des éléments iso-paramétriques volumiques n’ayant que deux nœuds
suivant la direction de l’épaisseur (voir
Figure 5). Ainsi, ces éléments respectent l’hypothèse de sections droites classiquement
admise pour les coques. À cela est ajoutée une modification du principe des puissances
virtuelles afin de négliger l’énergie engendrée par la déformation normale transverse. Cette
modification impose, en particulier, l’utilisation d’une relation de comportement matériau
prenant en compte l’hypothèse des contraintes planes ( σ zz = 0 ).
⎡σ x ⎤ ⎡εx ⎤
⎢σ ⎥ ⎢ε ⎥
⎢ y⎥ ⎢ y⎥
σ = ⎢τ xy ⎥ = C ⋅ ε = C ⋅ ⎢ε xy ⎥
⎢ ⎥ ⎢ ⎥
⎢τ xz ⎥ ⎢ε xz ⎥
⎢τ yz ⎥ ⎢ε yz ⎥
⎣ ⎦ ⎣ ⎦
où C est une matrice de comportement du matériau (5x5) anisotrope en général. Dans le cas
d’un matériau isotrope, cette matrice s’écrit :
21
Etudes bibliographiques
⎡ 1 ν 0 0 0 ⎤
⎢ 1 0 0 0 ⎥⎥
⎢
⎢ 1 −ν ⎥
⎢ 0 0 ⎥
E 2
C= ⎢ ⎥
(
1 −ν 2 ) ⎢ 1 −ν
0 ⎥⎥
⎢ 2k
⎢ ⎥
⎢ sym. 1 −ν ⎥
⎢⎣ 2k ⎥⎦
Il reste ensuite à choisir un degré d’interpolation suffisant dans les directions du plan
moyen pour éviter les phénomènes de blocage en cisaillement transverse. La qualité
principale de cette modélisation, outre son aspect volumique, est qu’elle s’appuie seulement
sur des degrés de liberté de déplacement. Aucun degré de liberté de rotation n’étant introduit,
le passage de l’analyse linéaire à non-linéaire et la connexion aux éléments volumiques 3D
deviennent des opérations simples. Son inconvénient est la modification du principe de
puissances virtuelles, ce qui revient à modifier la relation de comportement de matériau.
Basés sur cette modélisation, un ensemble d’éléments ont été développés par plusieurs
auteurs. Nous pouvons mentionner ici les travaux de Brendel et Ramm [27], Parisch [64],
Hughes et Liu [47].
Afin d’éviter de modifier le comportement du matériau, Kim et al. [49] ont proposé un
élément iso-paramétrique hexaédrique à 20 nœuds (voir Figure 6). La cinématique, interpolée
par l’intermédiaire des degrés de liberté de déplacement (DDL) aux nœuds, est associée à une
loi matériau tridimensionnelle classique pour le calcul des contraintes. Afin de réduire les
blocages en cisaillement transverse, l’évaluation du principe de puissances virtuelles est
effectuée par une intégration numérique réduite 2x2x3 points de Gauss. Pourtant, cette
modélisation reste chère en temps de calcul par rapport aux éléments coques car elle fait
intervenir 9 DDL le long d’une fibre épaisseur, alors que les éléments coques n’en font
intervenir que 6. Cela incite à plutôt s’orienter vers des éléments de faible degré.
22
Etudes bibliographiques
Avant d’aborder les modélisations de coque en formulation mixte, nous rappelons ci-
après les hypothèses cinématiques de Reissner–Mindlin. Afin de rester le plus généraliste
possible, on se limitera à la modélisation capable de prendre en compte le cisaillement
transverse. Cette théorie s’appuie sur l’hypothèse des sections droites qui consiste à supposer
qu’une droite normale à la surface moyenne de la coque reste droite au cours de la
transformation. Cette droite subit donc seulement une rotation sans élongation (voir Figure 7).
no n
no β
Q
Qo
UP Fibre moyenne
Po P après déformation
Z Configuration
Configuration Fibre moyenne après déformation
initiale avant déformation
X
Figure 7. Hypothèse cinématique de Reissner–Mindlin
G
où β est un vecteur de rotation de la fibre droite, Po représente la projection d’un point Qo sur
la surface moyenne de la configuration initiale, P représente la projection du point Q sur la
surface moyenne de la configuration déformée et z est la distance du point Qo par rapport à la
surface moyenne.
n ⎧β x ⎫ n ⎧ β xI ⎫
w = ∑ N I wI et ⎨ ⎬ = ∑ NI ⎨ ⎬ (25)
I =1 ⎩ β y ⎭ I =1 ⎩ β yI ⎭
23
Etudes bibliographiques
N1 = 1 − ξ − η ; N 2 = ξ ; N 3 = η où ξ ∈ [ 0,1] ; η ∈ [ 0,1 − ξ ] et n = 3
η
3
(0,1)
1
(0,0)
2 ξ
(1,0)
Figure 8. Élément de référence de coque triangulaire à 3 nœuds
Batoz et al. [14], Batoz et Ben Tahar [15] ont montré qu’un élément mixte simple,
triangulaire à 3 nœuds avec approximations linéaires en déplacement w et rotations β x , β y et
efforts tranchants Tx , Ty constants, conduit à un blocage sévère en cisaillement transverse (CT).
Afin d’éliminer ce blocage, plusieurs auteurs ont utilisé la technique nommée des
fonctions « bulle ». Le principe de cette technique est d’enrichir l’interpolation de la formule
(25) en introduisant des variables généralisées. Par exemple, l’approximation de l’élément
coque triangulaire à 3 nœuds ci-dessus peut être enrichie en introduisant une variable α de la
façon suivante :
3 3 ⎧βx ⎫ n ⎧ β xI ⎫
w = ∑ N I wI + ∑ N J α J et ⎨ ⎬ = ∑ NI ⎨ ⎬ (26)
I =1 J =1 ⎩ β y ⎭ I =1 ⎩ β yI ⎭
où α J sont les déplacements des points sur les trois côtés du triangle et N J = (1 − ξ − η ) ξη .
Le terme N J représente une fonction « bulle » qui est nulle sur les trois côtés du triangle.
Plus généralement, ces variables généralisées peuvent être introduites dans l’interpolation des
déplacements ou des rotations et l’ordre des fonctions « bulle » peut être quadratique, cubique
ou d’ordre 4. En se basant sur cette technique, une famille d’éléments finis de plaque ou
coque a vu le jour. Nous citons en particulier :
• L’élément de Pinsky et Jasti [67] basé sur un modèle mixte général dans lequel
l’approche des variables cinématiques est associée à un ensemble de fonctions
« bulle » indépendantes exprimées en fonction des paramètres généralisés que les
auteurs éliminent par condensation statique au niveau local. L’élément « 4-node
bubble » est défini par une interpolation bilinéaire des variables cinématiques w,
24
Etudes bibliographiques
Pourtant, tous les éléments mixtes précédents ont un inconvénient majeur qui est le
temps de calcul élevé causé par l’introduction des fonctions « bulle ». Ils nécessitent plus de
points d’intégration numérique au calcul des matrices élémentaires. De plus, le temps de
calcul devient encore plus grand lorsque l’on traite des problèmes non-linéaires…
L’augmentation des variables généralisées associées aux fonctions « bulle » est également un
inconvénient, lié cette fois-ci aux opérations d’inversion de matrices lors du processus de
condensation statique.
Sans passer par l’utilisation des fonctions « bulle », une autre famille d’éléments de
plaque ou coque de type mixte-hybride a été développée. Nous citons ici l’élément
quadrilatéral MiQ4 présenté par Ayad et al. [8]. Cet élément utilise la même interpolation de
l’équation (25) pour les variables cinématiques. En revanche, les variables mécaniques {M }
et {T } sont définis par :
{M } = [ PM ]{α M } ;
T
M = Mx My M xy
⎡ p 0 0 ⎤
[ PM ] = ⎢⎢ 0 p
⎥
0 ⎥ ; p = 1 ξ η ξη
⎢⎣ 0 0 p ⎥⎦
et par :
25
Etudes bibliographiques
⎧T ⎫ + M xy , y ⎫
⎧M ⎡ p1 0 p2 ⎤
{T } = ⎨Tx ⎬ = ⎨ M x, x ⎬ = [ PT ]{α M } ; [ PT ] = ⎢ ⎥ ;
⎩ y ⎭ ⎩ xy , x + M y , y ⎭ ⎣ 0 p2 p1 ⎦
p1 = 0 j11 j12 η j11 + ξ j12 ; p2 = 0 j21 j22 η j21 + ξ j22
où jmn sont des termes de la matrice Jacobienne inverse et {α M } représentent des paramètres
indéterminés.
⎧γ ⎫ ⎧γ ⎫
{γ } = ⎨ γ xz ⎬ = [ j ] ⎨γ ξ z ⎬ = [ j ][ A]{γ k }
⎩ yz ⎭ ⎩ ηz ⎭
Ensuite, les quantités {γ k } sont projetées sur les deux axes locaux ξ ,η de l’élément en
utilisant les hypothèses de Mindlin sous forme discrète :
⎛ ⎛ ∂w ⎞ ⎞
∫ ⎜ γ S − ⎜ βS +
côtés ⎝ ⎝
⎟ ⎟ ds = 0 ; s est l'axe ξ ou η
∂s ⎠ ⎠
(27)
À partir des relations (27), les quantités {γ k } peuvent être déterminées en fonction des
variables cinématiques w, β x , β y aux nœuds. Les paramètres {α M } , éliminés au niveau
local par condensation statique, sont aussi donnés en fonctions des variables w, β x , β y aux
nœuds.
Cette approche a été utilisée également pour les modèles en déplacements de Bathe et
Dvorkin [11], Chapelle et Bathe [30], Cheung et Chen [32]. Les modèles MiSP3 et MiSP4 ne
présentent pas de blocage de CT et passent bien les patch-tests standards. Ils peuvent être
interprétés comme une amélioration du modèle mixte standard.
Pourtant, les modélisations mixtes présentées ci-dessus basées sur l’approche en
contrainte plane négligent, donc, la variation de l’épaisseur des éléments au cours des
26
Etudes bibliographiques
De plus, les formulations de coques mixtes ont un point commun : les éléments sont
modélisés par leurs surfaces moyennes. Or, dans la réalité, les structures minces (plaques,
coques) et les structures volumiques tridimensionnelles coexistent fréquemment, et ces deux
types d’éléments doivent pouvoir être utilisés simultanément. Cela limite l’utilisation des
modèles mixtes surtout dans les structures dont l’épaisseur varie.
Des éléments effectifs à la fois pour des structures minces et des structures 3D,
simplifieraient considérablement la modélisation de telles structures, et éviteraient deux
procédures supplémentaires : la définition arbitraire de zones de séparation (ex. zones
structurales avec parties 3D) pour assurer la continuité structurale et la connexion des
différents types d’éléments (ex. coque avec 3D). De plus, les éléments 3D ont plusieurs
avantages : le non recours à la cinématique élaborée et complexe de coques ; l’utilisation des
lois de comportement générales tridimensionnelles ; le calcul direct des variations
d’épaisseurs ; la facilité de traitement des grandes rotations avec adaptation simple des
configurations, la connexion naturelle aux autres éléments 3D, et le traitement naturel des
conditions de contact sur les deux faces de la structure. Par conséquent, beaucoup d’efforts
ont été investis pour développer de tel type d’éléments finis nommés « solide–coques » que
nous allons présenter dans les prochains paragraphes.
Pour donner un contexte historique des développements des éléments finis de type
solide–coques, nous pouvons commencer par des travaux de Wilson et al. [88] qui ont
proposé un élément iso-paramétrique hexaédrique à 8 nœuds et de 3 degrés de liberté en
déplacement par nœud. Cet élément est associé à une loi de comportement modifiée de façon
à prendre en compte les contraintes planes. Malheureusement, le faible degré d’interpolation
du modèle est très sensible au phénomène de blocage. Quand les structures modélisées sont en
flexion dominante ou quand le rapport épaisseur sur largueur tend vers zéro, le blocage en CT
se manifeste de façon plus évidente. Afin de surmonter ces verrouillages, Wilson et al. [88]
ont introduit la méthode appelée « modes incompatibles ». Grâce à cette technique, le
comportement de l’élément hexaédrique à 8 nœuds a été nettement amélioré et surtout dans
les problèmes à flexion dominante. Mais, cet élément ne passait pas le patch test pour une
forme géométrique autre que le parallélogramme. Une version modifiée de cet élément, qui
permet de passer le patch test, a été ensuite introduite par Taylor et al. [82]. Plus récemment,
Simo et Rifai [78] ont montré que l’introduction des modes incompatibles pouvait être
justifiée dans le contexte des méthodes mixtes de type « Assumed Strain ».
27
Etudes bibliographiques
Une approche alternative pour contrôler les modes de hourglass, basée sur le principe
variationnel de Hu–Washizu, a été utilisée par Belytschko et al. [23] pour formuler deux
éléments à intégration réduite quadrilatéral et hexaédrique. Leur élément quadrangle à quatre
nœuds présente des blocages pour les matériaux incompressibles mais possède un bon
comportement dans les problèmes de flexion dans le plan. Ce concept a été poursuivi par Liu
et al. [56], Wang et Belytschko [87], Bachrach [10] et bien d’autres afin de résoudre les
problèmes impliquant des matériaux incompressibles.
Au cours des années 1990, les travaux de Simo et Armero [76], Simo et al. [77], Simo et
Rifai [78] ont introduit une nouvelle méthode nommée « Enhanced Assumed Strain » (EAS).
Le principe de cette méthode consiste à enrichir le champ de déformations en rajoutant un
champ de variables qui produit des modes de déformations supplémentaires. Ces champs de
variables sont ensuite projetés sur un sous-espace approprié afin d’éviter les différents types
de blocages.
28
Etudes bibliographiques
Domissy [36], Cho et al. [33], Lemosse [55], Sze et Yao [81], Doll et al. [35], Hauptmann et
al. [41], Hauptmann et Schweizerhof [42], Hauptmann et al. [43], Miehe [59], Klinkel et al.
[50], Klinkel et Wagner [51], Vu-Quoc et Tan [84], Chen et Wu [31], Kim et al. [48], Alves
de Sousa et al. [4], Reese [70], Abed-Meraim et Combescure [1], Legay et Combescure [52]
et bien d’autres. Les deux approches ont été largement utilisées et évaluées à travers plusieurs
applications structurales dans les travaux de Belytschko et Bindeman [22], Zhu et Cescotto
[90], Wriggers et Reese [89], Klinkel et Wagner [51], Wall et al. [86], Reese et al. [71], Puso
[69], pour n’en citer que quelques uns.
29
Chapitre II
Éléments coques volumiques SHB en linéaire
Dans ce chapitre, nous allons présenter les modélisations éléments finis coques
volumiques SHB adaptées à l’étude des structures minces dans le cas linéaire.
Dans la deuxième section, nous ferons une extension de cette famille d’éléments finis
coques volumiques SHB en quadratique. Nous présenterons donc les modélisations de
l’élément prismatique SHB15 à quinze nœuds et de l’élément hexaédrique SHB20 à vingt
nœuds. Les performances de ces éléments finis seront également mises en évidence à travers
une série de cas tests numériques dans le cas linéaire en comparant avec d’autres éléments
finis 3D quadratiques.
Formulation de l’élément SHB6 en linéaire
L’élément SHB6 est formulé dans les axes locaux du plan moyen. Les directions x, y, z
(ou encore x1 , x2 , x3 ) sont respectivement parallèles aux axes ξ ,η , ζ . La Figure 9 représente
la géométrie d’un élément de référence SHB6 élémentaire et ses points d’intégration.
ξ
6 5
5 ς
4
4 O 3
2
3 2
1
31
Formulation de l’élément SHB6 en linéaire
6
xi = xiI N I (ξ ,η , ζ ) = ∑ xiI N I (ξ ,η , ζ ) (28)
I =1
Les mêmes fonctions de forme sont utilisées pour définir le champ de déplacement de
l’élément ui en termes des déplacements nodaux uiI :
6
ui = uiI N I (ξ ,η , ζ ) = ∑ uiI N I (ξ ,η , ζ ) (29)
I =1
ui , j = uiI N I , j (30)
1
ε ij =
2
( ui, j + u j ,i ) (31)
32
Formulation de l’élément SHB6 en linéaire
⎛ (1 − ξ )η ⎞
⎜ ⎟
⎜ (1 − ξ ) ζ ⎟
1 ⎜ (1 − ξ )(1 − η − ζ ) ⎟
N (ξ ,η , ζ ) = ⎜ ⎟
2⎜ (1 + ξ )η ⎟ (32)
⎜ (1 + ξ ) ζ ⎟
⎜ ⎟
⎜ (1 + ξ )(1 − η − ζ ) ⎟
⎝ ⎠
ζ = [ 0,1] ; η = [ 0,1 − ζ ] ; ξ = [ −1,1]
L’origine du repère est confondue avec le coin droit du triangle du plan médian de
l’élément (celui qui possède un angle droit, voir Figure 9).
En évaluant l’équation (33) aux nœuds de l’élément, on arrive aux trois systèmes de six
équations suivants :
Dans l’équation précédente et dans toute la suite, les caractères soulignés ou gras
désignent des tenseurs d’ordre au moins un (vecteurs, tenseurs de contraintes, de
déformation…).
⎧⎪d Ti = ( ui1 , ui 2 , ui 3 , ui 4 , ui 5 , ui 6 )
⎨ T (35)
⎪⎩ x i = ( xi1 , xi 2 , xi 3 , xi 4 , xi 5 , xi 6 )
33
Formulation de l’élément SHB6 en linéaire
⎧ST = (1, 1, 1, 1, 1, 1)
⎪⎪ T
⎨h1 = ( −1, 0, 0, 1, 0, 0 ) (36)
⎪ T
⎪⎩h 2 = ( 0, −1, 0, 0, 1, 0 )
∂N
b i = N ,i ( 0 ) = i = 1, 2, 3 (37)
∂xi ξ =η =ζ =0
En fait, ce calcul est pratique car nous verrons par la suite qu’il est nécessaire de
déterminer la quantité suivante :
⎛ T 2
⎞
⎜
⎝
b j + ∑
α =1
hα , j γ αT ⎟ = N , j ( ξ ,η , ζ
⎠
T
)
b Tj = N T, j (0, 0, 0) = csteT
⎛ ∂ N1 ∂ N2 ∂ N3 ∂ N4 ∂ N5 ∂ N6 ⎞
bTi = NT,i ( 0 ) = ⎜ ⎟
⎝ ∂xi ∂xi ∂xi ∂xi ∂xi ∂xi ⎠ ξ =η =ζ =0
∂ N I ⎛ ∂ N I ∂ξ ∂ N I ∂η ∂ N I ∂ζ ⎞ ⎛ ∂ NI ∂ NI ∂ NI ⎞
=⎜ + + ⎟⎟ = ⎜ j1 j + j2 j + j3 j ⎟
∂x j ⎜⎝ ∂ξ ∂x j ∂η ∂x j ∂ζ ∂x j ⎠ ⎝ ∂ξ ∂η ∂ζ ⎠
avec I = 1, 2,..., 6 et j = 1, 2,3
⎛ ∂ξ ∂ξ ∂ξ ⎞
⎜ ⎟
⎜ ∂x ∂y ∂z ⎟
⎛ j11 j12 j13 ⎞
⎜ ∂η ∂η ∂η ⎟ ⎜ ⎟
F −1 = ⎜ ⎟ = j21 j22 j23 ⎟
⎜ ∂x ∂y ∂z ⎟ ⎜⎜
j31 j32 j33 ⎠⎟
⎜ ∂ζ ∂ζ ∂ζ ⎟ ⎝
⎜ ∂x ∂y ∂z ⎟⎠
⎝
En ξ = η = ζ = 0, nous avons :
34
Formulation de l’élément SHB6 en linéaire
∂ N1 1 ∂ N1 1 1 ∂ N1
=− η =0 = (1 − ξ ) = =0
∂ξ 2 ∂η 2 2 ∂ζ
∂ N2 1 ∂ N2 ∂ N2 1 1
=− ζ =0 =0 = (1 − ξ ) =
∂ξ 2 ∂η ∂ζ 2 2
∂ N3 1 1 ∂ N3 1 1 ∂ N3 1 1
= − (1 − η − ζ ) = − = − (1 − ξ ) = − = − (1 − ξ ) = −
∂ξ 2 2 ∂η 2 2 ∂ζ 2 2
∂ N4 1 ∂ N4 1 1 ∂N 4
= η =0 = (1+ξ ) = =0
∂ξ 2 ∂η 2 2 ∂ζ
∂N 5 1 ∂N 5 ∂N 5 1 1
= ζ =0 =0 = (1+ξ ) =
∂ξ 2 ∂η ∂ζ 2 2
∂N 6 1 1 ∂N 6 1 1 ∂N 6 1 1
= (1 − η − ζ ) = = − (1+ξ ) = − = − (1+ξ ) = −
∂ξ 2 2 ∂η 2 2 ∂ζ 2 2
De plus, on peut vérifier par des considérations algébriques que les conditions
d’orthogonalité suivantes sont satisfaites :
⎧bTi ⋅ hα = 0
⎪ T
⎪b i ⋅ S = 0
⎪ T i, j = 1,2,3
⎨bi ⋅ x j = δ ij ; avec (38)
⎪ T α , β = 1,2
⎪hα ⋅ S = 0
⎪hT ⋅ h = 2δ
⎩ α β αβ
A ce stade, on peut déterminer les constantes inconnues qui interviennent dans l’écriture
(34) du champ de déplacement en multipliant scalairement l’équation (34) par b Tj et hαT ,
respectivement, et en utilisant les relations d’orthogonalité (38). On obtient ainsi :
35
Formulation de l’élément SHB6 en linéaire
bTj ⋅ di = a0i bTj ⋅ S + a1i bTj ⋅ x1 + a2i bTj ⋅ x 2 + a3i bTj ⋅ x3 + c1i bTj ⋅ h1 + c2i bTj ⋅ h 2
= 0 + a1i bTj ⋅ x1 + a2i bTj ⋅ x 2 + a3i bTj ⋅ x3 + 0 = a ji
h1T ⋅ di = a0i h1T ⋅ S + a1i h1T ⋅ x1 + a2i h1T ⋅ x 2 + a3i h1T ⋅ x3 + c1i h1T ⋅ h1 + c2i h1T ⋅ h 2
( )( )
= 0 + h1T ⋅ x j bTj ⋅ di + 2c1i
hT2 ⋅ di = a0i hT2 ⋅ S + a1i hT2 ⋅ x1 + a2i hT2 ⋅ x 2 + a3i hT2 ⋅ x3 + c1i hT2 ⋅ h1 + c2i hT2 ⋅ h 2
( )( )
= 0 + hT2 ⋅ x j bTj ⋅ di + 2c2i
⎧⎪a ji = bTj ⋅ di
⎨ (39)
⎪⎩cα i = γα ⋅ di
T
1⎡ ⎛ 3 T ⎞⎤
γα =
2 ⎢⎣
h − (
⎢ α ⎜ ∑ hα ⋅ x j b j ⎟ ⎥) α = 1, 2 (40)
⎝ j =1 ⎠ ⎥⎦
⎛ 2
⎞
( )
ui , j = ⎜ bTj + ∑ hα , j γαT ⎟ ⋅ di = bTj + hα , j γαT ⋅ di (42)
⎝ α =1 ⎠
∇s (u ) = B ⋅ d (43)
où :
36
Formulation de l’élément SHB6 en linéaire
⎡ ux,x ⎤
⎢ u ⎥
⎢ y,y ⎥ ⎡ d1 ⎤
⎢ u z,z ⎥ ,
∇ (u ) = ⎢
s
⎥
d = ⎢⎢d 2 ⎥⎥ (44)
⎢u x , y + u y , x ⎥ ⎣⎢ d3 ⎦⎥
⎢ ux,z + uz,x ⎥
⎢ ⎥
⎣⎢ u y , z + u z , y ⎦⎥
⎡ b1T + hα , x1 γ αT 0 0 ⎤
⎢ ⎥
⎢ 0 b + hα , x2 γ α
T
2
T
0 ⎥
⎢ T ⎥
⎢ 0 0 b 3 + hα , x3 γ α ⎥
T
B=⎢ T ⎥
(45)
b + hα , x2 γ αT b1T + hα , x1 γ αT 0
⎢ 2 ⎥
⎢ b T3 + hα , x γ αT 0 b1T + hα , x1 γ αT ⎥
⎢ 3
⎥
⎢⎣ 0 b + hα , x3 γ αT
T
3 b T2 + hα , x2 γ αT ⎥⎦
Cette écriture de l’opérateur gradient discrétisé utilisant les formules de Hallquist [40]
est très commode car les vecteurs γα , qui interviennent dans l’expression de B , vérifient les
conditions d’orthogonalité suivantes :
γ αT ⋅ x j = 0 , γ αT ⋅ h β = δ αβ (46)
Ceci permet de manipuler séparément chacun des modes de déformation pour obtenir
simplement la forme du champ de déformation postulée. Notons qu’un élément basé sur la
formulation (45) est convergent lorsqu’il est évalué exactement. Cependant, l’évaluation de
cet opérateur B , donné en (45), en chacun des points d’intégration rend cet élément coûteux
en temps de calcul pour les applications pratiques, et une forme simplifiée de cet élément est
souvent préférée.
δπ ( u , ε , σ ) = ∫ δ ε T ⋅ σ d V + δ ∫ σ T ⋅ ( ∇ s u − ε ) d V − δ d T ⋅ f ext = 0 (47)
Ve Ve
37
Formulation de l’élément SHB6 en linéaire
δπ ( u , ε ) = ∫ δ ε T ⋅ σ d V − δ d T ⋅ f ext = 0 (48)
Ve
6
u ( x, t ) = ∑ N I ( x ) d I ( t ) (49)
I =1
Ceci conduit à :
∇ s u ( x, t ) = B ( x ) ⋅ d ( t ) (50)
ε ( x , t ) = B ( x ) ⋅ d ( t ) (51)
⎛ T ⎞
δ d T ⋅ ⎜ ∫ B ⋅ σ dV − f ext ⎟ = 0 (52)
⎜V ⎟
⎝ e ⎠
Dans l'équation ci-dessus, il est bien précisé que la contrainte σ est calculée par la loi
constitutive à partir du taux de déformation postulée ε . Pour les problèmes non linéaires, σ
peut aussi être une fonction intégrale du taux de déformation postulée et des autres variables
internes :
38
Formulation de l’élément SHB6 en linéaire
σ = F ( ε , α,...) (55)
où α représente les variables internes. La formulation ainsi obtenue est valable pour des
problèmes incluant les deux types de non linéarités : géométriques et matériaux. Dans le cas
de problèmes linéaires, on a :
σ = C⋅ ε = C⋅B ⋅d (56)
La matrice de comportement élastique C , dans le cas d’un matériau isotrope, est choisie
comme suit :
⎡ E Eν ⎤
⎢1 −ν 2 0 0 0 0 ⎥
1 −ν 2
⎢ ⎥
⎢ Eν E
0 0 0 0 ⎥
⎢1 −ν 2 1 −ν 2 ⎥
⎢ ⎥
⎢ 0 0 E 0 0 0 ⎥
C=⎢ 0 0 0
E
0 0 ⎥
⎥
⎢
⎢ 2 (1 + ν ) ⎥
⎢ E ⎥
⎢ 0 0 0 0 0 ⎥
⎢ 2 (1 + ν ) ⎥
⎢ E ⎥
⎢ 0 0 0 0 0 ⎥
⎢⎣ 2 (1 + ν ) ⎥⎦
Dans cette matrice, E est le module d’Young et ν est le coefficient de Poisson. Cette loi
est spécifique aux éléments SHB. Elle ressemble à celle que l’on aurait dans le cas de
l’hypothèse des contraintes planes, mise à part le terme (3,3). Nous pouvons noter que ce
choix entraîne un comportement anisotrope artificiel. Ce choix permet de satisfaire tous les
tests sans introduire de blocage.
f int = K e ⋅ d (57)
où :
Ke = ∫ B ⋅C ⋅ B dV
T
Ve
(58)
Ke = ∫ B ⋅C ⋅B dV
T
Ve
(59)
39
Formulation de l’élément SHB6 en linéaire
η ζ ξ w(ξ,η,ζ)
P(1) 1/3 1/3 -0.906179845938664 0.236926885056189
P(2) 1/3 1/3 -0.538469310105683 0.478628670499366
P(3) 1/3 1/3 0 0.568888888888889
P(4) 1/3 1/3 0.538469310105683 0.478628670499366
P(5) 1/3 1/3 0.906179845938664 0.236926885056189
5
K e = ∑ w ( Pj ) J ( Pj )BT ( Pj ) ⋅ C ⋅ B ( Pj )
j =1
Les modes de « hourglass » sont des modes cinématiques qui sont dus à la sous-
intégration et sont associés à une énergie nulle alors qu’ils induisent une déformation non
nulle. Cette anomalie s’explique par la différence, qu’induit la sous-intégration, entre le noyau
de l’opérateur de rigidité discrétisé et celui continu. Commençons d’abord par remarquer que
l’opérateur gradient discrétisé sous-intégré (i.e., associé à cinq points d’intégration situés en
1
ζ = η = , ξ = ξ Gi ) prend la forme de l’équation (45) avec α = 1,2 .
3
5
K e = ∑ w ( Pj ) J ( Pj ) BT ( Pj ) ⋅ C ⋅ B ( Pj ) (60)
j =1
∇ su = B ⋅ d = 0 (61)
40
Formulation de l’élément SHB6 en linéaire
Nous allons chercher désormais quels sont les modes de déformations qui donnent une
1
énergie de déformation nulle. L’énergie de déformation s’écrit w(ε ) = ∫ ε T ⋅ C ⋅ ε d V , et
2V
comme ε = B ⋅ d , nous avons donc :
⎡1 ⎤
w(ε ) = d T ⋅ ⎢ ∫ B T ⋅ C ⋅ B dV ⎥ ⋅ d
⎣2 V ⎦
1 T
w(ε ) = d ⋅Ke ⋅d
2
Ainsi rechercher les modes de déformations à énergie nulle, c’est rechercher le noyau de
Ke :
K e ⋅ X = 0 ⇒ B (ξ Gj ) ⋅ X = 0 ∀ ξ Gj
En d’autres termes, les modes de hourglass sont les vecteurs X tels que :
Il est naturel de retrouver dans le noyau de la rigidité K e les modes associés aux
mouvements de corps rigides. Pour un élément tridimensionnel tel que le prisme à 6 nœuds,
ces mouvements rigidifiant sont composés de trois translations et de trois rotations. Ainsi le
noyau de l’opérateur continu de rigidité est de dimension six et se réduit aux seuls modes
suivants :
⎛S 0 0⎞ ⎛y z 0⎞
⎜ ⎟ ⎜ ⎟
⎜0 S 0⎟ et ⎜ -x 0 z⎟ (63)
⎜0 0 S⎟ ⎜0 -x -y ⎟⎠
⎝ ⎠ ⎝
On vérifie aisément que chacun des six vecteurs colonnes ci-dessus satisfait à l’équation
(62) et appartient donc au noyau de K e . Il suffit, pour le voir, d’utiliser l’expression (45) de
B et les conditions d’orthogonalité (38) et (46). Les trois premiers vecteurs colonnes
correspondent aux translations selon les axes Ox, Oy et Oz, respectivement. Les trois autres
vecteurs sont relatifs aux rotations autour des axes Oz, Oy et Ox.
Nous recherchons désormais, outre les modes rigides précédents, des modes qui
annulent également l’opérateur gradient discrétisé donné en (45). Prenons une base de dix-
huit vecteurs suivants :
41
Formulation de l’élément SHB6 en linéaire
⎡S 0 0 x 0 0 y 0 0 z 0 0 h1 0 0 h2 0 0⎤
e i = ⎢⎢ 0 S 0 0 x 0 0 y 0 0 z 0 0 h1 0 0 h2 0 ⎥⎥
(a)
⎣⎢ 0 0 S 0 0 x 0 0 y 0 0 z 0 0 h1 0 0 h 2 ⎦⎥
i = 1 à 18
On peut montrer que les vecteurs ci-dessus sont linéairement indépendants dans
l’espace de dimension dix-huit. Pour cela, nous étudions un système de dix-huit d’équations
suivantes dont ai sont des inconnus:
18
∑a e
i =1
i i =0 (b)
⎧ T
) ⋅ ⎛⎜ ∑ a e ⎞⎟ = a
18
(
⎪ bi 0 0 i i 3i +1 = 0
⎪ ⎝ i =1 ⎠
⎪⎪ ⎛ 18 ⎞
( T
⎨ 0 bi 0 ) ⋅ ⎜ ∑ ai ei ⎟ = a3i + 2 = 0 , i = 1, 2,3 (c)
⎪ ⎝ i =1 ⎠
⎪ ⎛ 18
⎞
(⎪ 0 0 bTi )⋅ ⎜ ∑ ai ei ⎟ = a3i +3 = 0
⎩⎪ ⎝ i =1 ⎠
⎧ T ⎛ 18 ⎞
(
⎪ hα 0 0 ) ⋅ ⎜ ∑ ai ei ⎟ = 2a3α +10 = 0
⎪ ⎝ i =1 ⎠
⎪⎪ ⎛ 18
⎞
( T
⎨ 0 hα 0 ) ⋅ ⎜ ∑ ai ei ⎟ = 2a3α +11 = 0 , α = 1, 2 (d)
⎪ ⎝ i =1 ⎠
⎪ ⎛ 18
⎞
(
⎪ 0 0 hαT
⎪⎩
)⋅ ⎜ ∑ ai ei ⎟ = 2a3α +12 = 0
⎝ i =1 ⎠
a1 = a2 = a3 = 0 (e)
Cela montre que les vecteurs ei sont linéaires indépendants, et ils forment, donc, une
base pour l’espace vectoriel de dimension dix-huit.
42
Formulation de l’élément SHB6 en linéaire
18
X = ∑ ci ei (f)
i =1
⎛S⎞ ⎛ 0⎞ ⎛ 0⎞ ⎛ −y ⎞ ⎛ −z ⎞ ⎛ 0⎞
⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟
X = c1 ⎜ 0 ⎟ + c2 ⎜ S ⎟ + c3 ⎜ 0 ⎟ + c5 ⎜ x ⎟ + c6 ⎜ 0 ⎟ + c9 ⎜ − z ⎟
⎜ 0⎟ ⎜ 0⎟ ⎜S⎟ ⎜ 0⎟ ⎜ x⎟ ⎜ y⎟
⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠
Cela veut dire qu’il n’y a pas d’autres modes que les modes rigides qui annulent
l’opérateur gradient discrétisé donné en (45). Autrement dit, l’élément SHB6 ne présente pas
de modes de hourglass.
La première étape consiste à se placer dans le repère local de l’élément défini par le
repère (ξ ,η , ζ ) décrit dans la Figure 9. Les déformations seront désormais calculées dans ce
repère. L’opérateur gradient discrétisé B sera projeté sur un sous espace approprié afin
d'éviter les différents problèmes de blocage. Cette méthode est variationnellement cohérente
avec le principe de Hu–Washizu si l'interpolation de la contrainte est choisie judicieusement
(Simo et Hughes [75]). Cependant, il est très difficile de sélectionner de manière générale et
43
Formulation de l’élément SHB6 en linéaire
Nous présentons ici un choix efficace et acceptable. L’opérateur B est tout d’abord
séparé en deux parties B1 et B 2 . La matrice B1 contient les gradients dans le plan moyen de
la coque et la déformation normale, tandis que B 2 contient les gradients associés aux
déformations de cisaillement transverse.
⎡ b1T + hα , x1 γ αT 0 0 ⎤ ⎡ 0 0 0 ⎤
⎢ ⎥ ⎢ ⎥
⎢ 0 b + hα , x2 γ α
T
2
T
0 ⎥ ⎢ 0 0 0 ⎥
⎢ T ⎥ ⎢ 0 0 0 ⎥
0 0 b 3 + hα , x3 γ α ⎥
T
B=⎢ + ⎢ ⎥
⎢b + h γ T
T
b1T + hα , x1 γ αT 0 ⎥ ⎢ 0 0 0 ⎥
α , x2 α
⎢ 2 ⎥ ⎢ b 3 + hα , x γ αT
T
0 b1T + hα , x1 γ αT ⎥
⎢ 0 0 0 ⎥ ⎢ 3
⎥
⎢ ⎥ ⎢⎣ 0 b T3 + hα , x3 γ αT b T2 + hα , x2 γ αT ⎥⎦
⎣⎢ 0 0 0 ⎦⎥
= B1 + B 2
⎡ 0 0 0 ⎤
⎢ 0 0 0 ⎥
⎢ ⎥
⎢ 0 0 0 ⎥
B2 = c ⎢ ⎥
⎢ 0 0 0 ⎥
⎢ b T3 + hα , x γ αT 0 b1 + hα , x1 γ αT
T ⎥
⎢ 3
⎥
⎢⎣ 0 b T3 + hα , x3 γ αT b T2 + hα , x2 γ αT ⎥⎦
Les matrices K e1 , K e 2 , K e 3 et K e 4 sont intégrées avec les cinq points de Gauss définis
précédemment. La décomposition additive donnée plus haut, B = B1 + B 2 , pour l’opérateur
gradient discrétisé, fait que les termes croisés K e 2 et K e3 s’annulent. Suite aux nombreux
tests numériques, il a été choisi de caractériser la matrice B 2 par le coefficient : c = 0, 45 , qui
joue le rôle ici d’un facteur de réduction du cisaillement.
44
Formulation de l’élément SHB6 en linéaire
Ce choix donne à l’élément un bon comportement dans les cas tests de référence. Il est
clair que cette stratégie, comme celle mise en place pour les éléments coques volumiques
quadratiques, n’est adaptée qu’au comportement quasi isotrope du matériau choisi.
( K + µK σ ) ⋅ u = 0 ⇒ K ⋅ u = λK σ ⋅ u
avec λ = − µ , et µ est le coefficient multiplicateur du chargement.
3
eijQ ( δu , ∆u ) = ∑ δ uk ,i ∆uk , j = δ uk ,i ∆uk , j
k =1
δ uT ⋅ K σ ⋅ ∆u = ∫ σ : eQ (δ u, ∆u ) d Ω = ∫ σ : ∇δ uT ⋅∇∆u d Ω
Ω0 Ω0
Afin d’exprimer cette matrice dans l’espace discrétisé, introduisons les opérateurs
gradient quadratique discrétisés BQ tels que :
⎡ e11Q ⎤ ⎡ δ uT ⋅ B11 Q
⋅ ∆u ⎤
⎢ Q ⎥ ⎢ T Q ⎥
⎢ e22 ⎥ ⎢δ u ⋅ B 22 ⋅ ∆u ⎥
⎢ e Q ⎥ ⎢ δ u T ⋅ B Q ⋅ ∆u ⎥
eQ ( δu , ∆u ) = ⎢ Q 33 Q ⎥ = ⎢ T 33 ⎥
⎢ e12 + e21 ⎥ ⎢δ u ⋅ B12 ⋅ ∆u ⎥
Q
⎢ eQ + eQ ⎥ ⎢δ uT ⋅ BQ ⋅ ∆u ⎥
⎢ 13 31
⎥ ⎢ T 13 ⎥
⎣⎢ 23 32 ⎦⎥ ⎣⎢δ u ⋅ B 23 ⋅ ∆u ⎦⎥
+
Q Q Q
e e
Les différents termes BQij sont donnés par les équations suivantes :
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
⎢ B1B1T 0 0 ⎥ ⎢ B 2 B T2 0 0 ⎥ ⎢ B3B3T 0 0 ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
Q
B11 = ⎢
⎢ 0 B1B1T 0 ⎥
⎥ ; BQ22 = ⎢
⎢ 0 B 2 B T2 0 ⎥
⎥
Q
; B33 = ⎢
⎢ 0 B3B3T 0 ⎥
⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
T⎥ T⎥ T⎥
⎢
⎣⎢
0 0 B1B 1 ⎦⎥
⎢
⎣⎢
0 0 B2B 2 ⎦⎥
⎢
⎣⎢
0 0 B 3B 3 ⎦⎥
45
Formulation de l’élément SHB6 en linéaire
⎡ ⎤
B1B T2 + B 2B1T
⎢ 0 0 ⎥
⎢ ⎥
Q
B12 = ⎢
⎢ 0 B1B 2 + B 2B1T
T
0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B1B T2 + B 2 B T⎥
1 ⎦⎥
⎡ ⎤
⎢ B1B3T + B3B1T 0 0 ⎥
⎢ ⎥
Q
B13 =c 2 ⎢
⎢ 0 B1B3T + B3B1T 0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B1B + B3B
T
3
T⎥
1 ⎦⎥
⎡ ⎤
B 2 B3T + B3B T2
⎢ 0 0 ⎥
⎢ ⎥
BQ23 = c 2 ⎢
⎢ 0 B 2 B3 + B3B T2
T
0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B 2 B3T + B3B T⎥
2 ⎦⎥
(
Bi = bi + hα ,i γα )
Q
Remarque : Nous devons multiplier les matrices B13 et BQ23 par le coefficient
c 2 = 0, 452 = 0, 2025 car l’élément SHB6 est projeté par la technique « Assumed local strain
method » voir section 1.1.5.
k σ (ξ j ) = σ 11 (ξ j ) B11
Q
(ξ j ) + σ 22 (ξ j ) BQ22 (ξ j ) + σ 33 (ξ j ) BQ33 (ξ j )
+ σ 12 (ξ j ) B12
Q
(ξ j ) + σ13 (ξ j ) B13Q (ξ j ) + σ 23 (ξ j ) BQ23 (ξ j )
Par intégration sur les points de Gauss de l’élément, la matrice de rigidité géométrique
s’obtient par la formule :
5
K σ = ∑ ω (ξ j ) J (ξ j ) k σ (ξ j )
j =1
Les forces de pression suiveuses sont présentes dans la matrice tangente via la matrice
K P , car les forces externes suiveuses dépendent du déplacement. Les forces de pression
suiveuses s’écrivent :
∫ p nT ⋅ u dS = ∫
T
p det[ F (u)] nT0 ⋅ F (u) −1 dS0 = p F0 − p K P ⋅ U
∂Ω ∂Ω0
46
Formulation de l’élément SHB6 en linéaire
F (u) = I + ∇u
La formulation précédente conduit à une matrice non-symétrique. On sait que l’on peut
néanmoins utiliser une formulation symétrique si les forces extérieures dues à la pression
dérivent d’un potentiel. C’est le cas si les forces de pression ne travaillent pas sur la frontière
du domaine modélisé. On considère donc que la partie symétrique de la matrice suffit. La
matrice symétrisée prend la forme suivante :
⎡ 0
T T
b 2 n1 − b 1 n 2 b 3 n1 − b 1 n 3 ⎤
T T
⎢ ⎥
⎢ b 2 n1 − b 1 n 2 b 3 n1 − b 1 n 3 ⎥
T T T T
0
⎢ ⎥
b 2 n1 − b 1 n 2 b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0
⎢T T ⎥
b 3 n 2 − b 2 n 3 ⎥
T T
⎢b 1 n 2 − b 2 n1 0
K P = S 0 ⎢b 1T n 2 − b T2 n1 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢ ⎥
⎢b T n − b T n 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢ 1T 2 2 1
⎥
⎢ b 1 n 3 − b 3 n1 b 2 n 3 − b T3 n 2
T T
0 ⎥
⎢ T ⎥
⎢ b 1 n 3 − b 3 n1 b 2 n 3 − b 3 n 2
T T T
0 ⎥
⎢T T T T ⎥
⎣⎢ b 1 n 3 − b 3 n1 b 2 n 3 − b 3 n 2 0 ⎦⎥
C’est une matrice ( 9 × 9 ) , qu’il faut multiplier par les déplacements des 3 nœuds de la
face sur laquelle on applique une pression.
47
Validation de l’élément SHB6 en linéaire
Pour valider le nouvel élément SHB6, nous l’avons testé sur un ensemble de cas tests
élastiques. La réalisation de divers cas tests nous permet de voir l’influence du choix d’une
projection de l’opérateur gradient discrétisé sur les résultats de problèmes significatifs en
mécanique. Pour chaque cas test, le résultat obtenu est comparé, d’une part, à la solution de
référence, et d’autre part, à la solution donnée par l’élément fini volumique usuel existant
PRI6 et à l’élément SHB6 initial sans projection. Tous les maillages n’ont qu’un seul élément
dans l’épaisseur. Dans tous les tableaux apparaissent le nombre de découpage dans chacune
des directions, le nombre total d’éléments étant le double (car chaque rectangle est découpé en
deux triangles). L’élément projeté sera dénommé SHB 6 .
Les données de géométrie, matériau et chargement du test sont listées dans le Tableau 1.
Longueur L = 50 m
Largeur l=4m
Épaisseur e=1m
Module d’Young E = 6,825x107 Pa
Coefficient de Poisson ν = 0,3
G
Chargement P =4N
PL3 l h 3 l e3
flèche = u z = avec I= =
3EI 12 12
PL3 4 × 4 × 503
flèche = u z = = = 7,326 × 10−3 m
3EI 6,825 × 10 × 4 × 1
7 3
48
Validation de l’élément SHB6 en linéaire
z
x y
P
P L e
L e
A l
A
l
a) Géométrie et le chargement de la poutre b) Configuration initiale et déformée (échelle 3000:1)
Figure 10. Géométrie, chargement et conditions aux limites pour le test de la poutre en
flexion simple ; un exemple de maillage (12x2x1)x2
1,2
1
Uz/Uref du point A
0,8
0,6
0,4 SHB6
initial
0,2 PRI6
0
SHB6
0 100 200 300 400
Nombre d'éléments
Graphe 1. Déplacement normalisé du point A de la poutre en flexion simple
49
Validation de l’élément SHB6 en linéaire
La plaque est encastrée parfaitement à l’extrémité gauche. Le chargement est une force
P appliquée à l'extrémité droite de la poutre ce qui entraîne une répartition de contraintes
parabolique sur la hauteur et constante sur l’épaisseur. Le problème est un cas de contraintes
planes. Par ailleurs, on s’attend à ce que les déplacements ne soient pas « infiniment petits »
devant les dimensions de la structure (voir Figure 11a).
Longueur L = 48 mm
Hauteur H = 12 mm
Épaisseur h = 1 mm
Module d’Young E = 3x1010 Pa
Coefficient de Poisson ν = 0,25
G
Chargement P = 40 N
y
y z x
z x
h
P
h
A
P
L
A H
L b) Configuration initiale
H et déformée (échelle 50:1)
a) Géométrie et chargement
La solution de référence de ce cas test est calculée par la méthode de la fonction d’Airy.
Nous n’utilisons pas les formules de RDM.
50
Validation de l’élément SHB6 en linéaire
C 3
φ = Axy + xy + By + Dy 3 .
6
⎡ ⎤
⎢0 0 0 ⎥
⎢ ⎥
φ = ⎢0 0 0 ⎥
⎢ C 3 ⎥
⎢0 0 Axy + xy + By + Dy 3 ⎥
⎣ 6 ⎦
⎛⎡ C 2 ⎤⎞ ⎡ C ⎤
⎜ ⎢0 0 Ax + xy + B + 3Dy 2 ⎥ ⎟ ⎢Cxy + 6 Dy − A − y 2 0⎥
⎜⎢ 2 2
⎥⎟ ⎢ ⎥
( ( ))
σ = rot rot φ = rot ⎜⎜ ⎢ 0 0
⎢
C 3
− Ay − y
6 ⎥⎟ ⎢
C
⎥ ⎟ = ⎢− A − y2
2
0 0⎥
⎥
⎜⎢ ⎥⎟ ⎢
⎜ ⎢0 0 0 ⎥⎟ ⎢ 0 0 0 ⎥⎥
⎜⎢ ⎟
⎝⎣ ⎦⎥ ⎠ ⎣⎢ ⎥⎦
H
Sur la face libre y = + on a :
2
⎡ H H CH 2 ⎤
⎢ Cx + 6 D − A − 0⎥
⎢ 2 2 8 ⎥ ⎡0⎤ ⎡0 ⎤
⎢ CH 2
⎥
T = σ ⋅y = ⎢ −A− 0 0 ⎥ ⋅ ⎢⎢1 ⎥⎥ = ⎢⎢0 ⎥⎥
8
⎢ ⎥ ⎢0⎥ ⎢0 ⎥
⎢ 0 0 0⎥ ⎣ ⎦ ⎣ ⎦
⎢ ⎥
⎣ ⎦
Ce qui conduit à :
CH 2 CH 2
− A− =0⇒ A=−
8 8
Sur la face x = L :
⎡ C 2 ⎤
⎢CLy + 6 Dy − A − 2 y 0 ⎥
⎢ ⎥ ⎡1 ⎤ ⎡ 0 ⎤
C 2
⎢
T = σ ⋅x = −A− y 0 0 ⎥ ⋅ ⎢⎢ 0 ⎥⎥ = ⎢⎢σ xy ⎥⎥
⎢ 2 ⎥
⎢ ⎢0⎥ ⎢ 0 ⎥
⎢ 0 0 0 ⎥⎥ ⎣ ⎦ ⎣ ⎦
⎣⎢ ⎦⎥
Ce qui entraîne :
CL
D=−
6
Il nous faut calculer :
51
Validation de l’élément SHB6 en linéaire
H
h H h H 2
2 2 2 2
⎛ C 2⎞ ⎡ C 3⎤ ChH 3
P= ∫∫
h H
σ xy d s = ∫h ∫H ⎜⎝ − A − 2 y ⎟⎠ d y d z = h ⎣⎢− Ay − 6 y ⎦⎥ =
12
− − − − H
2 2 2 2 −
2
12 P 2 PL 3P
Il vient donc que nous avons : C = 3
ce qui entraîne D = − 3
et A = − .
hH hH 2hH
⎡ 12 Py ( x − L ) 6P ⎛ H 2 ⎞ ⎤
⎢ 3 ⎜
− y2 ⎟ 0⎥
⎢ 2 hH 3 2 hH ⎝ 4 ⎠ ⎥
⎢ 6P ⎛ H 2 ⎞ ⎥
σ=⎢ 3 ⎜
− y2 ⎟ 0 0⎥
⎢ 2 hH ⎝ 4 ⎠ ⎥
⎢ ⎥
⎢ 0 0 0⎥
⎢ ⎥
⎣ ⎦
12 Py( x − L ) 6P ⎛ H 2 ⎞
σ xx = 3
; σ yy = 0 ; σ xy = 3 ⎜
− y2 ⎟
2hH 2hH ⎝ 4 ⎠
⎛ 1 +ν ⎞ ν
ε=⎜ ⎟ σ − sI avec s = tr ( σ )
⎝ E ⎠ E
D’où :
⎡ 12 Py ( x − L ) ⎛ 1 +ν ⎞ 6 P ⎛ H
2
⎞ ⎤
⎢ 3 ⎜ ⎟ 3 ⎜ − y2 ⎟ 0 ⎥
⎢ 2 EhH ⎝ E ⎠ 2 hH ⎝ 4 ⎠ ⎥
⎢ 1 +ν 12 Py ( x − L ) ⎥
⎞ 6P ⎛ H ⎞
2
⎛
ε = ⎢⎜ ⎟ 3 ⎜
− y2 ⎟ −ν 0 ⎥
⎢ ⎝ E ⎠ 2 hH ⎝ 4 ⎠ 2 EhH 3 ⎥
⎢ ⎥
⎢ 12 Py ( x − L ) ⎥
⎢ 0 0 −ν
⎣ 2 EhH 3 ⎥⎦
52
Validation de l’élément SHB6 en linéaire
1,1
1
Uy/Ureference du point A
0,9
SHB6 initial
0,8
PRI6
0,7 SHB6
Reference
0,6
0,5
0 200 400 600 800 1000
Nombre d'éléments
Graphe 2. Déplacement normalisé du point A de la plaque en flexion
et cisaillement dans son plan
53
Validation de l’élément SHB6 en linéaire
Longueur L = 12 m
Largeur l = 1,1 m
Épaisseur h = 0,32 m
Module d’Young E = 29x106 Pa
Coefficient de Poisson ν = 0,22
Chargement vertical P =1 N
wA = 0,542 × 10−02 m
z z
y y
x x
P P
L L
A l A l
a) Géométrie et chargement
de la poutre vrillée h b) Géométrie initiale et déformée h
(échelle de déformation 500:1)
Figure 12. Poutre vrillée soumise à un effort tranchant ; un exemple de maillage (12x4x1)x2
Les résultats sont ceux du nœud A. On observe les déplacements de ce nœud suivant z
pour le chargement vertical P. Les points de Gauss sont placés à travers l’épaisseur de la
poutre.
Avec le chargement vertical P, les résultats obtenus sont donnés dans le Tableau 6.
54
Validation de l’élément SHB6 en linéaire
1,2
Uz/Uref du point A
0,8
0,6
Référence
0,4 SHB6
SHB6 initial
0,2 PRI6
0
0 100 200 300 400 500 600
Nombre d'éléments
Ces résultats sont représentés sous forme graphique sur le Graphe 3. Nous nous rendons
compte que l’élément SHB6 initial verrouille un peu lors du passage de ce cas test. L’élément
SHB6 est le plus performant par rapport aux SHB6 initial et PRI6, puisque sa courbe de
convergence est située au-dessus des deux autres quel que soit le nombre d’éléments utilisé
(voir Graphe 3).
Les conditions aux limites sont théoriquement libres. Toutefois pour éviter les
mouvements de corps rigides, nous bloquons le déplacement suivant z en bloquant un des
pôles. Par ailleurs le problème ayant deux plans de symétries, nous ne traiterons qu’un quart
G G
de la géométrie de l’hémisphère. Le problème stipule que 4F = 4 N , donc F = 1 N (voir
Figure 13a).
Rayon moyen R = 10 m
Épaisseur h = 0,04 m
Module d’Young E = 6,825x107 Pa
Coefficient de Poisson ν = 0,3
Chargement F =1 N
55
Validation de l’élément SHB6 en linéaire
sym
sym
x
0
1 A F = −1 z
F= 0 0
0 free
a) Géométrie et chargement b) Géométrie initiale et déformée
(échelle de déformation 1000:1)
Les résultats sont ceux du nœud A. On observe les déplacements de ce nœud suivant
Ox. Les résultats obtenus sont représentés dans le Tableau 8.
Nous constatons que les trois éléments ont un blocage sévère. Bien que l’élément SHB6
converge très lentement vers la solution de référence, il reste cependant meilleur que le SHB6
initial et bien meilleur que le PRI6.
Le problème précédent a été repris en mélangeant les éléments SHB6 et SHB8PS. Les
éléments SHB6 sont successivement placés au sommet, (loin du point d’application des
efforts), et aux trois coins du quart d’hémisphère. Pour le premier cas, l’angle de maillage en
SHB8PS est de 75°. Pour les autres cas, les SHB8PS sont disposés sur 60°, 75° ou 80°. Pour
56
Validation de l’élément SHB6 en linéaire
le premier cas où les SHB6 sont mis au sommet de l’hémisphère (voir Figure 14), on constate
une convergence très similaire à celle de la version améliorée de l’élément SHB8PS.
z
x y
y
sym
sym
x
z
0
1 A F = −1 75°
F= 0 0
0 free
a) Géométrie et chargement b) Géométrie initiale et déformée
(échelle de déformation 20:1)
Les résultats obtenus pour le premier cas sont représentés dans le Tableau 9.
z
x y
y
sym
sym
x
z
0
1 A F = −1 60°
F= 0 0
0 free
a) Géométrie et chargement b) Géométrie initiale et déformée
(échelle de déformation 50:1)
57
Validation de l’élément SHB6 en linéaire
sym
sym
x y
0 z
1 A F = −1 75°
F= 0 0
0 free
c) Géométrie et chargement d) Géométrie initiale et déformée
(échelle de déformation 10:1)
sym
sym
0 z y
1 A F = −1 80°
F= 0 0
0 free
e) Géométrie et chargement f) Géométrie initiale et déformée
(échelle de déformation 10:1)
Figure 15. Maillage mixte de l’hémisphère pincé (SHB6 aux 3 sommets)
Pour les autres cas, on observe que la convergence est assez rapide lorsque l’angle de
maillage en SHB8PS est supérieur ou égal à 75° (voir Tableau 10 et Graphe 4). Lorsque les
SHB6 occupent une large part du maillage, on revient à la convergence lente des SHB6 seuls.
Ce test montre aussi que la convergence lente observée est probablement due à la très pauvre
base des éléments finis triangulaire linéaires. Il montre également que l’élément SHB6 se
mélange naturellement aux éléments coques massifs à 8 nœuds SHB8PS. Il permet, donc, de
faire des calculs sur toutes les géométries de pièces, ce qui n’est pas possible avec uniquement
des éléments hexaédriques.
58
Validation de l’élément SHB6 en linéaire
1,2
Ux/Ureference du point A
0,8
0,6
alpha = 60°
0,4
alpha = 75°
0,2 alpha = 80°
Reference
0
0 1000 2000 3000 4000 5000 6000 7000 8000
Nombre d'éléments
Graphe 4. Déplacement normalisé du point A suivant Ox ;
maillage mixte SHB8PS et SHB6 aux trois sommets de l’hémisphère pincé
Les deux extrémités de la coque sont fermées par un diaphragme infiniment rigide. Un
huitième du cylindre est maillé grâce aux symétries du problème. La coque est soumise à un
JG
effort ponctuel au point A, P = 1 (Figure 16).
Les données de géométrie, chargement et de matériau pour ce cas test sont listées dans
le Tableau 11.
Longueur L = 600
Rayon moyen R = 300
Épaisseur h=3
Module d’Young E = 3x106
Coefficient de Poisson ν = 0,3
Chargement P =1
59
Validation de l’élément SHB6 en linéaire
Ri
m gi
P/4 sy d
sym P/4 di
ap
hr
ag
A A m
Rigid
h
h
d
iaphr
sym
z
sy
x
m
agm
y
z
y
R
2
x
L/
m
R
sy
sym L/2
b) Géométrie initiale et déformée
a) Géométrie et chargement (échelle de déformation 2x107:1)
60
Validation de l’élément SHB6 en linéaire
1,2
Uz/Uref du point A
0,8
0,6
PRI6
0,4
SHB6
0,2
initial
SHB6
0
0 2500 5000 7500 10000 12500 15000 17500
Nombre d'éléments
Graphe 5. Déplacement normalisé du point A suivant Oz
de la coque cylindrique pincée avec diaphragmes
L’élément SHB6 converge vers la solution de référence nettement mieux que le SHB6
initial et le PRI6 (voir Graphe 5).
Le contour bas circulaire de la plaque est encastré, il permet donc de supprimer les
mouvements de corps rigides. La plaque est soumise à un effort ponctuel au centre P = 1 N .
On remarquera, par ailleurs la symétrie du problème. On ne modélisera qu'un quart de la
plaque et on considérera les conditions de symétries. La charge appliquée à ce quart de plaque
circulaire sera égale au quart de la charge totale (voir Figure 17).
clamped
z x
y
clamped
h
sym
m
h
sy
R A
sym
m
sy
R P/4
A
a) Géométrie b) Géométrie initiale et déformée
et chargement P/4
(échelle de déformation 105:1)
Figure 17. Géométrie, chargement et déformée de la plaque circulaire soumise à une force
ponctuelle ; un exemple de maillage 3x(4x4x1)x2
61
Validation de l’élément SHB6 en linéaire
Les données de géométrie, de chargement et de matériau ce cas test sont listées dans le
Tableau 13.
Tableau 13. Données de géométrie, chargement et de matériau
de la plaque circulaire soumise à un effort ponctuel
Rayon moyen R = 10 m
Épaisseur h = 0,5 m
Module d’Young E = 10x106 Pa
Coefficient de Poisson ν = 0,25
Chargement P =1 N
Sur ce cas test, l’élément SHB6 se comporte de façon très satisfaisante. Il converge
mieux vers la solution de référence que les deux autres éléments SHB6 initial et PRI6 (voir
Graphe 6).
62
Validation de l’élément SHB6 en linéaire
1,2
Uz/Uref du point A
0,8
0,6
Référence
0,4 SHB6 initial
PRI6
0,2
SHB6
0
0 100 200 300 400 500
Nombre d'éléments
Graphe 6. Déplacement normalisé du point A de la plaque circulaire
soumise à un effort ponctuel
Nous cherchons à déterminer les différentes fréquences d’une poutre dont une extrémité
est encastrée, comprises entre 0 et 150 Hz (voir Figure 18).
Les données géométriques et matériau du problème sont regroupées dans le Tableau 15.
Longueur L=1m
Largeur b = 0,1 m
Épaisseur h = 0,01 m
Module d’Young E = 2,1x1011 Pa
Coefficient de Poisson ν = 0,3
Densité volumique ρ = 7800 kg/m3
La solution analytique des fréquences de vibration en flexion de la poutre est telle que :
λk EI
ω k
= avec cos λk cosh λk +1 = 0
l² ρS
Les solutions analytiques des fréquences de vibration de flexion obtenues dans les deux
directions classées dans l’ordre croissant nous donnent le tableau suivant :
63
Validation de l’élément SHB6 en linéaire
SHB6
f1ref=8.4 f2ref=52.5 f3ref=83.8 f4ref=147.1
Maillage
f1 f1/f1ref f2 f2/f2ref f3 f3/f3ref f4 f4/f4ref
(30x3x1)x2 = 180 9,559 1,14 59,61 1,14 97,06 1,16 165,60 1,13
(40x4x1)x2 = 320 9,184 1,09 57,15 1,09 91,80 1,10 159,20 1,08
(80x8x1)x2 = 1280 8,652 1,03 54,10 1,03 85,87 1,02 151,00 1,03
(120x12x1)x2 = 2880 8,538 1,02 53,40 1,02 84,57 1,01 149,20 1,01
PRI6
f1ref=8.4 f2ref=52.5 f3ref=83.8 f4ref=147.1
Maillage
f1 f1/f1ref f2 f2/f2ref f3 f3/f3ref f4 f4/f4ref
(30x3x1)x2 = 180 19,2 2,29 96,40 1,84 120,40 1,44 337,90 2,30
(40x4x1)x2 = 320 15,8 1,88 91,30 1,74 99,10 1,18 277,70 1,89
(80x8x1)x2 = 1280 11,2 1,33 70,30 1,34 85,60 1,02 196,90 1,34
(120x12x1)x2 = 2880 10 1,19 62,80 1,20 84,40 1,01 175,90 1,20
SHB6 initial
f1ref=8.4 f2ref=52.5 f3ref=83.8 f4ref=147.1
Maillage
f1 f1/f1ref f2 f2/f2ref f3 f3/f3ref f4 f4/f4ref
(30x3x1)x2 = 180 12,3 1,46 76,80 1,46 97,10 1,16 213,90 1,45
(40x4x1)x2 = 320 11,1 1,32 69,04 1,32 91,80 1,10 192,80 1,31
(80x8x1)x2 = 1280 9,28 1,10 58,10 1,11 85,90 1,03 162,70 1,11
(120x12x1)x2 = 2880 8,848 1,05 55,40 1,06 84,50 1,01 155,10 1,05
1er mode de flexion suivant l’épaisseur 2nd mode de flexion suivant l’épaisseur
z
z
y
y
x x
L b L h
h b
64
Validation de l’élément SHB6 en linéaire
1er mode de flexion suivant la largeur 3ème mode de flexion suivant l’épaisseur
z z
y y
x x
h h
L b L b
Ce test représente un calcul de stabilité d’une enveloppe cylindrique mince libre à ses
extrémités soumises à une pression externe. On calcule les charges critiques conduisant au
flambement élastique d’Euler. La matrice de rigidité géométrique utilisée dans la résolution
du problème aux valeurs propres est celle qui est due aux contraintes initiales.
La charge critique et le mode propre obtenus sont comparés à une solution de référence
analytique.
Les données géométriques et matériau du cas test sont listées dans le Tableau 17.
Hauteur L=2m
Rayon moyen R=2m
Épaisseur e = 0,02 m
Module d’Young E = 2x1011 Pa
Coefficient de Poisson ν = 0,3
Les conditions aux limites et chargement du cas test sont (voir Figure 19) :
* Pression uniformément répartie de Pcr = 1 Pa sur la partie cylindrique.
* Conditions de symétrie : sur AB : Uz = 0 ; sur BC : Ux = 0 ; sur DA : Uy = 0.
65
Validation de l’élément SHB6 en linéaire
z
y
x
R Pcr C
sym
L/2
sym
B
e
A sym
Figure 19. Géométrie, chargement et conditions aux limites du cylindre sous pression
externe ; un exemple de maillage (20x30x1)x2
66
Validation de l’élément SHB6 en linéaire
11
10
Référence
9
8 PRI6 (mode 2)
7 SHB6 initial (mode 2)
P/Pcr
6
SHB6 (mode 2)
5
4
3
2
1
0
1200 1700 2200 2700 3200
Nombre d'éléments
z z z
y y y
x x x
R C R C R C
D D D
sym
sym
sym
B B
B
e e e
A sym A sym A sym
Les résultats montrent que l’élément SHB6 converge bien vers la solution de référence.
67
2. Éléments coques volumiques SHB15 et SHB20
Dans ce paragraphe, nous allons présenter les modélisations des éléments finis coques
volumiques quadratiques SHB15 et SHB20.
L’élément SHB15 est un prisme à quinze nœuds purement tridimensionnel avec trois
degrés de liberté en déplacement à chaque nœud, et il a également une direction privilégiée
appelée « épaisseur » qui est normale au plan moyen du prisme. L’intégration numérique
réduite utilisant 3x5 points de Gauss est employée. L’intégration dans le plan utilise 3 points
et celle à travers l’épaisseur s’appuie sur cinq points de Gauss (voir Figure 21).
L’élément SHB20 est un hexaèdre à vingt nœuds purement tridimensionnel avec trois
degrés de liberté en déplacement à chaque nœud, et il a aussi une direction privilégiée appelée
« épaisseur » qui est normale au plan moyen de l’hexaèdre. L’intégration numérique réduite
utilisant 4x5 points de Gauss est employée. L’intégration dans le plan utilise 4 points et celle
à travers l’épaisseur s’appuie sur cinq points de Gauss (voir Figure 22)
On montre que ces éléments finis n’ont pas besoin de stabilisation et, pour l’instant,
nous n’avons pas utilisé de projections. Afin d’évaluer leurs performances, une série de cas
tests linéaires standards de la littérature leurs sera appliquée. Des résultats relatifs à la rapidité
de convergence seront aussi donnés avec des comparaisons avec d’autres éléments finis 3D.
L’élément SHB15 est formulé dans les axes locaux du plan moyen. Les directions x, y, z
(ou encore x1 , x2 , x3 ) sont respectivement parallèles aux axes ξ ,η , ζ . La Figure 21 représente
la géométrie d’un élément de référence SHB15 et ses points d’intégration.
10 15 14
10
11 13
7 5
9
8 9 η
15
4 7
12 14 O
3 6
13 5
12 2
11 1 1 6
8
ξ 2 4
3
15
xi = xiI N I (ξ ,η , ζ ) = ∑ xiI N I (ξ ,η , ζ ) (64)
I =1
Les mêmes fonctions de forme sont utilisées pour définir le champ de déplacement de
l’élément ui en termes des déplacements nodaux uiI :
15
ui = uiI N I (ξ ,η , ζ ) = ∑ uiI N I (ξ ,η , ζ ) (65)
I =1
ui , j = uiI N I , j (66)
1
ε ij =
2
( ui, j + u j ,i ) (67)
69
Formulation des éléments SHB15 et SHB20 en linéaire
⎛ (ζ − 1) (1 − ξ − η )(ζ + 2ξ + 2η ) / 2 ⎞
⎜ ⎟
⎜ 2ξ (1 − ξ − η )(1 − ζ ) ⎟
⎜ ξ (1 − ζ ) (2ξ − 2 − ζ ) / 2 ⎟
⎜ ⎟
⎜ 2ξη (1 − ζ ) ⎟
⎜ η (2η − 2 − ζ )(1 − ζ ) / 2 ⎟
⎜ ⎟
⎜ 2η (1 − ξ − η ) (1 − ζ ) ⎟
⎜ (1 − ζ ) (1 − ξ − η )
2 ⎟
⎜ ⎟
N (ξ ,η , ζ ) = ⎜ (1 − ζ 2 )ξ ⎟
⎜ ⎟
⎜ (1 − ζ )η
2
⎟
⎜ ( −ζ − 1)(1 − ξ − η )( 2ξ + 2η − ζ ) / 2 ⎟
⎜ ⎟
⎜ 2ξ (1 − ξ − η )(1 + ζ ) ⎟
⎜ ξ (1 + ζ ) (2ξ − 2 + ζ ) / 2 ⎟
⎜ ⎟
⎜ 2ξη (1 + ζ ) ⎟
⎜ η ( 2η − 2 + ζ )(1 + ζ ) / 2 ⎟
⎜ ⎟
⎜ 2 η (1 − ξ − η )(1 + ζ ) ⎟ (68)
⎝ ⎠
ξ ∈ [ 0,1] ; η ∈ [ 0,1 − ξ ] ; ζ ∈ [ −1,1].
L’origine du repère est confondue avec le coin droit du triangle du plan médian de
l’élément (celui qui possède un angle droit, voir Figure 21).
En évaluant l’équation (69) aux nœuds de l’élément, on arrive aux trois systèmes de
quinze équations suivants :
⎧d i = a0 i S + a1i x1 + a2 i x 2 + a3i x 3 + c1i h1 + c2 i h 2 + c3i h 3 +
⎪
⎨ c4 i h 4 + c5 i h 5 + c6 i h 6 + c7 i h 7 + c8i h 8 + c9 i h 9 + c10 i h10 + c11i h11 (70)
⎪i = 1, 2, 3
⎩
70
Formulation des éléments SHB15 et SHB20 en linéaire
⎨ T (71)
⎪⎩ x i = ( xi1 , xi 2 , xi 3 , xi 4 , xi 5 , xi 6 , xi 7 , xi 8 , xi 9 , xi10 , xi11 , xi12 , xi13 , xi14 , xi15 )
⎧ ST = (1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)
⎪
⎪hT = ⎛ 0 − 1 −1 − 1 0 0 0 0 0 0 1 1 1 0 0 ⎞
⎪ 1 ⎜⎝ 2 2 2 2
⎟
⎠
⎪
⎪h T = ⎛ 0 0 0 − 1 − 1 − 1 0 0 0 0 0 0 1 1 1 ⎞
⎪ 2 ⎜⎝ 2 2 2
⎟
2⎠
⎪
⎪ T ⎛ 1 1 ⎞
⎪h3 = ⎜⎝ 0 0 0 4 0 0 0 0 0 0 0 0 4 0 0 ⎟⎠
⎪
⎪ T ⎛ 1 1 ⎞
⎪h 4 = ⎜ 0 0 0 − 4 0 0 0 0 0 0 0 0 4 0 0 ⎟
⎪ ⎝ ⎠
⎪ T ⎛ 1 1 1 1 ⎞
⎪h 5 = ⎜ 0 1 0 0 0 1 0 0 1 0 0⎟
⎪⎪ ⎝ 4 4 4 4 ⎠
⎨
⎪h T = ⎛ 0 0 0 1 1 1 0 0 1 0 0 0 1 1 1 ⎞
⎪ 6 ⎜⎝ 4 4 4
⎟
4⎠
⎪ T
⎪h 7 = (1 1 1 1 1 1 0 0 0 1 1 1 1 1 1)
⎪
⎪hT8 = ⎛⎜ 0 − 1 −1 − 1 0 0 0 0 0 0 1 1 1 0 0 ⎞⎟
⎪ ⎝ 4 4 4 4 ⎠
⎪
⎪h T = ⎛ 0 0 0 − 1 − 1 − 1 0 0 0 0 0 0 1 1 1 ⎞
⎪ 9 ⎜⎝ 4 4 4
⎟
4⎠
⎪ (72)
⎪h T = ⎛ 0 1 1 1 0 0 0 0 0 0 1 1 1 0 0 ⎞
⎪ 10 ⎜⎝ 2 2 2 2
⎟
⎠
⎪
⎪h T = ⎛ 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 ⎞
⎪⎪⎩ 11 ⎜⎝ 2 2 2 2⎠
⎟
∂N
bTi = N ,i ( 0 ) = i = 1, 2, 3 (73)
∂xi ξ =η =ζ = 0
En fait, ce calcul est pratique car nous verrons par la suite qu’il est nécessaire de
déterminer la quantité suivante :
⎛ T 11 ⎞
⎜ b j + ∑ hα , j γ α ⎟ = N , j (ξ ,η , ζ )
T T
⎝ α =1 ⎠
71
Formulation des éléments SHB15 et SHB20 en linéaire
b Tj = N T, j (0, 0, 0) = csteT
⎛ ∂ N1 ∂ N2 ∂ N3 ∂ N15 ⎞
bTi = NT,i ( 0 ) = ⎜ . . . ⎟
⎝ ∂xi ∂xi ∂xi ∂xi ⎠ ξ =η =ζ =0
∂ N I ⎛ ∂ N I ∂ξ ∂ N I ∂η ∂ N I ∂ζ ⎞ ⎛ ∂ NI ∂ NI ∂ NI ⎞
=⎜ + + ⎟⎟ = ⎜ j1 j + j2 j + j3 j ⎟
∂x j ⎜⎝ ∂ξ ∂x j ∂η ∂x j ∂ζ ∂x j ⎠ ⎝ ∂ξ ∂η ∂ζ ⎠
avec I = 1, 2,...,15 et j = 1, 2,3
⎛ ∂ξ ∂ξ ∂ξ ⎞
⎜ ⎟
⎜ ∂x ∂y ∂z ⎟
⎛ j11 j12 j13 ⎞
⎜ ∂η ∂η ∂η ⎟ ⎜ ⎟
F −1 = ⎜ ⎟ = j21 j22 j23 ⎟
⎜ ∂x ∂y ∂z ⎟ ⎜⎜
j31 j32 j33 ⎠⎟
⎜ ∂ζ ∂ζ ∂ζ ⎟ ⎝
⎜ ⎟
⎝ ∂x ∂y ∂z ⎠
72
Formulation des éléments SHB15 et SHB20 en linéaire
(1 − ζ ) (1 − ζ ) (1 − ξ − η )
N1,ξ = ( 4ξ + 4η + ζ − 2) N1,η = ( 4ξ + 4η + ζ − 2) N1,ζ = ( 2ξ + 2η + 2ζ − 1)
2 2 2
N 2,ξ = 2 (1 − ζ )(1 − 2ξ − η ) N 2,η = −2 (1 − ζ ) ξ N 2,ζ = −2ξ (1 − ξ − η )
(1 − ζ ) ξ
N 3,ξ = ( −2 + 4ξ − ζ ) N 3,η = 0 N 3,ζ = (1 − 2ξ + 2ζ )
2 2
N 4,ξ = 2 (1 − ζ )η N 4,η = 2 (1 − ζ ) ξ N 4,ζ = −2ξη
(1 − ζ ) η
N 5,ξ = 0 N 5,η = ( −2 + 4η − ζ ) N 5,ζ = (1 − 2η + 2ζ )
2 2
N 6,ξ = −2 (1 − ζ )η N 6,η = 2 (1 − ζ )(1 − ξ − 2η ) N 6,ζ = −2η (1 − ξ − η )
N 7,ξ = −1 + ζ N 7,η = −1 + ζ N 7,ζ = −2ζ (1 − ξ − η )
2 2
(1 + ζ ) (1 + ζ ) (1 − ξ − η )
N10,ξ = ( 4ξ + 4η − ζ − 2) N10,η = ( 4ξ + 4η − ζ -2 ) N10,ζ = ( −2ξ − 2η + 2ζ + 1)
2 2 2
N11,ξ = 2 (1 + ζ )(1 − 2ξ − η ) N11,η = −2 (1 + ζ ) ξ N11,ζ = 2ξ (1 − ξ − η )
(1 + ζ ) ξ
N12,ξ = ( −2 + 4ξ + ζ ) N12,η = 0 N12,ζ = ( −1 + 2ξ + 2ζ )
2 2
N13,ξ = 2 (1 + ζ )η N13,η = 2 (1 + ζ ) ξ N13,ζ = 2ξη
(1 + ζ ) η
N14,ξ = 0 N14,η = ( −2 + 4η + ζ ) N14,ζ = ( −1 + 2η + 2ζ )
2 2
N15,ξ = −2 (1 + ζ )η N15,η = 2 (1 + ζ )(1 − 2η − ξ ) N15,ζ = 2η (1 − ξ − η )
⎡ ⎤
⎢ −1 2 −1 0 0 0 −1 1 0 − 1 2 − 1 0 0 0 ⎥
⎢ ⎥
bTi = ( j1i j2i j3i ) ⎢ −1 0 0 0 −1 2 −1 0 1 − 1 0 0 0 − 1 2 ⎥
⎢ 1 1 ⎥
⎢− 0 0 0 0 0 0 0 0 0 0 0 0 0⎥
⎣ 2 2 ⎦
De plus, on peut vérifier par des considérations algébriques que les conditions
d’orthogonalité suivantes sont satisfaites :
73
Formulation des éléments SHB15 et SHB20 en linéaire
⎧b Ti ⋅ hα = 0; b Ti ⋅ S = 0; b Ti ⋅ x j = δ ij
⎪
⎪ h T ⋅ S = 0; h T ⋅ S = 0; h T ⋅ S = 1 ; h T ⋅ S = 0; h T5 ⋅ S = 4; h T6 ⋅ S = 4;
⎪ 1 2 3
2
4
⎪ T
⎪h 7 ⋅ S = 12; h 8 ⋅ S = 0 h 9 ⋅ S = 0 h10 ⋅ S = 4 ⋅S = 4
T T T T
h11
⎪ ⎡ 1 1 5 1 ⎤
⎪ ⎢ 3 −2 0 4 0 0 0
2 4
0 0⎥
⎪ ⎢ ⎥
⎪ ⎢− 1 1 1 5
⎪ 3 0 0 0 0 0 0⎥
⎢ 2 4 4 2 ⎥
⎪ ⎢
⎪ ⎢ 0 1 1 1 1 1 1 ⎥⎥
⎪ 0 0 0 0
⎢ 8 8 8 2 4 4⎥
⎪ ⎢ 1 ⎥
⎪ 1 1 1 1
⎢ 0 0 0 0 0 0⎥
⎪ ⎢ 4 4 8 8 8 ⎥
⎪ ⎢ 1 13 1 5 1⎥
⎪ ⎢ 0 0 0 3 0 0 ⎥
⎪⎪ 8 4 8 2 4⎥
⎨ ⎢
⎪h T ⋅ h = ⎢ 0 0
1
0
1 13
3 0 0
1 5⎥
⎪ α β ⎢ 8 8 4 4 2⎥
⎪ ⎢ ⎥
⎢ 0 1
⎪ 0 0 3 3 12 0 0 4 4⎥
⎪ ⎢ 2 ⎥
⎪ ⎢ ⎥
⎪ ⎢ 5 1
0
1
0 0 0
9 1
0 0⎥
⎪ ⎢ 2 4 8 4 8 ⎥
⎪ ⎢ 1 5 1 1 9 ⎥
⎪ ⎢ 0 0 0 0 0 0⎥
⎪ ⎢ 4 2 8 8 4 ⎥ (74)
⎪ ⎢ 1 5 1 1⎥
⎪ ⎢ 0 0 0 4 0 0 3 ⎥
⎢ 4 2 4 2⎥
⎪
⎪ ⎢ 1 1 5 1 ⎥
⎢ 0 0 0 4 0 0 3⎥
⎪ ⎢⎣ 4 4 2 2 ⎥⎦
⎪
⎪ α , β = 1, 2,...,11 ; i, j = 1, 2, 3
⎪⎪⎩ δ désigne le symbol de Kronecker
À ce stade, on peut déterminer les constantes inconnues qui interviennent dans l’écriture
(69) du champ de déplacement en multipliant scalairement l’équation (70) par b Tj , S T et hαT ,
respectivement, et en utilisant les relations d’orthogonalité (74). On obtient ainsi :
⎪
⎪⎩ = 0 + a1i b j ⋅ x1 + a2 i b j ⋅ x 2 + a3i b j ⋅ x 3 + 0 = a ji
T T T
74
Formulation des éléments SHB15 et SHB20 en linéaire
⎧
⎪S T ⋅ d = a S T ⋅ S + a S T ⋅ x + a S T ⋅ x + a S T ⋅ x + c S T ⋅ h + c S T ⋅ h + ... + c S T ⋅ h
⎪ i 0i 1i 1 2i 2 3i 3 1i 1 2i 2 11i 11
⎪ 1
⎨ = 15 a0 i + ( S ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0 c2 i + c3 i + 0 c4 i + 4 c5 i + 4 c6 i + 12 c7 i + 0 c8 i + 0 c9 i + 4 c10 i + 4 c11i
T T
⎪ 2
⎪
⎩
1 T
( ( ) ) 1 ⎛1
⎪ a0 i = 15 S − S ⋅ x j ⋅ b j ⋅ d i − 15 ⎜ 2 c3i + 4c5 i + 4c6 i + 12c7 i + 4c10 i + 4c11i ⎟
T T
⎝
⎞
⎠
⎧h1T ⋅ di = a0i h1T ⋅ S + a1i h1T ⋅ x1 + a2i h1T ⋅ x 2 + a3i h1T ⋅ x3 + c1i h1T ⋅ h1 + c2i h1T ⋅ h 2 + ... + c11i h1T ⋅ h11
⎪
⎨ 1 1 5 1
⎩
( )(T
)
⎪= h1 ⋅ x j b j ⋅ di + 3c1i − c2i + 0c3i + c4i + 0c5i + 0c6i + 0c7 i + c8i + c9i + 0c10i + 0c11i
T
2 4 2 4
2 4 4 2
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 3 i 0i 3 1i 3 1 2i 3 2 3i 3 3 1i 3 1 2i 3 2 11i 3 11
⎪ 1 1 1 1 1 1 1
⎨ = a0 i + ( h 3 .x j )( b j .d i ) + 0 c1i + 0 c2 i + c3i + 0 c4 i + c5 i + c6 i + c7 i + 0 c8 i + 0 c9 i + c10 i + c11i
T T
⎪ 2 8 8 8 2 4 4
⎪⎛ T 1 T ⎞ ⎡⎛ T 1 T ⎞ ⎤ T 13 1 1 1 7 7
⎪⎜ h 3 − S ⎟ ⋅ d i = ⎢ ⎜ h 3 − S ⎟ ⋅ x j ⎥ ( b j ⋅ d i ) + c3 i − c5 i − c6 i + c7 i + c10 i + c11i
⎩⎝ 30 ⎠ ⎣⎝ 30 ⎠ ⎦ 120 120 120 10 60 60
⎩ 4 4 8 8 8
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 5 i 0i 5 1i 5 1 2i 5 2 3i 5 3 1i 5 1 2i 5 2 11i 5 11
⎪ 1 13 1 5 1
⎨ = 4 a0 i + ( h 5 ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0c2 i + c3i + 0 c4 i + c5 i + c6 i + 3c7 i + 0 c8 i + 0c9 i + c10 i + c11i
T T
⎪ 8 4 8 2 4
⎪⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T 1 131 113 1 43 49
⎪⎜ h 5 − S ⎟ ⋅ d i = ⎢ ⎜ h 5 − S ⎟ ⋅ x j ⎥ ( b j ⋅ d i ) − c3 i + c5 i − c 6 i − c7 i + c10 i − c11i
⎩⎝ 15 ⎠ ⎣⎝ 15 ⎠ ⎦ 120 60 120 5 30 60
75
Formulation des éléments SHB15 et SHB20 en linéaire
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 6 i 0i 6 1i 6 1 2i 6 2 3i 6 3 1i 6 1 2i 6 2 11i 6 11
⎪ 1 1 13 1 5
⎨ = 4 a0 i + ( h 6 ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0c2 i + c3i + 0 c4 i + c5 i + c6 i + 3c7 i + 0 c8 i + 0c9 i + c10 i + c11i
T T
⎪ 8 8 4 4 2
⎪⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T 1 113 131 1 49 43
⎪⎜ h 6 − S ⎟ ⋅ d i = ⎢ ⎜ h 6 − S ⎟ ⋅ x j ⎥ ( b j ⋅ d i ) − c3 i − c5 i + c 6 i − c7 i − c10 i + c11i
⎩⎝ 15 ⎠ ⎣⎝ 15 ⎠ ⎦ 120 120 60 5 60 30
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 7 i 0i 7 1i 7 1 2i 7 2 3i 7 3 1i 7 1 2i 7 2 11i 7 11
⎪ 1
⎨ = 12 a0 i + ( h 7 ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0c2 i + c3 i + 0 c4 i + 3c5 i + 3c6 i + 12 c7 i + 0 c8 i + 0 c9 i + 4 c10 i + 4 c11i
T T
⎪ 2
⎪⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T 1 1 1 12 4 4
( )
⎪⎜ h 7 − S ⎟ ⋅ d i = ⎢ ⎜ h 7 − S ⎟ ⋅ x j ⎥ b j ⋅ d i + c3i − c5 i − c6 i + c7 i + c10 i + c11i
5 ⎠ 5 ⎠ 10 5 5 5 5 5
⎩⎝ ⎣⎝ ⎦
2 4 8 4 8
⎩ 4 2 8 8 4
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 10 i 0 i 10 1i 10 1 2 i 10 2 3 i 10 3 1i 10 1 2 i 10 2 11i 10 11
⎪ 1 5 1 1
⎨ = 4 a0 i + ( h10 ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0 c2 i + c3 i + 0 c4 i + c5 i + c6 i + 4 c7 i + 0c8 i + 0 c9 i + 3c10 i + c11i
T T
⎪ 4 2 4 2
⎪⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T 7 43 49 4 29 17
⎪⎜ h10 − 15 S ⎟ ⋅ d i = ⎢ ⎜ h10 − 15 S ⎟ ⋅ x j ⎥ ( b j ⋅ d i ) + 60 c3i + 30 c5 i − 60 c6 i + 5 c7 i + 15 c10 i − 30 c11i
⎩⎝ ⎠ ⎣⎝ ⎠ ⎦
⎧
⎪h T ⋅ d = a h T ⋅ S + a h T ⋅ x + a h T ⋅ x + a h T ⋅ x + c h T ⋅ h + c h T ⋅ h + ... + c h T ⋅ h
⎪ 11 i 0 i 11 1i 11 1 2 i 11 2 3 i 11 3 1i 11 1 2 i 11 2 11i 11 11
⎪ 1 1 5 1
⎨ = 4 a0 i + ( h11 ⋅ x j )( b j ⋅ d i ) + 0 c1i + 0 c2 i + c3i + 0 c4 i + c5 i + c6 i + 4 c7 i + 0 c8 i + 0 c9 i + c10 i + 3c11i
T T
⎪ 4 4 2 2
⎪⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T 7 49 43 4 17 29
⎪⎜ h11 − 15 S ⎟ ⋅ d i = ⎢ ⎜ h11 − 15 S ⎟ ⋅ x j ⎥ ( b j ⋅ d i ) + 60 c3 i − 60 c5 i + 30 c6 i + 5 c7 i − 30 c10 i + 15 c11i
⎩⎝ ⎠ ⎣⎝ ⎠ ⎦
76
Formulation des éléments SHB15 et SHB20 en linéaire
⎡ 1 1 5 1 ⎤ ⎡
⎢
( ( ) )
h1T − h1T ⋅ x j bTj ⋅ d i ⎤
⎥
⎢3 0 0 0 0 0 0 ⎥
⎢
⎢1
2 4
1
2
1
4
5
⎥ ⎢
⎢ ( ( ) )
hT2 − hT2 ⋅ x j bTj ⋅ d i ⎥
⎥
3 0 0 0 0 0 0 ⎥ ⎢ ⎥
⎢2 4 4 2 ⎥ ⎢ ⎛ ⎛ hT − 1 ST ⎞ − ⎡⎛ hT − 1 S T ⎞ ⋅ x ⎤ bT ⎞ ⋅ d ⎥
⎢ ⎥
⎢0 13 1 1 1 7 7 ⎥ ⎢ ⎜⎝ ⎜⎝ 3 30 ⎟⎠ ⎢⎣⎜⎝ 3 30 ⎟⎠ j ⎥⎦ j ⎟⎠ i ⎥
0 0 − − 0 0
⎢ 60 ⎥ ⎡ c1i ⎤ ⎢ ⎥
⎢1 1
120
1
120 120 10
1 1
60
⎥ c⎢ ⎥ ⎢ h (
T
4 − h (
T
4 ⋅ x j b ) )
T
j ⋅ d i
⎥
⎢ 0 0 0 0 0 0 ⎥ ⎢ 2i ⎥ ⎢ ⎥
⎢4 4 8 8 8 ⎢
⎥ 3i c ⎥ ⎢ ⎛ ⎛ 4 ⎞ ⎡ ⎛ 4 ⎞ ⎤ ⎞ ⎥
⎢ ⎥ ⎢ ⎜ ⎜ h T
− S T
−
⎟ ⎜ h T
− S T
⎟ ⋅ x b T
⎟ ⋅ d ⎥
⎢ 1 131 113 1 43 49 ⎥ 15 ⎠ ⎢⎣⎝ 15 ⎠ ⎥⎦ ⎠
5 5 j j i
⎢0 0 − 0 − − 0 0 − ⎥ ⎢ c4i ⎥ ⎢ ⎝ ⎝ ⎥
⎢ 120 60 120 5 30 60 ⎥ ⎢ c ⎥ ⎢ ⎥
⎛⎛ 4 ⎞ ⎡⎛ 4 ⎞ ⎤ ⎞
43 ⎥ ⎢ ⎥ ⎢ ⎜ ⎜ hT6 − ST ⎟ − ⎢⎜ hT6 − ST ⎟ ⋅ x j ⎥ bTj ⎟ ⋅ d i ⎥
5i
⎢ 1 113 131 1 49
⎢0 0 − 0 − − 0 0 − ⎢ c6i ⎥ = ⎢ 15 ⎠ ⎣⎝ 15 ⎠ ⎦ ⎠ ⎥
120 120 60 5 60 30 ⎥ ⎢ ⎥ ⎢ ⎝ ⎝ ⎥
⎢ ⎥ c
⎢0 1 1 1 12 4 4 ⎥ ⎢ 7 i ⎥ ⎢ ⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T ⎞ ⎥
0 0 − − 0 0 ⎢ c ⎥ ⎢ ⎜ ⎜ h 7 − S ⎟ − ⎢⎜ h 7 − S ⎟ ⋅ x j ⎥ b j ⎟ ⋅ d i ⎥
⎢ 10 5 5 5 5 5 ⎥ 8i ⎝ 5 ⎠ ⎣⎝ 5 ⎠ ⎦ ⎠
⎢ ⎥ ⎢ c ⎥ ⎢⎢ ⎝ ⎥
⎥
⎢5
⎢2
1
4
0
1
8
0 0 0
9
4
1
8
0 0 ⎥ ⎢ 9i ⎥
⎥ ⎢ c10i ⎥ ⎢
( (
T
) )
h8 − h8 ⋅ x j b j ⋅ d i
T T
⎥
⎢ ⎥ ⎢ ⎥
⎢1
⎢
5
0
1
0 0 0
1 9
0 0 ⎥
⎥ ⎢c ⎥ ⎢
⎣ 11i ⎦
⎢
h (
T
9 − h (
T
9 ⋅ x j b ) )
T
j ⋅ d i ⎥
⎥
⎢4 2 8 8 4 ⎥ ⎢ ⎥
⎢ ⎛ ⎛ 4 ⎞ ⎛ ⎛ 4 ⎞ ⎞ ⎞
7 43 49 4 29 17 ⎥ ⎢⎜ ⎜ h10 − S ⎟ − ⎜ ⎜ h10 − S ⎟ ⋅ x j ⎟ b j ⎟ ⋅ di ⎥
T T T T T
⎢0 0 0 − 0 0 − ⎥ 15 ⎠ ⎝ ⎝ 15 ⎠
⎢ 60 30 60 5 15 30 ⎥ ⎢⎝ ⎝ ⎠ ⎠ ⎥
⎢ 29 ⎥ ⎢ ⎛⎛ T 4 T ⎞ ⎛⎛ T 4 T ⎞ ⎥
7 49 43 4 17
⎢ ⎞ T⎞ ⎥
0 0 0 − 0 0 − ⎜ h − S − h − S ⋅ x b ⎟ ⋅ d
⎢⎢⎣ 60 60 30 5 30 15 ⎥⎥⎦ ⎜ 11
⎢⎢⎣ ⎝ ⎝
⎟ ⎜ ⎜ 11
15 ⎠ ⎝ ⎝ 15 ⎠
⎟ j⎟ j
⎠ ⎠ ⎥⎥⎦
i
ou encore :
⎡
⎢
( h1T − ( h1T ⋅ x j ) bTj ) ⋅ di ⎤
⎥
⎡17 ⎤
⎢2
⎢
0 0 −8 0 0 0 −9 0 0 0 ⎥⎢
⎥⎢
( 2 2 j j) i
h T
− ( h T
⋅ x ) b T
⋅ d ⎥
⎥
⎢0 17 ⎥ ⎢ ⎛ ⎛ 1 ⎞ ⎡ ⎛ 1 ⎞ ⎤ ⎞ ⎥
0 −8 0 0 0 0 −9 0 0
⎥ ⎢ ⎜ ⎜ h 3 − 30 S ⎟ − ⎢⎜ h 3 − 30 S ⎟ ⋅ x j ⎥ b j ⎟ ⋅ d i ⎥
T T T T T
⎢ 2
⎡ c1i ⎤ ⎢ ⎢⎝⎝ ⎠ ⎣⎝ ⎠ ⎦ ⎠ ⎥
256 36 36 58 58 ⎥ ⎢
⎢c ⎥ ⎢ 0
⎢ 2i ⎥ ⎢
0
17
0
17 17
2 0 0 −
17
− ⎥
17 ⎥ ⎢
( h 4 − ( h 4 ⋅ x j ) b j ) ⋅ di
T T T ⎥
⎥
⎢ c3i ⎥ ⎢ −8 −8 0 24 0 0 0 8 8 0 0 ⎥ ⎢ ⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T ⎞ ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎜ ⎜ h 5 − S ⎟ − ⎢⎜ h 5 − S ⎟ ⋅ x j ⎥ b j ⎟ ⋅ d i ⎥
⎢ c4i ⎥ ⎢ 36 316 146 324 171 ⎥ ⎢ ⎝ ⎝ 15 ⎠ ⎣⎝ 15 ⎠ ⎦ ⎠ ⎥
⎢ c5i ⎥ ⎢ 0 0 0 1 0 0 − − ⎢ ⎥
⎥
187 ⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎞
⎥ ⎢ ⎜ ⎜ h 6 − S ⎟ − ⎢⎜ h 6 − S ⎟ ⋅ x j ⎥⎤ bTj ⎟ ⋅ d i ⎥
17 187 187 187
⎢ ⎥ ⎢
⎢ c6i ⎥ = ⎢ 0 0
36
0
146 316
1 0 0 −
171
−
324 ⎥ ⎢ ⎝ ⎝ 15 ⎠ ⎣⎝ 15 ⎠ ⎦ ⎠ ⎥
⎢c ⎥ ⎢ 17 187 187 187 187 ⎥ ⎢ ⎥
⎢ 7i ⎥ ⎢ ⎥ ⎢ ⎛ ⎛ hT − 4 S T ⎞ − ⎡⎛ hT − 4 S T ⎞ ⋅ x ⎤ b T ⎞ ⋅ d ⎥
3 3 3 ⎜ ⎜ 7 5 ⎟ ⎢⎜ 7 5 ⎟ j ⎥ j ⎟ i ⎥
⎢ c8i ⎥ ⎢ 0 0 2 0 1 1 0 0 − − ⎥ ⎢ ⎝⎝ ⎠ ⎣⎝ ⎠ ⎦ ⎠
⎢c ⎥ ⎢ 2 2 2 ⎥⎢ ⎥
⎢ 9 i ⎥ ⎢ −9 0 0 8 0 0 0 10 0 0 0 ⎥⎢ ( h8 − ( h8 ⋅ x j ) b j ) ⋅ d i
T T T
⎥
⎢ c10i ⎥ ⎢ ⎥⎢ ⎥
⎢ ⎥ ⎢0
⎢⎣ c11i ⎥⎦ ⎢
−9 0 8 0 0 0 0 10 0 0 ⎥⎢
⎢
( hT9 − ( hT9 ⋅ x j ) bTj ) ⋅ d i ⎥
⎥
58 324 171 3 505 585 ⎥
⎢0 0 − 0 − − − 0 0 ⎥ ⎢⎛ ⎛ T 4 T ⎞ ⎛ ⎛ T 4 T ⎞ ⎞ T ⎞ ⎥
⎢ 17 187 187 2 187 374 ⎥ ⎢⎜ ⎜ h10 − S ⎟ − ⎜ ⎜ h10 − S ⎟ ⋅ x j ⎟ b j ⎟ ⋅ d i ⎥
⎝⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎠ ⎥
⎢ 58 171 324 3 585 505 ⎥ ⎢
⎢⎢ 0 0 − 0 − − − 0 0 ⎥ ⎢⎛ ⎛ T 4 T ⎞ ⎛ ⎛ T 4 T ⎞ ⎞ T ⎞ ⎥
⎣ 17 187 187 2 374 187 ⎥⎦ ⎢ ⎜ ⎜ h11 − S ⎟ − ⎜ ⎜ h11 − S ⎟ ⋅ x j ⎟ b j ⎟ ⋅ d i ⎥
⎣⎢⎢ ⎝ ⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎠ ⎦⎥⎥
En posant :
77
Formulation des éléments SHB15 et SHB20 en linéaire
⎡17 ⎤
⎢2 0 0 −8 0 0 0 −9 0 0 0 ⎥
⎢ ⎥
⎢ 0 17 0 −8 0 0 0 0 −9 0 0 ⎥
⎢ 2 ⎥
⎢ 256 36 36 58 58 ⎥⎥
⎢0 0 0 2 0 0 − −
⎢ 17 17 17 17 17 ⎥
⎢ −8 −8 0 24 0 0 0 8 8 0 0 ⎥
⎢ ⎥
⎢ 36 316 146 324 171 ⎥
⎢0 0 0 1 0 0 − −
17 187 187 187 187 ⎥
⎢ ⎥
⎡⎣nαβ ⎤⎦ = ⎢ 36 146 316 171 324 ⎥
0 0 0 1 0 0 − −
⎢ 17 187 187 187 187 ⎥
⎢ 3 3 3 ⎥
⎢0 0 2 0 1 1 0 0 − − ⎥
⎢ 2 2 2 ⎥
⎢ −9 0 0 8 0 0 0 10 0 0 0 ⎥
⎢ ⎥
⎢ 0 −9 0 8 0 0 0 0 10 0 0 ⎥
⎢ 58 324 171 3 505 585 ⎥
⎢0 0 − 0 − − − 0 0 ⎥
⎢ 17 187 187 2 187 374 ⎥
⎢ 58 171 324 3 585 505 ⎥
⎢0 0 − 0 − − − 0 0 ⎥
⎣ 17 187 187 2 374 187 ⎦
α , β = 1, 2,...,11
( ( ) ) ( (
γαT = nα 1 h1T − h1T ⋅ x j bTj + nα 2 hT2 − hT2 ⋅ x j bTj ) )
⎡⎛ ⎞ ⎛⎛ ⎞ ⎤
⎣⎝
1 1
30 ⎠
⎞
⎠ ⎦
( (
+ nα 3 ⎢⎜ hT3 − ST ⎟ − ⎜ ⎜ hT3 − ST ⎟ ⋅ x j ⎟ bTj ⎥ + nα 4 hT4 − hT4 ⋅ x j bTj
30 ⎠ ⎝ ⎝
) )
⎡⎛ 4 ⎞ ⎛⎛ 4 ⎞ ⎞ ⎤ ⎛⎛ 4 ⎞ ⎛⎛ 4 ⎞ ⎞ ⎞
+ nα 5 ⎢⎜ hT5 − ST ⎟ − ⎜ ⎜ hT5 − ST ⎟ ⋅ x j ⎟ bTj ⎥ + nα 6 ⎜ ⎜ hT6 − ST ⎟ − ⎜ ⎜ hT6 − ST ⎟ ⋅ x j ⎟ bTj ⎟
⎣⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎦ ⎝⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎠
⎛⎛ 4 ⎞ ⎛⎛ ⎞ ⎞
⎝⎝
4 ⎞
5 ⎠ ⎠ ⎠
( ( ) ) ( (
+ nα 7 ⎜ ⎜ hT7 − ST ⎟ − ⎜ ⎜ hT7 − ST ⎟ ⋅ x j ⎟ bTj ⎟ + nα 8 hT8 − hT8 ⋅ x j bTj + nα 9 hT9 − hT9 ⋅ x j bTj +
5 ⎠ ⎝⎝
) )
⎛⎛ T 4 T ⎞ ⎛⎛ T 4 T ⎞ ⎞ ⎞ ⎛⎛ T 4 T ⎞ ⎛⎛ T 4 T ⎞ ⎞ ⎞
nα 10 ⎜ ⎜ h10 − S ⎟ − ⎜ ⎜ h10 − S ⎟ ⋅ x j ⎟ bTj ⎟ + nα 11 ⎜ ⎜ h11 − S ⎟ − ⎜ ⎜ h11 − S ⎟ ⋅ x j ⎟ bTj ⎟
⎝⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎠ ⎝⎝ 15 ⎠ ⎝ ⎝ 15 ⎠ ⎠ ⎠
78
Formulation des éléments SHB15 et SHB20 en linéaire
ui = a0 i + (b1T x1 + b T2 x2 + b T3 x3 + γ 1T h1 + γ T2 h2 + γ T3 h3 + γ T4 h4 + γ T5 h5
+ γ T6 h6 + γ T7 h7 + γ T8 h8 + γ T9 h9 + γ 10
T
h10 + γ 11
T
h11 ) ⋅ d i
⎛ 11
⎞
(
ui , j = ⎜ bTj + ∑ hα , j γαT ⎟ ⋅ di = bTj + hα , j γαT ⋅ di ) (77)
⎝ α =1 ⎠
∇s (u ) = B ⋅ d (78)
où :
⎡ u x,x ⎤
⎢ u ⎥
⎢ y,y ⎥ ⎡ d1 ⎤
⎢
∇ s (u ) = ⎢
u z,z
⎥ ,
⎥
d = ⎢⎢d 2 ⎥⎥ (79)
⎢u x , y + u y , x ⎥ ⎢⎣ d3 ⎥⎦
⎢ u x,z + u z,x ⎥
⎢ ⎥
⎣⎢ u y , z + u z , y ⎦⎥
⎡ b1T + hα , x1 γ αT 0 0 ⎤
⎢ ⎥
⎢ 0 b T2 + hα , x2 γ αT 0 ⎥
⎢ T ⎥
⎢ 0 0 b 3 + hα , x3 γ α ⎥
T
B=⎢ T ⎥
(80)
b + hα , x2 γ αT b + hα , x1 γ αT
T
0
⎢ 2 1
⎥
⎢ b T3 + hα , x γ αT 0 b1T + hα , x1 γ αT ⎥
⎢ 3
⎥
⎢⎣ 0 b T3 + hα , x3 γ αT b T2 + hα , x2 γ αT ⎥⎦
Cette écriture de l’opérateur gradient discrétisé utilisant les formules de Hallquist [40]
est très pratique car les vecteurs γα , qui interviennent dans l’expression de B , vérifient les
conditions d’orthogonalité suivantes :
γ αT ⋅ x j = 0 , γ αT ⋅ h β = δ αβ (81)
79
Formulation des éléments SHB15 et SHB20 en linéaire
L’élément SHB20 est formulé dans les axes locaux du plan moyen. Les directions x, y, z
(ou encore x1 , x2 , x3 ) sont respectivement parallèles aux axes ξ ,η , ζ . La Figure 22 représente
la géométrie d’un élément de référence SHB20 et ses points d’intégration.
5
ς 20 8
5 20
17 19
4 19
16
η
6 13 3
18
10 18 15 7
2 17
9 14
1 12 16 4
14 8
1 13 15
7 12
9
6
11
11
2
3
ξ 10
80
Formulation des éléments SHB15 et SHB20 en linéaire
Les mêmes fonctions de forme sont utilisées pour définir le champ de déplacement de
l’élément ui en termes des déplacements nodaux uiI :
20
ui = uiI N I (ξ ,η , ζ ) = ∑ uiI N I (ξ ,η , ζ ) (83)
I =1
ui , j = uiI N I , j (84)
1 1
N1 = (1 − ξ )(1 − η )(1 − ζ )( −2 − ξ − η − ζ ) N6 = (1 + ξ )(1 − η )(1 + ζ )( −2 + ξ − η + ζ )
8 8
1 1
N2 = (1 + ξ )(1 − η )(1 − ζ )( −2 + ξ − η − ζ ) N 7 = (1 + ξ )(1 + η )(1 + ζ )( −2 + ξ + η + ζ )
8 8
1 1
N3 = (1 + ξ )(1 + η )(1 − ζ )( −2 + ξ + η − ζ ) N8 = (1 − ξ )(1 + η )(1 + ζ )( −2 − ξ + η + ζ )
8 8
1 1
N4 = (1 − ξ )(1 + η )(1 − ζ )( −2 − ξ + η − ζ ) N 9 = (1 − ξ 2 ) (1 − η )(1 − ζ )
8 4
1 1
N5 = (1 − ξ )(1 − η )(1 + ζ )( −2 − ξ − η + ζ ) N10 = (1 − η 2 ) (1 + ξ )(1 − ζ )
8 4
1 1
N11 = (1 − ξ 2 ) (1 + η )(1 − ζ ) N16 = (1 − ζ 2 ) (1 − ξ )(1 + η )
4 4
1 1
N12 = (1 − η 2 ) (1 − ξ )(1 − ζ ) N17 = (1 − ξ 2 ) (1 − η )(1 + ζ )
4 4
1 1
N13 = (1 − ζ 2 ) (1 − ξ )(1 − η ) N18 = (1 − η 2 ) (1 + ξ )(1 + ζ ) (86)
4 4
1 1
N14 = (1 − ζ 2 ) (1 + ξ )(1 − η ) N19 = (1 − ξ 2 ) (1 + η )(1 + ζ )
4 4
1 1
N15 = (1 − ζ 2 ) (1 + ξ )(1 + η ) N 20 = (1 − η 2 ) (1 − ξ )(1 + ζ )
4 4
ξ ∈ [ −1;1] ; η ∈ [ −1;1] ; ζ ∈ [ −1,1].
81
Formulation des éléments SHB15 et SHB20 en linéaire
Ces fonctions de forme transforment un cube régulier dans l’espace (ξ, η, ζ) en un cube
quelconque dans l’espace ( x1 , x2 , x3 ) ou bien ( x, y, z ) . En combinant les équations (82), (83)
et (86), on arrive à développer le champ de déplacement comme un terme constant, des termes
linéaires en xi , et des termes faisant intervenir les fonctions hα :
En évaluant l’équation (87) aux nœuds de l’élément, on arrive aux trois systèmes de
vingt équations suivants :
82
Formulation des éléments SHB15 et SHB20 en linéaire
⎧ST = (1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)
⎪ T
⎪h1 = (1 −1 −1 1 −1 1 1 −1 0 −1 0 1 0 0 0 0 0 1 0 −1)
⎪hT2 = (1 1 −1 −1 −1 −1 1 1 1 0 −1 0 0 0 0 0 −1 0 1 0 )
⎪ T
⎪h3 = (1 −1 1 −1 1 −1 1 −1 0 0 0 0 1 −1 1 −1 0 0 0 0 )
⎪hT = (1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 0 1)
⎪ T4
⎪h5 = (1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 0)
⎪h T = 1 1 1 1 1 1 1 1 1 1
⎪ 6 (
1 1 0 0 0 0 1 1 1 1)
⎪hT7 = ( −1 1 −1 1 1 −1 1 −1 0 0 0 0 0 0 0 0 0 0 0 0)
⎪ T
⎨h 8 = ( −1 − 1 1 1 − 1 − 1 1 1 0 0 0 0 −1 −1 1 1 0 0 0 0 )
⎪hT = ( −1 −1 −1 −1 1 1 1 1 0 −1 0 −1 0 0 0 0 0 1 0 1)
⎪ T9
⎪h10 = ( −1 1 1 −1 −1 1 1 −1 0 0 0 0 −1 1 1 −1 0 0 0 0 )
⎪ T
⎪h11 = ( −1 −1 −1 −1 1 1 1 1 −1 0 −1 0 0 0 0 0 1 0 1 0 )
⎪h12
T
= ( −1 1 1 − 1 − 1 1 1 − 1 0 1 0 −1 0 0 0 0 0 1 0 −1)
⎪ T
⎪h13 = ( −1 −1 1 1 −1 −1 1 1 −1 0 1 0 0 0 0 0 −1 0 1 0 )
⎪h14
T
= (1 1 −1 −1 −1 −1 1 1 0 0 0 0 0 0 0 0 0 0 0 0)
⎪ T
⎪h15 = (1 −1 −1 1 −1 1 1 −1 0 0 0 0 0 0 0 0 0 0 0 0) (90)
⎪hT = (1 −1 1 −1 1 −1 1 −1 0 0 0 0 0 0 0 0 0 0 0 0)
⎪⎩ 16
∂N
bTi = N ,i ( 0 ) = i = 1, 2, 3 (91)
∂xi ξ =η =ζ = 0
En fait, ce calcul est pratique car nous verrons par la suite qu’il est nécessaire de
déterminer la quantité suivante :
⎛ T 16 T ⎞
⎜ b j + ∑ hα , j γ α ⎟ = N , j (ξ ,η , ζ )
T
⎝ α =1 ⎠
b Tj = N T, j (0, 0, 0) = csteT
⎛ ∂ N1 ∂ N2 ∂ N3 ∂ N 20 ⎞
bTi = NT,i ( 0 ) = ⎜ . . . ⎟
⎝ ∂xi ∂xi ∂xi ∂xi ⎠ ξ =η =ζ =0
83
Formulation des éléments SHB15 et SHB20 en linéaire
∂ N I ⎛ ∂ N I ∂ξ ∂ N I ∂η ∂ N I ∂ζ ⎞ ⎛ ∂ NI ∂ NI ∂ NI ⎞
=⎜ + + ⎟⎟ = ⎜ j1 j + j2 j + j3 j ⎟
∂x j ⎜⎝ ∂ξ ∂x j ∂η ∂x j ∂ζ ∂x j ⎠ ⎝ ∂ξ ∂η ∂ζ ⎠
avec I = 1, 2,..., 20 et j = 1, 2,3
⎛ ∂ξ ∂ξ ∂ξ ⎞
⎜ ⎟
⎜ ∂x ∂y ∂z ⎟
⎛ j11 j12 j13 ⎞
⎜ ∂η ∂η ∂η ⎟ ⎜ ⎟
F −1 = ⎜ ⎟ = j21 j22 j23 ⎟
⎜ ∂x ∂y ∂z ⎟ ⎜⎜
j31 j32 j33 ⎠⎟
⎜ ∂ζ ∂ζ ∂ζ ⎟ ⎝
⎜ ⎟
⎝ ∂x ∂y ∂z ⎠
( )
N10,ξ = 0, 25 1 − η 2 (1 − ζ ) N10,η = −0,5η (1 + ξ )(1 − ζ )
N11,ξ = −0,5ξ (1 + η )(1 − ζ ) ( )
N11,η = 0, 25 1 − ξ 2 (1 − ζ )
( )
N12,ξ = −0, 25 1 − η 2 (1 − ζ ) N12,η = −0,5η (1 − ξ )(1 − ζ )
N13,ξ = −0, 25 (1 − η ) (1 − ζ )
2
(
N13,η = −0, 25 (1 − ξ ) 1 − ζ 2 )
N14,ξ = 0, 25 (1 − η ) (1 − ζ )
2
N14,η = −0, 25 (1 + ξ ) (1 − ζ )
2
N15,ξ = 0, 25 (1 + η ) (1 − ζ )
2
N15,η = 0, 25 (1 + ξ ) (1 − ζ )
2
N16,ξ = −0, 25 (1 + η ) (1 − ζ )
2
N16,η = 0, 25 (1 − ξ ) (1 − ζ )
2
( )
N18,ξ = 0, 25 1 − η 2 (1 + ζ ) N18,η = −0,5η (1 + ξ )(1 + ζ )
N19,ξ = −0,5ξ (1 + η )(1 + ζ ) ( )
N19,η = 0, 25 1 − ξ 2 (1 + ζ )
( )
N 20,ξ = −0, 25 1 − η 2 (1 + ζ ) N 20,η = −0,5η (1 − ξ )(1 + ζ )
84
Formulation des éléments SHB15 et SHB20 en linéaire
⎡1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1⎤
⎢8 − − − − 0 0 − − − 0 0 − ⎥
8 8 8 8 8 8 8 4 4 4 4 4 4 4 4
⎢ ⎥
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
bTi = ( j1i j2i j3i ) ⎢ − − − − − 0 0 − − − 0 0 ⎥
⎢8 8 8 8 8 8 8 8 4 4 4 4 4 4 4 4 ⎥
⎢ ⎥
⎢1 1 1 1 1 1 1 1 1 1 1 1
− − − − − − − − 0 0 0 0
1 1 1 1 ⎥
⎣⎢ 8 8 8 8 8 8 8 8 4 4 4 4 4 4 4 4 ⎦⎥
De plus, on peut vérifier par des considérations algébriques que les conditions
d’orthogonalité suivantes sont satisfaites :
⎡12 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0⎤
⎢0 12 0 0 0 0 0 0 0 0 0 0 0 8 0 0⎥
⎢0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 8⎥
⎢0 0 0 16 12 12 0 0 0 0 0 0 0 0 0 0⎥
⎢0 0 0 12 16 12 0 0 0 0 0 0 0 0 0 0⎥
⎢0 0 0 12 12 16 0 0 0 0 0 0 0 0 0 0⎥
⎢0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0⎥
⎢ 0⎥
hα ⋅ h β = ⎢ 0
T 0 0 0 0 0 0 12 0 0 0 0 8 0 0
0 0 0 0 0 0 0 0 12 0 8 0 0 0 0 0⎥
⎢0 0 0 0 0 0 0 0 0 12 0 8 0 0 0 0⎥
⎢0 0 0 0 0 0 0 0 8 0 12 0 0 0 0 0⎥
⎢0 0 0 0 0 0 0 0 0 8 0 12 0 0 0 0⎥
⎢0 0 0 0 0 0 0 8 0 0 0 0 12 0 0 0⎥
⎢0 8 0 0 0 0 0 0 0 0 0 0 0 8 0 0 ⎥⎥
⎢
⎢8 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0⎥
⎣⎢ 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 8 ⎦⎥
α , β = 1, 2,...,16
85
Formulation des éléments SHB15 et SHB20 en linéaire
b Ti ⋅ hα = 0 ; b Ti ⋅ S = 0 ; b Ti ⋅ x j = δ ij ; h1T ⋅ S = 0 ; h T2 ⋅ S = 0 ; h T3 ⋅ S = 0 ; h T4 ⋅ S = 16 ;
h T5 ⋅ S = 16 ; h T6 ⋅ S = 16 ; h T7 ⋅ S = 0 ; h T8 ⋅ S = 0 ; h T9 ⋅ S = 0 ; h10
T
⋅S = 0 ;
h11 ⋅ S = 0 ; h12 ⋅ S = 0 ; h13 ⋅ S = 0 ; h14 ⋅ S = 0 ; h15 ⋅ S = 0 ; h16 ⋅ S = 0 ;
T T T T T T (92)
α = 1, 2,...,16 ; i , j = 1, 2, 3
δ désigne le symbol de Kronecker
À ce stade, on peut déterminer les constantes inconnues qui interviennent dans l’écriture
(87) du champ de déplacement en multipliant scalairement l’équation (88) par b Tj , S T et hαT ,
respectivement, et en utilisant les relations d’orthogonalité (92). On obtient ainsi :
⎧
⎪S T ⋅ d i = a0 i S T ⋅ S + a1i S T ⋅ x1 + a2 i S T ⋅ x 2 + a3i S T ⋅ x 3 + c1i S T ⋅ h1 + c2 i S T ⋅ h 2 + ... + c16 i S T ⋅ h16
⎪⎪
( )( )
⎨ = 20 a0 i + S ⋅ x j b j ⋅ d i + 0c1i + 0c2 i + 0c3i + 16c4 i + 16c5 i + 16c6 i + 0c7 i + 0c8 i + ... + 0c16 i
T T
⎪
⎪⎩ 0 i 20 ( ( j ) j i )
⎪ a = 1 S T − S T ⋅ x ⋅ b T .d − 1 (16c + 16c + 16c )
20
4i 5i 6i
⎧⎪h1T ⋅ di = a0i h1T ⋅ S + a1i h1T ⋅ x1 + a2i h1T ⋅ x 2 + a3i h1T ⋅ x3 + c1i h1T ⋅ h1 + c2i h1T ⋅ h 2 + ... + c16i h1T ⋅ h16
⎨
( )( )
⎪⎩= 0 + h1 ⋅ x j b j ⋅ di + 12c1i + 0c2i + 0c3i + ... + 8c15i + 0c16i
T T
⎧
⎪ T
⎪h 4 ⋅ d i = a0 i h 4 ⋅ S + a1i h 4 ⋅ x1 + a2 i h 4 ⋅ x 2 + a3i h 4 ⋅ x 3 + c1i h 4 ⋅ h1 + c2 i h 4 ⋅ h 2 + ... + c16 i h 4 ⋅ h16
T T T T T T T
⎪
( )( )
⎨ = 16 a0 i + h 4 ⋅ x j b j ⋅ d i + 0 c1i + ... + 16c4 i + 12c5 i + 12c6 i + ... + 0 c16 i
T T
⎪
⎪⎛ h T − 4 S T ⎞ ⋅ d = ⎡ ⎛ h T − 4 S T ⎞ ⋅ x ⎤ b T ⋅ d + 16 c − 4 c − 4 c
⎪⎜⎝ 4 5 ⎟⎠ i ⎢⎣ ⎜⎝ 4 5 ⎟⎠ j ⎥⎦ j i ( )
5
4i
5
5i
5
6i
⎩
⎧
⎪ T
⎪h 5 ⋅ d i = a0 i h 5 ⋅ S + a1i h 5 ⋅ x1 + a2 i h 5 ⋅ x 2 + a3i h 5 ⋅ x 3 + c1i h 5 ⋅ h1 + c2 i h 5 ⋅ h 2 + ... + c16 i h 5 ⋅ h16
T T T T T T T
⎪
( )( )
⎨ = 16 a0 i + h 5 ⋅ x j b j ⋅ d i + 0 c1i + 0 c2 i + 0c3i + 12c4 i + 16c5 i + 12c6 i + 0 c7 i + ... + 0c16 i
T T
⎪
⎪⎛ h T − 4 S T ⎞ ⋅ d = ⎡ ⎛ h T − 4 S T ⎞ ⋅ x ⎤ b T ⋅ d − 4 c + 16 c − 4 c
⎪⎜⎝ 5 5 ⎟⎠ i ⎢⎣ ⎜⎝ 5 5 ⎟⎠ j ⎥⎦ j i ( 5
) 4i
5
5i
5
6i
⎩
86
Formulation des éléments SHB15 et SHB20 en linéaire
⎧
⎪ T
⎪h 6 ⋅ d i = a0 i h 6 ⋅ S + a1i h 6 ⋅ x1 + a2 i h 6 ⋅ x 2 + a3i h 6 ⋅ x 3 + c1i h 6 ⋅ h1 + c2 i h 6 ⋅ h 2 + ... + c16 i h 6 ⋅ h16
T T T T T T T
⎪
( )( )
⎨ = 16 a0 i + h 6 ⋅ x j b j ⋅ d i + 0 c1i + 0 c2 i + 0c3i + 12c4 i + 12c5 i + 16c6 i + ... + 0c16 i
T T
⎪
⎪⎛ h T − 4 S T ⎞ ⋅ d = ⎡ ⎛ h T − 4 S T ⎞ ⋅ x ⎤ b T ⋅ d − 4 c − 4 c + 16 c
⎪⎜⎝ 6 5 ⎟⎠ i ⎢⎣ ⎜⎝ 6 5 ⎟⎠ j ⎥⎦ j i ( 5
) 4i
5
5i
5
6i
⎩
⎧⎪h T7 ⋅ d i = a0 i h T7 ⋅ S + a1i h T7 ⋅ x1 + a2 i h T7 ⋅ x 2 + a3i h T7 ⋅ x 3 + c1i h T7 ⋅ h1 + c2 i h T7 ⋅ h 2 + ... + c16 i h T7 ⋅ h16
⎨
( )( T
)
⎪⎩ = 0 + h 7 ⋅ x j b j ⋅ d i + 0c1i + 0 c2 i + ... + 8c7 i + 0c8 i + ... + 0 c16 i
T
⎧⎪h10 T
⋅ d i = a0 i h10
T
⋅ S + a1i h10
T
⋅ x1 + a2 i h10T
⋅ x 2 + a3i h10
T
⋅ x 3 + c1i h10T
⋅ h1 + c2 i h10
T
⋅ h 2 ... + c16 i h10
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h10 ⋅ x j b j ⋅ d i + 0c1i + 0c2 i + ... + 12c10 i + 0c11i + 8c12 i + 0c13i + ... + 0c16 i
T T
⎧⎪h11 T
⋅ d i = a0 i h11
T
⋅ S + a1i h11
T
⋅ x1 + a2 i h11
T
⋅ x 2 + a3i h11
T
⋅ x 3 + c1i h11T
⋅ h1 + c2 i h11
T
⋅ h 2 ... + c16 i h11
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h11 ⋅ x j b j ⋅ d i + 0c1i + 0c2 i + ... + 8c9 i + 0c10 i + 12c11i + 0 c12 i + ... + 0c16 i
T T
⎧⎪h12 T
⋅ d i = a0 i h12
T
⋅ S + a1i h12
T
⋅ x1 + a2 i h12T
⋅ x 2 + a3i h12
T
⋅ x 3 + c1i h12T
⋅ h1 + c2 i h12
T
⋅ h 2 ... + c16 i h12
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h12 ⋅ x j b j ⋅ d i + 0c1i + 0c2 i + ... + 8c10 i + 0c11i + 12c12 i + 0c13i + ... + 0c16 i
T T
⎧⎪h13 T
⋅ d i = a0 i h13
T
⋅ S + a1i h13
T
⋅ x1 + a2 i h13T
⋅ x 2 + a3i h13
T
⋅ x 3 + c1i h13T
⋅ h1 + c2 i h13
T
⋅ h 2 ... + c16 i h13
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h13 ⋅ x j b j ⋅ d i + 0c1i + 0 c2 i + ... + 8c8 i + 0c9 i + .. + 12c13i + ... + 0c16 i
T T
⎧⎪h14 T
⋅ d i = a0 i h14
T
⋅ S + a1i h14
T
⋅ x1 + a2 i h14T
⋅ x 2 + a3i h14
T
⋅ x 3 + c1i h14T
⋅ h1 + c2 i h14
T
⋅ h 2 ... + c16 i h14
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h14 ⋅ x j b j ⋅ d i + 0c1i + 8c2 i + 0 c3i + ... + 8c14 i + 0c15 i + 0c16 i
T T
⎧⎪h15 T
⋅ d i = a0 i h15
T
⋅ S + a1i h15
T
⋅ x1 + a2 i h15T
⋅ x 2 + a3i h15
T
⋅ x 3 + c1i h15T
⋅ h1 + c2 i h15
T
⋅ h 2 ... + c16 i h15
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + d15 ⋅ x j b j ⋅ d i + 8c1i + 0 c2 i + ... + 8c15 i + 0c16 i
T T
⎧⎪h16 T
⋅ d i = a0 i h16
T
⋅ S + a1i h16
T
⋅ x1 + a2 i h16T
⋅ x 2 + a3i h16
T
⋅ x 3 + c1i h16T
⋅ h1 + c2 i h16
T
⋅ h 2 ... + c16 i h16
T
⋅ h16
⎨
( )( )
⎪⎩ = 0 + h16 ⋅ x j b j ⋅ d i + 0c1i + 0c2 i + 8c3i + 0c4 i + ... + 0 c16 i
T T
87
Formulation des éléments SHB15 et SHB20 en linéaire
⎡
⎢
( h1T − (h1T ⋅ x j )bTj ) ⋅ d i ⎤
⎥
⎢
⎢
( hT2 − (hT2 ⋅ x j )bTj ) ⋅ d i ⎥
⎥
⎡12
⎢0
0 0 0 0 0 0 0 0 0 0 0 0 0 8 0⎤ ⎢ ( hT3 − (hT3 ⋅ x j )bTj ) ⋅ d i ⎥
⎢ 12 0 0 0 0 0 0 0 0 0 0 0 8 0 0 ⎥⎥ ⎡ c ⎤ ⎢ ⎥
1i ⎢ ⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T ⎞ ⎥
⎢0 0 12 0 0 0 0 0 0 0 0 0 0 0 0 8 ⎥ ⎢ c ⎥ ⎢ ⎜ ⎜ h 4 − S ⎟ − ⎢⎜ h 4 − S ⎟ ⋅ x j ⎥ b j ⎟ ⋅ d i ⎥
⎢ ⎥ ⎢ 2 i ⎥ ⎢⎝ ⎝ 5 ⎠ ⎣⎝ 5 ⎠ ⎦ ⎠ ⎥
⎢0 16 4 4 ⎥ ⎢ c ⎥
0 0 − − 0 0 0 0 0 0 0 0 0 0 3i ⎢ ⎥
⎢ 5 5 5 ⎥ ⎢ c ⎥ ⎢ ⎛ ⎛ h T − 4 S T ⎞ − ⎡⎛ hT − 4 ST ⎞ ⋅ x ⎤ b T ⎞ ⋅ d ⎥
⎥ ⎢ 4i ⎥ ⎢ ⎝⎜ ⎜ ⎟ ⎜ ⎟ ⎟
⎢ 5 ⎠ ⎢⎣⎝ 5 ⎠ ⎥⎦ ⎠ ⎥
5 5 j j i
4 16 4 ⎝
⎢0 0 0 − − 0 0 0 0 0 0 0 0 0 0 ⎥ ⎢ c5i ⎥ ⎢ ⎥
⎢ 5 5 5 ⎥ ⎢ ⎥ ⎢ ⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T ⎞ ⎥
⎢ 4 4 16 ⎥ ⎢ c6i ⎥ ⎢⎜ ⎜ h 6 − S ⎟ − ⎢⎜ h 6 − S ⎟ ⋅ x j ⎥ b j ⎟ ⋅ di ⎥
⎢0 0 0 − − 0 0 0 0 0 0 0 0 0 0 ⎥ ⎢ c ⎥ ⎢⎝ ⎝ 5 ⎠ ⎣⎝ 5 ⎠ ⎦ ⎠ ⎥
⎢ 5 5 5 ⎥ ⎢ 7i
⎥
⎢0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0 ⎥ ⎢ c8i ⎥ = ⎢
⎢
( hT7 − (hT7 ⋅ x j )bTj ) ⋅ di ⎥
⎥
⎢ ⎥⎢ ⎥
⎢0 0 0 0 0 0 0 12 0 0 0 0 8 0 0 0 ⎥ ⎢ c9i ⎥ ⎢
⎢ ( 8 8 j j) i
h T
− (h T
⋅ x )b T
⋅ d ⎥
⎥
⎢0 0 ⎥ ⎢ c10i ⎥ ⎢
⎢
0 0 0 0 0 0 0 12 0 8 0 0 0 0
⎥ ⎢c ⎥ ⎢ ( 9 9 j j) i
h T
− (h T
⋅ x )b T
⋅ d ⎥
⎥
⎢0 0 0 0 0 0 0 0 0 12 0 8 0 0 0 0 ⎥ ⎢ 11i ⎥ ⎢
⎢0 0 0 0 0 0 0 0 8 0 12 0 0 0 0 0 ⎥ ⎢ c12i ⎥ ⎢
( 10 10 j j ) i
h T
− (h T
⋅ x )b T
⋅ d ⎥
⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢0 0 0 0 0 0 0 0 0 8 0 12 0 0 0 c
0 ⎥ ⎢ 13i ⎥ ⎢ ( 11 11 j j ) i
h T
− (h T
⋅ x )b T
⋅ d ⎥
⎢0 ⎢ ⎥
0 0 0 0 0 0 8 0 0 0 0 12 0 0 0 ⎥ ⎢ c14i ⎥ ⎢ ( 12 12 j j ) i
h T
− (h T
⋅ x )b T
⋅ d ⎥
⎢ ⎥ ⎢ ⎥
⎢0 8 0 0 0 0 0 0 0 0 0 0 0 8 0 0 ⎥ ⎢ c15i ⎥ ⎢ ⎥
⎢ ⎥ ⎢⎢ c16i ⎥⎥ ⎢
0⎥ ⎣ ⎦
( h T
13 − (h T
13 ⋅ x j )b T
j ) ⋅ d i ⎥
⎢8 0 0 0 0 0 0 0 0 0 0 0 0 0 8
⎢ ⎥
⎢⎣ 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 8 ⎥⎦ ⎢ ( h14 − (h14 ⋅ x j )b j ) ⋅ di
T T T
⎥
⎢ ⎥
⎢ ( h15 − (h15 ⋅ x j )b j ) ⋅ di
T T T
⎥
⎢ ⎥
⎢⎣ ( h16 − (h16 ⋅ x j )b j ) ⋅ di
T T T
⎥⎦
ou encore :
88
Formulation des éléments SHB15 et SHB20 en linéaire
⎡ 1 1 ⎤
⎢ 4 0 0 0000 0 0 0 0 0 0 0 −4 0 ⎥
⎢ 1 1 ⎥
⎢ 0 0 0000 0 0 0 0 0 0 − 0 0 ⎥
⎢ 4 4
1 1⎥ ( h1T − (h1T ⋅ x j )bTj ) ⋅ di
⎢ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − ⎥ ⎡⎢ ⎤
⎥
⎢
⎢ 0 0 0
4
3 1 1
4⎥
⎢ ( h 2 − (h 2 ⋅ x j )b j ) ⋅ di
T T T
⎥
0 0 0 0 0 0 0 0 0 0 ⎥⎢ ( h T
− (h T
⋅ x )b T
) ⋅ d ⎥
⎢ 888 ⎥ 3 3 j j i
⎢ 131 ⎥ ⎢ ⎥
0 0 0 0 0 0 0 0 0 0 ⎥ ⎢⎛⎜ ⎛⎜ hT4 − ST ⎞⎟ − ⎡⎛⎜ hT4 − ST ⎞⎟ ⋅ x j ⎤ bTj ⎞⎟ ⋅ d i ⎥
4 4
⎡ c1i ⎤ ⎢ 0 0 0
⎢ c2i ⎥ ⎢ 888 ⎥ ⎢⎝ ⎝ 5 ⎠ ⎢⎣⎝ 5 ⎠ ⎥⎦
⎠ ⎥
⎢ c3i ⎥ ⎢ 0 0 0 1 1 3 0 0 0 0 0 0 0 0 0 0 ⎥ ⎢⎛ ⎛ T 4 T ⎞ ⎡⎛ T 4 T ⎞ ⎤ T ⎞ ⎥
h − S − h − S ⋅x b ⋅d
⎢ c4i ⎥ ⎢ 888 ⎥ ⎢⎜⎝ ⎜⎝ 5 5 ⎟⎠ ⎢⎣⎜⎝ 5 5 ⎟⎠ j ⎥⎦ j ⎟⎠ i ⎥
⎢ c5i ⎥ ⎢ 1 ⎥ ⎢ ⎥
⎢ c ⎥ ⎢ 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0 ⎥ ⎢ ⎛ ⎛ h T − 4 S T ⎞ − ⎡⎛ h T − 4 S T ⎞ ⋅ x ⎤ b T ⎞ ⋅ d ⎥
⎜ ⎜ 6 5 ⎟ ⎢⎜ 6 5 ⎟ j ⎥ j ⎟ i ⎥
⎢ c6i ⎥ ⎢ 3 1 ⎥ ⎢⎝ ⎝ ⎠ ⎣⎝ ⎠ ⎦ ⎠
0 0 0 0 − 0 0 0 ⎥⎢ ⎥
⎢c ⎥ ⎢ 0 0 0 0 0 0 0
( )
7i
20 10 ⎢ h T
− (h T
⋅ x )b j ⋅ di
T
⎥
⎢ 8i ⎥ = ⎢ ⎥
7 7 j
En posant :
( ) ( ) (
γαT = nα 1 h1T − ( h1T ⋅ x j ) bTj + nα 2 hT2 − ( hT2 ⋅ x j ) bTj + nα 3 hT3 − ( hT3 ⋅ x j ) bTj )
⎡⎛ 4 ⎞ ⎛⎛ 4 ⎞ ⎞ ⎤ ⎡⎛ 4 ⎞ ⎛⎛ 4 ⎞ ⎞ ⎤
+ nα 4 ⎢⎜ hT4 − ST ⎟ − ⎜ ⎜ hT4 − ST ⎟ ⋅ x j ⎟ bTj ⎥ + nα 5 ⎢⎜ hT5 − ST ⎟ − ⎜ ⎜ hT5 − ST ⎟ ⋅ x j ⎟ bTj ⎥
⎣⎝ 5 ⎠ ⎝⎝ 5 ⎠ ⎠ ⎦ ⎣⎝ 5 ⎠ ⎝⎝ 5 ⎠ ⎠ ⎦
⎡⎛ 4 ⎞ ⎛⎛ ⎞ ⎤
⎣⎝ 5 ⎠ ⎝⎝
4 ⎞
5 ⎠ ⎠ ⎦
(
+ nα 6 ⎢⎜ hT6 − ST ⎟ − ⎜ ⎜ hT6 − ST ⎟ ⋅ x j ⎟ bTj ⎥ + nα 7 hT7 − ( hT7 ⋅ x j ) bTj )
( ) ( )
+ nα 8 hT8 − ( hT8 ⋅ x j ) bTj + nα 9 hT9 − ( hT9 ⋅ x j ) bTj + nα 10 h10
T
− ( h10
T
(
⋅ x j ) bTj )
+ nα 11 (h − (h
T
11
T
11 )
⋅ x j ) bTj + nα 12 h12
T
(
− ( h12
T
)
⋅ x j ) bTj + nα 13 h13
T
− ( h13
T
(
⋅ x j ) bTj )
+ nα 14 h14
T
(− ( h14
T
)
⋅ x j ) bTj + nα 15 h15
T
(
− ( h15
T
)
⋅ x j ) bTj + nα 16 h16
T
− ( h16
T
(
⋅ x j ) bTj )
89
Formulation des éléments SHB15 et SHB20 en linéaire
⎡ 1 1 ⎤
⎢ 4 0 0 0000 0 0 0 0 0 0 0 − 0 ⎥
4
⎢ 1 1 ⎥
⎢ 0 0 0000 0 0 0 0 0 0 − 0 0 ⎥
⎢ 4 4
1 1⎥
⎢ 0 0 0000 0 0 0 0 0 0 0 0 − ⎥
⎢ 4 4⎥
⎢ 0 0 0 3 1 1
0 0 0 0 0 0 0 0 0 0 ⎥
⎢ 8 8 8 ⎥
⎢ 1 3 1 ⎥
⎢ 0 0 0 8 8 8
0 0 0 0 0 0 0 0 0 0 ⎥
⎢ 1 1 3 ⎥
⎢ 0 0 0 0 0 0 0 0 0 0 0 0 0 ⎥
⎢ 8 8 8 ⎥
⎢ 0 0 0 1
0 00 0 0 0 0 0 0 0 0 0 ⎥
⎢ 8 ⎥
⎢ 3 1 ⎥
⎢ 0 0 0 0000 0 0 0 0 − 0 0 0 ⎥
⎡⎣nαβ ⎤⎦ = ⎢ 20 10
3 1 ⎥
⎢ 0 0 0 0000 0 0 −
0 0 0 0 0 ⎥
⎢ 20 10 ⎥
3 1
⎢ 0 0 0 000 0 0 0 0 − 0 0 0 0 ⎥
⎢ 20 10 ⎥
⎢ 0 0 0 1 3
000 0 0 − 0 0 0 0 0 0 ⎥
⎢ 10 20 ⎥
⎢ 1 3 ⎥
⎢ 0 0 0 000 0 0 0 − 0 0 0 0 0 ⎥
⎢ 10 20 ⎥
1 3
⎢ 0 0 0 000 0− 0 0 0 0 0 0 0 ⎥
⎢ 10 20 ⎥
⎢ 0 −1 0 000 0 0 0 0 0 0 0
1
0 0 ⎥
⎢ 4 8 ⎥
⎢ 1 1 ⎥
−
⎢ 4 0 0 0000 0 0 0 0 0 0 0 0 ⎥
8
⎢ 1 1 ⎥
⎢ 0 0 − 0000 0 0 0 0 0 0 0 0 ⎥
⎣ 4 8 ⎦
α , β = 1, 2,...,16
90
Formulation des éléments SHB15 et SHB20 en linéaire
⎛ 16
⎞
ui , j = ⎜ bTj + ∑ hα , j γαT ⎟ ⋅ di = ( bTj + hα , j γαT ) ⋅ di (95)
⎝ α =1 ⎠
∇s (u ) = B ⋅ d (96)
où:
⎡ ux,x ⎤
⎢ u ⎥
⎢ y,y ⎥ ⎡ d1 ⎤
⎢
∇ s (u ) = ⎢
u z,z
⎥ ,
⎥
d = ⎢⎢d 2 ⎥⎥ (97)
⎢u x , y + u y , x ⎥ ⎢⎣ d3 ⎥⎦
⎢ ux,z + uz,x ⎥
⎢ ⎥
⎣⎢ u y , z + u z , y ⎦⎥
⎡ b1T + hα , x1 γ αT 0 0 ⎤
⎢ ⎥
⎢ 0 b T2 + hα , x2 γ αT 0 ⎥
⎢ T ⎥
⎢ 0 0 b 3 + hα , x3 γ α ⎥
T
B=⎢ T ⎥
(98)
b + hα , x2 γ αT b + hα , x1 γ αT
T
0
⎢ 2 1
⎥
⎢ b T3 + hα , x γ αT 0 b1T + hα , x1 γ αT ⎥
⎢ 3
⎥
⎢⎣ 0 b T3 + hα , x3 γ αT b T2 + hα , x2 γ αT ⎥⎦
Cette écriture de l’opérateur gradient discrétisé utilisant les formules de Hallquist [40]
est très commode car les vecteurs γα , qui interviennent dans l’expression de B , vérifient les
conditions d’orthogonalité suivantes:
γ αT ⋅ x j = 0 , γ αT ⋅ h β = δ αβ (99)
91
Formulation des éléments SHB15 et SHB20 en linéaire
δπ ( u , ε , σ ) = ∫ δ ε T ⋅ σ d V + δ ∫ σ T ⋅ ( ∇ s ( u ) − ε ) d V − δ d T ⋅ f ext = 0 (100)
Ve Ve
δπ ( u , ε ) = ∫ δ ε T ⋅ σ d V − δ d T ⋅ f ext = 0 (101)
Ve
∇ s ( u ( x, t ) ) = B ( x ) ⋅ d ( t ) (103)
ε ( x , t ) = B ( x ) ⋅ d ( t ) (104)
92
Formulation des éléments SHB15 et SHB20 en linéaire
⎛ T ⎞
δ d T ⋅ ⎜ ∫ B ⋅ σ dV − f ext ⎟ = 0 (105)
⎜V ⎟
⎝ e ⎠
T
f int = ∫ B ⋅ σ dV (107)
Ve
Dans l'équation ci-dessus, il est bien précisé que la contrainte σ est calculée par la loi
constitutive à partir du taux de déformation postulée ε . Pour les problèmes non linéaires, σ
peut aussi être une fonction intégrale du taux de déformation postulée et des autres variables
internes :
σ = F ( ε , α,...) (108)
où α représente les variables internes. La formulation ainsi obtenue est valable pour des
problèmes incluant les deux types de non linéarités : géométriques et matériau. Dans le cas de
problèmes linéaires, on a :
La matrice de comportement élastique C , dans le cas d’un matériau isotrope, est choisie
comme suit :
⎡ E Eν ⎤
⎢1 −ν 2 1 −ν 2 0 0 0 0 ⎥
⎢ ⎥
⎢ Eν E
0 0 0 0 ⎥
⎢1 −ν 2 1 −ν 2 ⎥
⎢ ⎥
⎢ 0 0 E 0 0 0 ⎥
C=⎢ 0 0 0
E
0 0 ⎥
⎥
⎢
⎢ 2 (1 + ν ) ⎥
⎢ E ⎥
⎢ 0 0 0 0 0 ⎥
⎢ 2 (1 + ν ) ⎥
⎢ E ⎥
⎢ 0 0 0 0 0 ⎥
⎢⎣ 2 (1 + ν ) ⎥⎦
Dans cette matrice, E est le module d’Young et ν est le coefficient de Poisson. Cette loi
est spécifique aux éléments SHB. Elle ressemble à celle que l’on aurait dans le cas de
l’hypothèse des contraintes planes, mise à part le terme (3,3). Nous pouvons noter que ce
choix entraine un comportement anisotrope artificiel. Ce choix permet de satisfaire tous les
tests sans introduire de blocage.
93
Formulation des éléments SHB15 et SHB20 en linéaire
Ve
(111)
Dans une approche standard en déplacement, qui sera retenue dans la suite, le taux de
déformation postulée s'identifie à la partie symétrique du gradient de vitesse, ce qui revient à
remplacer B par B dans les expressions précédentes. On obtient donc simplement :
Ke = ∫ B ⋅C ⋅ B dV
T
Ve
(112)
Point de w ( ξ ,η ,ζ
Gauss
ξ η ζ )
P(1) 1/2 1/2 ζ G1 = − 0.906179845938664 0.23692688 5056189/6
P(2) 1/2 1/2 ζ G 2 = −0.538469310105683 0.478628670 499366/6
P(3) 1/2 1/2 ζ G3 = 0 0.568888888888889/6
P(4) 1/2 1/2 ζ G 4 = 0.538469310105683 0.478628670499366/6
P(5) 1/2 1/2 ζ G 5 = 0.906179845938664 0.23692688 5056189/6
P(6) 0 1/2 ζ G 6 = −0.906179845938664 0.23692688 5056189/6
P(7) 0 1/2 ζ G 7 = −0.538469310105683 0.478628670499366/6
P(8) 0 1/2 ζ G8 = 0 0.56888888 8888889/6
P(9) 0 1/2 ζ G 9 = 0.538469310105683 0.478628670499366/6
P(10) 0 1/2 ζ G10 = 0.906179845938664 0.23692688 5056189/6
P(11) 1/2 0 ζ G11 = − 0.906179845938664 0.23692688 5056189/6
P(12) 1/2 0 ζ G12 = − 0.538469310105683 0.478628670499366/6
P(13) 1/2 0 ζ G13 = 0 0.568888888888889/6
P(14) 1/2 0 ζ G14 = 0.538469310105683 0.478628670499366/6
P(15) 1/2 0 ζ G15 = 0.906179845938664 0.23692688 5056189/6
15
K e = ∑ w (ξGj ) J (ξGj )BT (ξGj ) ⋅ C ⋅ B (ξGj )
j =1
94
Formulation des éléments SHB15 et SHB20 en linéaire
Les coordonnées des points de Gauss et leurs poids pour l’élément SHB20 sont donnés
dans le tableau ci-dessous :
Point de w (ξ ,η ,ζ
Gauss
ξ η ζ )
P(1) −1 3 −1 3 ζ G1 = − 0.906179845938664 0.236926885056189
P(2) −1 3 −1 3 ζ G 2 = −0.538469310105683 0.478628670499366
P(3) −1 3 −1 3 ζ G3 = 0 0.568888888888889
P(4) −1 3 −1 3 ζ G 4 = 0.538469310105683 0.478628670499366
P(5) −1 3 −1 3 ζ G 5 = 0.906179845938664 0.236926885056189
P(6) 1 3 −1 3 ζ G 6 = −0.906179845938664 0.236926885056189
P(7) 1 3 −1 3 ζ G 7 = −0.538469310105683 0.478628670499366
P(8) 1 3 −1 3 ζ G8 = 0 0.568888888888889
P(9) 1 3 −1 3 ζ G 9 = 0.538469310105683 0.478628670499366
P(10) 1 3 −1 3 ζ G10 = 0.906179845938664 0.236926885056189
P(11) 1 3 1 3 ζ G11 = − 0.906179845938664 0.236926885056189
P(12) 1 3 1 3 ζ G12 = − 0.538469310105683 0.478628670499366
P(13) 1 3 1 3 ζ G13 = 0 0.568888888888889
P(14) 1 3 1 3 ζ G14 = 0.538469310105683 0.478628670499366
P(15) 1 3 1 3 ζ G15 = 0.906179845938664 0.236926885056189
P(16) −1 3 1 3 ζ G16 = − 0.906179845938664 0.236926885056189
P(17) −1 3 1 3 ζ G17 = − 0.538469310105683 0.478628670499366
P(18) −1 3 1 3 ζ G18 = 0 0.568888888888889
P(19) −1 3 1 3 ζ G19 = 0.538469310105683 0.478628670499366
P(20) −1 3 1 3 ζ G 20 = 0.906179845938664 0.236926885056189
20
K e = ∑ w ( PGj ) J ( PGj )BT ( PGj ) ⋅ C ⋅ B ( PGj )
j =1
( K + µK σ ) ⋅ u = 0 ⇒ K ⋅ u = λK σ ⋅ u
95
Formulation des éléments SHB15 et SHB20 en linéaire
3
eijQ ( δu , ∆u ) = ∑ δ uk ,i ∆uk , j = δ uk ,i ∆uk , j
k =1
δ uT ⋅ K σ ⋅ ∆u = ∫ σ : eQ (δ u, ∆u ) d Ω = ∫ σ : ∇δ uT ⋅∇∆u d Ω
Ω0 Ω0
Afin d’exprimer cette matrice dans l’espace discrétisé, introduisons les opérateurs
gradient quadratique discrétisés BQ sous forme matricielle tels que :
⎡ e11Q ⎤ ⎡ δ uT ⋅ B11 Q
⋅ ∆u ⎤
⎢ Q ⎥ ⎢ T Q ⎥
⎢ e22 ⎥ ⎢δ u ⋅ B 22 ⋅ ∆u ⎥
⎢ eQ ⎥ ⎢δ uT ⋅ BQ ⋅ ∆u ⎥
eQ ( δu , ∆u ) = ⎢ Q 33 Q ⎥ = ⎢ T 33 ⎥
⎢ e12 + e21 ⎥ ⎢δ u ⋅ B12 ⋅ ∆u ⎥
Q
⎢ eQ + eQ ⎥ ⎢δ uT ⋅ BQ ⋅ ∆u ⎥
⎢ 13 31
⎥ ⎢ T 13 ⎥
⎢⎣e23 + e32 ⎥⎦ ⎢⎣δ u ⋅ B 23 ⋅ ∆u ⎥⎦
Q Q Q
Les différents termes BQij sont donnés par les équations suivantes :
⎡
B1B1T 0
⎢ 0 ⎤⎥ ⎡
⎢B2B2
T
0 0 ⎤⎥ ⎡ T
⎢ B 3B 3 0 0 ⎤⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
Q
B11 = 0 B1B1T 0 ⎥⎥ ; BQ22 = ⎢⎢ 0 B 2B T2 0 ⎥⎥ ; B33
⎢
⎢
Q
= ⎢⎢ 0 B3B3T 0 ⎥⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢0
⎢⎣
0 B1B1T ⎥⎥⎦ ⎢
⎢⎣
0 0 B 2 B T2 ⎥⎥⎦ ⎢
⎢⎣
0 0 B3B3T ⎥⎥⎦
⎡ ⎤
B1B T2 + B 2B1T
⎢ 0 0 ⎥
⎢ ⎥
Q
B12 = ⎢
⎢ 0 B1B 2 + B 2B1T
T
0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B1B T2 + B 2 B1T ⎥⎦⎥
⎡ ⎤
B1B3T + B3B1T
⎢ 0 0 ⎥
⎢ ⎥
Q
B13 = ⎢
⎢ 0 B1B3 + B3B1T
T
0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B1B3 + B3B1 ⎦⎥
T T⎥
⎡ ⎤
B 2B3T + B3B T2
⎢ 0 0 ⎥
⎢ ⎥
BQ23 = ⎢
⎢ 0 B 2B3 + B3B T2
T
0 ⎥
⎥
⎢ ⎥
⎢
⎣⎢
0 0 B 2B3 + B3B 2 ⎦⎥
T T⎥
96
Formulation des éléments SHB15 et SHB20 en linéaire
(
Bi = bi + hα ,i γα )
Avec ces notations, la contribution à la matrice de rigidité géométrique, k σ , au point de
Gauss ξ j est donnée par :
15
K σ = ∑ ω ( PGj ) J ( PGj ) k σ ( PGj ) pour l'élément SHB15
j =1
20
K σ = ∑ ω ( PGj ) J ( PGj ) k σ ( PGj ) pour l'élément SHB20
j =1
Les forces de pression suiveuses sont présentes dans la matrice tangente via la matrice
K P , car les forces externes suiveuses dépendent du déplacement. Les forces de pression
suiveuses s’écrivent :
∫ p nT ⋅ u dS = ∫
T
p det[ F (u)] nT0 ⋅ F (u) −1 dS0 = p F0 − p K P ⋅ u
∂Ω ∂Ω0
F (u) = I + ∇ ( u )
en utilisant les notations :
T
n 0 = (n1 , n2 , n3 ) , normale à la surface extérieure de l’élément dans la
configuration de référence ;
b i , vecteur de dimension 6 (pour le SHB15) ou 8 (pour le SHB20), dérivée des
fonctions de forme aux 6 (pour le SHB15) ou 8 (pour le SHB20) nœuds de la
face de l’élément chargée en pression ;
S0 aire de la face chargée en pression.
La formulation précédente conduit à une matrice non-symétrique. On sait que l’on peut
néanmoins utiliser une formulation symétrique si les forces extérieures dues à la pression
dérivent d’un potentiel. C’est le cas si les forces de pression ne travaillent pas sur la frontière
du domaine modélisé. On considère donc que la partie symétrique de la matrice suffit. La
matrice symétrisée prend la forme suivante :
97
Formulation des éléments SHB15 et SHB20 en linéaire
⎡ b 2 n1 − b 1 n 2 b 3 n1 − b 1 n 3 ⎤
T T T T
0
⎢ ⎥
b n − b n b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1 1 2
⎢ ⎥
b n1 − b n 2 b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1
⎢ ⎥
b n − b n b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1 1 2
⎢ 0 b n1 − b n 2
T T
b 3 n1 − b 1 n 3 ⎥
T T
⎢ 2 1
⎥
⎢ 0 b n − b n
T T
b 3 n1 − b 1 n 3 ⎥
T T
⎢ T 2 1 1 2
⎥
⎢b 1 n 2 − b T2 n1 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢ T ⎥
⎢b 1 n 2 − b 2 n1 b 3 n 2 − b 2 n 3 ⎥
T T T
0
⎢T ⎥
b 1 n 2 − b 2 n1 b 3 n 2 − b 2 n 3 ⎥
T T T
K P = S 0 ⎢⎢ T
0
pour l'élément SHB15
T b 3 n 2 − b 2 n 3 ⎥⎥
T T
⎢b 1 n 2 − b 2 n1 0
⎢b T n − b T n 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢ 1 2 2 1
⎥
T
⎢b 1 n 2 − b 2 n1T 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢ T ⎥
⎢ b 1 n 3 − b 3 n1 b 2 n 3 − b 3 n 2
T T T
0 ⎥
⎢T ⎥
⎢ b 1 n 3 − b 3 n1 b n − b n
T T T
2 3 3 2 0 ⎥
⎢T T ⎥
b 2 n 3 − b 3 n 2
T T
⎢ b 1 n 3 − b 3 n1 0 ⎥
⎢ b T n − b T n b n − b n
T T
0 ⎥
⎢ 1 3 3 1 2 3 3 2 ⎥
⎢ b n − b n
T T
b n 3 − b n 2
T T
0 ⎥
⎢ 1T 3 3 1 2 3
⎥
⎢⎣ b 1 n 3 − b 3 n1
T
b n − b n
T T 0 ⎥⎦
2 3 3 2
et :
98
Formulation des éléments SHB15 et SHB20 en linéaire
⎡ b 2 n1 − b 1 n 2 b 3 n1 − b 1 n 3 ⎤
T T T T
0
⎢ ⎥
b n − b n b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1 1 2
⎢ ⎥
b n1 − b n 2 b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1
⎢ ⎥
b n − b n b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 0 2 1 1 2
⎢ 0 b n1 − b n 2
T T
b 3 n1 − b 1 n 3 ⎥
T T
⎢ 2 1
⎥
⎢ 0 b n − b n
T T
b 3 n1 − b 1 n 3 ⎥
T T
⎢ 2 1 1 2
⎥
⎢ 0 b n1 − b n 2
T T
b 3 n1 − b 1 n 3 ⎥
T T
2 1
⎢ 0 ⎥
b n − b n b 3 n1 − b 1 n 3 ⎥
T T T T
⎢ 2 1 1 2
⎢T T ⎥
b 3 n 2 − b 2 n 3 ⎥
T T
⎢b 1 n 2 − b 2 n1 0
⎢T T ⎥
b 3 n 2 − b 2 n 3 ⎥
T T
⎢b 1 n 2 − b 2 n1 0
⎢b T n − b T n b 3 n 2 − b 2 n 3 ⎥
T T
⎢ 1 2 2 1 0 ⎥
⎢b 1T n 2 − b T2 n1 0 b 3 n 2 − b 2 n 3 ⎥
T T
⎢T T 0 ⎥
b 3 n 2 − b 2 n 3 ⎥
T T
⎢b 1 n 2 − b 2 n1
⎢b T n − b T n 0
b 3 n 2 − b 2 n 3 ⎥⎥
T T
⎢ 1 2 2 1
⎢ b n − b n
T T
b 2 n 3 − b 3 n 2
T T
0 ⎥
⎢ 1 3 3 1
⎥
⎢ b 1 n 3 − b 3 n1
T T
b n − b n
2
T
3
T
3 2 0 ⎥
⎢ T ⎥
⎢ b 1 n 3 − b 3 n1 b n 3 − b n 2
T T T
2 3 0 ⎥
⎢T ⎥
⎢ b 1 n 3 − b 3 n1 b n − b n
T T T
2 3 3 2 0 ⎥
⎢T T ⎥
b n 3 − b n 2
T T
⎢ b 1 n 3 − b 3 n1 2 3 0 ⎥
⎢ b T n − b T n b n − b n
T T
0 ⎥
⎢ 1 3 3 1 2 3 3 2
⎥
⎢ b 1 n 3 − b 3 n1
T T
b 2 n 3 − b 3 n 2
T T 0 ⎥
⎢ T 0 ⎥
⎢⎣ b 1 n 3 − b 3 n1 b n − b n
T T T
2 3 3 2 ⎥⎦
⎢ ⎥
C’est une matrice (18 ×18 ) (pour le SHB15) et ( 24 × 24 ) (pour le SHB20), qu’il faut
multiplier par les déplacements des 6 (pour le SHB15) ou 8 (pour le SHB20) nœuds de la face
sur laquelle on applique une pression.
99
2.3. Validation des éléments coques volumiques SHB15 et SHB20 en
linéaire
Comme l’élément SHB6, nous allons évaluer la performance des éléments SHB15 et
SHB20 sur un ensemble de cas tests linéaires. Pour chaque cas test, le résultat obtenu est
comparé, d’une part, à la solution de référence, et d’autre part, à la solution donnée par les
éléments finis volumiques usuels existants PENTA15 et HEXA20. Tous les maillages n’ont
qu’un élément dans l’épaisseur. Dans tous les tableaux apparaissent le nombre de découpage
dans chacune des directions et le nombre d’éléments total. Comme les données et les
méthodes de calcul des résultats de référence de chaque cas test ont été déjà détaillées dans la
section de validation de l’élément SHB6, nous ne les rappellerons plus ici. Nous allons
présenter directement des résultats numériques de simulation et des remarques pour chaque
modélisation SHB15 ou SHB20.
Nous reprenons les mêmes données du test 1.2.1. Ce cas test permet de vérifier le
comportement de l'élément en flexion simple. En effet dans ce cas test, la flexion est
dominante devant le cisaillement.
PENTA15 SHB15
Maillages
Uz Uz/Uref Uz Uz/Uref
2x1x2 ; 4 6,555E-03 0,895 6,642E-03 0,907
4x2x2 ; 16 7,064E-03 0,965 7,067E-03 0,965
6x2x2 ; 24 7,172E-03 0,980 7,173E-03 0,980
6x4x2 ; 48 7,188E-03 0,982 7,194E-03 0,983
12x4x2 ; 96 7,260E-03 0,992 7,260E-03 0,992
HEXA20 SHB20
Maillages
Uz Uz/Uref Uz Uz/Uref
2x1 ; 2 6,589E-03 0,900 7,085E-03 0,968
4x2 ; 8 6,971E-03 0,952 7,201E-03 0,984
6x2 ; 12 7,196E-03 0,983 7,254E-03 0,991
6x4 ; 24 7,229E-03 0,988 7,271E-03 0,993
12x4 ; 48 7,266E-03 0,993 7,280E-03 0,995
Tous les éléments testés convergent vers la solution de référence. On note que les
modélisations SHB15 et SHB20 convergent un peu plus rapidement que leurs homologues
PENTA15 et HEXA20 (voir Graphe 7).
100
Validation des éléments SHB15 et SHB20 en linéaire
1,02
1
Uz/Uref du point A
0,98
SHB15
0,96
PENTA15
0,94
HEXA20
0,92 Référence
0,9 SHB20
0,88
0 20 40 60 80 100
Nombre des éléments
Graphe 7. Déplacement normalisé du point A de la poutre en flexion simple
Nous reprenons les mêmes données du test 1.2.2. Ce cas test permet de vérifier le
comportement de l'élément en cisaillement. En effet dans ce cas test, le cisaillement est
dominant. Par ailleurs ce cas test a été repris par de nombreux auteurs de la littérature et
notamment [21].
PENTA15 SHB15
Maillages
Uy Uy/Uref Uy Uy/Uref
2x1x2 ; 4 -3,309E-04 0,926 -3,376E-04 0,945
2x2x2 ; 8 -3,368E-04 0,943 -3,425E-04 0,959
3x3x2 ; 18 -3,465E-04 0,970 -3,515E-04 0,984
6x4x2 ; 48 -3,538E-04 0,990 -3,563E-04 0,997
12x4x2 ; 96 -3,567E-04 0,998 -3,582E-04 1,003
HEXA20 SHB20
Maillages
Uy Uy/Uref Uy Uy/Uref
2x1 ; 2 -3,348E-04 0,937 -3,453E-04 0,966
2x2 ; 4 -3,408E-04 0,954 -3,528E-04 0,987
3x2 ; 6 -3,482E-04 0,975 -3,552E-04 0,994
6x4 ; 24 -3,545E-04 0,992 -3,585E-04 1,003
12x4 ; 48 -3,571E-04 0,999 -3,591E-04 1,005
101
Validation des éléments SHB15 et SHB20 en linéaire
Nous trouvons que les quatre éléments convergent bien vers la solution de référence.
1,02
1
Uy/Uref du point A
0,98
0,96 SHB15
0,94
PENTA15
SHB20
0,92
Référence
0,9
HEXA20
0,88
0,86
0 10 20 30 40 50 60 70 80 90 100
Nombre d'éléments
Graphe 8. Déplacement normalisé du point A de la plaque en flexion
et cisaillement dans son plan
Nous reprenons les mêmes données du test 1.2.3. Ce cas test permet de vérifier le
comportement de l'élément en flexion et cisaillement. Ce cas test a été repris par de nombreux
travaux de la littérature ([1], [2], [52]).
PENTA15 SHB15
Maillages
Uz Uz/Uref Uz Uz/Uref
2x2x2 ; 8 3,790E-03 0,699 3,887E-03 0,717
3x2x2 ; 12 4,707E-03 0,868 4,763E-03 0,879
4x2x2 ; 16 5,126E-03 0,946 5,153E-03 0,951
5x2x2 ; 20 5,315E-03 0,981 5,329E-03 0,983
6x2x2 ; 24 5,403E-03 0,997 5,411E-03 0,998
HEXA20 SHB20
Maillages
Uz Uz/Uref Uz Uz/Uref
4x2 ; 8 5,272E-03 0,973 5,675E-03 1,047
6x2 ; 12 5,417E-03 0,999 5,534E-03 1,021
8x2 ; 16 5,439E-03 1,003 5,486E-03 1,012
10x2 ; 20 5,440E-03 1,004 5,465E-03 1,008
12x2 ; 24 5,438E-03 1,003 5,452E-03 1,006
102
Validation des éléments SHB15 et SHB20 en linéaire
1,1
Uz/Uref du point A
0,9 Référence
SHB15
0,8 PENTA15
0,7
SHB20
HEXA20
0,6
8 10 12 14 16 18 20 22 24
Nombre d'éléments
Nous constatons que les quatre éléments convergent très bien vers la solution de
référence. Les éléments hexaèdres convergent plus rapidement que les éléments pentaèdres.
Nous reprenons les mêmes données du test 1.2.4. Ce cas test permet de vérifier le
comportement de l'élément en flexion et cisaillement. Ce cas test a été traité par de nombreux
auteurs ([1], [2], [52], [56], [57]). On observe les déplacements du nœud A suivant la
direction Ox. Les résultats obtenus sont donnés dans le Tableau 22.
PENTA15 SHB15
Maillages
Ux Ux/Uref Ux Ux/Uref
3x(5x5x2) ; 150 2,115E-02 0,229 1,557E-02 0,168
3x(10x10x2) ; 600 7,345E-02 0,795 6,544E-02 0,708
3x(15x15x2) ; 1350 8,814E-02 0,954 8,631E-02 0,934
3x(20x20x2) ; 2400 9,069E-02 0,981 8,999E-02 0,974
3x(25x25x2) ; 3750 9,179E-02 0,993 9,162E-02 0,992
HEXA20 SHB20
Maillages
Ux Ux/Uref Ux Ux/Uref
3x(1x1) ; 3 9,928E-05 0,001 1,058E-02 0,115
3x(2x2) ; 12 1,287E-03 0,014 5,213E-02 0,564
3x(3x3) ; 27 5,984E-03 0,065 8,564E-02 0,927
3x(4x4) ; 48 1,635E-02 0,177 9,099E-02 0,985
3x(5x5) ; 75 3,161E-02 0,342 9,209E-02 0,997
103
Validation des éléments SHB15 et SHB20 en linéaire
1,2
Uz/Uref du point A
0,8
0,6
SHB15
0,4 PENTA15
0,2 Référence
0
0 1000 2000 3000 4000
Nombre d'éléments
1,2
1
Ux/Uref du point A
0,8 SHB20
0,6
HEXA20
Référence
0,4
0,2
0
0 20 40 60 80
Nombre d'éléments
Graphe 10. Déplacement normalisé du point A suivant Ox de l’hémisphère pincé
Nous reprenons les mêmes données du test 1.2.6. Ce cas test permet de vérifier le
comportement de l'élément en flexion et cisaillement. Ce test a été étudié par plusieurs auteurs
([4], [48], [57], [79]).
Nous trouvons que les éléments SHB15 et PENTA15 convergent de la même manière
vers la solution de référence mais plus lentement que les éléments SHB20 et HEXA20. En
plus, l’élément SHB20 converge nettement mieux que le HEXA20 (voir Graphe 11).
104
Validation des éléments SHB15 et SHB20 en linéaire
PENTA15 SHB15
Maillages
Uz Uz/Uref Uz Uz/Uref
10x10x2 ; 200 1,141E-05 0,625 1,179E-05 0,646
15x15x2 ; 450 1,530E-05 0,838 1,667E-05 0,913
20x20x2 ; 800 1,677E-05 0,919 1,736E-05 0,951
25x25x2 ; 1250 1,745E-05 0,956 1,817E-05 0,996
30x30x2 ; 1800 1,781E-05 0,976 1,818E-05 0,996
HEXA20 SHB20
Maillages
Uz Uz/Uref Uz Uz/Uref
4x4 ; 16 2,554E-06 0,140 1,611E-05 0,883
6x6 ; 36 5,987E-06 0,328 1,753E-05 0,961
8x8 ; 64 9,539E-06 0,523 1,787E-05 0,979
10x10 ; 100 1,233E-05 0,675 1,805E-05 0,989
12x12 ; 144 1,419E-05 0,777 1,818E-05 0,996
1
Uz/Uref du point A
0,9
SHB15
0,8 PENTA15
Référence
0,7
0,6
200 700 1200 1700
Nombre d'éléments
0,8
Uz/Uref du point A
0,6
0,4 SHB20
HEXA20
0,2
Référence
0
0 50 100 150
Nombre d'éléments
105
Validation des éléments SHB15 et SHB20 en linéaire
Nous reprenons les mêmes données du test 1.2.7. Ce cas test permet de vérifier le
comportement de l'élément en flexion et cisaillement. Ce test a été repris dans de nombreux
travaux, notamment par [21].
PENTA15 SHB15
Maillages
Uz Uz/Uref Uz Uz/Uref
3x(2x2x2) ; 24 -2,332E-05 0,878 -2,518E-05 0,948
3x(5x5x2) ; 150 -2,552E-05 0,960 -2,609E-05 0,982
3x(8x8x2) ; 384 -2,579E-05 0,971 -2,621E-05 0,986
3x(11x11x2) ; 726 -2,602E-05 0,979 -2,636E-05 0,992
3x(14x14x2) ; 1176 -2,615E-05 0,984 -2,646E-05 0,996
HEXA20 SHB20
Maillages
Uz Uz/Uref Uz Uz/Uref
3x(2x2) ; 12 -2,427E-05 0,913 -2,543E-05 0,957
3x(5x5) ; 75 -2,556E-05 0,962 -2,584E-05 0,972
3x(8x8) ; 192 -2,585E-05 0,973 -2,609E-05 0,982
3x(11x11) ; 363 -2,606E-05 0,981 -2,630E-05 0,990
3x(14x14) ; 588 -2,623E-05 0,987 -2,647E-05 0,996
1,02
1
0,98
Uz/Uref du point A
0,96 SHB15
0,94 PENTA15
0,92 SHB20
0,9 HEXA20
0,88 Référence
0,86
0 200 400 600 800 1000 1200
Nombre d'éléments
Graphe 12. Déplacement normalisé du point A suivant Oz
de la plaque circulaire soumise à un effort ponctuel
Nous remarquons que les éléments SHB quadratiques convergent vers la solution de
référence un peu mieux que leurs homologues purement 3D (voir Graphe 12).
106
Validation des éléments SHB15 et SHB20 en linéaire
Nous reprenons les mêmes données du test 1.2.8. Ce test a pour objectif de tester
l’élément SHB15 et SHB20 en vibration. La matrice de masse de ces éléments est prise
comme celle d’un élément solide 3D.
Nous cherchons à déterminer les différentes fréquences d’une poutre dont une extrémité
est encastrée, comprises entre 0 et 150 Hz.
Nous trouvons que les éléments SHB15 et SHB20 donnent de très bons résultats.
107
Validation des éléments SHB15 et SHB20 en linéaire
Nous reprenons les mêmes données du test 1.2.9. Ce test représente un calcul de
stabilité d’une enveloppe cylindrique mince libre à ses extrémités soumises à une pression
externe. On calcule les charges critiques conduisant au flambement élastique d’Euler. La
matrice de rigidité géométrique utilisée dans la résolution du problème aux valeurs propres est
celle qui est due aux contraintes initiales.
Les résultats obtenus pour la pression critique de flambement de ce cylindre libre sous
pression externe sont donnés dans le Tableau 25.
Nous trouvons que les éléments SHB15 et SHB20 convergent bien vers la solution de
référence. L’élément SHB20 converge plus rapidement que l’élément SHB15.
Tableau 25. Pression critique de flambage du cylindre libre sous pression externe
108
Validation des éléments SHB15 et SHB20 en linéaire
Les déplacements obtenus sont comparés à la solution analytique élastique d’une poutre
en flexion. Ce test permet également de montrer les limites des éléments SHB en termes
d’élancement.
G
À l’extrémité droite, la poutre est soumise à un effort tranchant de P = 1 N dans la
direction verticale Oz. Ce chargement est réparti de façon adéquate aux nœuds de l’extrémité
droite. On simule un encastrement parfait à l'extrémité gauche de la poutre (voir Figure 23).
Les données de géométrie, matériau et de chargement du problème sont listées dans le
Tableau 26.
Longueur L = 1000 m
Largeur B = 100 m
Épaisseurs variables e = 100 à 0,2 m
Module d’Young E = 2x1011 Pa
Coefficient de Poisson ν = 0,3
G
Chargement P =1 N
Nous trouvons que les éléments SHB20, SHB15 et SHB8PS convergent très bien vers la
solution de référence et ils peuvent arriver à un élancement (rapport largeur/épaisseur) de 500
pour un maillage grossier de 10 ou 20 éléments seulement. L’élément SHB6 admet une limite
d’élancement nettement plus faible que les autres éléments SHB. Son élancement peut
atteindre 2 seulement.
109
Validation des éléments SHB15 et SHB20 en linéaire
z
x y
L A e =100 m
e =10 m
e =1 m
e =0,5 m
e =0,2 m
Figure 23. Géométrie, chargement et conditions aux limites pour le test de la poutre en
flexion avec divers élancements ; un exemple de maillage (10x1x1) SHB20
Tableau 27. Déplacement du point A suivant Oz de la poutre en flexion simple avec divers
élancements
110
Validation des éléments SHB15 et SHB20 en linéaire
Dans cette analyse linéaire, la pression critique d’Euler est déterminée ainsi que le mode
de flambage correspondant. Cet état critique est associé à la pression minimale à laquelle la
matrice de rigidité globale devient singulière et, il est obtenu par la résolution des valeurs
propres du problème suivant :
( K e + λc K σ ) ⋅ Xc = 0 (113)
111
Validation des éléments SHB15 et SHB20 en linéaire
z y
sym
x
sy
m
a) Géométrie initiale
sym
Point A
b) 1er mode de flambage
hr
L
Lr
La
ha
sym
Figure 24. Géométrie, conditions aux limites et 1er mode de flambage d’un quart de coque
cylindrique avec raidisseur ; un exemple de maillage mixte de 260 éléments SHB8PS pour le
raidisseur et 360 éléments SHB6 pour la coque principale
112
Validation des éléments SHB15 et SHB20 en linéaire
a) SHB15 b) SHB20
Figure 25. Géométrie initiale et 1er mode de flambage de la coque cylindrique avec
raidisseur ; un exemple de maillage de : a) 300 éléments SHB15 ; b) 100 éléments SHB20
113
Chapitre III
Éléments coques volumiques SHB en non
linéaire
1. Non-linéarités géométriques
On traite ici le cas des grands déplacements, mais rotations faibles et petites
déformations. On adopte pour cela une formulation lagrangienne mise à jour.
En non linéaire, nous cherchons à écrire l’équilibre entre forces internes et forces
externes à la fin de l’incrément de charge (repéré par l’indice 2) :
F2int = F2extr
F2int = ∫B ⋅ σ 2 dV
T
2
Ω2
Remarques importantes :
⎡ b1T + hα , x1 γ αT 0 0 ⎤
⎢ ⎥
⎢ 0 b T2 + hα , x2 γ αT 0 ⎥
⎢ ⎥
⎢ 0 0 b T3 + hα , x3 γ αT ⎥
B 2 = ⎢ bT + h γ T b1T + hα , x1 γ αT 0 ⎥ , avec c = 0,45
α , x2 α
⎢ 2 ⎥
⎢ (
⎢c bT + h γ T
3 α , x3 α ) 0 (c b1T + hα , x1 γ αT ) ⎥
⎥
⎢
⎣⎢
0 (
c b T3 + hα , x3 γ αT ) c (b T
2 + hα , x2 γ αT ) ⎥
⎦⎥
⎡ b1T + hα , x1 γ αT 0 0 ⎤
⎢ ⎥
⎢ 0 b T2 + hα , x2 γ αT 0 ⎥
⎢ ⎥
⎢ 0 0 b T3 + hα , x3 γ αT ⎥
B 2 = ⎢ bT + h γ T b1T + hα , x1 γ αT 0 ⎥
α , x2 α
⎢ 2 ⎥
(
⎢ bT + h γ T
⎢ 3 α , x3 α ) 0 (b T
1
T
+ hα , x1 γ α ) ⎥
⎥
⎢ ⎥
⎢⎣ 0 (b T
3 + hα , x3 γ αT ) (b T
2 + hα , x2 γ αT ) ⎥⎦
1
∆E =
2
( ∇1 ( ∆u ) + ∇1T ( ∆u ) )
L’opérateur gradient est calculé sur la géométrie de début de pas. Cette écriture de la
déformation est limitée aux petites rotations (<5 degrés).
On peut sans difficulté étendre la formulation aux grandes rotations en incluant dans la
déformation les termes de second ordre :
1
∆E =
2
( ∇1 ( ∆u ) + ∇1T ( ∆u ) + ∇1T ( ∆u ) ⋅∇1 ( ∆u ) )
∆π = C ' : ∆E
où C ' est la matrice de Hooke. Remarquons que pour les éléments SHB, cette matrice est une
matrice orthotrope transverse qui s’écrit dans les axes du lamina :
⎡λ + 2µ µ 0 0 0 0⎤
⎢ µ λ + 2µ 0 0 0 0 ⎥⎥
⎢
⎢ 0 0 E 0 0 0⎥
C' = ⎢ ⎥
⎢ 0 0 0 µ 0 0⎥
⎢ 0 0 0 0 µ 0⎥
⎢ ⎥
⎣⎢ 0 0 0 0 0 µ ⎦⎥
115
Formulation des éléments SHB en non linéaire
⎧ π 2 = σ1 + ∆π
⎪
⎪ 1
⎨σ 2 = FT ⋅ π 2 ⋅ F
⎪ det ( F )
⎪ F = I + ∇ ( ∆u )
⎩ 1
La combinaison des quatre dernières équations avec l’expression des forces internes
donne la formulation de l’élément en grandes déformations en lagrangien mis à jour.
Remarquons que cette formulation lagrangienne mise à jour est totalement équivalente à
la formulation lagrangienne totale pour laquelle les forces internes s’écrivent :
∫ (B + B )
T
F2int = NL
(u) ⋅ π 2 dV
0
Ω0
Dans ce cas toutes les intégrations sont faites sur la géométrie initiale Ω 0 la contrainte
π 2 utilisée est la contrainte de Piola–Kirchhoff II. Cette dernière méthode est probablement
préférable lorsque le maillage se déforme significativement et permet donc de traiter les
problèmes en grandes déformations mais nécessite le développement de l’opérateur B NL ( u ) .
1
∆E =
2
(
∇ 0 ( ∆u ) + ∇T0 ( ∆u ) + ∇T0 ( ∆u ) ⋅∇ 0 ( ∆u ) )
La combinaison des deux équations précédentes donne la formulation de l’élément en
grandes déformations en comportement linéaire matériau.
Dans le cas des petits déplacements, on confond géométrie en début et fin de pas,
contrainte de Cauchy et de Piola–Kirchhoff II, de plus on utilise l’expression linéaire des
déformations.
2. Non-linéarités matériaux
Comme les éléments SHB ont un comportement 3D orthotrope particulier, par
conséquent, ils ne sont pas utilisables pour un matériau de comportement non-linéaire
quelconque disponible dans le code Aster. Afin d’éviter cette limite d’utilisation, nous
proposons une méthode de construction particulière de la matrice tangente CT qui permet de
connecter ces éléments SHB à toutes les lois de comportement disponibles dans le code Aster.
Nous voulons déterminer la matrice CT de la relation suivante :
π 2 = σ1 + CT : ( ∆ε − ∆ε p )
où :
116
Formulation des éléments SHB en non linéaire
⎡π 2 xx ⎤ ⎡σ 1xx ⎤ ⎡ ∆ε xx − ∆ε xxp ⎤
⎢π ⎥ ⎢σ ⎥ ⎢ p ⎥
⎢ 2 yy ⎥ ⎢ 1 yy ⎥ ⎢ ∆ε yy − ∆ε yy ⎥
⎢π 2 zz ⎥ ⎢σ 1zz ⎥ ⎢ ∆ε zz − ∆ε zzp ⎥
π2 = ⎢ ⎥=⎢ ⎥+C
T
⎢ p ⎥
π σ
⎢ 2 xy ⎥ ⎢ 1xy ⎥ ⎢ ∆ε xy − ∆ε xy ⎥
⎢π 2 xz ⎥ ⎢σ 1xz ⎥ ⎢ ∆ε − ∆ε p ⎥
⎢ ⎥ ⎢ ⎥ ⎢ xz xz
⎥
⎢⎣π 2 yz ⎥⎦ ⎢⎣σ 1 yz ⎥⎦ ⎢⎣ ∆ε yz − ∆ε yz ⎥⎦
p
Nous supposons d’abord que l’élément est en état de contrainte plane dans le repère
local de chaque point d’intégration de Gauss et les déformations hors plan sont élastiques.
Cela entraîne alors immédiatement que les déformations totales hors plan sont égales aux
déformations élastiques. Nous avons donc :
Ensuite, les contraintes hors plan sont calculées de façon élastique. Nous avons donc :
⎡Cxxxx
CPT CPT
Cxxyy CPT
0 Cxxxy 0 0⎤
⎢ CPT CPT CPT ⎥
⎢C yyxx C yyyy 0 C yyxy 0 0⎥
⎢ 0 0 E 0 0 0⎥
CT = ⎢ CPT CPT CPT ⎥
⎢Cxyxx Cxyyy 0 Cxyxy 0 0⎥
⎢ 0 0 0 0 µ 0⎥
⎢ ⎥
⎣⎢ 0 0 0 0 0 µ ⎦⎥
π 2 = σ1 + CT : ( ∆ε − ∆ε p )
où :
117
Formulation des éléments SHB en non linéaire
⎢π ⎥ ⎢σ ⎥ ⎢ CPT ⎥ ⎢ p ⎥
⎢ 2 yy ⎥ ⎢ 1 yy ⎥ ⎢C yyxx
CPT
C yyyy CPT
0 C yyxy 0 0 ⎥ ⎢ ∆ε yy − ∆ε yy ⎥
⎢π 2 zz ⎥ ⎢σ 1zz ⎥ ⎢ 0 0 E 0 0 0 ⎥ ⎢ ∆ε zz ⎥
⎢ ⎥=⎢ ⎥ + ⎢ CPT ⎥⋅⎢ p ⎥
⎢π 2 xy ⎥ ⎢σ 1xy ⎥ ⎢Cxyxx 0 ⎥ ⎢ ∆ε xy − ∆ε xy ⎥
CPT CPT
Cxyyy 0 Cxyxy 0
⎢π 2 xz ⎥ ⎢σ 1xz ⎥ ⎢ 0 0 0 0 µ 0 ⎥ ⎢ ∆ε xz ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢⎣π 2 yz ⎥⎦ ⎢⎣σ 1 yz ⎥⎦ ⎢⎣ 0 0 0 0 0 µ ⎥⎦ ⎢⎣ ∆ε yz ⎥⎦
Nous devons trouver la contrainte à la fin du pas qui vérifie l’équilibre. Cette équation
est résolue dès que l’incrément de déformation plastique est connu. Cette déformation se
détermine en imposant à la contrainte finale d’être plastiquement admissible. Cette méthode
est complètement similaire à la méthode tridimensionnelle usuelle sauf qu’il n’y a pas de
solution explicite à ce problème si on utilise par exemple l’approximation efficace du retour
radial pour calculer la solution. Nous avons choisi de résoudre ce problème non linéaire par
une méthode de Newton.
Cette méthode permet ainsi de connecter simplement les éléments SHB à toutes les lois
de comportement disponibles dans le Code Aster.
118
3. Cas tests non linéaires géométriques et matériaux
Dans cette section, nous allons tester les éléments SHB sur un ensemble de cas tests non
linéaires géométriques et/ou matériaux. Pour chaque cas test, le résultat obtenu est comparé, à
la solution de référence. Tous les maillages n’ont qu’un élément dans l’épaisseur, sauf
nécessité particulière motivée par l’application de conditions aux limites appropriées dans le
cas de flambage par point limite. Dans tous les tableaux apparaissent le nombre de découpage
dans chacune des directions et le nombre total d’éléments.
Le cas test de la poutre console est un cas test élastique non-linéaire géométrique
(grands déplacements). C’est un cas test répandu dans la littérature scientifique sur les
éléments finis et qui a été traité notamment dans la référence [80].
La poutre est encastrée à une extrémité. Sur l’autre extrémité libre une force concentrée
P est appliquée suivant la direction z. Cette force est appliquée par incréments successifs
pour aller de de 0 à Pmax :
EI
P0 = =1 N et Pmax = 4P0
L2
Longueur L = 10 m
Largeur b=1m
Épaisseur h = 0,1 m
Module d’Young E = 1,2x106 Pa
Coefficient de Poisson ν=0
Nous avons testé les trois éléments SHB6, SHB15 et SHB20 sur plusieurs maillages
différents. Nous avons retenu trois maillages qui donnent de bonnes convergences par rapport
à la solution de référence (moins de 1% d’erreur). Ces résultats sont représentés dans le
Tableau 32 et sur le Graphe 13.
Validation des éléments SHB en non linéaire
Tableau 31. Solution de référence pour le cas test de la poutre console en non-linéaire
géométrique
Tableau 32. Résultats obtenus pour le cas test de la poutre console en non-linéaire
géométrique
Nous trouvons que les trois éléments SHB convergent bien vers la solution de référence.
L’élément SHB6 converge le moins rapidement. Il nous faut un maillage de 2000 éléments
SHB6 pour avoir une convergence avec moins de 1% d’erreur.
120
Validation des éléments SHB en non linéaire
1,2 Référence - Ux
SHB6 -Ux
SHB15 -Ux
1
SHB20 -Ux
Référence Uz
0,8 SHB6 Uz
SHB15 Uz
P/Pmax 0,6 SHB20 Uz
0,4
0,2
0
0 1 2 3 4 5 6 7
Déplacement de l'extrémité suivant Uz ou -Ux
Graphe 13. Résultats du cas test de la poutre console en non-linéaire géométrique
P P
z z
y A b y A b
x x
h h
L L
1er cas : le matériau utilisé est élastique. Donc, nous avons un problème de
•
flambage élastique (par point limite) non-linéaire géométrique ;
• 2ème cas : le matériau utilisé est élasto-plastique. Donc, nous avons un
problème non-linéaire géométrique et matériau.
Nous présentons d’abord le 1er cas. Les données géométriques et matériau sont
représentées dans le Tableau 33 et Figure 27a.
121
Validation des éléments SHB en non linéaire
Ce cas test a été traité par Sze et al. [80] ou Klinkel et Wagner [51]. Comme le
problème est symétrique, nous étudions seulement un quart de la structure. Le panneau est en
appui simple sur le côté BC (fibre neutre du panneau), libre sur la côté CD et soumis à une
force ponctuelle P au point A suivant la direction verticale Oz (voir Figure 27a). Pour pouvoir
tenir compte des bonnes conditions aux limites du problème, nous maillons le panneau avec
deux couches d’éléments dans l’épaisseur. Nous étudions ici trois maillages :
libre sy
m
C P
sim
pl
ea
pp A
ui
sym
L z
B y
h x
122
Validation des éléments SHB en non linéaire
3500
Référence Klinkel
3000
SHB6
2500 SHB15
Force P (N)
2000 SHB20
1500
1000
500
0
0 5 10 15 20 25 30 35
Déplacement suivant Oz du point A (mm)
Graphe 14. Déplacement du point A du panneau cylindrique épais soumis à une force
ponctuelle
Les résultats obtenus du déplacement du point A suivant Oz pour les trois maillages
sont donnés dans le Tableau 34, et sur le Graphe 14. Nous trouvons que les trois éléments
SHB donnent des résultats assez proches de ceux de la référence jusqu’au moment où le
déplacement du point A atteint Uz = 16mm ensuite, le chargement P est un peu sous-estimé.
Déplacement Référence
SHB6 SHB15 SHB20
du point A de Klinkel
Uz (mm) P (N) P (N) P (N) P (N)
2 706 651,97 674,49 698,049
4 1273 1189,92 1229,64 1263,31
6 1707 1616,83 1668,66 1701,51
8 2007 1930,47 1989,55 2012,2
10 2160 2119,03 2180,48 2183,92
12 2129 2151,3 2208,97 2183,07
14 1827 1948,47 1989,13 1922,24
16 1180 1316,44 1282,1 1203,39
18 677 527,37 414,75 440,53
20 592 339,53 240,64 277,26
21,78 690,9 458,32 374,8 413,39
22,875 885 612,62 539,53 582,19
24,049 1161,3 839,38 778,51 829,63
25,293 1527,9 1147,28 1100,5 1167,06
26,601 1993,2 1534,22 1503,03 1593,78
27,964 2564,7 2002,76 1988,04 2114,07
28,663 2892,9 2275,99 2269,37 2419,56
29,374 3205,5 2569,49 2570,7 2748,76
123
Validation des éléments SHB en non linéaire
Maintenant, nous traitons le 2ème cas. Nous gardons les mêmes données géométriques et
les mêmes conditions aux limites du 1er cas. Nous n’avons pas de solution de référence
existante. Donc, nous allons établir une solution de référence en faisant un calcul obtenu avec
un maillage de (20x20) éléments coques S4R5 du logiciel Abaqus. Nous utilisons un matériau
élasto-plastique de von Mises à écrouissage isotrope linéaire dont les paramètres sont définis
dans le Tableau 35. Nous calculons le déplacement vertical suivant Oz du point d’application
A en fonction de la force ponctuelle P (P varie de 0 à Pmax = 3000N). Ces résultats obtenus
sont pris comme des résultats de référence (voir le Tableau 37 et le Graphe 16).
Nous utilisons le même matériau élasto-plastique que nous avons utilisé pour établir la
solution de référence. Mais ses paramètres sont définis par une courbe de traction qui est
donnée dans le Tableau 36 et le Graphe 15.
124
Validation des éléments SHB en non linéaire
15
Sigma (Mpa)
10
0
0 0,5 1 1,5 2
Epsilon
3000
Reference
2500
SHB15
2000 SHB20
SHB8
P (N)
1500
SHB6
1000
500
0
0 5 10 15 20 25 30 35 40
Déplacement -Uz
Graphe 16. Déplacement du point A du panneau cylindrique épais soumis à une force
ponctuelle
Les résultats obtenus du déplacement du point A suivant Oz pour les quatre maillages
sont donnés dans le Tableau 37 et le Graphe 16. Nous trouvons que les quatre éléments SHB
donnent des résultats très proches de ceux de référence.
125
Validation des éléments SHB en non linéaire
126
Validation des éléments SHB en non linéaire
Ce cas test a été traité par plusieurs auteurs et en particulier dans Sze et al. [80]. Ce test
représente un calcul de stabilité d’un panneau cylindrique simplement supporté soumis à un
effort concentré en son centre. Le comportement du panneau change complètement et montre
nettement des points de retour en charge et en déplacement « snap-through/snap-back ». Dans
ce cas un pilotage en déplacement diverge et un pilotage en longueur d’arc doit être choisi. Il
permet de valider les modélisations SHB dans le domaine quasi-statique non-linéaire
géométrique en présence de fortes instabilités.
127
Validation des éléments SHB en non linéaire
D
libre sy
m
P
C
sim
ple
ap A
pu
i sym
L
B h z
y
x
0,8
0,6 Référence
SHB6
SHB8PS
0,4
P/Pmax
0,2
0
0 5 10 15 20 25 30 35
-0,2
Déplacement du point A suivant -Uz (mm)
Graphe 17. Déplacement du point A du panneau cylindrique mince soumis à une force
ponctuelle
128
Validation des éléments SHB en non linéaire
0,8
0,6
Référence
0,4 SHB15
P/Pmax SHB20
0,2
0
0 5 10 15 20 25 30 35
-0,2
Déplacement du point A suivant -Uz (mm)
Graphe 18. Déplacement du point A du panneau cylindrique mince soumis à une force
ponctuelle
129
Validation des éléments SHB en non linéaire
Les résultats obtenus du déplacement du point A suivant Oz pour les quatre maillages
sont représentés sur les Graphe 17 et Graphe 18Graphe 20. Nous trouvons que les quatre
éléments SHB donnent des résultats bien cohérents avec ceux de référence.
Ce cas test a été traité par Sze et al. [80]. Comme le problème est symétrique, nous
étudions seulement un huitième de la structure. Par des conditions de symétrie, les
déplacements des surfaces de symétrie sont nuls suivant la direction de leurs normales :
⎛ L ⎞
• Face AB : v ⎜ x, , z ⎟ = 0 ∀x, z
⎝ 2 ⎠
• Face AD : u ( 0, y, z ) = 0 ∀y, z
• Face BC : w ( x, y, 0 ) = 0 ∀x, y
La coque est libre sur la côté CD et soumise à une force ponctuelle P au point A suivant
la direction verticale Oz (voir Figure 29a). La force P est appliquée par incréments successifs
pour aller de 0 à Pmax = 40000. Nous étudions ici trois maillages :
130
Validation des éléments SHB en non linéaire
Pmax/4 = 10000.
libre D
h
P/4
sym
A
m
sy
C z
x
L/2 y
s ym
B
R
Les résultats obtenus sont donnés dans le Tableau 41 et sur les Graphe 19, Graphe 20,
Graphe 21. Nous trouvons que trois éléments SHB donnent des résultats bien cohérents avec
ceux de référence.
131
Validation des éléments SHB en non linéaire
Tableau 41. Résultats obtenus des déplacements des points A, B et C de la coque cylindrique
soumise à une force ponctuelle
132
Validation des éléments SHB en non linéaire
1,1
1
0,9 Wa reference
0,8 -Ub reference
0,7 -Uc reference
P/Pmax
0,6
Wa SHB6
0,5
-Ub SHB6
0,4
0,3
-Uc SHB6
0,2
0,1
0
0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5
Déplacements
1,1
1
0,9 Wa reference
0,8 -Ub reference
0,7 -Uc reference
P/Pmax
0,6
Wa SHB15
0,5
-Ub SHB15
0,4
-Uc SHB15
0,3
0,2
0,1
0
0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5
Déplacements
Graphe 20. Déplacement du point A, B et C de la coque cylindrique soumise à une force
ponctuelle ; un exemple de maillage SHB15 de (25x25x1)x2
133
Validation des éléments SHB en non linéaire
1,1
1
0,9 Wa reference
0,8 -Ub reference
0,7 -Uc reference
P/Pmax
0,6
Wa SHB20
0,5
-Ub SHB20
0,4
0,3
-Uc SHB20
0,2
0,1
0
0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5
Déplacements
134
Conclusions et perspectives
Conclusions et perspectives
La présente étude entre dans le cadre général du développement et de la validation
d’une famille d’éléments finis de solide–coques construite à partir d’une formulation
purement 3D avec trois degrés de liberté en déplacement à chaque nœud. Nous avons
développé 3 éléments finis solide–coques SHB6, SHB15 et SHB20. Ces éléments ont été bien
implantés dans le code des éléments finis implicite ASTER.
L’élément SHB6 est un prisme à six nœuds et il a une direction privilégiée appelée
« épaisseur » qui est normale au plan moyen du prisme. L’intégration numérique réduite est
utilisée pour économiser du temps de calcul. L’intégration à travers l’épaisseur s’appuie sur
cinq points de Gauss. L’intégration dans les deux autres directions est faite avec un seul point
de Gauss situé au centre de gravité de l’élément. Une analyse des modes de hourglass pour
cet élément montre que cet élément ne présente pas de mode de parasite à énergie nulle. La
stabilisation n’est donc pas nécessaire. Afin d’améliorer sa rapidité de convergence et de
réduire les blocages en membrane ou en cisaillement, nous avons appliqué une projection de
type « Assumed Strain » à l’opérateur gradient de l’élément. L’élément SHB6 projeté est
dénommé SHB6. Afin d’évaluer ses performances, une série de cas tests linéaires et non-
linéaires standards dans la littérature lui a été appliquée. À travers ces exemples numériques,
l’élément SHB6 montre des améliorations significatives par rapport à l’élément
tridimensionnel à six nœuds nommé PRI6. En plus, ce type d’élément se connecte
naturellement à l’élément solide–coque à huit-nœuds SHB8PS. Cela permet de modéliser
facilement des structures à géométrie complexe. Surtout dans le cas où l’on ne peut pas
utiliser l’élément hexaédrique SHB8PS seul pour modéliser la structure. Pourtant, la
convergence de l’élément SHB6 reste lente dans le cas test d’une coque hémisphérique
pincée, qui est probablement dû au fait que les déformations des éléments triangulaires
linéaires sont constantes. De même, sa limite d’élancement (rapport largeur/épaisseur de la
structure) reste faible autour de 2. Il est donc, conseillé de mailler la structure avec l’élément
solide–coque hexaédrique à huit-nœuds SHB8PS partout où c’est possible, et de n’utiliser
l’élément SHB6 que pour mailler le reste.
L’élément SHB20 est un hexaèdre à vingt nœuds, et il a aussi une direction privilégiée
appelée « épaisseur » qui est normale au plan moyen de l’hexaèdre. L’intégration numérique
réduite est utilisée. Cette intégration utilise vingt points de Gauss au total : quatre dans le plan
de l’élément et cinq à travers l’épaisseur.
135
(voir 2.3.4) et la coque cylindrique pincée avec diaphragmes (voir 2.3.5). En plus, les
éléments SHB15 et SHB20 peuvent se connecter naturellement ; et leurs limites d’élancement
sont élevées jusqu’à 500. Nous pouvons donc les utiliser pour modéliser des structures minces
et avec des formes géométriques complexes.
Cette famille d’éléments finis solide–coques SHB (SHB6, SHB8PS, SHB15 et SHB20)
est actuellement connectée à toutes les lois de comportement disponibles du code ASTER.
Dans la suite de ce travail, il sera intéressant d’évaluer la performance de cette famille en
termes de convergence et de temps de calcul par rapport aux autres éléments 3D et aux autres
éléments coques classiques pour une série plus complète de cas tests non-linéaires.
136
Bibliographies
Bibliographie
[1] Abed-Meraim F, Combescure A. SHB8PS a new intelligent assumed strain continuum
mechanics shell element for impact analysis on a rotating body. First M.I.T. Conference on
Computational Fluid and Solid Mechanics, U.S.A., 2001.
[2] Abed-Meraim F, Combescure A. SHB8PS, a new adaptive, assumed strain continuum
mechanics shell element for impact analysis. Computers and Structures, 80:791–803, 2002.
[3] Ahmad S. Analysis of thick and thin shell structures by curved finite elements.
International Journal for Numerical Methods in Engineering, 2:419–451, 1970.
[4] Alves de Sousa R.J, Cardoso R.P.R, Fontes Valente R.A, Yoon J.W, Gracio J.J, Natal
Jorge R.M. A new one-point quadrature enhanced assumed strain (EAS) solid–shell element
with multiple integration points along thickness: Part I – Geometrically linear applications.
International Journal for Numerical Methods in Engineering, 62:952–977, 2005.
[5] Alves de Sousa R.J, Cardoso R.P.R, Fontes Valente R.A, Yoon J.W, Gracio J.J, Natal
Jorge R.M. A new one-point quadrature enhanced assumed strain (EAS) solid–shell element
with multiple integration points along thickness: Part II – Nonlinear applications.
International Journal for Numerical Methods in Engineering, 67:160–188, 2006.
[6] Alves de Sousa R.J, Yoon J.W, Cardoso R.P.R, Fontes Valente R.A, Gracio J.J. On
the use of a reduced enhanced solid–shell (RESS) element for sheet forming simulations.
International Journal of Plasticity, 23:490–515, 2007.
[7] Ashwell D.G, Sabir A.B. A new cylindrical shell finite element based on simple
independent strain function. Int. J. Mech. Science, 14:171–185, 1972.
[8] Ayad R, Batoz J.L, Dhatt G. Un élément quadrangulaire de plaque basé sur une
formulation mixte-hybride avec projection en cisaillement. Revue Européenne des Eléments
Finis, 4:415–440, 1995.
[9] Ayad R, Dhatt G, Batoz J.L. A new hybrid-mixed variational approach for Reissner–
Mindlin plates: The MiSP model. International Journal for Numerical Methods in
Engineering, 42:1149–1179, 1998.
[10] Bachrach W.E. An efficient formulation of hexahedral elements with high accuracy
for bending and incompressibility. Computers and Structures, 26:453–467, 1987.
[11] Bathe K.J, Dvorkin E.N. Four-node plate bending element based on Mindlin/Reissner
plate theory and a mixed interpolation. International Journal for Numerical Methods in
Engineering, 21:367–383, 1985.
[12] Bathe K.J, Dvorkin E.N. Formulation of general shell elements – The use of mixed
interpolation of tensorial components. International Journal for Numerical Methods in
Engineering, 22:697–722, 1986.
[13] Bathe K.J, Ho L.W. A simple and effective element for analysis of general shell
structures. Computers and Structures, 13:673–681, 1981.
[14] Batoz J.L, Bathe K.J, Ho L.W. A study of three node triangular plate bending
elements. International Journal for Numerical Methods in Engineering, 15:1771–1812, 1980.
[15] Batoz J.L, Ben Tahar M. Evaluation of a new thin plate quadrilateral element.
International Journal for Numerical Methods in Engineering, 18:1655–1678, 1982.
137
[16] Batoz J.L, Dhatt G. Modélisation des structures par éléments finis. Vol 1, Solides
élastiques, Editions Hermès, 1990.
[17] Batoz J.L, Dhatt G. Modélisation des structures par éléments finis. Vol 2, Poutres et
plaques, Editions Hermès, 1990.
[18] Batoz J.L, Dhatt G. Modélisation des structures par éléments finis. Vol 3, Coques,
Editions Hermès, 1990.
[19] Batoz J.L, Dhatt G. Revue et bilan des éléments de plaque de type Kirchhoff discret.
Calcul des Structures et I.A, Vol 2, Fouet et al. Eds. Pluralis, 137–160, 1988.
[20] Belytschko T. Correction of article by Flanagan D.P, Belytschko T. International
Journal for Numerical Methods in Engineering, 19:467–468, 1983.
[21] Belytschko T, Bindeman L.P. Assumed strain stabilization of the 4-node quadrilateral
with 1-point quadrature for nonlinear problems. Computer Methods in Applied Mechanics and
Engineering, 88:311–340, 1991.
[22] Belytschko T, Bindeman L.P. Assumed strain stabilization of eight node hexahedral
element. Computer Methods in Applied Mechanics and Engineering, 105:225–260, 1993.
[23] Belytschko T, Ong J.S.J, Liu W.K, Kennedy J.M. Hourglass control in linear and
nonlinear problems. Computer Methods in Applied Mechanics and Engineering, 43:251-276,
1984.
[24] Betsch P, Gruttmann F, Stein E. A 4-node finite shell element for the implementation
of general hyperelastic 3D-elasticity at finite strains. Computer Methods in Applied
Mechanics and Engineering, 130:57–79, 1996.
[25] Bischoff M, Ramm E. Shear deformable shell elements for large strains and rotations.
International Journal for Numerical Methods in Engineering, 40:4427–4449, 1997.
[26] Brank B, Korelc J, Ibrahimbegovic A. Nonlinear shell problem formulation
accounting for through-the-thickness stretching and its finite element implementation.
Computers and Structures, 80:699–717, 2002.
[27] Brendel B, Ramm E. Linear and nonlinear stability analysis of cylindrical shells.
Internat. Conf. on Engineering Application of Finite Element Method, Hovik, Norway, May
9-11, 1979.
[28] Buechter N, Ramm E, Roehl D. Three-dimensional extension of non-linear shell
formulation based on the enhanced assumed strain concept. International Journal for
Numerical Methods in Engineering, 37:2551–2568, 1994.
[29] Cardoso R.P.R, Yoon J.W. One point quadrature shell element with through-thickness
stretch. Computer Methods in Applied Mechanics and Engineering, 194:1161–1199, 2005.
[30] Chapelle D, Bathe K.J. Fundamental considerations for the finite element analysis of
shell structures. Computers and Structures, 66:19–36, 1998.
[31] Chen Y.I, Wu G.Y. A mixed 8-node hexahedral element based on the Hu–Washizu
principle and the field extrapolation technique. Structural Engineering and Mechanics,
17:113–140, 2004.
[32] Cheung Y, Chen W.J. Refined hybrid method for plane iso-parametric element using
an orthogonal approach. Computers and Structures, 42:683–694, 1992.
[33] Cho C, Park H.C, Lee S.W. Stability analysis using a geometrically nonlinear assumed
strain solid shell element model. Finite Elements in Analysis and Design, 29:121–135, 1998.
[34] Dhatt G, Touzot G. Une présentation de la méthode des éléments finis. Edition des
Universités de Laval, Québec, 1986.
[35] Doll S, Schweizerhof K, Hauptmann R, Freischlager C. On volumetric locking of low-
order solid and solid–shell elements for finite elastoviscoplastic deformations and selective
reduced integration. Engineering Computations, 17:874–902, 2000.
[36] Domissy E. Formulation et évaluation d’éléments finis volumiques modifiés pour
l’analyse linéaire et non linéaire des coques. Ph.D. Thesis, UT Compiègne, France, 1997.
[37] Fish J, Belytschko T. Elements with embedded localization zones for large
deformation problems. Computers and Structures, 30:247–256, 1988.
[38] Flanagan D.P, Belytschko T. A uniform strain hexahedron and quadrilateral with
orthogonal hourglass control. International Journal for Numerical Methods in Engineering,
17:679–706, 1981.
[39] Hsiao K.M. Non-linear analysis of general shell structures by flat triangular shell
element. Computers and Structures, 25:665–75, 1987.
[40] Hallquist J.0. Theoretical manual for DYNA3D. UC1D-19401 Lawrence Livermore
National Lab., University of California, 1983.
[41] Hauptmann R, Doll S, Harnau M, Schweizerhof K. Solid–shell elements with linear
and quadratic shape functions at large deformations with nearly incompressible materials.
Computers and Structures, 79:1671–1685, 2001.
[42] Hauptmann R, Schweizerhof K. A systematic development of ‘solid–shell’ element
formulations for linear and non-linear analyses employing only displacement degrees of
freedom. International Journal for Numerical Methods in Engineering, 42:49–69, 1998.
[43] Hauptmann R, Schweizerhof K, Doll S. Extension of the ‘solid–shell’ concept for
application to large elastic and large elastoplastic deformations. International Journal for
Numerical Methods in Engineering, 49:1121–1131, 2000.
[44] Hermann L.R. A bending analysis for plates. Proc. 1st Conf. on Matrix Meth. In
structural Mech. Ohio, page 577–602, October 1965.
[45] Hughes T.J.R. Generalization of selective integration procedures to anisotropic and
nonlinear media. International Journal for Numerical Methods in Engineering, 15:1413–
1418, 1980.
[46] Hughes T.J.R, Cohen M, Haroun M. Reduced and selective integration techniques in
finite element analysis of plates. Nuclear Engineering Design, 46:203–222, 1978.
[47] Hughes T.J.R, Liu W.K. Nonlinear finite element analysis of shells: Part I. Three-
dimensional shells. Computer Methods in Applied Mechanics and Engineering, 26:331–362,
1981.
[48] Kim K.D, Liu G.Z, Han S.C. A resultant 8-node solid–shell element for geometrically
nonlinear analysis. Computational Mechanics, 35:315–331, 2005.
[49] Kim Y.H, Jones R.F, Lee S.W. Study of 20-node solid element. Communication in
Applied Numerical Methods in Engineering, 6:197–205, 1990.
139
Bibliographies
140
[67] Pinsky P.M, Jasti R.V. A mixed finite element formulation for Reissner–Mindlin
plates based on the use of bubble functions. International Journal for Numerical Methods in
Engineering, 28:1677–1702, 1989.
[68] Prathap B.G, Bhashyam G.R. Reduced integration and the shear flexible beam
element. International Journal for Numerical Methods in Engineering, 18:195–210, 1982.
[69] Puso M.A. A highly efficient enhanced assumed strain physically stabilized
hexahedral element. International Journal for Numerical Methods in Engineering, 49:1029–
1064, 2000.
[70] Reese S. A large deformation solid–shell concept based on reduced integration with
hourglass stabilization. International Journal for Numerical Methods in Engineering, 69:
1671–1716, 2007.
[71] Reese S, Wriggers P, Reddy B.D. A new locking-free brick element technique for
large deformation problems in elasticity. Computers and Structures, 75:291–304, 2000.
[72] Sabir A.B, Lock A.C. A curved cylindrical shell finite element. Int. J. Mech. Science,
14:125–135, 1972.
[73] Saleeb A.F, Chang T.Y. An efficient quadrilateral element for plate bending analysis.
International Journal for Numerical Methods in Engineering, 24:1123–1155, 1987.
[74] Saleeb A.F, Chang T.Y, Yingyeunyong S. A mixed formulation of C° - linear
triangular plate–shell element – the role of edge shear constraints. International Journal for
Numerical Methods in Engineering, 26:1101–1128, 1988.
[75] Simo J.C, Hughes T.J.R. On the variational foundations of assumed strain methods. J.
Appl. Mech. ASME, 53:51–54, 1986.
[76] Simo J.C, Armero F. Geometrically non-linear enhanced strain mixed methods and the
method of incompatible modes. International Journal for Numerical Methods in Engineering,
33:1413–1449, 1992.
[77] Simo J.C, Armero F, Taylor R.L. Improved versions of assumed enhanced strain tri-
linear elements for 3D finite deformation problems. Computer Methods in Applied Mechanics
and Engineering, 110:359–386, 1993.
[78] Simo J.C, Rifai M.S. A class of mixed assumed strain methods and the method of
incompatible modes. International Journal for Numerical Methods in Engineering, 29:1595–
1638, 1990.
[79] Stolarski H, Belytschko T, Carpenter N. A simple triangular curved shell element.
Engineering Computation, 1:210–218, 1984.
[80] Sze K.Y, Liu X.H, Lo S.H. Popular benchmark problems for geometric nonlinear
analysis of shells. Finite Elements in Analysis and Design, 40:1551–1569, 2004.
[81] Sze K.Y, Yao L.Q. A hybrid stress ANS solid–shell element and its generalization for
smart structure modelling. Part I-solid–shell element formulation. International Journal for
Numerical Methods in Engineering, 48:545–564, 2000.
[82] Taylor R.L, Beresford P.J, Wilson E.L. A nonconforming element for stress analysis.
International Journal for Numerical Methods in Engineering, 10:1211–1219, 1976.
[83] Timoshenko S.P, Gere J.M. Théorie de la stabilité élastique, deuxième édition
DUNOD, 500, 1966.
141
Bibliographies
[84] Vu-Quoc L, Tan X.G. Optimal solid–shells for non-linear analysis of multilayer
composites. Part I: statics. Computer Methods in Applied Mechanics and Engineering,
192:975–1016, 2003.
[85] Vu-Quoc L, Tan X.G. Optimal solid–shells for non-linear analysis of multilayer
composites. Part II: dynamics. Computer Methods in Applied Mechanics and Engineering,
192:1017–1059, 2003.
[86] Wall W.A, Bischoff M, Ramm E. A deformation dependent stabilization technique,
exemplified by EAS elements at large strains. Computer Methods in Applied Mechanics and
Engineering, 188:859–871, 2000.
[87] Wang X.J, Belytschko T. An efficient flexurally super convergent hexahedral element,
Engineering Computation, 4:281–288, 1987.
[88] Wilson E.L, Taylor R.L, Doherty W.P, Ghaboussi J. Incompatible displacement
models. In Numerical and Computer Methods in Structural Mechanics, (Eds. S.J. Fenves et
al.), 43–57, Academic Press, New York, 1973.
[89] Wriggers P, Reese S. A note on enhanced strain methods for large deformations.
Computer Methods in Applied Mechanics and Engineering, 135:201–209, 1996.
[90] Zhu Y.Y, Cescotto S. Unified and mixed formulation of the 8-node hexahedral
elements by assumed strain method. Computer Methods in Applied Mechanics and
Engineering, 129:177–209, 1996.
[91] Zienkiewicz O.C, Taylor R.L. The finite element method, Fourth Edition, Basic
formulation and linear problems. Vol 1, Editions Mac Graw-Hill Book Company, 1998.
[92] Zienkiewicz O.C, Taylor R.L. The finite element method, Fourth Edition, Solid and
fluid mechanics dynamics and non-linearity. Vol 2, Editions Mac Graw-Hill Book Company,
1998.
[93] Zienkiewicz O.C, Lefebvre D. A robust triangular plate bending element of the
Reissner–Mindlin type. International Journal for Numerical Methods in Engineering,
26:1169–1184, 1988.
[94] Zienkiewicz O.C, Taylor R.L, Too J.M. Reduced integration technique in general
analysis of plates and shells. International Journal for Numerical Methods in Engineering,
3:275–290, 1971.
142
DÉVELOPPEMENT ET VALIDATION D’ÉLÉMENTS FINIS DE TYPE COQUES
VOLUMIQUES SOUS-INTÉGRES STABILISÉS UTILISABLES POUR DES
PROBLÈMES A CINÉMATIQUE ET COMPORTEMENT NON-LINÉAIRES
RÉSUMÉ: Le recours aux logiciels de simulation basés sur la méthode des éléments finis
devenant de plus en plus systématique dans les différents secteurs de l’industrie, l’efficacité et la
précision de ces derniers deviennent des propriétés déterminantes. Dans les situations les plus
courantes, les structures minces nécessitent une analyse précise et efficace, rendue possible par
les éléments coques. En présence de structures dans lesquelles coexistent des parties minces et
des zones plus épaisses, l’utilisation de ces éléments est encore plus cruciale. Ce travail est une
contribution au développement et à la validation d’éléments finis solide-coques. Les déplacements
nodaux sont les seuls degrés de liberté et ils sont munis d’un ensemble de points d’intégration
distribués le long d’une direction préférentielle, désignée comme “l’épaisseur”. Une intégration
réduite dans le plan moyen est utilisée. La loi élastique 3D est modifiée pour s’approcher de la
situation coque et atténuer les verrouillages. Grâce à ses formulations particulières, ces éléments
solide-coques se connectent naturellement aux éléments 3D et présentent une bonne performance
dans des applications de structures minces et pour des problèmes dominés par la flexion. Il s’agit
des trois nouveaux éléments isoparamétriques SHB6, SHB15, et SHB20. L’analyse détaillée d’une
potentielle déficience du rang de la matrice de raideur a révélé que ces derniers ne possèdent pas
de modes à énergie nulle et qu’aucune stabilisation n’est donc nécessaire. Néanmoins, nous
proposons des modifications basées sur la méthode bien connue “Assumed Strain”, pour
l’opérateur gradient discrétisé de l’élément SHB6, dans le but d’améliorer sa vitesse de
convergence. Pour illustrer les capacités de ces éléments, ses performances sont évaluées sur un
ensemble de cas tests en configurations linéaire ou non-linéaire, communément utilisés dans la
littérature pour tester les éléments finis de type coques. En particulier, il est montré que le nouvel
élément SHB6 joue un rôle très utile en tant que complément à l’élément hexaèdre SHB8PS, ce
qui nous permet ainsi de mailler des géométries arbitraires.
ABSTRACT: The appeal to the software of simulation based on the finites elements method
becoming more and more systematic in the various industrial sectors, the efficiency and the
precision of these last ones become determining properties. In the most current situations, the thin
structures require a precise and effective analysis, made possible by shell elements. In the
presence of structures in which coexist thin parts and thicker zones, the use of these elements is
even more crucial. This work is a contribution to the development and to the validation of finites
elements solid-shells. The nodal displacements are the only degrees of freedom and they are
provided of a set of integration’s points distributed along a preferential direction called as “the
thickness”. A reduced integration in the average plan is used. The 3D elastic behaviour is modified
to approach the situation of shells and reduce the locking. Because of its particular formulations,
these elements solid- shells connect naturally in the 3D elements and present a good performance
in applications of thin structures and for problems dominated by the flexion. It is about three new
elements iso-parametric SHB6, SHB15 and SHB20. The detailed analysis of potential deficiency of
the kernel of the stiffness matrix revealed that these last ones don’t have modes associated with a
zero-energy and the stabilization is thus not necessary. Nevertheless, we propose modifications
based on the method “Assumed Strain”, for the operator discretized gradient of the element SHB6,
with the purpose to improve its speed of convergence. To illustrate the capacities of these
elements, its performances are evaluated on varied patch-tests in linear or non-linear
configurations, which are used in the literature to test the finites elements of type shells. In
particular, it is shown that the new element SHB6 plays a very useful role as complement in the
hexahedral element SHB8PS, what so allows us to mesh arbitrary geometries.