1. Introduction

Il y a plusieurs manières de définir ce qu’est une plaque ou une coque. Nous allons adopter la définition suivante: «une plaque (ou une coque) est une structure mécanique, qui «vue de loin» est une surface».

En pratique, cela signifiera qu’une de ses dimensions (son épaisseur) est négligeable par rapport aux deux autres. De manière plus fondamentale, du point de vue mathématique, seule la courbure est un paramètre permettant d’appliquer ou non une théorie de coque, son épaisseur n’est pas une hypothèse fondatrice de la théorie.

La distinction entre plaque (« plate» ) et coque (« shell» ) repose sur l’idée que la plaque est une structure initialement plane alors que la coque est une structure fondamentalement courbe. Pour représenter un composant courbe (cylindre par exemple), les éléments de plaque procèdent donc par «facettisation» de la surface, ce qui introduit à la fois une erreur géométrique (approximation de la courbure) mais aussi mécanique puisque, pour les comportements isotropes, cela implique que les plaques ne présentent pas de couplage entre les déformations de membrane et de flexion.

De manière générale, une structure mince a une propriété remarquable: sa tenue mécanique (son équilibre) dépend essentiellement de sa géométrie. Il y a une très forte différence de rigidité entre une sollicitation dans sa surface moyenne (membrane) et normale à sa surface moyenne (flexion essentiellement). Une structure mince est en effet beaucoup plus rigide en membrane qu’en flexion (il est beaucoup plus facile de déformer une feuille en la pliant qu’en l’allongeant).

En pratique, cette propriété va induire beaucoup de problèmes de modélisation, une série de «verrouillages» numériques, c’est-à-dire une mauvaise représentation des différentes situations mécaniques qu’on ne peut pas toujours vaincre en raffinant le maillage. Le problème est fondamental et nous allons passer en revue les différentes situations présentant ces verrouillages et les moyens de les supprimer.

1.1. Limitation des éléments isoparamétriques classiques

Les éléments isoparamétriques 3D classiques présentent une série de verrouillages numériques, dus soit à la mauvaise représentation des différents modes de déformation et de contraintes représentatif des structures minces,soit à des propriétés du matériau ou du comportement (comme l’incompressibilité ou la plasticité).

Pour illustrer les situations de verrouillage venant de l’hypothèse des structure mince, nous allons considérer le calcul d’une simple plaque carrée en flexion, représentée par une seule couche d’éléments de type cube à 8 nœuds (avec donc une approximation tri-linéaire des déplacements).

	Propriétés de la structure Coté: \(1m\) Épaisseur: \(0,1m\) Module de Young: \(210\mathit{GPa}\) Coefficient de Poisson: \(0,3\)
Figure 1.1-1: Exemple de calcul de structure mince en 3D

Le déplacement vertical maximal que l’on devrait obtenir, issu de la théorie des plaques minces (résultat qu’on obtient avec les éléments DKT par exemple) est de 1,9 mm.

Sur le maillage 3D:

Nombre d’éléments	Déplacement
Théorie des plaques	\(1,9\mathit{mm}\)
3D – 25	\(9,7\times {10}^{-4}\mathit{mm}\)
3D – 2250	\(8,2\times {10}^{-3}\mathit{mm}\)
3D – 75000	\(6,3\times {10}^{-2}\mathit{mm}\)

On voit qu’on a un déplacement très faible (et faux) et qu’on ne converge que très lentement vers le bon déplacement. C’est ce qu’on appelle un «verrouillage».

1.1.1. Le verrouillage en cisaillement

La plaque précédente est considérée en flexion pure. Dans cette configuration, la solution d’un modèle poutre (c’est la même chose en plaque, on considère juste une indépendance par rapport à la coordonnée \(y\) pour simplifier l’exposé) est la suivante:

(1.1)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{\mathit{poutre}}=\frac{M}{\mathit{EI}}z\\ {\epsilon}_{zz}^{\mathit{poutre}}=-\frac{M\nu }{\mathit{EI}}z=-\nu {\epsilon}_{xx}^{\mathit{poutre}}\\ {\gamma}_{xz}^{\mathit{poutre}}=0\end{array}\end{split}\]

Si on suppose que \(\nu =0\) dans un premier temps:

(1.2)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{\mathit{poutre}}=\frac{M}{\mathit{EI}}z\\ {\epsilon}_{zz}^{\mathit{poutre}}=0\\ {\gamma}_{xz}^{\mathit{poutre}}=0\end{array}\end{split}\]

Considérons maintenant l’approximation linéaire en 3D. En flexion et en considérant l’indépendance par rapport à \(y\) et des fonctions linéaires complètes 3D dans le plan \((x,z)\) , les déplacements s’écrivent:

(1.3)\[\begin{split}\left\lbrace \begin{array}{c}{u}^{3D}\\ {w}^{3D}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{C}_{1}+{C}_{3}x+{C}_{5}z+{C}_{7}xz\\ {C}_{2}+{C}_{4}x+{C}_{6}z+{C}_{8}xz\end{array}\right\rbrace\end{split}\]

On calcule ainsi les trois déformations:

(1.4)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{3D}={C}_{3}+{C}_{7}z\\ {\epsilon}_{zz}^{3D}={C}_{6}+{C}_{8}x\\ {\gamma}_{xz}^{3D}={C}_{4}+{C}_{5}+{C}_{7}x+{C}_{8}z\end{array}\end{split}\]

En fait, les huit coefficients \({C}_{i}\) correspondent des modes de déformations particuliers comme résumé dans le tableau suivant.

	Rigide \(x\)	Rigide \(z\)	Membrane \(x\)	Membrane \(z\)	Cisaillement \(x\)	Cisaillement \(z\)	Flexion \(x\)	Flexion \(z\)
\({u}^{3D}\)	\(C\mathrm{₁}\)	0	\(C\mathrm{₃}x\)	0	0	\(C\mathrm{₅}z\)	\(C\mathrm{₇}xz\)	0
\({w}^{3D}\)	0	\(C\mathrm{₂}\)	0	\(C\mathrm{₆}z\)	\(C\mathrm{₄}x\)	0	0	\(C\mathrm{₈}xz\)
\({\epsilon}_{xx}^{3D}\)	0	0	\(C\mathrm{₃}\)	0	0	0	\(C\mathrm{₇}z\)	0
\({\epsilon}_{zz}^{3D}\)	0	0	0	\(C\mathrm{₆}\)	0	0	0	\(C\mathrm{₈}x\)
\({\gamma}_{xz}^{3D}\)	0	0	0	0	\(C\mathrm{₄}\)	\(C\mathrm{₅}\)	\(C\mathrm{₇}x\)	\(C\mathrm{₈}z\)

Les deux premières colonnes correspondent à des déformations nulles, ce sont doncdes mouvements de corps rigide (suivant \(x\) ou suivant \(z\) ). Les deux suivantes ont des déformations \({\epsilon}_{xx}^{3D}\) et \({\epsilon}_{zz}^{3D}\) constantes, ce sont des déformations de membrane pures, les deux suivantes correspondant à une déformation \({\gamma}_{xz}^{3D}\) non-nulle (cisaillement pur). Enfin les deux dernières correspondent à une déformation de flexion (avec cisaillement).

Si on considère le cas de la flexion pure (sans composante de membrane ou de cisaillement) les déformations tridimensionnelles sont:

(1.5)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{3D}={C}_{7}z\\ {\epsilon}_{zz}^{3D}=0\\ {\gamma}_{xz}^{3D}={C}_{7}x\end{array}\end{split}\]

Par rapport à la solution poutre (voir () ), on voit que la composante de cisaillement \({\gamma}_{xz}^{3D}\) n’est pas nulle en 3D, c’est ce qu’on appelle le verrouillage en cisaillement .

1.1.2. Le verrouillage en pincement

Il ne s’agit pas, à proprement parler, d’un verrouillage car la plupart des modèles de plaque que l’on utilise en calcul (Kirchhoff, Mindlin) font l’hypothèse de l’inextensibilité de le fibre transverse (l’épaisseur) \({\epsilon}_{zz}=0\) . Mais dans certains modèles de plaque/coque (modèle de Cosserat par exemple), cette hypothèse n’est pas retenue.

La conjonction de cette hypothèse avec celle usuellement pratiquée des contraintes planes conduit à une inconsistance théorique dans les formulations plaques. En effet l’élasticité de Hooke ne permet pas, en théorie, d’avoir simultanément \({\sigma}_{zz}=0\) et \({\epsilon}_{zz}=0\) . Néanmoins, cette contradiction apparente disparaît de manière indirecte par l’application de la formulation variationnelle, et, du point de vue mécanique, par le fait que l’hypothèse des contraintes planes est «excessive»: on peut négliger les contraintes \({\sigma}_{.z}\) par rapport aux autres. Le fait de les prendre nulles est une facilité pratique.

Pour les théories de plaque les plus générales, on rappelle qu’on a:

(1.6)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{\mathit{poutre}}=\frac{M}{\mathit{EI}}z\\ {\epsilon}_{zz}^{\mathit{poutre}}=-\frac{M\nu }{\mathit{EI}}z=-\nu {\epsilon}_{xx}^{\mathit{poutre}}\\ {\gamma}_{xz}^{\mathit{poutre}}=0\end{array}\end{split}\]

Et on voit dans l’équation () que l’expression de \({\epsilon}_{zz}^{\mathit{poutre}}\) n’est pas retrouvée en 3D et vaut:

(1.7)\[{\epsilon}_{zz}^{3D}={C}_{6}+{C}_{8}x\]

La déformation ne dépend que de \(x\) et elle est constante dans l’épaisseur (hypothèse de l’inextensibilité de la fibre normale). L’un des apports de l’élément COQUE_SOLIDE sera de corriger ce point et donc d’être une formulation plaque plus riche que les autres disponibles dans code_aster commeDKT, COQUE_3D, etc. .

1.1.3. Le verrouillage en membrane

Considérons maintenant une poutre courbe en flexion de longueur \(2l\) et de rayon \(R\) , on a:

(1.8)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}^{\mathit{poutre}}=\frac{\mathit{du}}{\mathit{ds}}+\frac{w}{R}\\ {\chi }^{\mathit{poutre}}=\frac{\mathit{du}}{\mathit{ds}}\frac{1}{R}-\frac{{d}^{2}w}{{\mathit{ds}}^{2}}\end{array}\end{split}\]

Compte-tenu des conditions de dérivabilité, on doit avoir les approximations suivantes:

(1.9)\[\begin{split}\lbrace \begin{array}{c}u={a}_{0}+{a}_{1}\xi \\ w={b}_{0}+{b}_{1}\xi +{b}_{2}{\xi}^{2}+{b}_{3}{\xi}^{3}\end{array}\text{avec}\xi =\frac{s}{l}\end{split}\]

Et donc:

(1.10)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}^{\mathit{poutre}}=\left(\frac{{a}_{1}}{l}+\frac{{b}_{0}}{R}\right)+\frac{{b}_{1}}{R}\xi +\frac{{b}_{2}}{R}{\xi}^{2}+\frac{{b}_{3}}{R}{\xi}^{3}\\ {\chi }^{\mathit{poutre}}=\left(\frac{{a}_{1}}{R}l-\frac{2{b}_{2}}{{l}^{2}}\right)-\frac{6{b}_{3}}{{l}^{2}}\xi \end{array}\end{split}\]

En flexion pure, on a évidemment la déformation de membrane qui est nulle:

(1.11)\[{\epsilon}^{\mathit{poutre}}=\left(\frac{{a}_{1}}{l}+\frac{{b}_{0}}{R}\right)+\frac{{b}_{1}}{R}\xi +\frac{{b}_{2}}{R}{\xi}^{2}+\frac{{b}_{3}}{R}{\xi}^{3}=0\]

Ce qui donne deux conditions:

(1.12)\[ \begin{align}\begin{aligned}\left(\frac{{a}_{1}}{l}+\frac{{b}_{0}}{R}\right)=0\\\textrm{et}\\\frac{{b}_{1}}{R}=\frac{{b}_{2}}{R}=\frac{{b}_{3}}{R}=0\end{aligned}\end{align} \]

La première condition est assez facile à réaliser par un choix judicieux des constantes \({a}_{1}\) et \({b}_{0}\) . La seconde implique que \(w\) soit une constante et donc:

(1.13)\[{\epsilon}^{\mathit{poutre}}=\frac{\mathit{du}}{\mathit{ds}}+\frac{w}{R}\ne 0\]

Ce qui veut dire que même en flexion pure, on a une contribution de membrane qu’on ne peut annuler. Il y a donc un verrouillage en membrane .

1.1.4. Le verrouillage en trapèze

Considérons un trapèze et son interpolation dans un espace de référence (carré de coté \(2\) ).

L’espace paramétrique a pour coordonnées \((\xi ,\zeta )\) . Le changement de coordonnées s’écrit:

\(\lbrace \begin{array}{c}x=\Lambda \xi (1-\alpha \zeta )\\ z=\zeta \end{array}\text{avec}\alpha \Lambda =\tan(\delta )\) et \(\lbrace \begin{array}{c}\xi =\frac{x}{\Lambda (1-\alpha z)}\\ \zeta =z\end{array}\)

Le jacobien de changement de coordonnées s’écrit:

(1.14)\[ \begin{align}\begin{aligned}\begin{split}\lbrace \begin{array}{c}\frac{\partial \xi }{\partial x}=\frac{1}{\Lambda (1-\alpha \zeta )}\\ \frac{\partial \xi }{\partial z}=\frac{\alpha \xi }{(1-\alpha \zeta )}\end{array}\end{split}\\\textrm{et}\\\begin{split}\lbrace \begin{array}{c}\frac{\partial \zeta }{\partial x}=0\\ \frac{\partial \zeta }{\partial z}=1\end{array}\end{split}\end{aligned}\end{align} \]

On remarquera que le jacobien n’est pas constant. Si on suppose pour simplifier que \(\frac{M}{\mathit{EI}}=1\) et \(\nu =0\) , la solution en déplacement d’une poutre en flexion pure s’écrit:

(1.15)\[\begin{split}\lbrace \begin{array}{c}{u}^{\mathit{poutre}}=xz\\ {w}^{\mathit{poutre}}=-\frac{1}{2}{x}^{2}\end{array}\end{split}\]

Soit dans l’espace paramétrique:

(1.16)\[\begin{split}\lbrace \begin{array}{c}{u}^{\mathit{poutre}}=\Lambda (\zeta \xi -\alpha \xi {\zeta}^{2})\\ {w}^{\mathit{poutre}}=-\frac{1}{2}{\Lambda}^{2}\left({\xi}^{2}+{\alpha}^{2}{\xi}^{2}{\zeta}^{2}-2\alpha \zeta {\xi}^{2}\right)\end{array}\end{split}\]

[Bib2] a introduit la notion d’aliasing des fonctions d’interpolation en éléments finis. Le principe est le suivant: dans un élément fini, on ne peut pas approximer toutes les fonctions, tout dépendra de l’ordre des polynômes utilisés dans ces fonctions. Par exemple, sur un triangle avec des fonctions linéaires, on ne peut pas représenter des fonctions quadratiques. Dans la représentation paramétrique habituelle, en 2D, on a le tableau suivant:

Élément	TRIA3	QUAD4	TRIA6	QUAD8
Fonction
\({\xi}^{2}\)	\(\xi\)	\(1\)	\({\xi}^{2}\)	\({\xi}^{2}\)
\(\xi \eta\)	\(0\)	\(\xi \eta\)	\(\xi \eta\)	\(\xi \eta\)
\({\eta}^{2}\)	\(\eta\)	\(1\)	\({\eta}^{2}\)	\({\eta}^{2}\)
\({\xi}^{3}\)	\(\xi\)	\(\xi\)	\(\frac{1}{2}\left(3{\xi}^{2}-\xi \right)\)	\(\xi\)
\({\xi}^{2}\eta\)	\(0\)	\(\eta\)	\(\frac{1}{2}\xi \eta\)	\({\xi}^{2}\eta\)
\(\xi {\eta}^{2}\)	\(0\)	\(\xi\)	\(\frac{1}{2}\xi \eta\)	\(\xi {\eta}^{2}\)
\({\eta}^{3}\)	\(\eta\)	\(\eta\)	\(\frac{1}{2}\left(3{\eta}^{2}-\eta \right)\)	\(\eta\)
\({\xi}^{4}\)	\(\xi\)	\(1\)	\(\frac{1}{4}\left(7{\xi}^{2}-3\xi \right)\)	\({\xi}^{2}\)
\({\xi}^{3}\eta\)	\(0\)	\(\xi \eta\)	\(\frac{\xi \eta }{4}\)	\(\xi \eta\)
\({\xi}^{2}{\eta}^{2}\)	\(0\)	\(1\)	\(\frac{\xi \eta }{4}\)	\({\xi}^{2}+{\eta}^{2}-1\)
\({\xi}^{2}\eta\)	\(0\)	\(\xi \eta\)	\(\frac{\xi \eta }{4}\)	\(\xi \eta\)
\({\eta}^{4}\)	\(\eta\)	\(1\)	\(\frac{1}{4}\left(7{\eta}^{2}-3\eta \right)\)	\({\eta}^{2}\)

Ce tableau se lit de la manière suivante: pour une fonction en \({\eta}^{3}\) , un élément TRIA3, QUAD4 ou QUAD8 représentera, au mieux, le cas \(\eta\) et pour le TRIA6, on pourra représenter par \(\frac{1}{2}\left(3{\eta}^{2}-\eta \right)\) . L’aliasing transforment donc les déplacements () ainsi pour le QUAD4:

(1.17)\[\begin{split}\lbrace \begin{array}{c}{u}^{\mathit{poutre}}=\Lambda (\zeta \xi -\alpha \xi {\zeta}^{2})\to {u}^{\mathit{EF}}=\Lambda (\zeta \xi -\alpha \xi )\\ {w}^{\mathit{poutre}}=-\frac{1}{2}{\Lambda}^{2}\left({\xi}^{2}+{\alpha}^{2}{\xi}^{2}{\zeta}^{2}-2\alpha \zeta {\xi}^{2}\right)\to {w}^{\mathit{EF}}=-\frac{1}{2}{\Lambda}^{2}\left(1+{\alpha}^{2}-2\alpha \zeta \right)\end{array}\end{split}\]

Ce qui donne les déformations suivantes:

(1.18)\[ \begin{align}\begin{aligned}{\epsilon}_{xx}=\frac{\partial u}{\partial x}=\frac{\partial u}{\partial \xi }\frac{\partial \xi }{\partial x}+\frac{\partial u}{\partial \zeta }\frac{\partial \zeta }{\partial x}=\frac{\zeta -\alpha }{1-\alpha \zeta }\\{\epsilon}_{zz}=\frac{\partial w}{\partial z}=\frac{\partial w}{\partial \xi }\frac{\partial \xi }{\partial z}+\frac{\partial w}{\partial \zeta }\frac{\partial \zeta }{\partial z}=\alpha {\Lambda}^{2}\\{\gamma}_{xz}=\frac{1}{2}\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)=\Lambda \xi \left(1+\frac{\alpha (\zeta -\alpha )}{1-\alpha \zeta }\right)\end{aligned}\end{align} \]

Or la solution analytique (voir () ) nous donne:

(1.19)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{xx}^{\mathit{poutre}}=z=\zeta \\ {\epsilon}_{zz}^{\mathit{poutre}}=0\\ {\gamma}_{xz}^{\mathit{poutre}}=0\end{array}\end{split}\]

Pour la contribution \({\epsilon}_{xx}\) , la solution élément fini est exacte pour \(\alpha =0\) (quadrangle parfait) et l’erreur reste petite pour un petit \(\alpha\) . L’erreur sur le cisaillement \({\gamma}_{xz}\) peut être annulée si on estime le terme en \((0,0)\) . Par contre le terme en \({\epsilon}_{zz}\) a une erreur plus grande que \({\gamma}_{xz}\) dès lors que \(\alpha \Lambda =\tan\delta >1\) . On a donc un verrouillage en trapèze de l’élément linéaire dès lors qu’on a des arêtes qui s’éloignent du cas géométrique du rectangle (ou de l’hexaèdre pour le cas 3D).

1.1.5. Le verrouillagevolumique (incompressibilité)

Ce type de verrouillage n’est pas spécifique aux structures minces et apparaît systématiquement lorsque le matériau ou le comportement devient incompressible, ce qui correspond en général à deux situations:

• quand le coefficient de Poisson tend vers \(0,5\) , ce qui est le cas des matériaux hyperélastiques comme les élastomères par exemple;

• quand le comportement est plastique, la plasticité des métaux se faisant à volume constant

Il s’agit d’un verrouillage volumique ou en incompressibilité. Pour le cas des éléments «massifs», il existe plusieurs stratégies pour vaincre ce verrouillage, voir [R3.06.08] ou [R3.06.14].

1.1.6. Méthodes de correction des verrouillages pour COQUE_SOLIDE

Pour vaincre ces verrouillages, on propose trois ingrédients:

• Pour supprimer le verrouillage volumique et de membrane, on va utiliser une méthode d’ intégration réduite (RI) qui consiste à utiliser un schéma d’intégration moins riche que celui utilise classiquement pour des briques à 8 nœuds. L’intégration réduite va faire apparaître des modes parasites à énergie nulle (ou modes «sabliers») qu’on va éliminer par une méthode de stabilisation;

• Pour supprimer le verrouillage volumique, on propose en outre une approche basée sur les modes incompatibles dans sa formulation générique dite « EAS » (Enhanced Assumed Strain) qui permettra d’enrichir l’approximation pour corriger le problème du pincement;

• Pour supprimer le verrouillage en cisaillement transverse et trapézoïdal, on utilisera la méthode du champ assumé ( ANS «Assumed Natural Strain»).

La méthode EAS (Enhanced Assumed Strain), dont le formalisme est basé sur la formulation à trois champs de Hu-Washizu, a été initialement proposée par [Bib2, Bib3] en petites déformations, puis étendue aux grandes déformations par [Bib4,Bib5] avant d’être appliqué aux problèmes thermomécaniques par [Bib6]. Cette méthode est une généralisation de la méthode des modes incompatibles. Elle consiste à réduire l’ordre des polynômes pour certaines composantes des déformations.

Cette méthode permet de réduire significativement le verrouillage volumique. Des travaux récents [Bib8], [Bib9] ont montré qu’un seul paramètre suffisait, ce qui permet en plus d’être cohérent avec une intégration réduite. Néanmoins, dans certaines configurations (chargement en compression), des instabilités apparaissent, de plus il est nécessaire dans le cas des éléments plaques de vaincre aussi les verrouillages en cisaillement transverse. On va combiner la méthode EAS avec la méthode ANS et la prise en compte du pincement.

La technique ANS («Assumed Natural Strain») consiste à remplacer la déformation de cisaillement transverse par une déformation de cisaillement postulée obtenue par interpolation dans la base naturelle de l’élément (d’où le qualificatif de «naturel») . En plus de faciliter la gestion du verrouillage en cisaillement, l’interpolation dans la base naturelle donne des éléments fonctionnant avec un maillage plus ou moins déformé.