# Élément de coque 3D basé sur une formulation de déformation améliorée postulée

Philippe Jetteur, Philippe Pasquet

► To cite this version:

Philippe Jetteur, Philippe Pasquet. Élément de coque 3D basé sur une formulation de déformation améliorée postulée. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422179

HAL Id: hal-01422179

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

Submitted on 23 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Élément de coque 3D basé sur une formulation de déformation améliorée postulée

Ph. Jetteur, Ph. Pasquet

SAMTECH

Rue des Chasseurs Ardennais 8, B-4031Liège, Belgique

Philippe.Jetteur@samcef.com; Philippe.Pasquet@samcef.com

http://www.samcef.com

**Résumé** : Dans le but de simplifier le maillage en partant directement d'une géométrie 3D (ie sans extraire la surface moyenne), un nouvel élément de coque 3D robuste a été développé dans SAMCEF [1].

L'élément est de forme hexaédrique à 8 noeuds. Pour améliorer la qualité des résultats, nous ajoutons des degrés de libertés (ddl) internes proposés par Simo et al. [2]. Nous utilisons la « méthode de déformation améliorée postulée », Enhanced Assumed Strain (EAS) en anglais. Un traitement spécial des cisaillements transversaux est exécuté pour passer avec succès le test de rapiéçage (flexion constante) et éviter le verrouillage en cisaillement. La formulation est basée sur les travaux de Bathe et Dvorkin [3].

Comme les maillages industriels comportent souvent des triangles, nous avons aussi développé un élément prismatique à 6 nœuds sans ddl internes.

Nous introduisons des hypothèses de contraintes planes dans la loi de comportement pour être capable de représenter de grands rapports longueur/épaisseur : pour un élément carré soumis à une pression transverse, le rapport longueur sur épaisseur peut atteindre 100 000.

L'élément est développé pour les analyses linéaires et non linéaires. Il peut être utilisé en configuration mono ou multi-couche.

Il est d'abord validé sur des tests classiques [4] puis comparé à d'autres éléments de coques [5], [6] en termes de taux de convergence.

**Mots-Clés** : Coque, Elément 3D, Flambement, Formulation EAS

## 1. INTRODUCTION

La formulation consiste à ajouter des modes internes pour enrichir les déformations ou le gradient de déplacement. Elle concerne les briques à 8 nœuds. Elle a d'abord été introduite dans les éléments volumiques existant. Un élément spécialement conçu pour les coques, que nous décrivons et utilisons dans la suite, a ensuite été développé. Une version triangulaire (prisme à 6 nœuds) a été également développée. Le quadrangle et le triangle sont libres de verrouillage en cisaillement.

## 2. FORMULATION

Le point de départ est un principe variationnel à 3 champs où les positions déformées (x), le gradient de déplacement ($\tilde{F}$) et les contraintes de PK1 ($\tilde{P}$) sont discrétisés :

$$ \Pi(x, \tilde{F}, \tilde{P}) = \int [W(\tilde{F}) + \tilde{P}(\text{grad } x - \tilde{F})] dV + \Pi_{\text{ext}} $$

où le dernier terme est associé aux forces extérieures. Le gradient de déplacement est décomposé en un terme qui dépend des déplacements ($F_{com}$) et un terme amélioré ($F_{enh}$), « enhanced » : $\tilde{F} = F_{com} + F_{enh} = \text{grad } x + F_{enh}$. Le potentiel devient :

$$ \Pi(x, \tilde{F}, \tilde{P}) = \int [W(\tilde{F}) - \tilde{P}F_{enh}] dV + \Pi_{\text{ext}} $$
On choisit $\tilde{P}$ constant sur l'élément et $F_{enh}$ tel que $\int F_{enh}dV=0$ (condition d'orthogonalité), le potentiel se résume à :

$$ \Pi(x, \tilde{F}, \tilde{P}) = \int [W(\tilde{F})] dV + \Pi_{ext} $$

On prend : $F_{enh} = F_0 F'(\xi, \eta, \delta)$, $F_0 = F_{com}|_{\xi=\eta=\delta=0}$ où $F'$ est défini par :

$$ F' = \frac{j_0}{j_\xi} J_0 F''(\xi, \eta, \delta) J_0^{-1} \qquad J_0 = \begin{vmatrix} \frac{\partial X}{\partial \xi} & \frac{\partial X}{\partial \eta} & \frac{\partial X}{\partial \delta} \\ \frac{\partial Y}{\partial \xi} & \frac{\partial Y}{\partial \eta} & \frac{\partial Y}{\partial \delta} \\ \frac{\partial Z}{\partial \xi} & \frac{\partial Z}{\partial \eta} & \frac{\partial Z}{\partial \delta} \end{vmatrix} \quad j_0 = \det(J_0) $$

où $F''$ dépend des coordonnées intrinsèques et $J_0$ représente la matrice jacobienne au centre de l'élément. Si l'intégrale de $F''$ sur le volume unitaire est nulle, la condition d'orthogonalité sera vérifiée. Cette définition de $F'$ est proposée dans [2].

$F''$ est choisi de telle manière à compléter les champs manquant dans la partie conforme des déformations. Les 9 premiers termes permettent d'améliorer le comportement volumétrique de l'élément, les 15 suivants le comportement en cisaillement (verrouillage en volume et verrouillage en cisaillement). On prend :

$$ F'' = \begin{vmatrix} \alpha_1 \xi + \alpha_4 \xi \eta + \alpha_7 \xi \delta & \alpha_{13} \eta + \alpha_{19} \eta \xi + \alpha_{22} \eta \delta & \alpha_{14} \delta + \alpha_{20} \delta \xi + \alpha_{23} \delta \eta \\ \alpha_{10} \xi + \alpha_{16} \xi \delta & \alpha_2 \eta + \alpha_5 \eta \xi + \alpha_8 \eta \delta & \alpha_{15} \delta + \alpha_{21} \delta \xi + \alpha_{24} \delta \eta \\ \alpha_{11} \xi + \alpha_{17} \xi \eta & \alpha_{12} \eta + \alpha_{18} \eta \xi & \alpha_3 \delta + \alpha_6 \delta \xi + \alpha_9 \delta \eta \end{vmatrix} $$

## 2.1 Version introduite dans les volumes de SAMCEF

Une fois que le gradient de déplacement ou les déformations de Green sont définies, le développement de l'élément est classique. Il a été introduit dans les volumes mono et multicouches (types 8 et 11 dans SAMCEF), a priori les 24 modes sont gardés.

Dans SAMCEF, les deux versions (gradient de déplacement ou déformations) ont été introduites. Seule la version sur les déformations est accessible pour l'utilisateur.

## 2.2 Version coque, approche 1

Les développements décrits à la section précédente concernent le volume. Si on analyse des problèmes de coque avec cet élément, il donne des bons résultats sans devoir faire des hypothèses d'état plan de contrainte dans la loi constitutive.

Néanmoins, comme le montrent les exemples, si la coque devient trop mince, des problèmes numériques apparaissent lors de la résolution du système d'équations. Ces problèmes numériques ne sont pas dus aux termes EAS, ils apparaissent aussi avec l'élément standard.

Pour fixer les idées, on considère un cas 2D avec un élément rectangulaire de dimension $L$ selon $X$ et $t$ selon $Z$, de modules de Young $E_x$ et $E_z$ selon les directions $X$ et $Z$, de coefficient de Poisson nul et de module de cisaillement $G$. Pour avoir un ordre de grandeur des différents termes de raideurs, on calcule la raideur associée au premier nœud, en utilisant un schéma d'intégration à un point :

$$ U = \frac{1}{2} \int [E_x \varepsilon_x^2 + E_z \varepsilon_z^2 + G \gamma^2] dA \qquad \begin{aligned} \varepsilon_x &= u_1/2L \\ \varepsilon_z &= w_1/2t \\ \gamma &= u_1/2t + w_1/2L \end{aligned} $$

La matrice de raideur associée au premier nœud vaut :

$$ K = \frac{Lt}{4L^2} \begin{vmatrix} E_x + \alpha^2 G & \alpha G \\ \alpha G & G + \alpha^2 E_z \end{vmatrix} \qquad \alpha = \frac{L}{t} $$
Les termes en $\alpha^2$ peuvent créer des problèmes numériques quand l’élément devient mince.

Dans le terme (1,1), cela conduit à du verrouillage en cisaillement, éliminé par une intégration adéquate ou des termes supplémentaires. Une possibilité, non utilisée ici, est aussi d'introduire la notion de « residual bending flexibility » (introduite par Mac Neal) et de réduire le module de cisaillement quand l'élément est mince, étant entendu que, si l'élément est mince, le cisaillement a peu d'importance. Le terme (2,2) est associé à la contrainte transverse qui n'est pas considérée dans les problèmes de coques. Dans l'élément, on prend une loi constitutive en état plan de contrainte et on prend selon la direction normale un module de Young divisé par $max(1, \alpha^2)$, ce rapport est calculé automatiquement par le programme.

La deuxième modification due à l'hypothèse coque est le choix des modes EAS. La déformation transversale étant découplée des autres composantes, il est inutile de l'enrichir. Le cisaillement transverse peut être enrichi uniquement avec des termes constants sur l'épaisseur. Dans le cas des coques, on se limite à 12 modes EAS :

$$F'' = \begin{vmatrix} \alpha_1 \xi + \alpha_7 \xi \delta & \alpha_4 \eta + \alpha_{10} \eta \delta & 0 \\ \alpha_2 \xi + \alpha_8 \xi \delta & \alpha_5 \eta + \alpha_{11} \eta \delta & 0 \\ \alpha_3 \xi + \alpha_9 \xi \eta & \alpha_6 \eta + \alpha_{12} \eta \xi & 0 \end{vmatrix}$$

Cette modification est introduite dans les types 8 et 11 de Samcef, option coque.

## 2.3 Version coque, approche 2

### 2.3.1 Triangle

Quelques triangles apparaissent souvent dans les maillages industriels. Dans le cas d'un élément triangulaire du premier degré, les états membranaires et flexionnels sont constants sur la surface de l'élément. Il paraît compliqué de les enrichir par des modes EAS, qui ne seraient pas couplés aux déformations déduites des déplacements et ne serviraient à rien.

Si rien n'est fait, ces éléments triangulaires du premier degré sont beaucoup trop rigides en cisaillement transverse (verrouillage en cisaillement). La procédure adoptée est :

- calcul du cisaillement transverse au milieu des 3 cotés, on garde la composante le long du coté

- en chaque nœud, calcul du cisaillement dans des axes xy locaux à l'élément, à partir des 2 cisaillements des cotés arrivant à ce nœud

- le cisaillement en surface est pris comme la moyenne du cisaillement aux 3 nœuds

Les résultats obtenus avec cette procédure sont bons comme le montrent les exemples traités.

### 2.3.2 Quadrangle

Comme la procédure pour passer du cisaillement des cotés au cisaillement en volume a été introduite pour le triangle, on en tire parti pour le quadrangle. On adopte la procédure proposée par [3] :

- on a une base covariante, dont les 2 premiers axes sont parallèles aux coordonnées intrinsèques et dont le troisième est la normale

- on calcule $\gamma_{\xi3}$ $\gamma_{\eta3}$ aux milieux des cotés de l'élément

$$\gamma_{\xi3} = \frac{1-\eta}{2} \gamma_{\xi3}^{I} + \frac{1+\eta}{2} \gamma_{\xi3}^{III}$$

- on interpole le cisaillement en surface :

$$\gamma_{\eta3} = \frac{1+\xi}{2} \gamma_{\eta3}^{II} + \frac{1-\xi}{2} \gamma_{\eta3}^{IV}$$

- on passe ensuite du cisaillement exprimé dans la base covariante au cisaillement exprimé dans la base cartésienne parallèle à la surface moyenne

Les 12 modes EAS sont gardés. De cette manière, l'élément passe exactement un patch test en flexion et la convergence (en nombre d'éléments) est améliorée quand l'élément n'est pas rectangulaire.
## 2.4 Intégration selon l'épaisseur

Dans le cas d'une coque multicouche, une première solution est de mettre des points sur chaque couche, c'est couteux en terme de temps de calcul. En pratique, par point d'intégration de surface, on met deux points d'intégration sur l'épaisseur, on passe des déformations en volume en déformations de coque plus un pincement (variation de l'épaisseur), on calcule les efforts normaux et moments de flexion et la matrice associée, on repasse ensuite aux contraintes en volume.

Si, par simplicité, on néglige la direction y, on a pour les déformations :

$$ 
\begin{vmatrix} 
\epsilon_0 \\ 
\chi \\ 
\gamma \\ 
\epsilon_z 
\end{vmatrix} 
= 
\begin{vmatrix} 
1/2 & 0 & 0 & 1/2 & 0 & 0 \\ 
1/2t^* & 0 & 0 & -1/2t^* & 0 & 0 \\ 
0 & 0 & 1/2 & 0 & 0 & 1/2 \\ 
0 & 1/2 & 0 & 0 & 1/2 & 0 
\end{vmatrix} 
\begin{vmatrix} 
\epsilon_x^1 \\ 
\epsilon_z^1 \\ 
\gamma^1 \\ 
\epsilon_x^2 \\ 
\epsilon_z^2 \\ 
\gamma^2 
\end{vmatrix} 
\quad 
\begin{array}{l} 
\epsilon_{coque} = C \epsilon_{volume} \\ 
t^* = \frac{t}{2\sqrt{3}} 
\end{array} 
$$

Dans une coque monocouche, les efforts de coques se calculent par :

$$ 
\begin{vmatrix} 
N \\ 
M \\ 
T \\ 
\sigma_z 
\end{vmatrix} 
= 
\begin{vmatrix} 
Et \\ 
Et^3/12 \\ 
Gt \\ 
\alpha^{-2} E 
\end{vmatrix} 
\begin{vmatrix} 
\epsilon_0 \\ 
\chi \\ 
\gamma \\ 
\epsilon_z 
\end{vmatrix} 
\quad 
\sigma_{coque} = H_{coque} \epsilon_{coque} 
$$

Si $dz (= t/2)$ est le poids associé à l'intégration selon l'épaisseur, les contraintes en volume se calculent par :

$$ \langle \epsilon_1^T \quad \epsilon_2^T \rangle \begin{vmatrix} \sigma_1 \\ \sigma_2 \end{vmatrix} dz = \epsilon_{coque}^T \sigma_{coque} \quad \begin{vmatrix} \sigma_1 \\ \sigma_2 \end{vmatrix} = \frac{1}{dz} C^T \sigma_{coque} $$

La raideur est calculée à l'aide d'une intégrale de surface en utilisant $H_{coque}$, les contraintes locales ne sont utilisées que pour calculer la raideur géométrique. Si les matériaux sont linéaires, $H_{coque}$, est calculé de manière analytique.

La procédure décrite est utilisée aussi bien en version mono que multicouche. En version monocouche, si le matériau est non-linéaire, il y a 5 points sur l'épaisseur pour intégrer la loi constitutive. Avec cette procédure, le coût du calcul de la matrice de raideur de l'élément ne dépend ni du nombre de couches ni du nombre de points sur l'épaisseur.

# 3. RÉSULTATS NUMÉRIQUES

Parmi les nombreux tests effectués, nous en présentons trois en linéaire et un en non linéaire.

## 3.1 Plaque carrée uniformément chargée sur appuis simples

On considère une plaque carrée d'épaisseur t et de coté L. Nous présentons les résultats pour les éléments quadrangulaires et triangulaires en fonction de l'élancement L/t, 8*8 éléments sont utilisés pour modéliser ¼ de plaque. Les résultats sont comparés avec ceux obtenus dans [6] et ceux obtenus par SAMCEF (en italiques) avec l'ancienne coque volumique. Le tableau suivant donne le rapport entre le déplacement trouvé et le déplacement théorique.

<table><thead><tr><th>L/t</th><th>SAMCEF</th><th>SAMCEF quad</th><th>SAMCEF tria</th><th>[6]</th></tr></thead><tbody><tr><td>10</td><td></td><td></td><td></td><td>1,004</td></tr><tr><td>100</td><td></td><td></td><td></td><td>1,008</td></tr><tr><td>200</td><td>0,919</td><td>1</td><td>0,993</td><td>1,004</td></tr><tr><td>1000</td><td></td><td></td><td></td><td></td></tr><tr><td>2000</td><td>0,102</td><td>1</td><td>1,005</td><td></td></tr></tbody></table>
<table><tr><td>20000</td><td>0,001</td><td>1</td><td>1,017</td></tr><tr><td>200000</td><td></td><td>1</td><td>1,016</td></tr></table>

On remarque que :

- le fait d'introduire l'hypothèse coque dans la formulation EAS permet d'augmenter l'élancement

- le triangle se comporte bien.

## 3.2 Poutre vrillée

Il s'agit d'un test classique ([4]). Les principales caractéristiques sont données dans le tableau suivant. Deux cas de charge sont étudiés (P1 et P2). Le maillage est constitué de N éléments sur la largeur de la poutre et de 6N éléments sur la longueur. On utilise des quadrangles et on compare les résultats à ceux de [6]. Les tableaux suivants donnent le rapport entre le déplacement trouvé et le déplacement théorique (au point A).

<table><tr><td>L</td><td>12</td></tr><tr><td>I</td><td>1,1</td></tr><tr><td>h</td><td>0,32</td></tr><tr><td>E</td><td>29 10<sup>6</sup></td></tr><tr><td>v</td><td>0,22</td></tr><tr><td>torsion</td><td>90°</td></tr></table>

A

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

<table><tr><td>P1//z</td><td>1</td></tr><tr><td>W<sub>max</sub></td><td>5,424 10<sup>-3</sup></td></tr></table>

<table><thead><tr><th>N</th><th>SAMCEF</th><th>[6]</th></tr></thead><tbody><tr><td>2</td><td>0,987</td><td>0,999</td></tr><tr><td>4</td><td>0,990</td><td>0,999</td></tr><tr><td>8</td><td>0,998</td><td></td></tr></tbody></table>

<table><tr><td>P2//y</td><td>1</td></tr><tr><td>V<sub>max</sub></td><td>1,754 10<sup>-3</sup></td></tr></table>

<table><thead><tr><th>N</th><th>SAMCEF</th><th>[6]</th></tr></thead><tbody><tr><td>2</td><td>0,979</td><td>0,994</td></tr><tr><td>4</td><td>0,990</td><td>0,998</td></tr><tr><td>8</td><td>0,998</td><td></td></tr></tbody></table>

## 3.3 Cylindre pincé – diaphragme rigide

Il s'agit d'un test classique de coque repris par de nombreux auteurs ([5], [6]) : il met en évidence la prépondérance du verrouillage en cisaillement sur le verrouillage en membrane. Les données du problème sont :

<table><tr><td>L</td><td>600</td></tr><tr><td>R</td><td>300</td></tr><tr><td>t</td><td>3</td></tr><tr><td>E</td><td>3 10<sup>6</sup></td></tr><tr><td>v</td><td>0,3</td></tr><tr><td>P</td><td>1</td></tr><tr><td>W<sub>max</sub></td><td>-18.248 10<sup>-6</sup></td></tr></table>

Diaphragme rigide

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

Le tableau suivant donne le rapport entre le déplacement trouvé et le déplacement théorique en fonction du maillage (NxN) de 1/8e de cylindre. Les résultats sont comparés à ceux obtenus dans [5] & [6] et ceux obtenus par SAMCEF (en italiques) avec le volume et les 24 modes EAS.

<table><thead><tr><th>N</th><th>SAMCEF</th><th>SAMCEF</th><th>[5]</th><th>[6]</th></tr></thead><tbody><tr><td>4</td><td>0,107</td><td>0,357</td><td>0,104</td><td>0,387</td></tr><tr><td>8</td><td>0,497</td><td>0,744</td><td>0,494</td><td>0,754</td></tr><tr><td>16</td><td>0,912</td><td>0,934</td><td>0,912</td><td>0,940</td></tr><tr><td>32</td><td>0,990</td><td>0,991</td><td>0,995</td><td>0,997</td></tr></tbody></table>

La convergence est assez lente mais est conforme aux résultats publiés.
## 3.4 Carré cylindrique encastré sous pression – non linéaire

La structure étudiée est un carré cylindrique. Le ¼ de structure est maillée avec 8*8 éléments version quadrangulaire et triangulaire. On compare la solution avec la coque Mindlin utilisant le même maillage.

La figure suivante représente le déplacement transversal du point central en fonction de la pression.

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

## 4. CONCLUSION

Les éléments volumiques sont enrichis par des modes EAS. La nouvelle coque se comporte bien, tant en option quadrangulaire que triangulaire, qu'en linéaire ou non, qu'en mono ou multicouche.

## RÉFÉRENCES

[1] http://www.samcef.com

[2] Simo J.C., Armero F. and Taylor R.L.: Improved versions of assumed enhanced strain tri-linear element for 3D finite deformation problem. *Comput. Methods Appl. Mech Engng*, vol 110. pp 359-386 (1993)

[3] Dvorkin E.N. and Bathe K. J.: A continuum mechanics based four-node shell element for general nonlinear analysis. *Engineering Computations*, Vol 1, pp 77-88 (1984)

[4] Mac Neal R.H. and Harder R.L.: A proposed standard set of problems to test finite element accuracy. *Finite Elements in Analysis and Design*, 1, pp 3-20 (1985)

[5] Alves de Sousa R.J., Natal Jorge R.M., Fontes Valente R.A., César de Sà J.M.A.: A new volumetric and shear locking-free 3D enhanced strain element. *Engineering Computations*, Vol. 20 N°7, pp. 896-925 (2003)

[6] Abed-Meraim F., Combescure A.: A physically stabilized and locking-free formulation of the (SHB8PS) solid-shell element. *Revue Européenne de Mécanique Numérique*, Vol. 16 N°8, pp 1037-1072 (2007)