# Stabilisation et projection de type “ Assumed Strain ” de l’élément de coque massif “ SHB8PS ” pour éliminer les problèmes de verrouillage

Farid Abed-Meraim, Alain Combescure

► To cite this version:

Farid Abed-Meraim, Alain Combescure. Stabilisation et projection de type "Assumed Strain" de l'élément de coque massif "SHB8PS" pour éliminer les problèmes de verrouillage. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. pp.583-588. hal-01812918

HAL Id: hal-01812918

https://hal.science/hal-01812918v1

Submitted on 11 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Stabilisation et projection de type « Assumed Strain » de l'élément de coque massif « SHB8PS » pour éliminer les problèmes de verrouillage

Farid Abed-Meraim* — Alain Combescure**

* Laboratoire de Physique et Mécanique des Matériaux, ENSAM CER de Metz
UMR CNRS 7554
4 rue Augustin Fresnel, 57078 Metz
farid.abed-meraim@metz.ensam.fr

** Laboratoire de Mécanique des Contacts et des Solides, INSA de Lyon
UMR CNRS 5514
18-20 rue des Sciences, 69621 Villeurbanne
alain.combescure@insa-lyon.fr

RÉSUMÉ. Dans cette étude, la formulation de l'élément « SHB8PS » est revisitée dans le but d'éliminer certains blocages persistants en membrane ou cisaillement transverse. Rappelons que cet élément est de type coque épaisse obtenue à partir d'une formulation purement tridimensionnelle. Il possède donc huit nœuds et cinq points d'intégration répartis selon la direction de l'épaisseur. Ainsi, il peut être utilisé pour modéliser des structures minces tout en prenant correctement en compte les différents phénomènes à travers l'épaisseur. Afin d'améliorer ses performances de calcul et d'éviter certains blocages, l'intégration réduite a été employée. Les modes de hourglass générés par la sous-intégration sont efficacement stabilisés et les modes de blocages persistants sont éliminés par une technique de projection pouvant se mettre sous le formalisme de la « Méthode de Déformation Postulée ».

ABSTRACT. In this work, the formulation of the finite element "SHB8PS" is revisited in order to eliminate some persistent membrane or transverse shear locking. Let us recall that the "SHB8PS" is a continuum mechanics shell element obtained from a purely three-dimensional formulation. It therefore possesses eight nodes and five integration points distributed according to the thickness direction. Thus, it can be used to model thin structures while correctly taking into account the various phenomena throughout the thickness. In order to improve its calculation performances and to avoid some blockages, reduced integration was employed. The hourglass modes generated by under-integration are efficiently stabilized and the persistent locking modes are eliminated using a projection technique entering under the "Assumed Strain Method" formalism.

MOTS-CLÉS: Coque massive SHB8PS, Hourglass, Blocage, Déformation Postulée, Projection.

KEYWORDS: SHB8PS solid shell, Hourglass, Locking, Assumed Strain Method, Projection.
# 1. Introduction

Plusieurs travaux ont été consacrés ces dernières années au développement d'éléments finis solides capables de modéliser des structures minces (Lemosse 2000), (Sze et al., 2000), (Abed-Meraim et al., 2001, 2002), (Legay et al., 2003), (Chen et al., 2004). L'intérêt de disposer d'éléments solides/coques est motivé par de nombreux problèmes industriels. Un élément de cette famille appelé « SHB8PS » a récemment été développé à partir d'une formulation purement tridimensionnelle (Abed-Meraim et al., 2001, 2002), (Legay et al., 2003). Les avantages de cet élément coque massif sont multiples, on peut citer :

* La capacité à modéliser le comportement des structures tridimensionnelles minces avec très peu d'éléments dans l'épaisseur tout en prenant correctement en compte les différents phénomènes à travers l'épaisseur (flexion, élasto-plasticité, etc.).

* Le maillage simplifié de géométries complexes où éléments coques et éléments solides doivent cohabiter (renforts ou brides par exemple) sans que se posent les problèmes classiques de raccordements de maillages réalisés avec différents types d'éléments.

* Les performances obtenues en termes de temps de calcul sont très intéressantes grâce, d'une part, à l'élancement important toléré (peu d'éléments de maillage) et d'autre part, à l'élimination des blocages en membrane et en cisaillement par des techniques appropriées.

* La formulation associée est relativement simple et attractive (géométrie hexaédrique, seulement trois degrés de liberté pour chacun des huit nœuds) ce qui évite le recours aux formulations complexes de coques faisant appel à des métriques spécifiques.

Dans cette étude, la formulation du SHB8PS est revisitée avec de nouvelles projections pour éliminer certains blocages persistants en membrane ou en cisaillement. Rappelons qu'ayant une géométrie volumique (huit nœuds), on cherche à donner au SHB8PS un comportement de type coque en modifiant la loi de comportement tridimensionnelle pour se rapprocher des hypothèses de contraintes planes, et en utilisant cinq points d'intégration disposés à travers une direction privilégiée que l'on appelle épaisseur.

Afin d'améliorer les performances de l'élément en temps de calcul et pour éviter certains blocages en membrane ou en cisaillement, l'intégration réduite a été employée. Pour éliminer les modes à énergie nulle dus à cette sous-intégration, une technique efficace de stabilisation et de contrôle des modes de *Hourglass* a été utilisée en suivant la démarche de (Belytschko et al., 1993). On montrera d'abord que ces modes cinématiques sont les vecteurs du noyau de la raideur autres que les modes rigides. En exhibant une base de l'espace vectoriel des déplacements discrétisés, on arrive à expliciter ces modes de *Hourglass* puis à les stabiliser efficacement.
Pour éliminer les blocages numériques en cisaillement transverse ou en membrane, l'opérateur gradient discrétisé est projeté sur un sous-espace approprié. Cette technique de projection peut se mettre sous le formalisme de la méthode de déformation postulée « Assumed Strain Method ». On montre également que cette démarche peut être justifiée dans le cadre du principe variationnel mixte de Hu-Washizu. Cette étape d'élimination du verrouillage est longue car elle consiste à manipuler chacun des termes de l'opérateur gradient discrétisé séparément pour éliminer ceux responsables du blocage en cisaillement transverse ou en membrane.

Rappelons que la première version du SHB8PS a été formulée en explicite et implantée dans le code de dynamique rapide EUROPLEXUS du CEA pour simuler les problèmes d'impacts (Abed-Meraim et al., 2001, 2002). Cette version explicite a également été utilisée pour simuler les phénomènes d'ingestion de volatiles par des turboréacteurs d'avions ainsi que d'autres situations accidentelles posées à la SNECMA. Une version implicite a ensuite été formulée et introduite dans le code quasi statique INCA pour faire des analyses de stabilité élasto-plastique (Legay et al., 2003). Plus récemment, et au vu du bon retour d'expérience, cette version implicite a été implantée dans le code ASTER d'EDF.

Notons que les formulations précédentes : explicite et implicite du SHB8PS, bien que bénéficiant d'une projection pour éliminer certains blocages, présentaient encore une convergence relativement lente dans le cas test de la sphère pincée. Cette persistance de blocages en cisaillement transverse ou en membrane dans certaines situations nous a amené à retravailler sur la formulation de cet élément. Nous-nous sommes consacré plus particulièrement aux techniques de projection afin d'aller beaucoup plus loin dans l'élimination des phénomènes de verrouillage. La nouvelle version que nous obtenons pour le SHB8PS est présentée dans cette étude et ses performances en termes de rapidité de convergence sont clairement démontrées. Sur un grand nombre de cas tests bien connus, cette nouvelle formulation du SHB8PS montre son efficacité et permet d'atteindre la solution de référence sans aucun phénomène de blocage.

## 2. Formulation de l'élément SHB8PS

### 2.1. Cinématique et interpolation

Le SHB8PS est un élément hexaédrique à huit nœuds, isoparamétrique et à interpolation linéaire. Il possède cinq points d'intégration disposés le long de la direction $\zeta$ des coordonnées locales. La figure 1 montre la géométrie de référence de l'élément ainsi que la répartition des cinq points d'intégration. Les coordonnées $x_i$, $i = 1, 2, 3$ d'un point de l'élément sont reliées aux coordonnées nodales $x_{ij}$ via les fonctions de formes classiques linéaires isoparamétriques $N_i$ ($i = 1, ..., 8$) par les formules :
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Géométrie de référence de l'élément et points d'intégration.</figcaption></figure>

$$x_i = x_{ii} N_i(\xi, \eta, \zeta) = \sum_{j=1}^{8} x_{ij} N_i(\xi, \eta, \zeta) \quad [1]$$

Dans la suite, et sauf mention contraire, on adoptera la convention de sommation pour les indices répétés. Les indices en minuscule *i* varient de un à trois et représentent les directions des coordonnées spatiales. Ceux en majuscule *I* varient de un à huit et correspondent aux nœuds de l'élément.

Avec cette convention, l'interpolation du champ de déplacement $u_i$ dans l'élément en termes de déplacements nodaux $u_{ij}$ est similaire :

$$u_i = u_{ii} N_i(\xi, \eta, \zeta) \quad [2]$$

## 2.2. Opérateur gradient discrétisé

L'interpolation du champ de déplacement (équation [2]) va nous permettre d'écrire les relations reliant le champ de déformation aux déplacements nodaux. La partie linéaire du tenseur de déformation s'écrit :

$$\varepsilon_{ij} = \frac{1}{2} (u_{i,j} + u_{j,i}) = \frac{1}{2} (u_{ii} N_{i,i} + u_{jj} N_{j,j}) \quad [3]$$

Pour poursuivre la formulation, on se donne les fonctions de forme tri-linéaires :

$$\begin{gathered} N_i(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi_i \xi)(1 + \eta_i \eta)(1 + \zeta_i \zeta) \\ \xi, \eta, \zeta \in [-1, 1], \quad I = 1, \dots, 8 \end{gathered} \quad [4]$$
Ces fonctions de forme transforment un cube unitaire dans l'espace $(\xi, \eta, \zeta)$ en un hexaèdre quelconque dans l'espace $(x_1, x_2, x_3)$. Combinant les équations [1], [2] et [4], on arrive à développer le champ de déplacement comme un terme constant, des termes linéaires en $x_i$ et des termes faisant intervenir les fonctions $h_\alpha$ :

$$ egin{cases} u_i = a_{0i} + a_{1i}x + a_{2i}y + a_{3i}z + c_{1i}h_1 + c_{2i}h_2 + c_{3i}h_3 + c_{4i}h_4 \\ i = 1, 2, 3 \\ h_1 = \eta\zeta, \quad h_2 = \zeta\xi, \quad h_3 = \xi\eta, \quad h_4 = \xi\eta\zeta \end{cases} \quad [5] $$

En évaluant l'équation ci-dessus aux nœuds de l'élément, on arrive aux trois systèmes de huit équations suivants :

$$ \begin{cases} \underline{d}_i = a_{0i}s + a_{1i}\underline{x}_1 + a_{2i}\underline{x}_2 + a_{3i}\underline{x}_3 + c_{1i}\underline{h}_1 + c_{2i}\underline{h}_2 + c_{3i}\underline{h}_3 + c_{4i}\underline{h}_4 \\ i = 1, 2, 3 \end{cases} \quad [6] $$

Dans l'équation ci-dessus, les vecteurs $\underline{d}_i$ et $\underline{x}_i$ désignent, respectivement, les déplacements et les coordonnées aux nœuds et sont donnés par :

$$ \begin{cases} \underline{d}_i^T = (u_{i1}, u_{i2}, u_{i3}, \dots, u_{i8}) \\ \underline{x}_i^T = (x_{i1}, x_{i2}, x_{i3}, \dots, x_{i8}) \end{cases} \quad [7] $$

Les vecteurs $\underline{s}$ et $\underline{h}_\alpha$ ($\alpha = 1, \dots, 4$) sont données quant à eux par :

$$ \begin{cases} \underline{s}^T = (1, 1, 1, 1, 1, 1, 1, 1) \\ \underline{h}_1^T = (1, 1, -1, -1, -1, -1, 1, 1) \\ \underline{h}_2^T = (1, -1, -1, 1, -1, 1, 1, -1) \\ \underline{h}_3^T = (1, -1, 1, -1, 1, -1, 1, -1) \\ \underline{h}_4^T = (-1, 1, -1, 1, 1, -1, 1, -1) \end{cases} \quad [8] $$
Les constantes inconnues $a_{ji}$ et $c_{\alpha i}$ de l'équation [5] sont obtenues par l'introduction des vecteurs de (Hallquist, 1983), ($\underline{b}_i$, $i=1,...,3$) définis par :

$$ \underline{b}_i^T = \underline{N}_{,i}(0,0,0) \quad i=1,2,3 \quad \text{Hallquist Form} \quad [9] $$

Les expressions explicites de ces dérivées des fonctions de forme évaluées à l'origine du référentiel $(\xi, \eta, \zeta)$ ont été données par (Belytschko et al. 1984). On démontre d'abord les relations d'orthogonalité suivantes impliquant ces vecteurs $\underline{b}_i$ et les vecteurs : $\underline{x}^T = (x_{i1}, x_{i2}, ..., x_{i8})$, $\underline{d}^T = (u_{i1}, u_{i2}, ..., u_{i8})$, $\underline{s}$, $\underline{h}_1$, $\underline{h}_2$, $\underline{h}_3$ et $\underline{h}_4$ :

$$ \left\{ \begin{array}{lll} \underline{b}^T \cdot \underline{h}_\alpha = 0 & , & \underline{b}^T \cdot \underline{s} = 0 & , \underline{b}^T \cdot \underline{x}_j = \delta_{ij} \\ \\ \underline{h}_\alpha^T \cdot \underline{s} = 0 & , & \underline{h}_\alpha^T \cdot \underline{h}_\beta = 8\delta_{\alpha\beta} & \\ \\ i, j = 1, \dots, 3 & \alpha, \beta = 1, \dots, 4 & & \end{array} \right. \quad [10] $$

Ces relations d'orthogonalité permettent de calculer les constantes $a_{ji}$ et $c_{\alpha i}$ en multipliant scalairement l'équation [6] par $\underline{b}_j^T$ et $\underline{h}_\alpha^T$, respectivement. On obtient :

$$ \left\{ \begin{array}{ll} a_{ji} = \underline{b}_j^T \cdot \underline{d}_i & , \quad c_{\alpha i} = \underline{\gamma}_\alpha^T \cdot \underline{d}_i \\ \\ \text{avec :} \quad \underline{\gamma}_\alpha = \frac{1}{8} \left[ \underline{h}_\alpha - \sum_{j=1}^3 (\underline{h}_\alpha^T \cdot \underline{x}_j) \underline{b}_j \right] & \end{array} \right. \quad [11] $$

Le champ de déplacement se met finalement sous la forme très pratique :

$$ u_i = a_{0i} + (x_1 \underline{b}_1^T + x_2 \underline{b}_2^T + x_3 \underline{b}_3^T + h_1 \underline{\gamma}_1^T + h_2 \underline{\gamma}_2^T + h_3 \underline{\gamma}_3^T + h_4 \underline{\gamma}_4^T) \cdot \underline{d}_i \quad [12] $$

En différenciant cette dernière formule par rapport à $x_j$, on obtient le gradient de déplacement :

$$ u_{i,j} = \left( \underline{b}_j^T + \sum_{\alpha=1}^{4} h_{\alpha,j} \underline{\gamma}_{\alpha}^T \right) \cdot \underline{d}_i = \left( \underline{b}_j^T + h_{\alpha,j} \underline{\gamma}_{\alpha}^T \right) \cdot \underline{d}_i \quad [13] $$
Ceci permet d'exprimer l'opérateur gradient discrétisé reliant le champ des déformations aux coordonnées nodales par :

$$ \underline{\nabla}_s(\underline{u}) = \underline{B} \cdot \underline{d} $$
$$ \text{où : } \underline{\nabla}_s(\underline{u}) = \begin{bmatrix} u_{x,x} \\ u_{y,y} \\ u_{z,z} \\ u_{x,y} + u_{y,x} \\ u_{y,z} + u_{z,y} \\ u_{x,z} + u_{z,x} \end{bmatrix} , \quad \underline{d} = \begin{bmatrix} \underline{d}_x \\ \underline{d}_y \\ \underline{d}_z \end{bmatrix} \qquad [14] $$

Cet opérateur gradient discrétisé prend alors la forme matricielle pratique :

$$ \underline{B} = \begin{bmatrix} \underline{b}_{\leftarrow}^T + h_{a,x} \underline{\gamma}_{\alpha}^T & \underline{0} & \underline{0} \\ \underline{0} & \underline{b}_{\rightarrow}^T + h_{a,y} \underline{\gamma}_{\alpha}^T & \underline{0} \\ \underline{0} & \underline{0} & \underline{b}_{\leftarrow}^T + h_{a,z} \underline{\gamma}_{\alpha}^T \\ \underline{b}_{\rightarrow}^T + h_{a,y} \underline{\gamma}_{\alpha}^T & \underline{b}_{\leftarrow}^T + h_{a,x} \underline{\gamma}_{\alpha}^T & \underline{0} \\ \underline{0} & \underline{b}_{\rightarrow}^T + h_{a,z} \underline{\gamma}_{\alpha}^T & \underline{b}_{\leftarrow}^T + h_{a,y} \underline{\gamma}_{\alpha}^T \\ \underline{b}_{\leftarrow}^T + h_{a,z} \underline{\gamma}_{\alpha}^T & \underline{0} & \underline{b}_{\rightarrow}^T + h_{a,x} \underline{\gamma}_{\alpha}^T \end{bmatrix} \qquad [15] $$

Il est à noter que l'écriture ci-dessus de l'opérateur gradient discrétisé est très avantageuse car elle permet de manipuler séparément chacun des modes non constants de la déformation pour construire simplement la forme du champ de déformation postulée. De plus, on peut montrer facilement que les vecteurs $\underline{\gamma}_{\alpha}$ intervenant dans l'expression de $\underline{B}$, vérifient les relations d'orthogonalité suivantes :

$$ \underline{\gamma}_{\alpha}^T \cdot \underline{x}_j = 0, \quad \underline{\gamma}_{\alpha}^T \cdot \underline{h}_{\beta} = \delta_{\alpha\beta} \qquad [16] $$

Les propriétés précédentes seront très utiles dans l'analyse de stabilité de type hourglass de l'élément SHB8PS.
## 2.3. Modes de hourglass pour l'élément SHB8PS

En suivant la démarche de (Belytschko et al., 1993), nous allons analyser les modes de hourglass du SHB8PS. On montrera que ces modes cinématiques sont dus pour le SHB8PS à la position choisie pour les points d'intégration (sur une même droite). Ces modes sont associés à une énergie nulle alors qu'ils induisent une déformation non nulle. Cette anomalie s'explique par la différence qui existe entre le noyau de l'opérateur de rigidité discrétisé et celui continu. Rappelons que pour donner au SHB8PS un comportement de type coque, nous avons modifié la loi de comportement tridimensionnelle pour se rapprocher des hypothèses de contraintes planes et nous utilisons cinq points d'intégration disposés à travers une direction privilégiée que l'on appelle épaisseur. Cette sous-intégration a pour but aussi d'améliorer les performances de l'élément en termes de temps de calcul et d'éviter certains blocages en cisaillement dans les problèmes à flexion dominante. Ainsi, la raideur élastique s'obtient par intégration de Gauss :

$$ \underline{K}_e = \int_{\Omega_e} \underline{B}^T \underline{C} \underline{B} d\Omega = \sum_{i=1}^{5} \omega(\zeta_i) J(\zeta_i) \underline{B}^T(\zeta_i) \underline{C} \underline{B}(\zeta_i) \quad [17] $$

où $J(\zeta_i)$ représente le jacobien de la transformation entre la configuration de référence unitaire et la configuration courante d'un hexaèdre arbitraire. Le tableau ci-dessous donne les coordonnées des cinq points d'intégration du SHB8PS ainsi que les poids associés qui sont les racines du polynôme de Gauss-Legendre :

<table><caption>**Tableau 1.** Coordonnées et poids des points d'intégration du SHB8PS.</caption><thead><tr><th></th><th>ξ</th><th>η</th><th>ζ</th><th>ω</th></tr></thead><tbody><tr><td>P(1)</td><td>0</td><td>0</td><td>-0,91</td><td>0,24</td></tr><tr><td>P(2)</td><td>0</td><td>0</td><td>-0,54</td><td>0,48</td></tr><tr><td>P(3)</td><td>0</td><td>0</td><td>0</td><td>0,57</td></tr><tr><td>P(4)</td><td>0</td><td>0</td><td>0,54</td><td>0,48</td></tr><tr><td>P(5)</td><td>0</td><td>0</td><td>0,91</td><td>0,24</td></tr></tbody></table>

En remarquant qu'aux cinq points de Gauss (I=1,...,5) ci-dessus de coordonnées $\xi_i = \eta_i = 0$, $\zeta_i \neq 0$, les termes $h_{a,i}$ ($\alpha = 3,4$; $i = 1,2,3$) sont nuls, l'opérateur $\underline{B}$
défini par l'équation [15] se réduit à $\underline{B}_{12}$ dans lequel les sommations sur l'indice $\alpha$ ne vont que de 1 à 2 :

$$ \underline{B}_{12} = \begin{bmatrix} \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} & \underline{0} & \underline{0} \\ \underline{0} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,y} \underline{\gamma}_{\alpha}^{T} & \underline{0} \\ \underline{0} & \underline{0} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} \\ \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,y} \underline{\gamma}_{\alpha}^{T} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} & \underline{0} \\ \underline{0} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1,2} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,y} \underline{\gamma}_{\alpha}^{T} \\ \underline{b}_{\alpha}^{T} + \sum_{\alpha=1}^{2} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} & \underline{0} & \underline{b}_{\alpha}^{T} + \sum_{\alpha=1,2} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} \end{bmatrix} \quad [18] $$

Pour étudier le noyau de la rigidité, on construit une base des déplacements discrétisés et on montre que l'intégration réduite a baissé le rang de la raideur discrétisée. En effet d'après l'équation [17], on voit que le rang de la raideur $K_e$ est lié à celui de l'opérateur $B$. En d'autre termes, on est amené à rechercher les modes $d$ à déformation nulle, c'est-à-dire vérifiant en chacun des points de Gauss :

$$ \nabla_{\underline{s}}(\underline{u}) = \underline{B}(\zeta_{\underline{s}}) \cdot \underline{d} = \underline{0} \quad [19] $$

En utilisant l'expression [18] de l'opérateur gradient discrétisé évalué aux points de Gauss et avec l'aide des relations d'orthogonalité [10] et [16], on arrive à explicitier complètement le noyau de la raideur. On y trouve naturellement les six modes rigides qui sont les seuls vecteurs du noyau d'une raideur intégrée exactement :

$$ \begin{pmatrix} \underline{s} \\ \underline{0} \\ \underline{0} \end{pmatrix}, \begin{pmatrix} \underline{0} \\ \underline{s} \\ \underline{0} \end{pmatrix}, \begin{pmatrix} \underline{0} \\ \underline{0} \\ \underline{s} \end{pmatrix}, \begin{pmatrix} \underline{y} \\ -\underline{x} \\ \underline{0} \end{pmatrix}, \begin{pmatrix} \underline{z} \\ \underline{0} \\ -\underline{x} \end{pmatrix}, \begin{pmatrix} \underline{0} \\ \underline{z} \\ -\underline{y} \end{pmatrix} \quad [20] $$

Les trois premiers vecteurs colonnes correspondent aux translations selon les axes $Ox$, $Oy$ et $Oz$, respectivement. Les trois derniers vecteurs sont relatifs aux rotations autour des axes $Oz$, $Oy$ et $Ox$, respectivement.
Outre ces six mouvements de corps rigide, on montre que les six vecteurs suivants sont dans le noyau de la raideur $K$ :

$$ \begin{pmatrix} h_3 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ h_3 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ h_3 \end{pmatrix}, \begin{pmatrix} h_4 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ h_4 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ h_4 \end{pmatrix} \qquad [21] $$

La figure 2 schématise les modes de *hourglass* selon l'axe *Ox* pour un hexaèdre ayant un seul point de Gauss situé à l'origine des coordonnées de référence. En faisant tourner les axes des coordonnées, on peut obtenir de la même manière les modes de *hourglass* selon l'axe *Oy* ou selon *Oz*.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Modes de hourglass selon Ox pour un hexaèdre à un seul point de Gauss.</figcaption></figure>

## 2.4. Stabilisation et contrôle des modes cinématiques

Pour contrôler les six modes de *hourglass* du SHB8PS donnés par l'équation [21], il faut rajouter à la raideur $K$ une rigidité de stabilisation. Pour cela, nous nous sommes inspirés de la stratégie proposée par (Belytschko et al., 1993). Il est utile de noter que nous séparons parfaitement cette étape de stabilisation de celle de la projection de type « Assumed Strain », dont l'objectif est d'éliminer les blocages. Cette dernière étape de projection sera traitée plus tard dans le paragraphe 2.5.
On commence d'abord par décomposer l'opérateur $\underline{B}$ en la somme suivante de deux termes :

$$ \underline{B} = \underline{B}_{12} + \underline{B}_{34} \quad [22] $$

Le premier terme de cette décomposition est donné par l'équation [18]. Le second terme $\underline{B}_{34}$ est celui qui s'annule aux points de Gauss. Il est donné par la matrice :

$$ \underline{B}_{34} = \begin{bmatrix} \sum_{\alpha=3}^{4} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} & \underline{0} & \underline{0} \\ \underline{0} & \sum_{\alpha=3}^{4} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} & \underline{0} \\ \underline{0} & \underline{0} & \sum_{\alpha=3}^{4} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} \\ \sum_{\alpha=3}^{4} h_{\alpha,y} \underline{\gamma}_{\alpha}^{T} & \sum_{\alpha=3}^{4} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} & \underline{0} \\ \underline{0} & \sum_{\alpha=3,4} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} & \sum_{\alpha=3}^{4} h_{\alpha,y} \underline{\gamma}_{\alpha}^{T} \\ \sum_{\alpha=3}^{4} h_{\alpha,z} \underline{\gamma}_{\alpha}^{T} & \underline{0} & \sum_{\alpha=3,4} h_{\alpha,x} \underline{\gamma}_{\alpha}^{T} \end{bmatrix} \quad [23] $$

Dans une approche standard en déplacement, la raideur et les forces internes sont données par les formules :

$$ \underline{K}_e = \int_{\Omega_e} \underline{B}^T \cdot \underline{C} \cdot \underline{B} \, d\Omega \quad [24] $$

$$ \underline{f}^{\text{int}} = \int_{\Omega_e} \underline{B}^T \cdot \underline{\sigma} \cdot \underline{d} \Omega $$

En utilisant la décomposition de l'opérateur $\underline{B}$ donnée par l'équation [22], la rigidité s'exprime alors comme :
$$
\begin{equation}
\begin{aligned}
\underline{K}_e = & \int_{\Omega_e} \underline{B}_{12}^T \cdot \underline{C} \cdot \underline{B}_{12} \, d\Omega + \int_{\Omega_e} \underline{B}_{12}^T \cdot \underline{C} \cdot \underline{B}_{34} \, d\Omega \\
& + \int_{\Omega_e} \underline{B}_{34}^T \cdot \underline{C} \cdot \underline{B}_{12} \, d\Omega + \int_{\Omega_e} \underline{B}_{34}^T \cdot \underline{C} \cdot \underline{B}_{34} \, d\Omega
\end{aligned}
\tag{25}
\end{equation}
$$

Pour simplifier les écritures, on utilise alors la notation :

$$
\underline{K}_e = \underline{K}_{12} + \underline{K}_{STAB} \qquad [26]
$$

Le premier terme $\underline{K}_{12}$ est le seul pris en compte lorsque la rigidité est évaluée aux points de Gauss précédemment définis. Il s'écrit :

$$
\underline{K}_{12} = \int_{\Omega_e} \underline{B}_{12}^T \cdot \underline{C} \cdot \underline{B}_{12} \, d\Omega = \sum_{I=1}^{5} \omega(\zeta_I) J(\zeta_I) \underline{B}_{12}^T(\zeta_I) \cdot \underline{C} \cdot \underline{B}_{12}(\zeta_I) \quad [27]
$$

Le deuxième terme $\underline{K}_{STAB}$ représente la rigidité de stabilisation car il est nul si évalué aux points de Gauss :

$$
\underline{K}_{STAB} = \int_{\Omega_e} \underline{B}_{12}^T \cdot \underline{C} \cdot \underline{B}_{34} \, d\Omega + \int_{\Omega_e} \underline{B}_{34}^T \cdot \underline{C} \cdot \underline{B}_{12} \, d\Omega + \int_{\Omega_e} \underline{B}_{34}^T \cdot \underline{C} \cdot \underline{B}_{34} \, d\Omega \quad [28]
$$

Pour les forces internes de l'élément, le même type de raisonnement conduit à l'écriture :

$$
\underline{f}^{\text{int}} = \underline{f}_{12}^{\text{int}} + \underline{f}^{\text{STAB}} \quad [29]
$$

Le premier terme $\underline{f}_{12}^{\text{int}}$ est le seul pris en compte lors de l'évaluation des forces internes au moyen des points de Gauss et s'écrit :

$$
\underline{f}_{12}^{\text{int}} = \int_{\Omega_e} \underline{B}_{12}^T \cdot \underline{\sigma} \, d\Omega = \sum_{I=1}^{5} \omega(\zeta_I) J(\zeta_I) \underline{B}_{12}^T(\zeta_I) \cdot \underline{\sigma}(\zeta_I) \quad [30]
$$
Le deuxième terme $f^{STAB}$ de l'équation [29] représente les forces de stabilisation qu'il convient de calculer en accord avec la raideur de stabilisation donnée en [28].

Pour prendre en compte les contributions de la raideur et des forces de stabilisation, qui sont nulles aux points de Gauss, et éviter ainsi les phénomènes liés aux modes de hourglass, on choisit de les calculer dans le repère co-rotationnel proposé par (Belytschko et al., 1993). Une étape intermédiaire à ce calcul consiste à projeter $\underline{\underline{B}}$ en un opérateur $\underline{\underline{B}}$ pour éliminer les problèmes de verrouillage.

## 2.5. Projection de type « Assumed Strain »

Pour éliminer les différents blocages en cisaillement transverse ou en membrane, l'opérateur gradient discrétisé est projeté sur un sous-espace approprié. Cette technique de projection peut se mettre sous le formalisme de la méthode de déformation postulée « Assumed Strain ». On montre aussi que cette démarche peut être justifiée dans le cadre du principe variationnel mixte statique non linéaire de Hu-Washizu (voir par exemple Korelec et al., 1996). En effet, ce principe variationnel s'écrit :

$$ \delta\pi(\underline{\underline{v}}, \underline{\underline{\varepsilon}}, \underline{\underline{\sigma}}) = \int_{\Omega_e} \delta\underline{\underline{\varepsilon}}^T \cdot \underline{\underline{\sigma}} d\Omega + \delta\int_{\Omega_e} \underline{\underline{\sigma}}^T \cdot (\underline{\underline{\nabla}}_s(\underline{\underline{v}}) - \underline{\underline{\varepsilon}}) d\Omega - \delta\underline{\underline{d}}^T \cdot \underline{\underline{f}}^{\text{ext}} = 0 \quad [31] $$

où $\delta$ représente une variation, $\underline{\underline{v}}$ le champ de vitesse, $\underline{\underline{\varepsilon}}$ est le taux de déformation postulée, $\underline{\underline{\sigma}}$ la contraint postulée, $\underline{\underline{\varepsilon}}$ la contrainte évaluée par la loi de comportement, $\underline{\underline{d}}$ les vitesses nodales, $\underline{\underline{f}}^{\text{ext}}$ les forces nodales externes et $\underline{\underline{\nabla}}_s(\underline{\underline{v}})$ la partie symétrique du gradient des vitesses. La formulation « Assumed strain » retenue pour construire le SHB8PS est une forme simplifiée du principe variationnel de Hu-Washizu comme elle a été décrite par (Simo et al., 1986). Dans cette forme simplifiée, la contrainte postulée est choisie orthogonale à la différence entre la partie symétrique du gradient des vitesses et le taux de déformation postulée. Ainsi, le second terme de l'équation [31] s'élimine et l'on obtient :

$$ \delta\pi(\underline{\underline{v}}, \underline{\underline{\varepsilon}}) = \int_{\Omega_e} \delta\underline{\underline{\varepsilon}}^T \cdot \underline{\underline{\sigma}} d\Omega - \delta\underline{\underline{d}}^T \cdot \underline{\underline{f}}^{\text{ext}} = 0 \quad [32] $$

Sous cette forme, le principe variationnel est indépendant de l'interpolation de la contrainte, puisque la contrainte postulée n'intervient plus et n'a donc plus besoin d'être définie. Les équations discrétisées nécessitent donc la seule interpolation de la vitesse et de la déformation postulée. En exprimant le champ « Assumed Strain » $\underline{\underline{\varepsilon}}$
en fonction d'un opérateur $\overline{\mathbf{B}}$ projeté à partir de l'opérateur gradient discrétisé classique $\underline{\mathbf{B}}$ défini en [14] et [15], il vient :

$$ \dot{\underline{\mathbf{\epsilon}}}(x,t) = \overline{\mathbf{B}}(x) \cdot \dot{\underline{\mathbf{d}}}(t) \qquad [33] $$

Remplaçant cette expression dans le principe variationnel [32], on obtient immédiatement les nouvelles formules de la raideur élastique et des forces internes :

$$ \underline{\mathbf{K}}_e = \int_{\Omega_e} \overline{\mathbf{B}}^T \cdot \underline{\mathbf{C}} \cdot \overline{\mathbf{B}} \, d\Omega, \quad \underline{\mathbf{f}}^{\text{int}} = \int_{\Omega_e} \overline{\mathbf{B}}^T(x) \cdot \underline{\mathbf{\sigma}} \, d\Omega \qquad [34] $$

Avant de définir l'opérateur projeté $\overline{\mathbf{B}}$, nous remplaçons dans toutes les formules précédentes la forme de Hallquist des vecteurs $\underline{\mathbf{b}}_i$ (équation [9]) par la forme moyenne $\hat{\underline{\mathbf{b}}}_i$ de (Flanagan-Belytschko, 1981) :

$$ \hat{\underline{\mathbf{b}}}_i^T = \frac{1}{\Omega_e} \int_{\Omega_e} \underline{\mathbf{N}}_{,i}(\underline{\mathbf{\xi}}, \underline{\mathbf{\eta}}, \underline{\mathbf{\zeta}}) d\Omega, \quad i = 1, 2, 3 \qquad \text{Formules moyennes} \qquad [35] $$

Ceci consiste aussi à remplacer les vecteurs $\underline{\mathbf{\gamma}}_\alpha$ par les vecteurs $\hat{\underline{\mathbf{\gamma}}}_\alpha$ définis par :

$$ \hat{\underline{\mathbf{\gamma}}}_\alpha = \frac{1}{8} \left[ \underline{\mathbf{h}}_\alpha - \sum_{j=1}^{3} (\underline{\mathbf{h}}_\alpha^T \cdot \underline{\mathbf{x}}_j) \hat{\underline{\mathbf{b}}}_j \right] \qquad [36] $$

Par conséquent, l'opérateur $\underline{\mathbf{B}}$ donné en [15] est remplacé par l'opérateur $\hat{\underline{\mathbf{B}}}$ défini par :
$$
\hat{\mathbf{B}} = \begin{bmatrix} \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,x} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,y} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,z} \hat{\mathbf{\gamma}}_{\alpha}^{T} \\ \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,y} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,x} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \mathbf{0} \\ \mathbf{0} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,z} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,y} \hat{\mathbf{\gamma}}_{\alpha}^{T} \\ \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,z} \hat{\mathbf{\gamma}}_{\alpha}^{T} & \mathbf{0} & \hat{\mathbf{b}}_{\alpha}^{T} + h_{\alpha,x} \hat{\mathbf{\gamma}}_{\alpha}^{T} \end{bmatrix} \quad [37]
$$

Toute la démarche précédente reste valable ainsi que les expressions de la
rigidité et des forces internes de stabilisation si l'on adopte la même décomposition
additive :

$$
\bar{\hat{\mathbf{B}}} = \bar{\mathbf{B}}_{12} + \bar{\mathbf{B}}_{34} \quad [38]
$$

Notons ici une première différence importante avec la précédente version du SHB8PS dans laquelle, les formules de Hallquist $\hat{b}_i$ ont été remplacées par les formules moyennes $\hat{b}_i$ de Flanagan-Belytschko seulement dans les termes de stabilisations $\hat{B}_{34}$ et donc $\underline{K}_{STAB}$.

Il est important de noter aussi que les deux formes $b_i$ et $\hat{b}_i$ ont été testées sur un grand nombre de cas tests, et que la forme moyenne de Flanagan-Belytschko s'est montrée plus efficace pour l'ensemble des cas. Cette supériorité en termes de convergence de la forme moyenne est encore plus marquée pour des maillages avec peu d'éléments distordus. Des remarques similaires avaient été constatées dans (Belytschko et al., 1993) avec un élément « assumed strain » solide à huit nœuds.

On peut à ce stade projeter l'opérateur $\hat{\underline{\mathbf{B}}}$ de l'équation [38] en un opérateur $\underline{\hat{\mathbf{B}}}$ tel que :

$$
\bar{\underline{\hat{\mathbf{B}}}} = \underline{\hat{\mathbf{B}}}_{12} + \underline{\hat{\mathbf{B}}}_{34} \qquad [39]
$$

Il est à noter que seul le deuxième terme $\underline{\hat{\mathbf{B}}}_{34}$ de l'équation [38] est projeté, le premier terme $\underline{\hat{\mathbf{B}}}_{12}$ reste inchangé et donné par l'équation [18] dans laquelle les vecteurs $\underline{\hat{b}}_i$ sont remplacés par les vecteurs $\hat{b}_i$. Ainsi, l'opérateur $\underline{\hat{\mathbf{B}}}_{34}$ est projeté en $\underline{\hat{\mathbf{B}}}_{34}$ donné par l'expression :
$$
\hat{\underline{\underline{B}}}_{34} = \begin{bmatrix} \sum_{\alpha=3}^{4} h_{\alpha,x} \hat{\gamma}_{\alpha}^{T} & \underline{0} & \underline{0} \\ \underline{0} & \sum_{\alpha=3}^{4} h_{\alpha,y} \hat{\gamma}_{\alpha}^{T} & \underline{0} \\ \underline{0} & \underline{0} & h_{3,z} \hat{\gamma}_{3}^{T} \\ \underline{0} & \underline{0} & \underline{0} \\ \underline{0} & \underline{0} & h_{4,x} \hat{\gamma}_{4}^{T} \\ \underline{0} & \underline{0} & \underline{0} \end{bmatrix} \quad [40]
$$

La raideur élastique est alors donnée par l'équation [26] représentant le somme
des deux contributions suivantes :

$$
\underline{\underline{K}}_{12} = \int_{\Omega_e} \hat{\underline{\underline{B}}}^{T}_{12} \cdot \underline{\underline{C}} \cdot \hat{\underline{\underline{B}}}_{12} \, d\Omega = \sum_{l=1}^{5} \omega(\zeta_l) J(\zeta_l) \hat{\underline{\underline{B}}}^{T}_{12}(\zeta_l) \cdot \underline{\underline{C}} \cdot \hat{\underline{\underline{B}}}_{12}(\zeta_l) \quad [41]
$$

$$
\underline{\underline{K}}_{\text{STAB}} = \int_{\Omega_e} \hat{\underline{\underline{B}}}^{T}_{12} \cdot \underline{\underline{C}} \cdot \hat{\underline{\underline{B}}}_{34} \, d\Omega + \int_{\Omega_e} \hat{\underline{\underline{B}}}^{T}_{34} \cdot \underline{\underline{C}} \cdot \hat{\underline{\underline{B}}}_{12} \, d\Omega + \int_{\Omega_e} \hat{\underline{\underline{B}}}^{T}_{34} \cdot \underline{\underline{C}} \cdot \hat{\underline{\underline{B}}}_{34} \, d\Omega \quad [42]
$$

Pour le calcul de cette rigidité de stabilisation on se place dans un repère co-
rotationnel (Belytschko et al., 1993). Le choix de ce repère tournant avec l’élément
et aligné avec le repère de référence est justifié ici par la rotation extraite de la
décomposition en produit du gradient de la transformation (Abed-Meraim et al.,
2001, 2002). Les avantages de cette formulation co-rotationnelle sont multiples :
simplification des formules donnant la matrice de stabilisation ci-dessus dont les
deux premiers termes croisés s’annulent, le traitement du blocage en cisaillement est
plus efficace dans ce repère local, cette formulation co-rotationnelle est mieux
adaptée aux non linéarités géométriques.

Nous allons donner brièvement les principales formules définissant le repère co-
rotationnel choisi. On commence par calculer les composantes des vecteurs
colonnes formant la matrice de rotation :

$$
a_{li} = \underline{\underline{\Lambda}}_{1}^{T} \cdot \underline{x}_{i} \quad , \quad a_{2i} = \underline{\underline{\Lambda}}_{2}^{T} \cdot \underline{x}_{i} \quad , \quad i = 1,2,3 \qquad [43]
$$

avec :
$$
\left\{
\begin{array}{l}
\underline{\Delta}_1^T = (-1, 1, 1, -1, -1, 1, 1, -1) \\
\underline{\Delta}_2^T = (-1, -1, 1, 1, -1, -1, 1, 1) \\
\underline{\Delta}_3^T = (-1, -1, -1, -1, 1, 1, 1, 1)
\end{array}
\right.
\qquad [44]
$$

Ensuite, on calcule le terme de correction $\underline{a}_c$ de sorte à vérifier la relation d'orthogonalité : $\underline{a}_1^T \cdot (\underline{a}_2 + \underline{a}_c) = 0$, soit :

$$
\underline{a}_c = - \frac{\underline{a}_1^T \cdot \underline{a}_2}{\underline{a}_1^T \cdot \underline{a}_1} \underline{a}_1 \quad [45]
$$

Le troisième vecteur de base $\underline{a}_3$ est alors obtenu par le produit vectoriel :

$$
\underline{a}_3 = \underline{a}_1 \wedge (\underline{a}_2 + \underline{a}_c) \quad [46]
$$

La matrice de rotation $\underline{R}$ qui permet la transformation du repère général au repère co-rotationnel est finalement donnée après normalisation par :

$$
R_{1i} = \frac{a_{1i}}{\|a_i\|} , \quad R_{2i} = \frac{a_{2i} + a_{ci}}{\|a_2 + a_c\|} , \quad R_{3i} = \frac{a_{3i}}{\|a_3\|} , \quad i = 1, 2, 3 \quad [47]
$$

On fera le passage dans ce repère local pour calculer les termes de stabilisation (forces internes et rigidité de stabilisation [42]) dont beaucoup de termes se simplifient dans un tel repère. En effet, dans ce repère local on a les formules :

$$
\left\{
\begin{array}{l}
\displaystyle \int_{\Omega_i} h_{i,j} \, d\Omega = 0 \\[2ex]
H_{li} = \int_{\Omega_i} (h_{j,i})^2 \, d\Omega = \int_{\Omega_i} (h_{k,i})^2 \, d\Omega = 3 \int_{\Omega_i} (h_{4,i})^2 \, d\Omega = \frac{1}{3} \frac{(\underline{\Delta}_i^T \cdot \underline{x}) (\underline{\Delta}_i^T \cdot \underline{x})}{(\underline{\Delta}_i^T \cdot \underline{x})} \quad [48] \\[2ex]
H_{ij} = \int_{\Omega_i} h_{i,j} h_{j,i} \, d\Omega = \frac{1}{3} \underline{\Delta}_k^T \cdot \underline{x}
\end{array}
\right.
$$
Dans les formules ci-dessus, il n'y a pas de sommation sur les indices répétés. De plus, les indices *i*, *j* et *k* qui interviennent dans les termes $H_{ii}$ et $H_{ij}$ sont distincts 2 à 2 et prennent les valeurs 1, 2, 3 avec toutes les permutations possibles. Utilisant ces formules explicites, la rigidité de stabilisation exprimée en [42] s'obtient complètement analytiquement dans ce repère local par :

$$ \underline{K}_{\text{STAB}} = \begin{bmatrix} \underline{k}_{11} & \underline{k}_{12} & \underline{k}_{13} \\ \underline{k}_{21} & \underline{k}_{22} & \underline{k}_{23} \\ \underline{k}_{31} & \underline{k}_{32} & \underline{k}_{33} \end{bmatrix} \qquad [49] $$

où les matrices $\underline{k}_{ij}$ de dimension 8 x 8 sont données par :

$$ \left\{ \begin{array}{l} \underline{k}_{11} = (\bar{\lambda} + 2\mu)H_{11} \left[ \hat{\gamma}_3 \hat{\gamma}_3^T + \frac{1}{3} \hat{\gamma}_4 \hat{\gamma}_4^T \right] \\ \underline{k}_{22} = (\bar{\lambda} + 2\mu)H_{22} \left[ \hat{\gamma}_3 \hat{\gamma}_3^T + \frac{1}{3} \hat{\gamma}_4 \hat{\gamma}_4^T \right] \\ \underline{k}_{33} = \mu H_{11} \frac{1}{3} \hat{\gamma}_4 \hat{\gamma}_4^T \\ \underline{k}_{ij} = 0 \quad , \quad i \neq j \end{array} \right. \qquad [50] $$

Il est à noter qu'une loi de comportement de type contraintes planes améliorée est adoptée pour l'élément SHB8PS. Cette loi spécifique est donnée par :

$$ \underline{C} = \begin{bmatrix} \bar{\lambda} + 2\mu & \bar{\lambda} & 0 & 0 & 0 & 0 \\ \bar{\lambda} & \bar{\lambda} + 2\mu & 0 & 0 & 0 & 0 \\ 0 & 0 & E & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{bmatrix} \qquad [51] $$ $$ \mu = \frac{E}{2(1+\nu)} \qquad \bar{\lambda} = \frac{E\nu}{1-\nu^2} $$

Ce choix de matrice de comportement nous évite, d'une part, les blocages rencontrés avec la loi tridimensionnelle complète, et d'autre part, permet de prendre
en compte l'énergie de déformation associée aux déformations normales à la surface moyenne du SHB8PS.

Pour le calcul des forces internes de l'élément, nous adoptons la même démarche (Abed-Meraim et al., 2001, 2002). La décomposition additive [39] ainsi que la projection [40] permettent de calculer la contribution des forces de stabilisation :

$$ \underline{f}^{\text{int}} = \sum_{i=1}^{5} \omega(\zeta_i) J(\zeta_i) \underline{\hat{B}}_{12}^{T}(\zeta_i) \underline{\sigma}(\zeta_i) + \underline{f}^{\text{STAB}} \quad [52] $$

avec :

$$ \underline{f}^{\text{STAB}} = \begin{pmatrix} \underline{f}_{1}^{\text{STAB}} \\ \underline{f}_{2}^{\text{STAB}} \\ \underline{f}_{3}^{\text{STAB}} \end{pmatrix} \quad [53] $$

et :

$$ \underline{f}_{i}^{\text{STAB}} = \sum_{\alpha=3}^{4} Q_{i\alpha} \hat{\gamma}_{\alpha} \quad , \quad i=1,2,3 \quad [54] $$

Les $Q_{i\alpha}$, dites contraintes généralisées, intervenant dans le calcul des forces de stabilisation, sont reliées aux $q_{i\alpha}$, dites déformations généralisées, par les relations suivantes incrémentales :

$$ \left\{ \begin{aligned} \dot{Q}_{13} &= (\bar{\lambda} + 2\mu)H_{11}\dot{q}_{13} \\ \dot{Q}_{14} &= \frac{1}{3}(\bar{\lambda} + 2\mu)H_{11}\dot{q}_{14} \\ \dot{Q}_{23} &= (\bar{\lambda} + 2\mu)H_{22}\dot{q}_{23} \\ \dot{Q}_{24} &= \frac{1}{3}(\bar{\lambda} + 2\mu)H_{22}\dot{q}_{24} \\ \dot{Q}_{33} &= 0 \\ \dot{Q}_{34} &= \frac{1}{3}\mu H_{11}\dot{q}_{34} \end{aligned} \right. \quad [55] $$

Les taux de déformations généralisées $\dot{q}_{i\alpha}$ sont donnés par :

$$ \dot{q}_{i\alpha} = \hat{\gamma}_{\alpha}^{T} \cdot \dot{d}_{i} \quad , \quad i=1,2,3 \quad , \quad \alpha=3,4 \quad [56] $$
Les expressions précédentes de la raideur et des forces de stabilisation sont valables pour un comportement élastique. Dans le cas d'un comportement élasto-plastique, le module d'Young *E* est remplacé par la moyenne du module tangent sur les cinq points de Gauss. Ce choix évite d'avoir une réponse trop raide comme c'est le cas avec une stabilisation *anti-hourglass* complètement élastique. Ceci conduit ainsi à une stabilisation qui s'adapte automatiquement à la situation physique de l'élément : élastique ou élasto-plastique.

### 3. Validation sur des cas tests

Pour valider cette nouvelle version du SHB8PS, nous l'avons testée, dans un premier temps, sur un ensemble de plus de dix cas tests élastiques. 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 la version précédente du SHB8PS (Legay *et al.*, 2003). Notons que plusieurs projections ont été formulées dans cette étude et testées sur l'ensemble des cas tests. La projection retenue et présentée ici est celle qui a montré la meilleure convergence et l'absence des phénomènes de blocages numériques en cisaillement transverse et en membrane. Cette projection a amélioré la version précédente du SHB8PS sur tous les cas tests et notamment sur celui de la sphère pincée où l'amélioration est très significative. Nous donnons ici les résultats relatifs à ce dernier cas test et à celui de la poutre vrillée avec les deux directions de sollicitations.

#### 3.1. Hémisphère pincé

Ce cas test est devenu très populaire et a été repris par de nombreux auteurs depuis (MacNeal *et al.*, 1985). Il est particulièrement sévère car les blocages en cisaillement transverse et en membrane sont importants et accentués par la géométrie du problème (caractère "distordu" et "gauche" des éléments). Ce problème a été étudié en détail par (Belytschko *et al.*, 1989) qui montrent que, comme tous les éléments sont courbes, ceci augmente l'intensité des blocages en membrane et en cisaillement. Ils montrent aussi que pour ce problème particulier, le blocage en membrane est plus sévère que le blocage en cisaillement. La figure 3 montre la géométrie, le chargement et les conditions aux limites pour ce problème. Le rayon est R=10, l'épaisseur t=0,04, le module d'Young E=6,825×10⁷, le coefficient de Poisson ν=0,3. En utilisant les plans de symétries (XZ) et (YZ), un quart de l'hémisphère seulement est maillé avec deux forces unitaires selon les directions axiales Ox et Oy. Les conditions aux limites sont libres, toutefois pour éviter les mouvements de corps rigides on bloque le déplacement selon Oz d'un point en plus des conditions relatives aux deux plans de symétrie. La solution de référence fournie dans (MacNeal *et al.*, 1985) donne un déplacement du point A dans la direction Ox égal à 0,0924 (voir figure 3).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Géométrie et chargement pour la coque sphérique pincée.</figcaption></figure>

Les résultats en termes de convergence vers la solution de référence sont exprimés par le tableau 2 qui donne le déplacement normalisé du point A dans la direction Ox. La nouvelle version du SHB8PS est comparée à l'ancienne et aux trois éléments HEX8, HEXDS et H8-ct-cp. L'élément HEX8 est l'élément standard volumique à huit nœuds et à intégration exacte (huit points de Gauss). L'élément HEXDS est un élément volumique à huit nœuds et quatre points de Gauss (Liu et al., 1998). L'élément H8-ct-cp est présenté dans le mémoire de thèse de Lemosse (Lemosse, 2000). Le tableau 2 montre que la nouvelle version du SHB8PS donne une excellente convergence et ne présente aucun blocage.

<table><caption>**Tableau 2.** Déplacement normalisé du point A de l'hémisphère.</caption><thead><tr><th></th><th>SHB8PS<br>ancienne<br>version</th><th>HEX8</th><th>HEXDS</th><th>H8-ct-cp</th><th>SHB8PS<br>nouvelle<br>version</th></tr><tr><th>Nombre<br>d'éléments</th><th>Ux/Uréf</th><th>Ux/Uréf</th><th>Ux/Uréf</th><th>Ux/Uréf</th><th>Ux/Uréf</th></tr></thead><tbody><tr><td>12</td><td>0,0629</td><td>0,0005</td><td></td><td>0,05</td><td>0,8645</td></tr><tr><td>27</td><td>0,0474</td><td>0,0011</td><td></td><td></td><td>1,0155</td></tr><tr><td>48</td><td>0,1660</td><td>0,0023</td><td>0,408</td><td>0,35</td><td>1,0098</td></tr><tr><td>75</td><td>0,2252</td><td>0,0030</td><td>0,512</td><td>0,58</td><td>1,0096</td></tr><tr><td>192</td><td>0,6332</td><td>0,0076</td><td>0,701</td><td>0,95</td><td>1,0008</td></tr><tr><td>363</td><td>0,8592</td><td>0,0140</td><td>0,800</td><td></td><td>1,0006</td></tr><tr><td>768</td><td>0,9651</td><td>0,0287</td><td></td><td></td><td>1,0006</td></tr><tr><td>1462</td><td>0,9910</td><td>0,0520</td><td></td><td></td><td>1,0009</td></tr></tbody></table>
## 3.2. Poutre vrillée

Ce cas test a été proposé par (MacNeal *et al.*, 1985) et a ensuite été repris par de nombreux auteurs pour tester les formulations solides ou coques dans le cas de maillages « tordus ». En effet, la géométrie vrillée d'un angle de 90° conduit à des éléments distordus rendant ainsi ce cas test assez sévère.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Géométrie et chargements pour la poutre vrillée.</figcaption></figure>

La figure 4 ci-dessus donne la géométrie de la poutre vrillée, les conditions aux limites et les deux types de chargements appliqués. L'extrémité gauche de cette poutre est encastrée et un effort tranchant P=1 est appliqué à son extrémité droite. Deux directions de chargement sont étudiées : P1 et P2. Les caractéristiques géométriques et matérielles pour ce cas test sont données dans le tableau 3.

<table><caption>**Tableau 3.** Caractéristiques géométriques et matérielles de la poutre.</caption><tr><td>Longueur</td><td>L</td><td>12</td></tr><tr><td>Largeur</td><td>l</td><td>1,1</td></tr><tr><td>Epaisseur</td><td>h</td><td>0,32</td></tr><tr><td>Module d'Young</td><td>E</td><td>29.10<sup>6</sup></td></tr><tr><td>Coefficient de Poisson</td><td>ν</td><td>0,22</td></tr><tr><td>Charge appliquée</td><td>P</td><td>1</td></tr></table>
### 3.2.1. Chargement P1

Ce premier type de chargement est dans le plan de la poutre et consiste à appliquer la charge P=1 dans la direction verticale Oz (voir figure 4). La solution de référence pour ce problème est donnée par (MacNeal et al., 1985). Pour ce chargement P1, la flèche du point A dans la direction Oz (voir figure 4) doit être égale à $5,42 \times 10^{-3}$.

<table><caption>**Tableau 4. Déplacement normalisé selon Oz du point A de la poutre vrillée.**</caption><thead><tr><th></th><th>HEX8</th><th>SHB8PS ancienne version</th><th>SHB8PS nouvelle version</th></tr></thead><tbody><tr><th>Nombre d'éléments</th><td>Uz/Uréf</td><td>Uz/Uréf</td><td>Uz/Uréf</td></tr><tr><th>2x1x1=2</th><td>0,008</td><td>1,035</td><td>1,041</td></tr><tr><th>4x1x1=4</th><td>0,031</td><td>0,995</td><td>1,014</td></tr><tr><th>12x2x1=24</th><td>0,206</td><td>0,997</td><td>0,999</td></tr><tr><th>24x4x1=96</th><td>0,490</td><td>0,998</td><td>0,999</td></tr></tbody></table>

Le tableau 4 ci-dessus donne les résultats pour ce cas test en termes de convergence en fonction du nombre d'éléments dans le maillage. Le déplacement normalisé du point A de la poutre dans la direction Oz du chargement est donné pour trois types d'éléments. En effet, on compare les performances de la nouvelle version du SHB8PS (dernière colonne du tableau 4) par rapport à l'ancienne version de cet élément, d'une part, et par rapport à l'élément HEX8. Ce dernier représente l'élément solide standard linéaire à huit nœuds et à intégration exacte (huit points de Gauss). On voit donc que les deux versions du SHB8PS présentent une excellente convergence (dès quatre éléments de maillage) alors que l'élément HEX8 exhibe un verrouillage numérique et une convergence beaucoup plus lente.

### 3.2.2. Chargement P2

Ce deuxième type de chargement est une sollicitation hors plan et revient à appliquer l'effort P=1 dans la direction horizontale Oy (voir figure 4). La solution de référence pour ce problème est donnée par (MacNeal et al., 1985). Pour ce chargement P2, le déplacement du point A dans la direction du chargement (voir figure 4) doit être égal à $1,75 \times 10^{-3}$.
<table><caption>**Tableau 5.** Déplacement normalisé selon Oy du point A de la poutre vrillée.</caption><thead><tr><th></th><th>HEX8</th><th>SHB8PS<br>ancienne<br>version</th><th>SHB8PS<br>nouvelle<br>version</th></tr><tr><th>Nombre<br>d'éléments</th><th>Uy/Uréf</th><th>Uy/Uréf</th><th>Uy/Uréf</th></tr></thead><tbody><tr><td>2x1x1=2</td><td>0,023</td><td>0,705</td><td>0,867</td></tr><tr><td>4x1x1=4</td><td>0,081</td><td>0,906</td><td>0,952</td></tr><tr><td>12x2x1=24</td><td>0,333</td><td>0,986</td><td>0,994</td></tr><tr><td>24x4x1=96</td><td>0,592</td><td>0,996</td><td>0,998</td></tr></tbody></table>

Dans le dernier tableau ci-dessus relatif au chargement hors du plan, on donne les résultats de convergence en fonction du nombre d'éléments dans le maillage. Cette fois c'est le déplacement normalisé du point A de la poutre dans la direction Oy du chargement qui est donné pour les trois types d'éléments. On peut ainsi comparer les performances de la nouvelle version du SHB8PS (dernière colonne du tableau 5) par rapport à l'ancienne version et par rapport à l'élément HEX8. On voit donc que les deux versions du SHB8PS présentent une excellente convergence (dès une dizaine d'éléments dans le maillage) alors que l'élément HEX8 montre un verrouillage numérique et une convergence plus lente. Pour les maillages à 2 et 4 éléments, les résultats sont moins bons que pour le chargement vertical. En effet, si l'on regarde la position des points de Gauss (figure 5), on s'aperçoit que dans la zone où le moment de flexion est maximal (à l'encastrement), il n'y a qu'un seul point de Gauss dans la direction de la sollicitation P2 alors qu'il y en a 5 dans la direction du chargement P1.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Position des points de Gauss lors du maillage de la poutre vrillée.</figcaption></figure>
## 4. Discussion et conclusions

Une nouvelle formulation de l'élément de coque massif « SHB8PS » a été effectuée et implantée dans le code de calcul par éléments finis non linéaire implicite *Stanlax-INCA*. Cette nouvelle version a ensuite été testée, dans un premier temps en linéaire, sur un ensemble de plus de dix cas tests élastiques. Rappelons que cet élément est sous-intégré pour améliorer ses performances de calculs et bénéfice d'une stabilisation efficace pour contrôler les modes parasites (*hourglass*). De plus, grâce à la projection utilisée dans cette nouvelle version du SHB8PS, nous allons beaucoup plus loin dans l'élimination des phénomènes de verrouillage. En effet, les performances de l'élément en termes d'efficacité et de rapidité de convergence sont clairement démontrées à travers de nombreux cas tests et en particulier celui de l'hémisphère pincé. Les différents cas tests effectués prouvent bien l'absence de blocages en membrane ou en cisaillement transverse. La validation de cette version du SHB8PS en non linéaire géométrique et matériau (flambage, plasticité ...) est en cours de réalisation. De même, une nouvelle version explicite de cet élément sera implantée dans un code de dynamique rapide pour la valider sur des problèmes d'impact ou de crash.

## 5. Bibliographie

* Abed-Meraim F., Combescure A., Stabilité des éléments finis sous-intégrés, rapport interne n° 247, janvier 2001, LMT de Cachan.

* 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 Comput. Fluid and Solid Mechanics*, 12-15 June 2001, U.S.A.

* Abed-Meraim F., Combescure A., « SHB8PS- a new adaptive, assumed-strain continuum mechanics shell element for impact analysis », *Computers & Structures*, vol. 80, 2002, p. 791-803.

* 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*, vol. 43, 1984, p. 251-276.

* Belytschko T, Wong, B.L., « Assumed strain stabilization procedure for the 9-node Lagrange shell element », *International Journal for Numerical Methods and Engineering*, vol. 28, 1989, p. 385-414.

* Belytschko T., Bindeman, L.P., « Assumed strain stabilization of the eight node hexahedral element », *Computer Methods in Applied Mechanics and Engineering*, vol. 105, 1993, p. 225-260.

* 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*, vol. 17, No. 1, 2004, p. 113-140.
Flanagan D.P., Belytschko T., « A uniform strain hexahedron and quadrilateral with orthogonal hourglass control », *International Journal for Numerical Methods and Engineering*, vol. 17, 1981, p. 679-706.

Hallquist J.O., *Theoretical manual for DYNA3D*, UC1D-19401 Lawrence Livemore National Lab., University of California, 1983.

Korelec J., Wriggers P., « Consistent gradient formulation for a stable enhanced strain method for large deformations », *Engineering Computations*, vol. 13, 1996, p. 103-123.

Legay A., Combescure A., « Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS », *International Journal for Numerical Methods and Engineering*, vol. 57, 2003, p. 1299-1322.

Lemosse D., Eléments finis isoparamétriques tridimensionnels pour l'étude des structures minces, Thèse de doctorat, Ecole Doctorale SPMI/INSA-Rouen, 2000.

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, 1998, p. 69-132.

MacNeal R.H., Harder R.L., « A proposed standard set of problems to test finite element accuracy », *Finite Elements in Analysis and Design*, vol. 1, 1985, p. 3-20.

Simo J.C., Hughes T.J.R., « On the variational foundations of assumed strain methods », *Journal of Applied Mechanics*, ASME, vol. 53, 1986, p. 51-54.

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 and Engineering*, vol. 48, 2000, p. 545-564.