2. Sous intégration du QUAD4#
Malgré tout le potentiel de l’élément quadrilatère isoparamétrique pour le calcul par éléments finis, l’application d’une formulation bilinéaire classique pour définir son champ de déplacement, aboutit à des résultats médiocres. Ceci s’explique par diverses raisons:
il présente une rigidité excessive («lock»), lors d’une sollicitation dont la partie flexion plane est importante.
la formulation bilinéaire classique du champ de déplacement est très sensible à la distorsion du maillage et présente un sévère «locking» lorsqu’on l’applique à un matériau incompressible.
Une solution à ces problèmes de blocage numérique consiste à calculer la matrice de rigidité par l’intermédiaire d’une intégration réduite. Le principe de cette méthode est de considérer lors du schéma d’intégration numérique moins de points d’intégrations qu’il n’en faut habituellement pour évaluer les matrices de rigidité exacte de l’élément. Partant de l’élément QUAD4 isoparamétrique [r3.03.02], on modifie le nombre de points d’intégration ainsi que le poids et les coordonnées de ces derniers, pour créer l’élément sous-intégré que nous nommerons dans ce document: QUAS4
QUAD4 |
QUAS4 |
|---|---|
Nombre de points d’intégration: 4 |
Nombre de points d’intégration: 1 |
Poids de chaque point: 1 |
Poids du point: 4 |
Coordonnées: \((x,y)=(\pm \frac{1}{\sqrt{3}},\pm \frac{1}{\sqrt{3}})\) |
Coordonnées: \((x,y)=(0,0)\) |
2.1. Formulation#
Au centre, l’opérateur gradient discrétisé utilisé pour calculer la matrice de rigidité est de la forme suivante:
\(B={B}_{c}=\left[\begin{array}{cc}{b}^{{t}_{x}}& 0\\ 0& {b}^{{t}_{y}}\\ {b}^{{t}_{y}}& {b}^{{t}_{x}}\end{array}\right]\) éq 2.1-1
Avec
\(\begin{array}{}{b}_{x}^{t}=N,x(0)=\frac{\partial N}{\partial x}\mid \xi =\eta =0\\ {b}_{y}^{t}=N,y(0)=\frac{\partial N}{\partial y}\mid \xi =\eta =0\end{array}\) éq 2.1-2
et \(N\) représente \(({N}_{1},{N}_{2},{N}_{3},{N}_{4})\) , le vecteur des fonctions de forme. Rappelons que les vecteurs \(b\) représentent la forme la plus simple de l’opérateur gradient discrétisé sous-intégré introduit par Hallquist [bib1] et qui est basé sur l’évaluation des dérivées des fonctions de forme isoparamétriques à l’origine du référentiel \((\xi ,\eta )\)
Soit:
\(\begin{array}{c}{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]=\text{Constant sur l'élément}\\ {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]=\text{Constant sur l'élément}\end{array}\) éq 2.1-3
\(A=(1/2).((\mathit{x2}-\mathit{x4}).(\mathit{y3}-\mathit{y1})+(\mathit{x3}-\mathit{x1})\mathrm{\ast }(\mathit{y4}-\mathit{y2}))\) : Aire de l’élément
La matrice de rigidité s’écrit:
\({K}_{e}={K}_{c}=A.{B}_{{c}^{T}}.C.{B}_{c}\) éq 2.1-4
\(C\) étant:
soit la matrice de comportement élastique pour les calculs en élasticité
soit la matrice tangente pour les calculs en plasticité. Notons qu’au cours de tels calculs, c’est l’intégration de la loi de comportement au point de Gauss (au centre dans notre cas) qui détermine la valeur des coefficients de \(C\) .
Finalement, les forces internes s’écrivent:
\({F}_{int}={K}_{e}.U={B}_{{c}^{T}}.{\sigma}_{c}\) éq 2.1-5
\(U\) : Vecteur des déplacements nodaux.
2.2. Echec du calcul dans Code_Aster : modes parasites#
Les calculs lancés sur le cas test n°1 avec un tel élément échouent. L’étape du calcul qui consiste à inverser \({K}_{e}\) pour déterminer les déplacement nodaux ne peut être franchie. En effet, au centre, la matrice de rigidité est singulière. En affichant le noyau de \({K}_{e}\) , nous constatons que sa dimension n’est plus trois mais cinq. Deux vecteurs se sont ajoutés au noyau et ont rendu l’inversion de la matrice de rigidité impossible. L’apparition de ces deux vecteurs supplémentaires dans le noyau de \({K}_{e}\) est directement liée au fait que nous choisissons le centre comme seul point d’intégration. Autrement dit, il existe deux champs de déplacements nodaux autres que les champs de déplacement correspondant aux mouvement de solide rigide annulant les forces internes. Ces modes représentent les modes en sablier du QUAD4 [Fig. 2.1]. Par la suite, une des étapes de stabilisation consistera à enrichir l’opérateur gradient discrétisé afin de rendre \({K}_{e}\) inversible.
Fig. 2.1 Modes en sablier du QUAD4#
2.3. Interprétation graphique du problème lié à la sous intégration#
Le problème de la sous-intégration du QUAD4 est lié à ses modes de déplacements nodaux en sablier (sollicitation en flexion). La [Fig. 2.2] nous montre que, sur de tels modes, les coordonnées du centre restent inchangées. Si ceci est en accord avec la théorie des poutres (i.e. en flexion pure, xx = 0 sur la fibre neutre), la théorie classique des éléments finis ne permet pas de différencier état déformé et non déformé d’un élément dans un tel cas. C’est pour cela que ces modes sont également appelés modes à énergie nulle, qui au niveau de Code_Aster , rendent les calculs irréalisables.
Fig. 2.2 Modes#
Ce problème a pour origine une valeur de la déformation au centre non significative de la déformation du QUAD4. Le prochain chapitre va décrire la démarche qui nous permettra de calculer la déformation du QUAD4 quels que soient les modes de déplacement de ces nœuds.