r3.06.10 Élément quadrangulaire à un point d’intégration, stabilisé par la méthode « Assumed Strain »#
Résumé:
Malgré tout le potentiel de l’élément quadrilatère isoparamétrique linéaire pour le calcul par éléments finis, l’application d’une formulation bilinéaire classique pour définir son champ de déplacement abouti à des résultats médiocres. Si la sous intégration de l’élément permet d’améliorer ses performances, elle fait cependant apparaître des modes parasites qui rendent les calculs instables. Le présent document reprend les principales étapes d’une méthode de stabilisation des calculs nommée «assumed strainmethod » et explique la manière dont elle a été implantée dans le code de calcul Code_Aster . De nombreux résultats s’appuyant sur l’élément stabilisé sont comparés et commentés afin de conclure sur les performances de la méthode.
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)\) |
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.
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. 136]. Par la suite, une des étapes de stabilisation consistera à enrichir l’opérateur gradient discrétisé afin de rendre \({K}_{e}\) inversible.
Fig. 136 Modes en sablier du QUAD4#
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. 137] 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. 137 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.
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 ).
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\)
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 [r3.06.10-cas-test-deux] 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 [(367)] 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.
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 38].
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 [(368)] et [(369)] 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.
Intégration de l’élément dans Code_Aster#
L’élément QUAS4, possède deux familles de points de Gauss. La première famille est constituée d’un point se trouvant au centre et dont le poids vaut quatre. Les calculs élémentaires se font naturellement avec la première famille. Les options calculées sont identiques à celle calculées par le QUAD4. On calcule donc au centre de l’élément les contraintes et variables internes, ainsi que le comportement tangent. Dans le cadre d’un post traitement, chaque élément présente des champs constants.
La deuxième famille est identique à celle du QUAD4 classique et comporte 4 points de Gauss. Elle est utilisée pour calculer la matrice de stabilisation. Pour stocker les forces de stabilisation, on ajoute au champ de contraintes au point de Gauss 6 composantes.
Cet élément est activé en choisissant les modélisations “C_PLAN_SI’ou “D_PLAN_SI’dans AFFE_MODELE pour des mailles QUAD4. Seuls des calculs en petites déformations sont possibles. Il reste une limitation importante: dans la version 7.2, deux types de stabilisation sont programmés: ASBQI et ASOI(1/2), mais ne sont pas accessibles par un mot-clé dans le fichier de commandes. Il faut pour les activer modifier le paramètre PROJ dans la routine NMAS2D.
Mise en évidence de l’apport de l’élément QUAS4#
Afin d’évaluer l’apport de l’élément QUAS4, nous l’avons utilisé pour effectuer les calculs des cas tests numéro un et deux (cf. [r3.06.10-stabilisation-quad-point-gauss]).
Cas test n°1 (SSLP106)#
Parmi les éléments figurant dans le [Tableau 38], seul l’élément ASOI(1/2) a pu être implanté avec succès. De plus les calculs ne convergent que pour le cas test de la poutre droite, modélisée par un maillage régulier. Les même calculs effectués sur un maillage volontairement déformé divergent [Fig. 138].
Fig. 138 Maillage déformé divergent#
Concernant les calculs effectués sur un maillage régulier, nous constatons une nette amélioration des résultats. La solution théorique est atteinte avec un maillage de soixante quatre éléments. Ceci nous permet raisonnablement de dire que le blocage de l’élément QUAD4 a disparu. Ce test correspond au test SSLP106 de la base de tests de code_aster .
Cas test n°2 (SSNP123)#
Rappel du problème: Avec une maille de type QUAD4, le calcul fait apparaître d’importantes oscillations de contraintes sur le maillage. Nous allons comparer les résultats issues du maillage QUAD4 à ceux du maillage QUAS4:
Calcul en élasticité dans le plan: isovaleurs de \({\sigma}_{yy}\)
QUAD4
\(\nu =0,4999\)
QUAS4
\(\nu =0,4999\)
Calculs en plasticité dans le plan: isovaleurs de \({\sigma}_{yy}\)
QUAD4
\(\nu =0,3\)
QUAS4
\(\nu =0,3\)
Chaque durée de calcul indiquée est une moyenne effectuée sur trois lancement.
Pour la visualisation des résultats, chaque QUAD4 a été divisé en quatre zones contenant chacune un point de Gauss. Les QUAS4 quant à eux ne sont pas redivisés.
Les tests effectués sont volontairement sévères (\(\nu\) proche de 0,5). En effet le but ici est d’atteindre les limites de l’élément QUAD4 afin de constater les performances du nouvel élément.
Commentaires#
Nous constatons à travers cette série de calculs que les résultats obtenus avec des éléments QUAD4 comportent de fortes oscillations de contraintes. Les résultats issus des calculs s’appuyant sur le QUAS4 ne présentent presque plus d’oscillations. Lorsqu’elles apparaissent, elles sont très localisées. Que ce soit pour le QUAD4 ou pour le QUAS4, les déplacement calculés au nœud \(A\) sont identiques. En effet, l’opérateur gradient discrétisé du QUAS4 est déduit du champ de déformation du QUAD4 (cf[éq 3.2-2]). Par conséquent, qu’elle soit maillée avec l’un ou l’autre de ces éléments, la structure présente la même rigidité.
Compte tenu de la sévérité du test et afin de nous rendre compte de la qualité des résultats fournis par le QUAS4, nous avons effectué des calculs s’appuyant sur des mailles quadratiques telles que le QUAD8, le QUAS8 (QUAD8 sous intégré avec 4 points de Gauss) ou encore le QUAD8 INCO_UPG, élément traitant parfaitement l’incompressibilité. La comparaison d’éléments à interpolation quadratique avec des éléments à interpolation linéaire a peu de sens. De tels tests ont été effectués en vue d’avoir une solution de référence et ils nous ont permis de porter un jugement purement qualitatif sur l’élément QUAS4.
Résultats de référence: calculs en élasticité dans le plan avec éléments quadratiques
QUAS8
\(\nu =0,4999\)
\(\mathit{dy}(A)=6,01E-2\)
Durée du calcul: \(3,9s\)
QUAD8
\(\nu =0,4999\)
\(\mathit{dy}(A)=6,0E-2\)
Durée du calcul: \(5,6s\)
QUAD8 INCO_UPG
\(\nu =0,4999\)
\(\mathit{dy}(A)=6,01E-2\)
Durée du calcul: \(5,5s\)
Calculs en plasticité avec éléments quadratiques
QUAD8 INCO_UPG
\(\nu =0,4999\)
\(\mathit{dy}(A)=1,01E-1\)
Durée du calcul: \(34,4s\)
QUAD8
\(\nu =0,3\)
\(\mathit{dy}(A)=1,01E-1\)
Durée du calcul: \(17s\)
QUAS8
\(\nu =0,3\)
\(\mathit{dy}(A)=1,01E-1\)
Durée du calcul: \(13,1s\)
Un raffinement du maillage pour chaque calcul à été effectué. Ce raffinement consiste à doubler le nombre de mailles présentes sur chaque arrête de l’entaille. Nous résumons les résultats dans les tableaux suivants:
|
Sol ref : Depl. Nœud A |
6,01E-02 |
Nombre de Mailles |
504 |
2016 |
8064 |
|||
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
|
QUAD8 INCO_UPG |
5,5 |
6,01E-02 |
18 |
6,01E-02 |
||
QUAD8 |
5,6 |
6,00E-02 |
15,5 |
6,01E-02 |
||
QUAD8 SI |
3,9 |
6,01E-02 |
10,2 |
6,01E-02 |
||
QUAD4 |
3,3 |
2,70E-02 |
8,4 |
3,60E-02 |
||
QUAS4 |
2,72 |
2,70E-02 |
5,77 |
3,60E-02 |
17,3 |
4,77E-02 |
Tableau 5.3-1 : Calculs en élasticité \(\nu =0,4999\)
|
Sol ref : Depl. Nœud A |
1,01E-01 |
Nombre de Mailles |
504 |
2016 |
8064 |
|||
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
|
QUAD8 |
17 |
1,01E-01 |
75,1 |
1,01E-01 |
||
QUAD8 SI |
13,1 |
1,01E-01 |
62,8 |
1,01E-01 |
||
QUAD4 |
6,6 |
9,32E-02 |
22 |
1,00E-01 |
||
QUAS4 |
6,5 |
9,18E-02 |
17,33 |
1,00E-01 |
133,2 |
1,01E-01 |
Tableau 5.3-2 : Calculs en plasticité \(\nu =0,3\)
|
Sol ref : Depl. Nœud A |
1,01E-01 |
Nombre de Mailles |
504 |
2016 |
8064 |
|||
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
Temps de calcul |
Depl Nœud A |
|
QUAD8 INCO_UPG |
34,4 |
1,01E-01 |
240,8 |
1,01E-01 |
||
QUAD8 |
21,7 |
8,84E-02 |
88,17 |
1,01E-01 |
||
QUAD8 SI |
20,6 |
1,01E-01 |
87,83 |
1,01E-01 |
||
QUAD4 |
5,1 |
2,06E-02 |
13,07 |
2,25E-02 |
||
QUAS4 |
4,28 |
2,06E-02 |
9,87 |
2,25E-02 |
43,75 |
2,79E-02 |
Tableau 5.3-3 : Calculs en plasticité \(\nu =0,4999\)
Par ailleurs, nous avons tracé la valeur de la contrainte \({\sigma}_{yy}\) le long d’un segment \(\mathit{AB}\) traversant l’entaille pour un maillage de base à 504 mailles.
Commentaires sur les résultats
Aspects qualitatifs:
Malgré la richesse du champ de déplacement de la maille QUAD8, la majorité des résultats présentent des oscillations.
Malgré la sous intégration du QUAS8, des oscillations apparaissent sur des mailles QUAS8 pour certains calculs en plasticité.
Quel que soit l’élément utilisé (excepté pour le QUAD8 INCO_UPG), la sous intégration reste indispensable, particulièrement pour les matériaux dont \(\nu\) approche 0,5.
Ces calculs nous amènent à établir le constat suivant: le QUAS4 reste stable vis à vis des oscillations quel que soit les paramètres de calcul utilisés.
Aspect quantitatif:
Nous constatons toujours que la valeur du déplacement du nœud A issue d’un maillage QUAS4 reste entre 3 et 5 fois plus faible que la solution de référence (convergence lente). Cette valeur reste néanmoins identique à celle calculée avec un maillage QUAD4.
Les profils tracés le long de l’entaille nous permettent d’affirmer que les valeurs de contraintes calculées avec un maillage QUAS4 sont de bien meilleure qualité que celles calculées avec un maillage QUAD4. En prenant les durées de calcul du QUAD4, la maille QUAS4 permet une réduction substantielle de la durée du calcul:
Gain de temps (temps de ref. = temps QUAD4) |
|||
Maillage |
504 |
2016 |
|
Elasticité |
\(\nu =0,3\) |
26% |
35% |
\(\nu =0,4999\) |
18% |
31% |
|
Plasticité |
\(\nu =0,3\) |
16% |
21% |
\(\nu =0,4999\) |
18% |
24% |
|
Remarquons au passage que les gains de temps semblent augmenter avec le nombre de mailles. En plasticité notamment, ce gain de temps dépend en grande majorité du nombre d’itérations de Newton nécessaires à la convergence du calcul au sein de chaque pas de temps.
Conclusions#
En élasticité comme en plasticité, que ce soit pour des matériaux dont \(\nu\) vaut \(0,3\) ou \(0,5\) , le nouvel élément QUAS4 reste stable et les résultats toujours réalistes, sans aucune oscillations de contraintes. Ceci, nous l’avons constaté, est loin d’être le cas pour le QUAD4. La stabilité de ce nouvel élément face à des cas tests aussi sévères que ceux présentés dans ce rapport est comparable à celle de l’élément quadratique QUAD8 sous intégré.
Par contre, cet élément a la convergence d’un élément linéaire en termes de nombre de degrés de liberté. Il faut donc mailler avec une finesse suffisante pour capter les gradients de contraintes de la solution cherchée. Ce raffinement nécessaire doit être mis en balance avec le gain de temps induit par la sous‑intégration.
Sur les exemples traités, le QUAS4 a permis un gain de temps de calcul significatif de l’ordre de 20% en moyenne pour des lois de comportements élastiques et élasto-plastiques. Notons que ces lois sont relativement peu coûteuses à intégrer. Des gains de temps beaucoup plus importants sont attendus pour des lois plus difficiles à intégrer.
Bibliographie#
HALLQUIST : Theoretical manual for DYNA3D. UC1D-19401 Lawrence Livermore National Lab., University of California, 1983.
BELYTSCHKO and W.E. BACHRACH : Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 43 (1986) 279-301.
BELYTSCHKO and L.P. BINDEMAN : Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems. Comput. Methods Appl. Mech. Eng., 88:311‑340, 1991
J.L. BATOZ et G. DHATT: Modélisation des structures par éléments finis volume 2, poutres et plaques. HERMES 1990.
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
7.2 |
N. TARDIEU, S. LIMOUZI EDF-R&D/AMA |
Texte initial |