4. Stabilisation

4.1. Principes

Dans le calcul des quantités nécessaires, il apparaît des termes intégraux venant de la formulation variationnelle (voir § 3 ). La matrice de rigidité de la structure s’écrira de la manière suivante pour un élément fini \(e\) donné:

(4.1)\[{\int}_{{\Omega}_{e}}\left(\epsilon :\delta \sigma \right)d\Omega \to {K}_{e}={\int}_{{\Omega}_{e}}\left({B}^{T}CB\right)d\Omega\]

Avec \(C\) la matrice d’élasticité de Hooke et \(B\) la matrice des déformations linéarisées. On doit donc utiliser un schéma de quadrature numérique pour évaluer l’intégrale. Ce schéma contiendra \({n}_{\mathit{pt}}\) points de coordonnées \({\xi}_{p=1,{n}_{\mathit{pt}}}\) , chacun d’un poids \({\omega}_{p=1,{n}_{\mathit{pt}}}\) tels que:

(4.2)\[{K}_{e}={\int}_{{\Omega}_{e}}\left({B}^{T}CB\right)d\Omega =\sum_{p=1}^{{n}_{\mathit{pt}}}{\omega}_{p}({\xi}_{p})J({\xi}_{p}){B({\xi}_{p})}^{T}CB({\xi}_{p})\]

\(J\) est le jacobien de transformationde l’élément entre l’espace paramétrique et l’espace cartésien (voir § 2.2 ). Le choix de ce schéma, sa précision en particulier, va dépendre de la nature de l’intégrande. Pour l’élément COQUE_SOLIDE, on va utiliser un schéma de type Lobatto à 5 points.

Ce schéma est trop «pauvre» et il va sous-estimer la rigidité de l’élément avec pour objectif de contrôler les verrouillages. Mais en procédant ainsi, on fait également apparaître des modes de déformations qu’on appelle «modes sabliers». L’élément risque de changer de forme sans dépenser d’énergie, ce qui ne correspond pas à la physique. Pour minimiser ce problème, on va stabiliser l’élément par une technologie adéquate.

Le but de la procédure de stabilisation est de corriger le défaut de rang de la matrice de rigidité provenant du schéma d’intégration réduit adopté. Le lecteur est renvoyé aux importants travaux sur ce sujet de [bib18], [bib19], [bib13], [bib14], [bib8], [bib9], [bib10]. Dans certaines formulations de coque-solides, la matrice jacobienne et son inverse ne sont évaluées qu’au centre de l’élément, voir par exemple [bib21], [bib22], [bib23]. Une telle approximation suppose que l’élément réel peut être représenté par son parallélépipède équivalent. Et c’est assez précis pour les maillages fins et peu déformés.

Cependant pour les maillages qui sont fortement déformés à l’état initial, cette hypothèse montre un manque de précision [bib10]. Afin de prendre en considération la forme réelle de l’élément et de stabiliser correctement la matrice de rigidité, une décomposition polynomiale de la matrice jacobienne inverse a étéconstruite explicitement.

4.1.1. Stabilisation des déformations

Par ailleurs, pour le passage entre des quantités comme la déformation entre les coordonnées paramétriques et les coordonnées réelles, donné par la matrice \(T\) , il convient d’utiliser une décomposition polynomiale cohérente avec celle choisie pour le jacobien en particulier. Il a été montré dans [bib10] qu’une très bonne précision est atteinte en utilisant simplement les termes constants et linéaires de \(T\) . Par conséquent, les termes bilinéaires sont ignorés:

(4.3)\[T\approx {T}^{0}+\xi {T}^{\xi}+\eta {T}^{\eta}+\zeta {T}^{\zeta}\]

On va utiliser la même forme d’approximation sur \({E}^{u}\) que sur \({\widehat{E}}^{u}\) :

(4.4)\[{E}^{u}=T{\widehat{E}}^{u}\approx {E}^{u,0}+\zeta {E}^{u,\zeta }+{\zeta}^{2}{E}^{u,\zeta \zeta }+\xi {E}^{u,\xi }+\eta {E}^{u,\eta }+\xi \eta {E}^{u,\xi \eta }+\xi \zeta {E}^{u,\xi \zeta }+\eta \zeta {E}^{u,\eta \zeta }\]

Que l’on va diviser en deux composantes:

(4.5)\[ \begin{align}\begin{aligned}{E}^{u}={E}^{u,\mathit{RI}}+{E}^{u,\mathit{STAB}}\\{E}^{u,\mathit{RI}}={E}^{u,0}+\zeta {E}^{u,\zeta }+{\zeta}^{2}{E}^{u,\zeta \zeta }\\{E}^{u,\mathit{STAB}}=\xi {E}^{u,\xi }+\eta {E}^{u,\eta }+\xi \eta {E}^{u,\xi \eta }+\xi \zeta {E}^{u,\xi \zeta }+\eta \zeta {E}^{u,\eta \zeta }\end{aligned}\end{align} \]

Avec les expressions explicites suivantes:

(4.6)\[ \begin{align}\begin{aligned}{E}^{u,0}={T}^{0}{\widehat{E}}^{u,0}\\{E}^{u,\zeta }={T}^{\zeta}{\widehat{E}}^{u,0}+{T}^{0}{\widehat{E}}^{u,\zeta }\\{E}^{u,\zeta \zeta }={T}^{\zeta}{\widehat{E}}^{u,\zeta }+{T}^{0}{\widehat{E}}^{u,\zeta \zeta }\\{E}^{u,\xi }={T}^{0}{\widehat{E}}^{u,\xi }+{T}^{\xi}{\widehat{E}}^{u,0}\\{E}^{u,\eta }={T}^{0}{\widehat{E}}^{u,\eta }+{T}^{\eta}{\widehat{E}}^{u,0}\\{E}^{u,\xi \eta }={T}^{0}{\widehat{E}}^{u,\xi \eta }+{T}^{\xi}{\widehat{E}}^{u,\eta }+{T}^{\eta}{\widehat{E}}^{u,\xi }\\{E}^{u,\eta \zeta }={T}^{0}{\widehat{E}}^{u,\eta \zeta }+{T}^{\zeta}{\widehat{E}}^{u,\eta }+{T}^{\eta}{\widehat{E}}^{u,\zeta }\\{E}^{u,\xi \zeta }={T}^{0}{\widehat{E}}^{u,\xi \zeta }+{T}^{\xi}{\widehat{E}}^{u,\zeta }+{T}^{\zeta}{\widehat{E}}^{u,\xi }\end{aligned}\end{align} \]

On utilisera la même décomposition pour le tenseur des déformations linéarisées et donc pour la matrice gradient:

(4.7)\[ \begin{align}\begin{aligned}\widehat{B}={\widehat{B}}^{\mathit{RI}}+{\widehat{B}}^{\mathit{STAB}}\\\textrm{et}\\B=T\widehat{B}\\{B}^{\mathit{RI}}={B}^{0}+\zeta {B}^{\zeta}+{\zeta}^{2}{B}^{\zeta \zeta }\\{B}^{\mathit{STAB}}=\xi {B}^{\xi}+\eta {B}^{\eta}+\xi \eta {B}^{\xi \eta }+\xi \zeta {B}^{\xi \zeta }+\eta \zeta {B}^{\eta \zeta }\end{aligned}\end{align} \]

Avec:

(4.8)\[ \begin{align}\begin{aligned}{B}^{0}={T}^{0}{\widehat{B}}^{0}\\{B}^{\zeta}={T}^{\zeta}{\widehat{B}}^{0}+{T}^{0}{\widehat{B}}^{\zeta}\\{B}^{\zeta \zeta }={T}^{\zeta}{\widehat{B}}^{\zeta}+{T}^{0}{\widehat{B}}^{\zeta \zeta }\\{B}^{\xi}={T}^{0}{\widehat{B}}^{\xi}+{T}^{\xi}{\widehat{B}}^{0}\\{B}^{\eta}={T}^{0}{\widehat{B}}^{\eta}+{T}^{\eta}{\widehat{B}}^{0}\\{B}^{\xi \eta }={T}^{0}{\widehat{B}}^{\xi \eta }+{T}^{\xi}{\widehat{B}}^{\eta}+{T}^{\eta}{\widehat{B}}^{\xi}\\{B}^{\eta \zeta }={T}^{0}{\widehat{B}}^{\eta \zeta }+{T}^{\zeta}{\widehat{B}}^{\eta}+{T}^{\eta}{\widehat{B}}^{\zeta}\\{B}^{\xi \zeta }={T}^{0}{\widehat{B}}^{\xi \zeta }+{T}^{\xi}{\widehat{B}}^{\zeta}+{T}^{\zeta}{\widehat{B}}^{\xi}\end{aligned}\end{align} \]

\({E}^{u,\mathit{STAB}}\) et \({B}^{\mathit{STAB}}\) représentent les termes annulés dans l’intégration réduite. Le processus de stabilisation consiste à restituer analytiquement ces termes dans la matrice de rigidité.

De plus, pour éliminer le verrouillage volumétrique, l’approche B-bar [bib24] est adoptée. Par conséquent, les modes en sablier des opérateurs de déformation et de déplacement de déformation stabilisée est divisée en composants volumétriques et déviatoriques:

(4.9)\[{B}^{\mathit{STAB}}={B}_{\mathit{dev}}^{\mathit{STAB}}\left(\xi ,\eta ,\zeta \right)+{B}_{\mathit{vol}}^{\mathit{STAB}}\left(0,0,0\right)\]

Il n’y a a pas de terme constant dans l’expression de \({B}^{\mathit{STAB}}\) , donc:

(4.10)\[{B}_{\mathit{vol}}^{\mathit{STAB}}\left(0,0,0\right)=0\]

Seule la partie déviatorique est conservée. On peut démontrer que la contribution \({B}^{\xi \eta }\) (et \({E}^{u,\xi \eta }\) ) est nulle dans l’écriture de la stabilisation.

4.1.2. Stabilisation des contraintes

On reprend dans l’expression de l’énergie de déformation de la structure la séparation entre les contributions venues de l’intégration réduite et de la stabilisation:

(4.11)\[W(E)={W}^{\mathit{RI}}({E}^{u,\mathit{RI}})+{W}^{\mathit{STAB}}({E}^{u,\mathit{STAB}})\]

Le tenseur de contraintes de seconde espèce de Piola-Kirchhoff étant la dérivée de cette énergie par rapport aux déformations de Green-Lagrange, on a:

(4.12)\[S=\frac{\partial W}{\partial E}=\frac{\partial {W}^{\mathit{RI}}}{\partial {E}^{u,\mathit{RI}}}\frac{\partial {E}^{u,\mathit{RI}}}{\partial E}+\frac{\partial {W}^{\mathit{STAB}}}{\partial {E}^{u,\mathit{STAB}}}\frac{\partial {E}^{u,\mathit{STAB}}}{\partial E}\]

On retrouve la séparation entre les deux contributions (intégration réduite et stabilisation) pour le tenseur des contraintes:

(4.13)\[S={S}^{\mathit{RI}}({E}^{u,\mathit{RI}})+{S}^{\mathit{STAB}}({E}^{u,\mathit{STAB}})\]

Le tenseur \({S}^{\mathit{RI}}\) vient directement de l’intégration de la loi de comportement tandis que le tenseur \({S}^{\mathit{STAB}}\) correspond aux termes de contraintes pour la stabilisation.

Le point clé de la stabilisation est de trouver un moyen optimal d’évaluer cette contrainte de stabilisation sans l’intégrer afin que le temps de calcul soit minimisé. Or, comme expliqué précédemment, pour éviter le verrouillage volumétrique, seule la partie déviatorique du tenseur de stabilisation de Piola Kirchhoff est conservée. Ce faisant, la stabilisation de la contrainte est donnée comme:

(4.14)\[{S}^{\mathit{STAB}}={C}^{\mathit{STAB}}:{E}^{u,\mathit{STAB}}\]

\({C}^{\mathit{STAB}}\) étant la partie déviatorique du comportement de Saint-Venant Kirchhoff:

(4.15)\[\begin{split}{C}^{\mathit{STAB}}={\mu}^{\mathit{STAB}}\left[\begin{array}{cccccc}\frac{4}{3}& -\frac{2}{3}& -\frac{2}{3}& 0& 0& 0\\ -\frac{2}{3}& \frac{4}{3}& -\frac{2}{3}& 0& 0& 0\\ -\frac{2}{3}& -\frac{2}{3}& \frac{4}{3}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\end{split}\]

Avec le module de cisaillement \({\mu}^{\mathit{STAB}}\) , qui, pour un matériau élastique isotrope vaut:

(4.16)\[{\mu}_{\mathit{elas}}^{\mathit{STAB}}=\frac{E}{2(1+\nu )}\]

Pour les matériaux non-élastiques, l’utilisation du module de cisaillement élastique conduit à une surestimation de la contrainte de stabilisation. Pour surmonter ce problème, le module sécant tel que défini dans [bib20] est utilisé:

(4.17)\[2{\mu}^{\mathit{STAB}}=\sqrt{\frac{{\pi}_{S}}{{\pi}_{E}}}\]

Avec:

(4.18)\[ \begin{align}\begin{aligned}{\pi}_{S}=\frac{1}{2}{S}_{\mathit{dev}}:{S}_{\mathit{dev}}\\\textrm{et}\\{\pi}_{E}=\frac{1}{2}{E}_{\mathit{dev}}^{u}:{E}_{\mathit{dev}}^{u}\end{aligned}\end{align} \]

Cette manière d’évaluer le paramètre de stabilisation est très pratiques car il n’y a pas besoin de linéarisation cohérente de la matrice ce stabilisation, ce qui est très intéressant en termes de temps de calcul. De plus, la définition proposée de \({C}^{\mathit{STAB}}\) assure une stabilisation efficace et robuste.

4.1.3. Stabilisation des opérateurs discrets

En reprenant les définitions précédentes, on calcule la contribution \({\mu}^{\mathit{STAB}}\) de la manière suivante. A chaque point de Gauss \(g\) on calcule \({{\pi}_{S}}^{g}\) et \({{\pi}_{E}}^{g}\) puis on considère les deux situations suivantes:

1. Si \(|{{\pi}_{E}}^{g}|<0,001\) ) alors \({\mu}^{\mathit{STAB},g}={K}_{\mathit{vol}}^{g}\) . C’est à dire que le coefficient de stabilisation au point de Guass courant \(g\) est égal à la partie volumique de la matrice de comportement matériau si la déformation déviatorique est inférieure à 0.1%;

2. Sinon \({\mu}^{\mathit{STAB},g}=\frac{1}{2}\sqrt{\frac{{\pi}_{S}^{g}}{{\pi}_{E}^{g}}}\) ;

Puis le coefficient de stabilisation est sommé sur chaque point de Gauss:

(4.19)\[{\mu}^{\mathit{STAB}}={\mu}^{\mathit{STAB},g}+\frac{{\omega}_{g}}{8}\]

La stabilisation des contraintes dans le cas des petites déformations, consiste à calculer les différentes contributions de la manière suivante (on rappelle que le terme \(\xi \eta\) ne contribue pas):

(4.20)\[ \begin{align}\begin{aligned}{\epsilon}^{\xi}={B}^{\xi}u\\{\epsilon}^{\eta}={B}^{\eta}u` :math:`{\epsilon}^{\xi \zeta }={B}^{\xi \zeta }u` :math:`{\epsilon}^{\eta \zeta ,\mathit{STAB}}={B}^{\eta \zeta }u\\{\sigma}^{\xi ,\mathit{STAB}}={C}^{\mathit{STAB}}{\epsilon}^{\xi}\\{\sigma}^{\eta ,\mathit{STAB}}={C}^{\mathit{STAB}}{\epsilon}^{\eta}` :math:`{\sigma}^{\xi \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{\epsilon}^{\xi \zeta }` :math:`{\sigma}^{\eta \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{\epsilon}^{\eta \zeta }\end{aligned}\end{align} \]

Et en grandes déformations:

(4.21)\[ \begin{align}\begin{aligned}{\sigma}^{\xi ,\mathit{STAB}}={C}^{\mathit{STAB}}{E}^{u,\xi }\\{\sigma}^{\eta ,\mathit{STAB}}={C}^{\mathit{STAB}}{E}^{u,\eta }` :math:`{\sigma}^{\xi \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{E}^{\xi \zeta }` :math:`{\sigma}^{\eta \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{E}^{\eta \zeta }\end{aligned}\end{align} \]

Le processus de linéarisation de Newton de l’équation d’équilibre produit deux opérateurs: la matrice tangente consistante et le vecteur des forces internes. Les deux vont donc devoir intégrer la stabilisation puisqu’ils dépendent des contraintes. Pour la matrice tangente, il y a deux contributions:

• la première concerne la stabilisation liée aux contraintes issues de l’intégration de la loi de comportement (matrice matérielle), c’est-à-dire la matrice \({K}_{m}^{\mathit{uu}}\) ;

• la seconde correspond au modèle de grande déformation, c’est la partiegéométrique qui vient de la non-linéarité cinématique, soit la matrice \({K}_{g}^{\mathit{uu}}\) .

Pour la stabilisation de la matrice tangente matérielle, on a quatre contributions:

(4.22)\[ \begin{align}\begin{aligned}{K}_{m}^{\mathit{STAB},\xi }={({B}^{\xi})}^{T}{C}^{\mathit{STAB}}({B}^{\xi})\\{K}_{m}^{\mathit{STAB},\eta }={({B}^{\eta})}^{T}{C}^{\mathit{STAB}}({B}^{\eta})\\{K}_{m}^{\mathit{STAB},\xi \zeta }={({B}^{\xi \zeta })}^{T}{C}^{\mathit{STAB}}({B}^{\xi \zeta })\\{K}_{m}^{\mathit{STAB},\eta \zeta }={({B}^{\eta \zeta })}^{T}{C}^{\mathit{STAB}}({B}^{\eta \zeta })\end{aligned}\end{align} \]

Ce qui nous donne la matrice globale de stabilisation matérielle avec le jacobien de changement de coordonnées ():

(4.23)\[{K}_{m}^{\mathit{STAB}}=\frac{8}{3}{J}_{0}\left({K}_{m}^{\mathit{STAB},\xi }+{K}_{m}^{\mathit{STAB},\eta }\right)+\frac{8}{9}{J}_{0}\left({K}_{m}^{\mathit{STAB},\xi \zeta }+{K}_{m}^{\mathit{STAB},\eta \zeta }\right)\]

Pour la partie géométrique, comme la contribution est le produit d’une matrice (géométrique) avec les contraintes, il faut aussi le faire par un produit pour la stabilisation, pour chaque nœud \(i\) et \(j\) :

(4.24)\[{K}_{g,ij}^{\mathit{STAB}}=\frac{8}{3}{J}_{0}\left({G}_{ij}^{u,\xi }{\sigma}^{\xi ,\mathit{STAB}}+{G}_{ij}^{u,\eta }{\sigma}^{\eta ,\mathit{STAB}}\right)+\frac{8}{9}{J}_{0}\left({G}_{ij}^{u,\eta \zeta }{\sigma}^{\eta \zeta ,\mathit{STAB}}+{G}_{ij}^{u,\xi \zeta }{\sigma}^{\xi \zeta ,\mathit{STAB}}\right)\]