# Développement d'un nouvel élément fini “ SHB6 ” de type solide-coque

Vuong-Dieu Trinh, Farid Abed-Meraim, Alain Combescure

► To cite this version:

Vuong-Dieu Trinh, Farid Abed-Meraim, Alain Combescure. Développement d'un nouvel élément fini "SHB6" de type solide-coque. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504291

HAL Id: hal-01504291

https://hal.science/hal-01504291

Submitted on 9 Apr 2017

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.

Public Domain
# Développement d'un nouvel élément fini
## « SHB6 » de type solide-coque

Formulation et évaluation à travers des cas tests

Vuong-Dieu Trinh* — 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
[vuong-dieu.trinh@metz.ensam.fr](mailto:vuong-dieu.trinh@metz.ensam.fr), [farid.abed-meraim@metz.ensam.fr](mailto: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](mailto:alain.combescure@insa-lyon.fr)

RÉSUMÉ. Cet article décrit le développement d'un nouvel élément fini prismatique SHB6 de type solide-coque, obtenu à partir d'une formulation purement tridimensionnelle. Cet élément possède six nœuds et cinq points d'intégration répartis selon la direction de l'épaisseur. L'objectif étant d'avoir des éléments à géométrie volumique capables de 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. On montre d'abord que cette sous-intégration ne génère pas de modes de hourglass. Ensuite, on met en évidence que l'élément SHB6, sans aucune modification ou projection de son opérateur gradient discrétisé, peut souffrir de certains verrouillages de type cisaillement transverse ou membrane.

ABSTRACT. This paper presents the development of a new solid-shell finite element "SHB6" derived from a purely three-dimensional formulation. It has six nodes as well as five integration points, all distributed along the "thickness" direction. The main goal of this research is to develop low-order solid elements that are able to model thin structures while correctly taking into account the various through-thickness phenomena. In order to improve its calculation performances and to prevent some locking phenomena, reduced integration was used. We demonstrate first that there are no hourglass modes generated by the reduced integration. On the other hand, we show that, without any modification or projection of its discrete gradient operator, the SHB6 element could suffer from some membrane and shear locking phenomena.

MOTS-CLÉS: Coque massive SHB6, Hourglass, Verrouillage en membrane et cisaillement.

KEYWORDS: SHB6 solid-shell, Hourglass, Shear and membrane locking.
# 1. Introduction

Un effort important a été consacré ces dernières années au développement d'éléments finis de type coques volumiques. L'objectif est de disposer d'éléments permettant de modéliser des structures minces, à faible coût, tout en facilitant la connexion avec les éléments volumiques traditionnels. Parmi ces travaux, on peut citer (Lemosse, 2000), (Sze et al., 2000), (Abed-Meraim et al., 2001, 2002), (Legay et al., 2003), (Chen et al., 2004). Un élément hexaédrique de cette famille appelé SHB8PS a été récemment développé, on se propose ici de le compléter par son équivalent prismatique SHB6. Afin d'améliorer l'efficacité de l'élément et de réduire certains blocages en membrane et cisaillement, l'intégration réduite est employée. Cette technique se retrouve aussi dans les travaux de (Zienkiewicz et al., 1971), (Pawsey et al., 1971), (Hughes et al., 1977), (Flanagan et al., 1981) et bien d'autres. Cependant, le recours à la sous-intégration peut conduire à des éléments instables caractérisés par la présence de modes à énergie nulle dits aussi modes de hourglass.

# 2. Formulation de l'élément coque volumique SHB6

## 2.1. Interpolation de l'élément

Pour cet élément linéaire isoparamétrique à six nœuds SHB6, les coordonnées $x_i$ et les déplacements $u_i$ sont reliés aux coordonnées et déplacements aux nœuds $x_{il}$ et $u_{il}$, respectivement, au moyen des fonctions de forme $N_i$ par les formules :

$$ \left\{ \begin{aligned} x_i &= x_{il} N_i(\xi, \eta, \zeta) = \sum_{I=1}^{6} x_{il} N_i(\xi, \eta, \zeta) \\ u_i &= u_{il} N_i(\xi, \eta, \zeta) = \sum_{I=1}^{6} u_{il} N_i(\xi, \eta, \zeta) \end{aligned} \right. , \quad i=1,2,3; \quad I=1,2,...,6 \quad [1] $$

où les fonctions de forme linéaires $N_i(\xi, \eta, \zeta)$ sont données par les expressions :

$$ \left\{ \begin{aligned} N_1 &= \frac{1}{2}(1-\zeta)(1-\eta-\xi), & N_2 &= \frac{1}{2}(1-\zeta)\xi, & N_3 &= \frac{1}{2}(1-\zeta)\eta, \\ N_4 &= \frac{1}{2}(1+\zeta)(1-\eta-\xi), & N_5 &= \frac{1}{2}(1+\zeta)\xi, & N_6 &= \frac{1}{2}(1+\zeta)\eta, \\ \zeta &\in [0,1], & \eta &\in [0,1-\zeta], & \zeta &\in [-1,1] \end{aligned} \right. \quad [2] $$

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

Pour obtenir une forme explicite de l'opérateur gradient discrétisé, nous allons développer le champ de déplacement en une somme constituée d'un terme constant,
de termes linéaires en $x_i$, et de termes faisant intervenir les fonctions $h_\alpha$ ($\alpha = 1, 2$). Ceci est rendu possible par la combinaison des équations [1] et [2] :

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

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

$$ \underline{d}_i = a_{0i} \underline{s} + a_{1i} \underline{x} + a_{2i} \underline{y} + a_{3i} \underline{z} + c_{1i} \underline{h}_1 + c_{2i} \underline{h}_2, \quad i=1,2,3 \quad [4] $$

où $\underline{s}^T = (1,1,1,1,1)$, $\underline{h}_1^T = (0,0,-1,0,0,1)$, $\underline{h}_2^T = (0,-1,0,0,1,0)$ et les vecteurs $\underline{d}_i$ et $\underline{x}_i$ ($i=1,2,3$) désignent, respectivement, les déplacements et les coordonnées aux nœuds. Les constantes inconnues $a_{ji}$ et $c_{\alpha i}$ de l'équation [4] sont obtenues par l'introduction des vecteurs $\underline{b}_i^T$ de (Hallquist, 1983), définis par :

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

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

Ceci permet d'exprimer l'opérateur gradient discrétisé sous la forme pratique :

$$ \underline{B} = \begin{bmatrix} \underline{b}_x^T + h_{\alpha,x} \underline{\gamma}_\alpha^T & \underline{0} & \underline{0} \\ \underline{0} & \underline{b}_y^T + h_{\alpha,y} \underline{\gamma}_\alpha^T & \underline{0} \\ \underline{0} & \underline{0} & \underline{b}_z^T + h_{\alpha,z} \underline{\gamma}_\alpha^T \\ \underline{b}_y^T + h_{\alpha,y} \underline{\gamma}_\alpha^T & \underline{b}_x^T + h_{\alpha,x} \underline{\gamma}_\alpha^T & \underline{0} \\ \underline{0} & \underline{b}_z^T + h_{\alpha,z} \underline{\gamma}_\alpha^T & \underline{b}_y^T + h_{\alpha,y} \underline{\gamma}_\alpha^T \\ \underline{b}_z^T + h_{\alpha,z} \underline{\gamma}_\alpha^T & \underline{0} & \underline{b}_x^T + h_{\alpha,x} \underline{\gamma}_\alpha^T \end{bmatrix} \quad [7] $$

Notons que le développement de cette expression de $\underline{B}$ est basé, entre autres, sur les relations d'orthogonalité suivantes que l'on démontre au préalable :

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

Dans une approche standard en déplacement, la matrice de rigidité est déterminée simplement par la formule suivante :
$$ \underline{K}_e = \int_{\Omega_e} \underline{B}^T \cdot \underline{C} \cdot \underline{B} \, d\Omega \quad [9] $$

## 2.3. Analyse des modes de hourglass pour l'élément SHB6

Les modes de hourglass sont des modes non-physiques qui sont générés par l'intégration réduite. Ils produisent une énergie nulle alors qu'ils devraient induire 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.

Analysons maintenant le noyau de la matrice de rigidité obtenue par sous-intégration. Dans ce cas linéaire, la rigidité est obtenue par intégration de Gauss :

$$ \underline{K}_e = \sum_{i=1}^{5} \omega(\zeta_i) J(\zeta_i) \underline{B}^T(\zeta_i) \cdot \underline{C} \cdot \underline{B}(\zeta_i) \quad [10] $$

On voit ainsi que l'examen du noyau de la rigidité sous-intégrée revient à l'étude du rang de la matrice $\underline{B}$. En d'autres termes, il suffit de rechercher les modes de déplacement $\underline{X}$ à déformation nulle, c'est-à-dire vérifiant :

$$ \underline{B}(\zeta_i) \cdot \underline{X} = \underline{0} \quad \forall \zeta_i \quad [11] $$

Il est naturel de retrouver dans le noyau de la rigidité les modes de mouvements de corps rigides (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 rigides suivants :

$$ \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 [12] $$

On vérifie aisément que chacun des six vecteurs colonnes ci-dessus satisfait à l'équation [11] et appartient donc au noyau de la matrice de rigidité. Nous allons montrer que, outre les modes rigides précédents, il n'y a pas de modes qui annulent l'opérateur gradient discrétisé [7]. Construisons la base suivante de dix-huit vecteurs :

$$ \begin{bmatrix} \underline{s} & \underline{0} & \underline{0} & \underline{y} & \underline{z} & \underline{0} & \underline{x} & \underline{0} & \underline{0} & \underline{0} & \underline{0} & \underline{0} & \underline{h}_1 & \underline{0} & \underline{0} & \underline{h}_2 & \underline{0} & \underline{0} \\ \underline{0} & \underline{s} & \underline{0} & -\underline{x} & \underline{0} & \underline{z} & \underline{0} & \underline{x} & \underline{0} & \underline{y} & \underline{0} & \underline{0} & \underline{0} & \underline{h}_1 & \underline{0} & \underline{0} & \underline{h}_2 & \underline{0} \\ \underline{0} & \underline{0} & \underline{s} & \underline{0} & -\underline{x} & -\underline{y} & \underline{0} & \underline{0} & \underline{x} & \underline{0} & \underline{y} & \underline{z} & \underline{0} & \underline{0} & \underline{h}_1 & \underline{0} & \underline{0} & \underline{h}_2 \end{bmatrix} $$
On peut montrer que les vecteurs ci-dessus sont linéairement indépendants dans l'espace vectoriel des déplacements discrétisés de dimension dix-huit. Comme les douze derniers vecteurs ne vérifient pas l'équation [11], on montre alors qu'il n'y a pas d'autres modes que les modes rigides qui annulent l'opérateur gradient discrétisé donné en [7]. Autrement dit, l'élément SHB6 ne présente pas de modes de hourglass.

### 3. Résultats numériques de cas tests

Pour évaluer ce nouvel élément SHB6, nous l'avons testé sur un ensemble de 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 l'élément volumique existant PRI6. Dans la suite, on présente un seul cas test pour illustrer. C'est le cas test de la poutre vrillée encastrée et soumise à un effort tranchant à l'autre extrémité.

Ce cas test est devenu une référence pour tester les éléments de coques. 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 en termes de verrouillages.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Géométrie et chargement pour la poutre vrillée.</figcaption></figure>

<table><caption>**Tableau 1.** Déplacement normalisé selon Oz du point A de la poutre vrillée.</caption><thead><tr><th></th><th>PRI6</th><th>SHB6</th></tr><tr><th>Nombre d'éléments</th><th>Uz/Uréf</th><th>Uz/Uréf</th></tr></thead><tbody><tr><td>12x8x1</td><td>0,202</td><td>0,470</td></tr><tr><td>24x8x1</td><td>0,485</td><td>0,777</td></tr><tr><td>30x16x1</td><td>0,591</td><td>0,828</td></tr><tr><td>48x16x1</td><td>0,764</td><td>0,923</td></tr><tr><td>96x32x1</td><td>0,894</td><td>0,978</td></tr></tbody></table>

La figure 1 ci-dessus donne la géométrie de la poutre vrillée, les conditions aux limites et le chargement appliqué. L'extrémité gauche de cette poutre est encastrée et un effort tranchant P = 1N est appliqué à son extrémité droite. Les caractéristiques géométriques et matériau pour ce cas test sont données comme suit : longueur L = 12m, largeur l = 1,1m, épaisseur h = 0,32m, module d'Young E = 29.10⁶ Pa, coefficient de Poisson v = 0,22. La solution de référence pour ce problème est donnée dans plusieurs références de la littérature. Pour ce chargement P, la flèche du point A dans la direction Oz est égale à 5,42×10⁻³ m.

Le tableau 1 ci-dessus donne les résultats normalisés pour ce cas test en termes de convergence en fonction du nombre d'éléments dans le maillage. On voit que l'élément SHB6 converge bien et plus rapidement que l'élément PRI6.
## 4. Discussions et conclusions

Un nouvel élément de coque massive SHB6 a été formulé et implanté dans le code de calcul INCA. La version actuelle se comporte de façon relativement satisfaisante, puisqu'elle converge assez bien dans plusieurs cas tests programmés. D'ailleurs, les résultats montrent qu'il est plus performant que l'élément solide PRI6. Toutefois, l'élément SHB6 montre encore des verrouillages (membrane et cisaillement transverse) dans certaines configurations. Il est donc nécessaire de continuer l'effort sur l'amélioration de son opérateur gradient discrétisé (projections, 'assumed strain method', etc.) pour arriver à réduire les verrouillages rencontrés.

## Remerciements

Ces travaux sont financés par une bourse Cifre d'EDF R&D et un contrat industriel avec EDF R&D et le CETIM. Les auteurs tiennent à remercier J.M. Proix et X. Desroches (EDF R&D) et M. Afzali et J.F. Billaud du CETIM.

## 5. Bibliographie

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 adaptative, assumed-strain continuum mechanics shell element for impact analysis », *Computers & Structures*, vol. 80, 2002, p. 791-803.

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.

Hughes T., Taylor R., Kanoknukulchai W., « A simple and efficient finite element for plate bending », *Int. Journal Numerical Methods Engng.*, vol. 11, 1977, p. 1529-1543.

Legay A., Combescure A., « Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS », *Int. Jal. Num. Methods Engng.*, 2003, 57: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.
Pawsey S., Clough R., « Improved numerical integration of thick shell finite elements », *International Journal for Numerical Methods and Engineering*, vol. 3, 1971, p. 5575-86.

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.

Zienkiewicz O., Taylor R., Too J., « Reduced integration technique in general analysis of plates and shells », *International Journal for Numerical Methods Engng.*, 1971, 3:275-90.