3. Stabilisation du QUAD4 à un point de Gauss#
La stabilisation du QUAS4 s’opère en deux étapes:
Enrichir l’opérateur gradient discrétisé \({B}_{c}\) et permettre ainsi de calculer l’énergie de déformation liée aux modes de déplacements en sablier ( hourglass modes );
Interpoler un champ de déformation / de contraintes permettant de rendre compte de la déformation / des contraintes sur l’ensemble de l’élément tout en intégrant le QUAD4 au centre ( assumed strain stabilization ).
3.1. Principe variationnel du problème#
Celui ci est tiré de la forme faible du principe variationnel de Hu-Washizu:
\(\delta \pi (\upsilon ,\stackrel{ˉ}{\dot{\varepsilon}},\stackrel{ˉ}{\sigma})={\int}_{{\Omega}_{e}}{\delta \stackrel{ˉ}{\dot{\varepsilon}}}^{T}.\sigma d\Omega +\delta {\int}_{{\Omega}_{e}}{\stackrel{ˉ}{\sigma}}^{T}.(\nabla {}_{s}\text{}\upsilon -\stackrel{ˉ}{\dot{\varepsilon}})d\Omega -\delta {\dot{d}}^{T}.{f}^{\text{ext}}=0\) éq 3.1-1
Avec
\(\stackrel{ˉ}{\dot{\varepsilon}}(x,t)=B(x).\dot{d}(t)\)
La stabilisation «Assumed strain» s’appuie sur le fait que la contrainte postulée est choisie orthogonale à la différence entre la partie symétrique du gradient de vitesse et le taux de déformation. Par conséquent, cela nous permet d’écrire:
\(\delta {\dot{d}}^{T}.{\int}_{{\Omega}_{e}}{B}^{T}.\sigma d\Omega -\delta {\dot{d}}^{T}.{f}^{\text{ext}}=0\)
Et donc que :
\({f}^{int}={\int}_{{\Omega}_{e}}{B}^{T}.\sigma (\stackrel{ˉ}{\dot{\varepsilon}})d\Omega\)
3.2. Enrichissement de l’opérateur gradient discrétisé#
Enrichir l’opérateur gradient discrétisé \({B}_{c}\) consiste à en faire un nouvel opérateur \(B\) en lui ajoutant une troisième composante \(\gamma\) qui, comme \({b}_{x}\) et \({b}_{y}\) , est un vecteur de \({\mathrm{\Re }}^{4}\) . Cependant, comme l’opérateur initial \({B}_{c}\) calcule correctement le gradient des champs de déplacement linéaires, la nouvelle composante \(\gamma\) doit être orthogonale à ces derniers. Les étapes du calcul de cet opérateur enrichi sont détaillées dans le paragraphe [Section 1.2] du document joint intitulé: «Compte rendu bibliographique».
Note: Dans ce compte rendu, le nouvel opérateur \(B\) relie le tenseur taux de déformation \(\dot{\varepsilon}\) et le vecteur des vitesses nodales \({\dot{d}}_{i}\) . Cette formulation nous permet dans la suite du compte rendu de raisonner en terme d’incréments de déplacement, et par conséquent de traiter des problèmes incrémentaux, effectués sur plusieurs pas de temps. Pour les calculs élastiques (solutions obtenue en un pas de temps) nous formulons le problème à partir des déplacements nodaux: l’opérateur \(B\) relie alors tenseur des déformations \(\varepsilon\) et déplacements nodaux \({u}_{i}\) .
Le nouvel opérateur \(B\) s’appuie sur l’expression du champ de déplacement du QUAD4 écrit par Belytschko et Bachrach [bib2] :
\({u}_{i}=({\Delta}^{T}+x.{b}_{x}^{T}+y.{b}_{y}^{T}+{h}_{\gamma}^{T}){u}_{i}\) éq 3.2-1
Et s’écrit:
Avec:
\(\begin{array}{c}\Delta =1/4\left[t-({t}^{T}.x){b}_{x}-({t}^{T}.y){b}_{y}\right],\gamma =1/4\left[h-({h}^{T}.x){b}_{x}-({h}^{T}.y){b}_{y}\right]\\ {b}_{x}^{t}=\frac{1}{\mathrm{2A}}\left[({y}_{2}-{y}_{4}),({y}_{3}-{y}_{1}),({y}_{4}-{y}_{2}),({y}_{1}-{y}_{3})\right]\\ {b}_{y}^{t}=\frac{1}{\mathrm{2A}}\left[({x}_{4}-{x}_{2}),({x}_{1}-{x}_{3}),({x}_{2}-{x}_{4}),({x}_{3}-{x}_{1})\right]\\ h=\xi \eta \end{array}\)
\(\xi ,\eta\) sont les coordonnées de référence. \({u}_{i}\) sont les vecteurs de déplacements nodaux et \(h\) sont les valeurs prises par la fonction \(h\) aux quatre nœuds.
\({d}^{T}=\left[{u}_{x},{u}_{y}\right]\)
Notons que \(\gamma\) peut s’exprimer directement en fonction des coordonnées nodales:
\(\gamma =\frac{1}{4}\left[\begin{array}{c}{x}_{2}({y}_{3}-{y}_{4})+{x}_{3}({y}_{4}-{y}_{2})+{x}_{4}({y}_{2}-{y}_{3})\\ {x}_{3}({y}_{1}-{y}_{4})+{x}_{4}({y}_{3}-{y}_{1})+{x}_{1}({y}_{4}-{y}_{3})\\ {x}_{4}({y}_{1}-{y}_{2})+{x}_{1}({y}_{2}-{y}_{4})+{x}_{2}({y}_{4}-{y}_{1})\\ {x}_{1}({y}_{3}-{y}_{2})+{x}_{2}({y}_{1}-{y}_{3})+{x}_{3}({y}_{2}-{y}_{1})\end{array}\right]\) éq 3.2-3
La forme [(3.1)] de l’opérateur \(B\) est équivalente à celle du QUAD4. Cependant cette écriture particulière de l’opérateur permet de différencier les termes à intégrer au centre et les termes de stabilisation. C’est uniquement en intervenant sur la valeur de ces termes de stabilisation que nous allons pouvoir améliorer les performances de l’élément.
3.3. Interpolation du champ de déformation :#
La forme générale d’un champ «assumed strain» au sein d’un QUAD4 est donnée par Belytschko et Bindeman [bib3] . Elle est de la forme suivante:
\({\epsilon}_{\text{assumed}\text{strain}}=\left[\begin{array}{c}{\epsilon}_{x(\text{centre})}+{q}_{x}{e}_{1}{h}_{,x}+{q}_{y}{e}_{2}{h}_{,y}\\ {\epsilon}_{y(\text{centre})}+{q}_{x}{e}_{2}{h}_{,x}+{q}_{y}{e}_{1}{h}_{,y}\\ {\mathrm{2\epsilon }}_{xy(\text{centre})}+{q}_{x}{e}_{3}{h}_{,y}+{q}_{y}{e}_{3}{h}_{,x}\end{array}\right]=\left[\begin{array}{c}{\epsilon}_{x(\text{centre})}+{\epsilon}_{x(\text{stab})}\\ {\epsilon}_{y(\text{centre})}+{\epsilon}_{y(\text{stab})}\\ {\mathrm{2\epsilon }}_{xy(\text{centre})}+{\mathrm{2\epsilon }}_{xy(\text{stab})}\end{array}\right]\) éq 3.3-1
Avec:
\(\begin{array}{c}{q}_{x}=\gamma .{u}_{x}\\ {q}_{y}=\gamma .{u}_{y}\end{array}\)
et \({e}_{1},{e}_{2},{e}_{3}\) qui varient selon les considération physiques de chaque auteur [Tableau 3.1].
Chaque triplet de valeurs caractérise un élément et donne lieu à une interpolation particulière de la déformation:
Elément |
\(\mathit{e1}\) |
\(\mathit{e2}\) |
\(\mathit{e3}\) |
QUAD4 |
1 |
0 |
1 |
ASMD |
1/2 |
-1/2 |
1 |
ASBQI |
1 |
\(-\nu\) |
0 |
ASOI |
1 |
-1 |
0 |
ASOI(1/2) |
1/2 |
-1/2 |
0 |
Ainsi nous pouvons déduire l’expression de l’opérateur gradient discrétisé découlant du champ de déformation supposé (\({\varepsilon}_{(\mathit{stab})}\) ) de l’élément. Nous notons ce nouvel opérateur \({B}_{n}\) .
\(\begin{array}{}\text{QUAD4}:{B}_{n}=\left[\begin{array}{cc}{h}_{,x}{\gamma}^{t}& 0\\ 0& {h}_{,y}{\gamma}^{t}\\ {h}_{,y}{\gamma}^{t}& {h}_{,x}{\gamma}^{t}\end{array}\right]\\ \text{ASBQI}:{B}_{n}=\left[\begin{array}{cc}{h}_{\text{,x}}{\text{γ}}^{t}& -\stackrel{ˉ}{\nu}{\text{h}}_{\text{,y}}{\text{γ}}^{t}\\ -\stackrel{ˉ}{\nu}{\text{h}}_{\text{,x}}{\text{γ}}^{t}& {h}_{\text{,y}}{\text{γ}}^{t}\\ 0& 0\end{array}\right]\\ \text{ASOI}(\text{1/2}):{B}_{n}=\left[\begin{array}{cc}\frac{1}{2}{h}_{,x}{\gamma}^{t}& -\frac{1}{2}{h}_{,y}{\gamma}^{t}\\ -\frac{1}{2}{h}_{,x}{\gamma}^{t}& \frac{1}{2}{h}_{,y}{\gamma}^{t}\\ 0& 0\end{array}\right]\end{array}\) éq 3.3-2
Cette écriture nous permet de tester aisément les différents éléments.
La matrice de rigidité s’écrit alors de la manière suivante:
\({K}_{e}={K}_{c}+{K}^{\text{stab}}\) éq 3.3-3
Avec:
Et
\(\begin{array}{}{K}_{\text{stab}}=\underset{{O}_{e}}{\int}{B}_{{c}^{T}}C{B}_{n}\text{d}O+\underset{{O}_{e}}{\int}{B}_{{n}^{T}}{\text{C B}}_{c}\text{d}O+\underset{{O}_{e}}{\int}{B}_{{n}^{T}}{\text{C B}}_{n}\text{d}O\\ =\sum_{i=1}^{4}\text{JAC}(i).({B}_{c}^{T}.C.{B}_{n}(i)+{B}_{n}^{T}(i).C.{B}_{c}+{B}_{n}^{T}(i).C.{B}_{n}(i))\end{array}\) éq 3.3-5
Enfin nous calculons les forces internes de la manière suivante:
Notes:
Bien que le calcul de \({K}_{\text{stab}}\) nécessite une somme sur les quatre points de Gauss, l’intégration de la loi de comportement qui détermine la valeur des termes de \(C\) , s’effectue au centre. Par ailleurs, les équations [(3.2)] et [(3.3)] nous montre que les calculs restent relativement volumineux. Une solution à ce problème (non encore implantée dans Code_Aster ) consiste à effectuer les calculs en se plaçant dans un repère tournant avec l’élément (cf.[§3.4] du compte rendu bibliographique). Ceci a pour avantage:
de supprimer le calcul des termes croisés: \({B}_{c}.C.{B}_{n}\text{et}{B}_{n}.C.{B}_{c}\) ;
un meilleurs traitement du blocage en cisaillement transverse;
une écriture de la loi de comportement mieux adaptée aux problèmes incluant des non‑linéarités géométriques.