r3.07.10 Éléments finis de coques « solides » – COQUE_SOLIDE#
Résumé :
Dans le but de compléter la bibliothèque d’éléments finis de plaque plans [R3.07.03] actuellement disponibles (DKT, DST, Q4G…), on se propose d’introduire un élément fini de coque solide [bib7]. Cette modélisation COQUE_SOLIDE (solid-shell element) permet d’effectuer des calculs de structures minces de formes quelconques.
Avertissement
Cette documentation est provisoire, les éléments ne sont pas suffisamment validés pour être utilisés dans les études AIP et l’élément COQUE_SOLIDE est hors du périmètre de qualification.
Table des matières
Cinématique#
Description des configurations#
On considère la configuration initiale \({\Omega}_{0}\) et la configuration courante \({\Omega}_{t}\) . Le vecteur \(X\) donne la position d’une particule dans la configuration initiale et \(x\) est sa position dans la configuration courante. Avec \(u\) le déplacement de la particule entre les deux configurations:
On considérera deux systèmes de coordonnées. Le système de cordonnées cartésien qui a trois composantes \((x,y,z)\) et dont la base est constituée par les trois vecteurs \({e}_{i}\) et le système de coordonnées paramétriques de référence de l’élément fini qui a aussi trois composantes \(\xi =(\xi ,\eta ,\zeta )\) .
Notations:
Une quantité exprimée dans la configuration initiale sera en majuscule ;
Une quantité exprimée dans la configuration courante sera en minuscule .
On va définir la base dans la configuration initiale par les trois vecteurs suivants:
Le vecteur contravariant correspondant:
De même, dans la configuration courante:
Jacobien de la transformation#
On peut exprimer la matrice jacobienne (de dimension \(3\times 3\) )de la transformation entre les la configuration initiale et l’espace paramétrique:
Et l’inverse de la matrice jacobienne:
Le déterminant de \(J\) permet de transformer les intégrales de l’espace des coordonnées dans la configuration initiale vers l’espace paramétrique de référence des éléments finis. C’est aussi une mesure du volume réel de l’élément fini:
Pour stabiliser efficacement l’élémentet donc pour supprimer les modes sabliersprovenantdel’intégration réduite, il faut pouvoir calculer les différents termes intégraux avec une précision suffisante.
Comme le schéma d’intégration numérique est volontairement très pauvre, la meilleure façon de faire est de travailler avec des intégrandes qui sont des polynômes autant que possible. Pour cette raison, il est intéressant de rechercher une approximation polynomiale des matrices jacobienne \(J\) et jacobienne inverse \(\overline{J}\) , qui permet de représentern’importe quelle forme de l’élément avec une précision suffisante.
La matrice jacobienne \(J\) se décompose en composantes constantes, linéaires et bi-linéaires:
Une façon possible d’obtenir une forme adéquate de la matrice jacobienne inverse pour une bonne efficacité de la stabilisation est de l’évaluer qu’au centre de l’élément (voir [Bib12] et [Bib11]). Dans ce cas, la matrice jacobienne ainsi que son inverse sont constantes à l’intérieur de l’élément. Cependant, cette méthode n’est robuste que pour des maillages très fins.
Pour cette raison, nous approchons le jacobien inverse au moyen d’un développement de Taylor par rapport au centre de l’élément, suivant les travaux de [Bib10]. Cette méthode nous permet d’avoir une forme polynomiale de la matrice jacobienne inverse pour faciliter l’intégration numérique comme recherchée mais aussi que les éléments qui s’écartent de la forme parallélépipédique soient pris en compte de manière plus réaliste.
A partir de la décomposition (), on peut approcher \(\overline{J}\) en ne gardant que le terme constant et les termes linéaires:
Le terme constant est très facile à calculer puisque c’est l’inverse du terme constant de la jacobienne. Pour les autres termes, il faut évaluer les dérivée. On procède en développant le terme \(\overline{J}J\) au centre de l’élément:
Et comme \(\overline{J}J=I\) , sa dérivée est nulle et on a donc pour chaque composante \({\xi}_{i}\) :
Et donc:
Il faut évaluer le terme. Pour cela, on part de () et on dérive. Pour la composante \(\xi\) par exemple:
On n’a gardé que les termes linéaires en \(\xi\) Il reste:
On rappelle que les termes linéaires du développement de \(J\) sont constants sur l’élément (pas de dépendance à la coordonnée paramétrique \(\xi\) ). Et donc:
Finalement, de (), on a l’estimation suivante de la jacobienne inverse:
Attention à bien distinguer \({J}^{0}\) qui est le terme constant (à l’ordre zéro) du développement limité de \(J\) par rapport au terme \({J}_{0}\) qui est la valeur exacte de \(J\) au centre de l’élément (en \(\xi =0\) ).
On aura également besoin de la version du jacobien sur la configuration courante, celle qu’on définit de la manière suivante:
On procède exactement de la même manière pour \(D\) que pour \(J\) . Sa décomposition dans le repère paramétrique:
Tenseur des petites déformations#
La partie linéaire du tenseur des déformations s’écrit en fonction des déplacements:
On va décomposer ce tenseur dans l’espace paramétrique de la manière suivante:
En utilisant une méthode d’intégration réduite, il est nécessaire de modifier l’évaluation de cette matrice pour stabiliser le problème (faire disparaître les modes «sabliers» à énergie nulle), on écrira donc:
La contribution correspondant à l’intégration réduite dans l’épaisseur est notée \({\widehat{B}}^{\mathit{RI}}\) et celle correspondant à la stabilisation est \({\widehat{B}}^{\mathit{STAB}}\) . En décomposant à l’ordre deux sur les composantes paramétriques, l’intégration réduite utilisant un schéma de points alignés suivant l’épaisseur \(\zeta\) , on a:
Et la stabilisation qui correspond aux modes de déformation dans les autres directions:
Tenseur de déformations de Green-Lagrange#
Le tenseur gradient de déformation vaut:
Le tenseur de Green-Lagrange «standard» (non-modifié):
Explicitement, les coordonnées dans le repère paramétrique:
En notation de Voigt, ce tenseur a six composantes:
Pour passer de l’un à l’autre, on utilisera la matrice \(T\) :
Avec la matrice de passage \(T\) :
et \(\overline{J}\) les composantes de l’inverse de la matrice jacobienne:
On utilisera l’approximation () pour \(\overline{J}\) .
La décomposition du tenseur des grandes déformations est de la même forme que celle des petites déformations ():
Et comme dans le cas des petites déformations, en utilisant une méthode d’intégration réduite, il est nécessaire de modifier l’évaluation du tenseur des déformations de Green-Lagrange pour stabiliser le problème (faire disparaître les modes «sabliers» à énergie nulle), on écrira donc:
La contribution correspondant à l’intégration réduite dans l’épaisseur est notée \({\widehat{E}}^{u,\mathit{RI}}\) et celle correspondant à la stabilisation est \({\widehat{E}}^{u,\mathit{STAB}}\) . En décomposant sur les composantes paramétriques, l’intégration réduite utilisant un schéma de points alignés suivant l’épaisseur, on a:
Et la stabilisation qui correspond aux modes de déformation dans les autres directions:
Le tenseur des déformations de Green-Lagrange dans le repère paramétrique s’écrira:
Attention à ne pas confondre les composantes du tenseur des déformations () avec ceux de sa décomposition (-).
Formulations variationnelles#
Équilibre statique linéaire et non-linéaire#
On va exprimer les équations d’équilibre statique linéaire et non-linéaire sous forme variationnelle. On se place dans le cadre de la théorie des grandes déformations. Le tenseur des déformations de Green-Lagrange \({E}^{u}\) vaut:
Avec \(F\) le tenseur gradient de déplacement classique (voir § 2.4 ) et \(g\) le tenseur métrique.
Afin de pouvoir introduire l’ensemble des méthodes de stabilisation et de correction des verrouillages, on va partir d’une formulation variationnelle mixte de Hu-Washizu [Bib2] et [Bib17], écrite sur la configuration initiale et en grandesdéformations:
Avec \(W(E)\) l’énergie de déformation de la structure et \({\pi}_{\text{ext}}(u)\) le travail des efforts extérieurs. L’application de la méthode EAS consiste à reparamétrer le tenseur des déformations de Green-Lagrange tel que:
Ce qui nous donne une nouvelle expression de () qui transforme la dépendance de \(E\) vers \({E}^{\mathit{EAS}}\) :
Il est nécessaire de trouver des approximations correctes de ces trois champs, ce qui peut être compliqué en pratique. Cette formulation à trois champs est donc généralement transformée en formulation à deux champs en imposant la condition d’orthogonalité suivante:
Ce qui permet d’éliminer la contrainte. Au niveau discret, cette condition peut mener à des problèmes de convergence à cause d’un espace des contraintes \({S}^{h}\) trop restreint (voir [Bib3]). Ce problème peut être corrigé en imposant au niveau discret l’inclusion au moins des fonctions constantes par morceaux dans cet espace. Dans notre cas, nous allons procéder autrement pour permettre la prise en compte du pincement. Il y a deux situations:
s’il n’y a pas pincement, c’est-à-dire s’il n’y a pas de pression appliquée sur la surface inférieure ou supérieure de l’élément, on imposera la condition d’orthogonalité ();
s’il y a pincement et qu’on considère une pression sur la face supérieure \({P}_{u}\) et sur la surface inférieure \({P}_{l}\) , on choisit la contrainte de Piola-Kirchhoff de seconde espèce telle que
Avec \(\zeta\) la cordonnée paramétrique dans l’épaisseur.
On a donc:
Ou encore:
La première variation de () permet d’obtenir la forme faible:
Le problème est non-linéaire. On va appliquer la méthode de Newton (voir [R5.03.01]), avec les notations de ce document, \(\delta {\pi}_{int}\) sont les forces internes. Il faut donc exprimer les différentes variations en fonction des inconnues du système.
Pour la déformation \({E}^{u}\) , la variation étant petite, on ne considère que la partie linéaire. Dans le système des coordonnées paramétriques de l’élément, on a:
Et pour la déformation EAS, sur l’élément, en retenant une approximation avec un seul paramètre \(\alpha\) :
On a donc:
Et la force interne:
L’application de la méthode de Newton va faire apparaître quatre matrices tangentes:
Avec les quatre matrices:
Pour \({K}^{\mathit{uu}}\) est la matrice «classique» (que l’on retrouve dans les éléments isoparamétriques standards par exemple), elle comprend elle-même deux parties:
la matrice tangente matérielle \({K}_{m}^{\mathit{uu}}\) ;
la matrice tangente géométrique \({K}_{g}^{\mathit{uu}}\) dans le cas des grandes déformations.
Soit:
Calcul de stabilité#
Ce paragraphe permet d’exposer l’expression du problème aux valeurs propres pour cvalculer la stabiltié d’une structure (flambement).La matrice géométrique intervient dans deux cas:
en non-linéaire quand on utilise des grandes déformations, le processus de linéarisation fait apparaître une contribution d’origine géométrique, voir ();
pour faire des calculs de stabilité (flambement), ce qui correspond dans code_aster à l’option RIGI_GEOM, voir [R7.05.01].
Dans les deux cas, la forme générale de cette matrice est la suivante:
On va exprimer l’opérateur \(G\) tell que:
Comme l’approximation de \(B\) et \(E\) ont utilisé une décomposition en polynômes, on gardera la même logique de décomposition pour \(G\) (et le même schéma d’intégration numérique de Lobatto). De même on aura le même changement de repère entre l’espace paramétrique de l’élément et l’espace réel:
La décomposition choisie est donc de la forme suivante:
Par ailleurs, on aura également une contribution venant de l’intégration réduite et une pour stabiliser:
On écrit les composantes de ce tenseur sous forme de vecteur:
Pour chaque couple de nœud \(i\times j\) , on a les expressions explicites de ces composantes. Pour la membrane:
Calcul de dynamique#
Dans ce paragraphe, on exprime la masse du système pour calculer les problèmes de dynamique linéaire et non-linéaires transitoires ainsi que le calcul modal.
Stabilisation#
Principes#
Dans le calcul des quantités nécessaires, il apparaît des termes intégraux venant de la formulation variationnelle (voir § 3 ). La matrice de rigidité de la structure s’écrira de la manière suivante pour un élément fini \(e\) donné:
Avec \(C\) la matrice d’élasticité de Hooke et \(B\) la matrice des déformations linéarisées. On doit donc utiliser un schéma de quadrature numérique pour évaluer l’intégrale. Ce schéma contiendra \({n}_{\mathit{pt}}\) points de coordonnées \({\xi}_{p=1,{n}_{\mathit{pt}}}\) , chacun d’un poids \({\omega}_{p=1,{n}_{\mathit{pt}}}\) tels que:
\(J\) est le jacobien de transformationde l’élément entre l’espace paramétrique et l’espace cartésien (voir § 2.2 ). Le choix de ce schéma, sa précision en particulier, va dépendre de la nature de l’intégrande. Pour l’élément COQUE_SOLIDE, on va utiliser un schéma de type Lobatto à 5 points.
Ce schéma est trop «pauvre» et il va sous-estimer la rigidité de l’élément avec pour objectif de contrôler les verrouillages. Mais en procédant ainsi, on fait également apparaître des modes de déformations qu’on appelle «modes sabliers». L’élément risque de changer de forme sans dépenser d’énergie, ce qui ne correspond pas à la physique. Pour minimiser ce problème, on va stabiliser l’élément par une technologie adéquate.
Le but de la procédure de stabilisation est de corriger le défaut de rang de la matrice de rigidité provenant du schéma d’intégration réduit adopté. Le lecteur est renvoyé aux importants travaux sur ce sujet de [bib18], [bib19], [bib13], [bib14], [bib8], [bib9], [bib10]. Dans certaines formulations de coque-solides, la matrice jacobienne et son inverse ne sont évaluées qu’au centre de l’élément, voir par exemple [bib21], [bib22], [bib23]. Une telle approximation suppose que l’élément réel peut être représenté par son parallélépipède équivalent. Et c’est assez précis pour les maillages fins et peu déformés.
Cependant pour les maillages qui sont fortement déformés à l’état initial, cette hypothèse montre un manque de précision [bib10]. Afin de prendre en considération la forme réelle de l’élément et de stabiliser correctement la matrice de rigidité, une décomposition polynomiale de la matrice jacobienne inverse a étéconstruite explicitement.
Stabilisation des déformations#
Par ailleurs, pour le passage entre des quantités comme la déformation entre les coordonnées paramétriques et les coordonnées réelles, donné par la matrice \(T\) , il convient d’utiliser une décomposition polynomiale cohérente avec celle choisie pour le jacobien en particulier. Il a été montré dans [bib10] qu’une très bonne précision est atteinte en utilisant simplement les termes constants et linéaires de \(T\) . Par conséquent, les termes bilinéaires sont ignorés:
On va utiliser la même forme d’approximation sur \({E}^{u}\) que sur \({\widehat{E}}^{u}\) :
Que l’on va diviser en deux composantes:
Avec les expressions explicites suivantes:
On utilisera la même décomposition pour le tenseur des déformations linéarisées et donc pour la matrice gradient:
Avec:
\({E}^{u,\mathit{STAB}}\) et \({B}^{\mathit{STAB}}\) représentent les termes annulés dans l’intégration réduite. Le processus de stabilisation consiste à restituer analytiquement ces termes dans la matrice de rigidité.
De plus, pour éliminer le verrouillage volumétrique, l’approche B-bar [bib24] est adoptée. Par conséquent, les modes en sablier des opérateurs de déformation et de déplacement de déformation stabilisée est divisée en composants volumétriques et déviatoriques:
Il n’y a a pas de terme constant dans l’expression de \({B}^{\mathit{STAB}}\) , donc:
Seule la partie déviatorique est conservée. On peut démontrer que la contribution \({B}^{\xi \eta }\) (et \({E}^{u,\xi \eta }\) ) est nulle dans l’écriture de la stabilisation.
Stabilisation des contraintes#
On reprend dans l’expression de l’énergie de déformation de la structure la séparation entre les contributions venues de l’intégration réduite et de la stabilisation:
Le tenseur de contraintes de seconde espèce de Piola-Kirchhoff étant la dérivée de cette énergie par rapport aux déformations de Green-Lagrange, on a:
On retrouve la séparation entre les deux contributions (intégration réduite et stabilisation) pour le tenseur des contraintes:
Le tenseur \({S}^{\mathit{RI}}\) vient directement de l’intégration de la loi de comportement tandis que le tenseur \({S}^{\mathit{STAB}}\) correspond aux termes de contraintes pour la stabilisation.
Le point clé de la stabilisation est de trouver un moyen optimal d’évaluer cette contrainte de stabilisation sans l’intégrer afin que le temps de calcul soit minimisé. Or, comme expliqué précédemment, pour éviter le verrouillage volumétrique, seule la partie déviatorique du tenseur de stabilisation de Piola Kirchhoff est conservée. Ce faisant, la stabilisation de la contrainte est donnée comme:
\({C}^{\mathit{STAB}}\) étant la partie déviatorique du comportement de Saint-Venant Kirchhoff:
Avec le module de cisaillement \({\mu}^{\mathit{STAB}}\) , qui, pour un matériau élastique isotrope vaut:
Pour les matériaux non-élastiques, l’utilisation du module de cisaillement élastique conduit à une surestimation de la contrainte de stabilisation. Pour surmonter ce problème, le module sécant tel que défini dans [bib20] est utilisé:
Avec:
Cette manière d’évaluer le paramètre de stabilisation est très pratiques car il n’y a pas besoin de linéarisation cohérente de la matrice ce stabilisation, ce qui est très intéressant en termes de temps de calcul. De plus, la définition proposée de \({C}^{\mathit{STAB}}\) assure une stabilisation efficace et robuste.
Stabilisation des opérateurs discrets#
En reprenant les définitions précédentes, on calcule la contribution \({\mu}^{\mathit{STAB}}\) de la manière suivante. A chaque point de Gauss \(g\) on calcule \({{\pi}_{S}}^{g}\) et \({{\pi}_{E}}^{g}\) puis on considère les deux situations suivantes:
Si \(|{{\pi}_{E}}^{g}|<0,001\) ) alors \({\mu}^{\mathit{STAB},g}={K}_{\mathit{vol}}^{g}\) . C’est à dire que le coefficient de stabilisation au point de Guass courant \(g\) est égal à la partie volumique de la matrice de comportement matériau si la déformation déviatorique est inférieure à 0.1%;
Sinon \({\mu}^{\mathit{STAB},g}=\frac{1}{2}\sqrt{\frac{{\pi}_{S}^{g}}{{\pi}_{E}^{g}}}\) ;
Puis le coefficient de stabilisation est sommé sur chaque point de Gauss:
La stabilisation des contraintes dans le cas des petites déformations, consiste à calculer les différentes contributions de la manière suivante (on rappelle que le terme \(\xi \eta\) ne contribue pas):
Et en grandes déformations:
Le processus de linéarisation de Newton de l’équation d’équilibre produit deux opérateurs: la matrice tangente consistante et le vecteur des forces internes. Les deux vont donc devoir intégrer la stabilisation puisqu’ils dépendent des contraintes. Pour la matrice tangente, il y a deux contributions:
la première concerne la stabilisation liée aux contraintes issues de l’intégration de la loi de comportement (matrice matérielle), c’est-à-dire la matrice \({K}_{m}^{\mathit{uu}}\) ;
la seconde correspond au modèle de grande déformation, c’est la partiegéométrique qui vient de la non-linéarité cinématique, soit la matrice \({K}_{g}^{\mathit{uu}}\) .
Pour la stabilisation de la matrice tangente matérielle, on a quatre contributions:
Ce qui nous donne la matrice globale de stabilisation matérielle avec le jacobien de changement de coordonnées ():
Pour la partie géométrique, comme la contribution est le produit d’une matrice (géométrique) avec les contraintes, il faut aussi le faire par un produit pour la stabilisation, pour chaque nœud \(i\) et \(j\) :
Élément hexaédrique à neuf nœuds (SB9)#
L’élément SB9 a neuf nœuds, un point d’intégration dans le plan et plusieurs points d’intégration le long de la direction de l’épaisseur de l’élément.
Par rapport aux autres briques «à coque solide» à huit nœuds, la présence d’un nœud supplémentaire a deux objectifs principaux:
Obtenir d’abord une composante de déformation normale linéaire qui, avec un comportement de déformation-contrainte constitutif tridimensionnel complet, permet d’obtenir des résultats similaires dans les cas de flexion que ceux obtenus avec la contrainte plane habituelle. Pour cela, le neuvième nœud joue le rôle d’un paramètre supplémentaire essentiel pour obtenir une interpolation quadratique du déplacement dans le sens de l’épaisseur.
De plus, ce DDL a une signification physique et, par exemple, une force équivalente à une pression normale peut être prescrite pour améliorer la contrainte normale lorsque la structure de la coque est moyennement épaisse.
Géométrie de la brique à huit nœuds#
La position initiale \(X\) d’une particule est interpolée par des fonctions de formes définies dans l’espace paramétrique de l’élément:
De même pour la position courante \(x\) :
avec \({X}^{i}\) et \({x}^{i}\) les valeurs nodales de la position dans l’espace paramétrique, pour les huit nœuds:
Nœud |
\(\xi\) |
\(\eta\) |
\(\zeta\) |
1 |
-1 |
-1 |
-1 |
2 |
+1 |
-1 |
-1 |
3 |
+1 |
+1 |
-1 |
4 |
-1 |
+1 |
-1 |
5 |
-1 |
-1 |
+1 |
6 |
+1 |
-1 |
+1 |
7 |
+1 |
+1 |
+1 |
8 |
-1 |
+1 |
+1 |
Les fonctions de forme \({N}_{i}\) pour chaque nœud \(i\) sont le produit cartésien des fonctions linéaires suivant les trois directions:
Ce qu’on résume par la combinaison linéaire suivante:
Avec les vecteursisoparamétriques suivants:
Des vecteurs auxiliaires sont définis à partir des combinaisons précédentes des vecteurs isoparamétriques:
Les dérivées des fonctions de forme (par rapport aux coordonnées paramétriques) s’écrivent comme:
Interpolation des déplacements#
L’élément est isoparamétrique et n’utilise que les déplacements aux nœuds (pas les rotations). Ces déplacements s’écriront donc avec les mêmes fonctions de forme que la géométrie:
Et les dérivées:
Calcul des jacobiens#
De (), la matrice jacobienne \(J\) se décompose en composantes constantes, linéaires et bi-linéaires:
Chaque composante est une matrice \(3\times 3\) et vaut:
En pratique, il y a 12 quantités non-nulles différentes (des vecteurs lignes), que l’on notera ainsi:
La stratégie utilisée dans le cas de \(J\) sera identique pour \(D\) . En particulier, de manière similaire à (-), on a:
Avec une dépendance au déplacement \(U\) au lieu de la position \(X\) . De même:
Calcul des petites déformations#
On rappelle la décomposition du tenseur des petites déformations dans l’espace paramétrique ():
Dans le plan (membrane), les composantes pour chaque nœud \(i\) sont les suivantes:
Avec les différents vecteurs lignes composantes \({J}_{k}^{{\xi}_{i}}\) données par (-) et \({a}_{m}^{i}\) la i-ème composante (scalaire) du nœud \(i\) extrait des vecteurs \({a}_{m=\mathrm{1,2,3}}\) définis par () et de même pour \({h}_{n}^{i}\) la i-ème composante (scalaire) du nœud \(i\) extrait des vecteurs \({h}_{n=\mathrm{1,2,3,4}}\) définis par ().
Par ailleurs, la matrice globale au niveau de l’élément sera obtenu par assemblage sur tous les nœuds de la maille sur le modèle suivant:
Une matrice \({\widehat{B}}_{m}^{\cdot}\) sur l’élément est donc de dimension \(3\times 24\) .
Remarques:
on n’a pas écrit l’expression explicite de la composante \({\widehat{B}}_{m,i}^{\xi \eta }\) car la procédure de stabilisation supprime ce terme;
la partie hors-membrane sera vue lorsque l’on appliquera la méthode ANS.
Calcul des grandes déformations#
On rappelle la décomposition du tenseur des petites déformations dans l’espace paramétrique de ():
Les termes de décomposition en membrane valent:
Méthode des déformations postulées (ANS)#
Comme nous l’avons vu dans le § 1.1.6 , une façon de réduire le blocage par cisaillement est d’utiliser la méthode ANS. La déformation de cisaillement transverse est évaluée dans les quatre points des bords médians de la surface médiane de la coque solide. Pour interpoler à travers le volume, certains auteurs proposent de multiplier la déformation de cisaillement dans le plan médian avec la fonction bien connue de Reissner. Ce type d’interpolation assure le respect de la condition statique de contrainte de cisaillement sur les faces supérieure et inférieure de l’élément.
Cependant dans les problèmes non-linéaires (matériau et grandes déformations), cette interpolation devient fausse et elle ne garantit pas la satisfaction de la condition statique sur les faces. De plus, une telle interpolation ne permet pas à l’élément de gérer naturellement les problèmes de contact avec frottement puisque la déformation en cisaillement est supposée nulle sur les faces supérieure et inférieure de l’élément. Il a également été montré qu’une telle interpolation conduit à une matrice de rigidité mal conditionnée. Pour éviter ces inconvénients nous avons appliqué la technique ANS suivant les travaux de [Bib13], [Bib8], [Bib9] et [Bib11].
Dans l’espace paramétrique, chacune des déformations de cisaillement \({\widehat{E}}^{u,\xi \zeta }\) et \({\widehat{\epsilon}}^{\eta \zeta }\) est interpolée en utilisant la déformation de cisaillement équivalente à partir de quatre points de liaison différents (figure ).
Méthode ANS en petites déformations#
On rappelle que la partie linéaire du tenseur des déformations s’écrit en fonction des déplacements:
Sur la matrice \(\widehat{B}\) , on écrit donc une interpolation spécifique pour les composantes de cisaillement et normales. Pour \({\widehat{\epsilon}}^{\xi \zeta }\) , la matrice \({\widehat{B}}^{\xi \zeta }\) utilise une interpolation sur les nœuds J, K, L et M:
Avec les fonctions de forme suivantes:
Pour \({\widehat{\epsilon}}^{\eta \zeta }\) , la matrice \({\widehat{B}}^{\eta \zeta }\) utilise une interpolation sur les nœuds E, F, G et H:
Avec les fonctions de forme suivantes:
Les coordonnées des nœuds dans l’espace paramétrique \(\lbrace {\xi}_{i},{\eta}_{i},{\zeta}_{i}\rbrace\) sont données dans le § 5.1 . Afin d’éliminer le verrouillagedans l’épaisseur pour le cas courbe, nous postulerons la déformation de pincement suggérée par [Bib15], [Bib16], [Bib10] telle que:
Avec les fonctions de forme suivantes:
On note de manière compacte les sommes suivantes:
Où \({\eta}_{A}\) désigne par exemple la cordonnée paramétrique \(\eta\) pour le nœud A et de manière équivalente pour les groupes de nœuds E,F,G,H et J,K,L,M. Les composantes pour chaque nœud \(i\) sont les suivantes pour la partie relative à l’intégration réduite:
Et pour la partie correspondant à la stabilisation:
Méthode ANS en grandes déformations#
Matrice géométrique#
La matrice géométrique intervient dans deux cas:
en non-linéaire quand on utilise des grandes déformations, le processus de linéarisation fait apparaître une contribution d’origine géométrique, voir ();
pour faire des calculs de stabilité (flambement), ce qui correspond dans code_aster à l’option RIGI_GEOM, voir [R7.05.01].
Dans les deux cas, la forme générale de cette matrice est la suivante:
On va exprimer l’opérateur \(G\) tell que:
Comme l’approximation de \(B\) et \(E\) ont utilisé une décomposition en polynômes, on gardera la même logique de décomposition pour \(G\) (et le même schéma d’intégration numérique de Lobatto). De même on aura le même changement de repère entre l’espace paramétrique de l’élément et l’espace réel:
La décomposition choisie est donc de la forme suivante:
Par ailleurs, on aura également une contribution venant de l’intégration réduite et une pour stabiliser:
On écrit les composantes de ce tenseur sous forme de vecteur:
Pour chaque couple de nœud \(i\times j\) , on a les expressions explicites de ces composantes. Pour la membrane:
Les cotnributions de cisaillement utilisent la méthode ANS.
antes paramétriques, l’intégration réduite utilisant un schéma de points alignés suivant l’épaisseur, on a:
Et la stabilisation qui correspond aux modes de déformation dans les autres directions:
Le tenseur des déformations de Green-Lagrange dans le repère p
Références#
McNeal R. – Finite elements: their design and performance – Marcel Dekker Eds – 1994.
Simo J. C., Fox D. D., Rifai M. S. – On a stress resultant geometrically exact shell model. Part III: Computational aspects of the nonlinear theory . Computer Methods in Applied Mechanics and Engineering, vol. 79, N°1, 1990, p. 21–70.
Simo J. C., Rifai M. – A class of mixed assumed strain methods and the method of incompatible modes . International journal for numerical methods in engineering, vol. 29, N°8, 1990, p. 1595–1638.
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, vol. 33, N°7, 1992, p. 1413–1449.
Simo J., Armero F., Taylor R. – Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems . Computer methods in applied mechanics and engineering, vol. 110, N°3-4, 1993, p. 359–386.
Adam L., Ponthot J.-P. – Thermomechanical modeling of metals at finite strains: first and mixed order finite elements . International Journal of Solids and Structures, vol. 42, N°21-22, 2005, p. 5615–5655.
Dia M., Hamila N., Abbas M., Gravouil A. – A nine nodes solid-shell finite element with enhanced pinching stress . Computational Mechanics, vol. 65, 2020, p. 1377–1395.
Alves de Sousa R. J., Cardoso R. P., Fontes Valente R. A., Yoon J.-W., Grácio 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, vol. 62, N°7, 2005, p. 952–977.
Alves de Sousa R. J., Cardoso R. P., Fontes Valente R. A., Yoon J.-W., Grácio 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, vol. 67, N°2, 2006, p. 160–188.
Schwarze M., Reese S. – A reduced integration solid-shell finite element based on the EAS and the ANS concept—Geometrically linear problems . International Journal for Numerical Methods in Engineering, vol. 80, N°10, 2009, p. 1322–1355.
Reese S. – A large deformation solid-shell concept based on reduced integration with hourglass stabilization . International Journal for Numerical Methods in Engineering, vol. 69, N°8, 2007, p.1671–1716.
Legay A., Combescure A. – Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS . International Journal for Numerical Methods in Engineering, vol. 57, N°9, 2003, p. 1299–1322.
Cardoso R. P., Yoon J. W., Mahardika M., Choudhry S., Alves de Sousa – Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements . International Journal for Numerical Methods in Engineering, vol. 75, N°2, 2008, p. 156–187.
Cardoso R. P., Yoon J. W. One point quadrature shell element with through-thickness stretch . Computer Methods in Applied Mechanics and Engineering, vol. 194, N°9, 2005, p. 1161–1199.
Bischoff M., Ramm E. – Shear deformable shell elements for large strains and rotations . International Journal for Numerical Methods in Engineering, vol. 40, N°23, 1997, p. 4427–4449.
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, vol. 130, N°1-2, 1996, p. 57–79.
Militello C., Felippa C. A. – A variational justification of the assumed natural strain formulation of finite elements — I. Variational principles . Computers & Structures, vol. 34, N°3, 1990, p. 431–438.
Liu W. K., Guo Y., Tang S., Belytschko T. – A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis . Computer Methods in Applied Mechanics and Engineering, vol. 154, N°1-2, 1998, p. 69–132.
Belytschko T., Tsay C.-S. – A stabilization procedure for the quadrilateral plate element with one-point quadrature . International Journal for Numerical Methods in Engineering, vol. 19, N°3, 1983, p. 405–419.
Belytschko T., Bindeman L. P. – Assumed strain stabilization of the eight node hexahedral element . Computer Methods in Applied Mechanics and Engineering, vol. 105, N°2, 1993, p. 225–260.
Bassa B., Sabourin F., Brunet M. – A new nine-node solid-shell finite element using complete 3D constitutive laws . International Journal for Numerical Methods in Engineering, vol. 92, N°7, 2012, p. 589–636.
Abed-Meraim F., Combescure A. – SHB8PS—-a new adaptative, assumed-strain continuum mechanics shell element for impact analysis . Computers & Structures, vol. 80, N°9, 2002, p. 791–803.
Abed-Meraim F., Combescure A. – An improved assumed strain solid–shell element formulation with physical stabilization for geometric non-linear applications and elastic–plastic stability analysis . International Journal for Numerical Methods in Engineering, vol. 80, N°13, 2009, p. 1640–1686.
Hughes T. J. – Generalization of selective integration procedures to anisotropic and nonlinear media . International Journal for Numerical Methods in Engineering, vol. 15, N°9, 1980, p. 1413–1418.