r3.06.11 Élément hexaédrique à un point d’intégration, stabilisé par la méthode « Assumed Strain »#
Résumé:
L’élément hexaédrique à 8 nœuds sous-intégré standard à 1 point d’intégration introduit des modes parasites associés à une énergie nulle (modes de sablier) et peut conduire à une singularité de la matrice de raideur globale pour certaines conditions aux limites. La déficience du rang de la matrice de raideur, due à la sous-intégration, doit donc être comblée en rajoutant à la rigidité élémentaire une matrice dite de stabilisation. C’est l’objet de la méthode ASM (Assumed Strain Method) développée ici.
La principale caractéristique de cette méthode est que l’opérateur gradient discrétisé \(B\) ne dérive pas nécessairement du champ de déplacement et des relations classiques reliant la déformation au déplacement. En effet, cette méthode ASM consiste à projeter l’opérateur gradient discrétisé sur un sous-espace approprié afin d’éviter les différents types de blocage.
Table des matières
Conclusion Bibliographie Description_des_versions_du_document__
Formulation de l’élément HEXA8 à un point d’intégration#
Champ de déplacement et opérateur gradient discrétisé#
Dans l’élément hexaédrique, les coordonnées spatiales \({x}_{i}\) sont reliées aux coordonnées nodales \({x}_{\text{iI}}\) au moyen des fonctions de forme isoparamétriques \({N}_{I}\) par les formules:
\({x}_{i}={x}_{\text{iI}}{N}_{I}(\xi ,\eta ,\zeta )=\sum_{I=1}^{8}{N}_{I}(\xi ,\eta ,\zeta ){x}_{\text{iI}}\) éq 1.1-1
On utilisera dans la suite la convention de sommation pour les indices répétés. Les indices en minuscule \(i\) varient de 1 à 3 et représentent les directions spatiales. Les indices en majuscule \(I\) varient de 1 à 8 et correspondent aux nœuds de l’élément.
Les mêmes fonctions de forme sont utilisées pour définir le champ de déplacement de l’élément \({u}_{i}\) en fonction des déplacements nodaux \({u}_{\text{iI}}\) :
\({u}_{i}={u}_{\text{iI}}{N}_{I}(\xi ,\eta ,\zeta )\)
Puisque les mêmes fonctions de forme s’appliquent aux coordonnées et aux déplacements, leur dérivée matérielle s’annule et le champ de vitesse peut être donné par :
\({v}_{i}={v}_{\text{iI}}{N}_{I}(\xi ,\eta ,\zeta )\) éq 1.1-2
L’interpolation précédente du champ de vitesse va permettre de définir le taux de déformation et d’écrire les relations reliant les déformations aux vitesses nodales. Le gradient \({v}_{i,j}\) du champ de vitesse s’écrit:
\({v}_{i,j}={v}_{\text{iI}}{N}_{I,j}\)
Par convention, une virgule précédant un indice en minuscule représente une différentiation par rapport aux coordonnées spatiales. Le tenseur taux de déformation \({D}_{ij}\) est donné ensuite par la partie symétrique du gradient de vitesse: \({D}_{ij}=\frac{1}{2}({v}_{i,j}+{v}_{j,i})\)
On se donne des fonctions de forme isoparamétriques trilinéaires \({N}_{I}(\xi ,\eta ,\zeta )\) définies par :
\({N}_{I}(\xi ,\eta ,\zeta )=\frac{1}{8}(1+{\xi}_{I}\xi )(1+{\eta}_{I}\eta )(1+{\zeta}_{I}\zeta )\)
\(\xi ,\eta ,\zeta \in [-1,1],I=1,....,8\) avec \({\xi}_{I},{\eta}_{I},{\zeta}_{I}\) valant 1 ou –1 éq 1.1-3
Ces fonctions de forme transforment un cube unitaire dans l’espace \((\xi ,\eta ,\zeta )\) en un hexaèdre quelconque dans l’espace \(({x}_{1},{x}_{2},{x}_{3})=(x,y,z)\) . En combinant les équations 1.1-1, 1.1-2, 1.1-3 , on arrive à développer le champ de vitesse comme la somme d’un terme constant, des termes linéaires en \({x}_{i}\) , et des termes faisant intervenir les fonctions \({h}_{\alpha}\) :
\({v}_{i}={a}_{\mathrm{0i}}+{a}_{\mathrm{1i}}{x}_{1}+{a}_{\mathrm{2i}}{x}_{2}+{a}_{\mathrm{3i}}{x}_{3}+{c}_{\mathrm{1i}}{h}_{1}+{c}_{\mathrm{2i}}{h}_{2}+{c}_{\mathrm{3i}}{h}_{3}+{c}_{\mathrm{4i}}{h}_{4}\) i = 1,2,3 éq 1.1-4
\({h}_{1}=\eta \zeta\) \({h}_{2}=\xi \zeta\) \({h}_{3}=\eta \xi\) \({h}_{4}=\xi \eta \zeta\)
En effet, l’équation 1.1-1 permet d’écrire \(\xi ,\eta ,\zeta\) en fonction seulement des \({x}_{i}\) , des \({h}_{\alpha}\) et d’un paramètre constant. En injectant ensuite ces dernières relations dans l’équation 1.1-2 , on retrouve les expressions 1.1-4 recherchées.
En évaluant l’équation 1.1-4 aux nœuds de l’élément, on arrive aux 3 systèmes de 8 équations suivants:
\({\dot{d}}_{i}={a}_{\mathrm{0i}}s+{a}_{\mathrm{1i}}{x}_{1}+{a}_{\mathrm{2i}}{x}_{2}+{a}_{\mathrm{3i}}{x}_{3}+{c}_{\mathrm{1i}}{h}_{1}+{c}_{\mathrm{2i}}{h}_{2}+{c}_{\mathrm{3i}}{h}_{3}+{c}_{\mathrm{4i}}{h}_{4}\) \(i=1,2,3\) éq 1.1-5
Dans l’équation précédente et dans la suite, les caractères gras désignent des tenseurs d’ordre au moins 1. Ainsi, les vecteurs \({\dot{d}}_{i}\) et \({x}_{i}\) représentent, respectivement, les vitesses et les coordonnées nodales et sont donnés par:
\({\dot{d}}_{i}^{t}=({v}_{\mathrm{i1}},{v}_{\mathrm{i2}},{v}_{\mathrm{i3}},...,{v}_{\mathrm{i8}})\)
\({x}_{i}^{t}=({x}_{\mathrm{i1}},{x}_{\mathrm{i2}},{x}_{\mathrm{i3}},...,{x}_{\mathrm{i8}})\)
Les vecteurs \(s\) et \({h}_{\alpha}(\alpha =1,....,4)\) sont donnés quant à eux par
\({s}^{t}=(1,1,1,1,1,1,1,1)\)
\({h}_{1}^{t}=(1,1,-1,-1,-1,-1,1,1)\)
\({h}_{2}^{t}=(1,-1,-1,1,-1,1,1,-1)\)
\({h}_{3}^{t}=(1,-1,1,-1,1,-1,1,-1)\)
\({h}_{4}^{t}=(-1,1,-1,1,1,-1,1,-1)\)
Pour arriver à une écriture avantageuse de l’opérateur gradient discrétisé \(\mathrm{B}\) , on introduit les 3 vecteurs \({b}_{i}\) définis par: \({b}_{i}^{T}={N}_{,i}(0)=\frac{\partial N}{\partial {x}_{i}}{\mid}_{\xi =\eta =\zeta =0}\) \(i=1,2,3\) éq 1.1-6
où \(N\) représente \(({N}_{1},{N}_{2},...,{N}_{8})\) . Ces vecteurs \({b}_{i}\) représentent la forme la plus simple de l’opérateur gradient discrétisé sous-intégré introduit par Hallquist 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 ,\zeta )\) . Ils peuvent être déterminés explicitement. De plus, on peut vérifier les conditions d’orthogonalité suivantes:
\({b}_{i}^{T}.{h}_{\alpha}=0\) , \({b}_{i}^{T}.s=0\) , \({b}_{i}^{T}.{x}_{j}={\delta}_{ij}\)
\({h}_{\alpha}^{T}.s=0\) , \({h}_{\alpha}^{T}.{h}_{\beta}=8{\delta}_{\alpha \beta }\) éq 1.1-7
\(i=1à3\) , \(\alpha ,\beta =1à4\)
où \(\delta\) désigne le symbole de Kronecker.
A ce stade, on peut déterminer les constantes inconnues qui interviennent dans l’écriture du champ de vitesse ( éq 1.1-4) en multipliant scalairement l’équation ( éq 1.1-5) par \({b}_{j}^{T}\) et \({\mathrm{h}}_{\alpha}^{\mathrm{T}}\) , respectivement, et en utilisant les relations d’orthogonalité ci-dessus. On obtient:
\({a}_{ij}={b}_{j}^{T}.{\dot{d}}_{i}\) , \({c}_{\alpha i}={\gamma}_{\alpha}^{T}.{\dot{d}}_{i}\) avec \({\gamma}_{\alpha}=\frac{1}{8}\left[{h}_{\alpha}-\sum_{j=1}^{3}({h}_{\alpha}^{T}.{x}_{j}){b}_{j}\right]\)
Le champ de vitesse se met finalement sous la forme très pratique suivante:
\({v}_{i}={a}_{\mathrm{0i}}+({x}_{1}{b}_{1}^{T}+{x}_{2}{b}_{2}^{T}+{x}_{3}{b}_{3}^{T}+{h}_{1}{\gamma}_{1}^{T}+{h}_{2}{\gamma}_{2}^{T}+{h}_{3}{\gamma}_{3}^{T}+{h}_{4}{\gamma}_{4}^{T}).{\dot{d}}_{i}\) \(i=1,2,3\) éq 1.1-8
En différentiant la formule ci-dessus par rapport à \({x}_{j}\) , on obtient le gradient de vitesse:
\({v}_{i,j}=({b}_{j}^{T}+\sum_{\alpha =1}^{4}{h}_{\alpha ,j}{\gamma}_{\alpha}^{T}).{\dot{d}}_{i}=({b}_{j}^{T}+{h}_{\alpha ,j}{\gamma}_{\alpha}^{T}).{\dot{d}}_{i}\)
L’opérateur gradient discrétisé reliant le tenseur taux de déformation au vecteur des vitesses nodales par:
\({\nabla}_{s}v=B.\dot{d}\)
où:
\({\nabla}_{s}v=\left[\begin{array}{c}{v}_{x,x}\\ {v}_{y,y}\\ {v}_{z,z}\\ {v}_{x,y}+{v}_{y,x}\\ {v}_{x,z}+{v}_{z,x}\\ {v}_{y,x}+{v}_{z,y}\end{array}\right]\) , \(\dot{d}=\left[\begin{array}{c}{\dot{d}}_{x}\\ {\dot{d}}_{y}\\ {\dot{d}}_{z}\end{array}\right]\)
prend alors la forme matricielle pratique:
\(B=\left[\begin{array}{ccc}{b}_{x}^{T}+{h}_{\alpha ,x}{\gamma}_{\alpha}^{T}& 0& 0\\ 0& {b}_{y}^{T}+{h}_{\alpha ,y}{\gamma}_{\alpha}^{T}& 0\\ 0& 0& {b}_{z}^{T}+{h}_{\alpha ,z}{\gamma}_{\alpha}^{T}\\ {b}_{y}^{T}+{h}_{\alpha ,y}{\gamma}_{\alpha}^{T}& {b}_{x}^{T}+{h}_{\alpha ,x}{\gamma}_{\alpha}^{T}& 0\\ {b}_{z}^{T}+{h}_{\alpha ,z}{\gamma}_{\alpha}^{T}& 0& {b}_{x}^{T}+{h}_{\alpha ,x}{\gamma}_{\alpha}^{T}\\ 0& {b}_{z}^{T}+{h}_{\alpha ,z}{\gamma}_{\alpha}^{T}& {b}_{y}^{T}+{h}_{\alpha ,y}{\gamma}_{\alpha}^{T}\end{array}\right]\) éq 1.1-9
Les vecteurs \({\gamma}_{\alpha}\) qui interviennent dans \(B\) vérifient les conditions d’orthogonalité suivantes:
\({\gamma}_{\alpha}^{T}.{x}_{j}=0\) , \({\gamma}_{\alpha}^{T}{h}_{\beta}={\delta}_{\alpha \beta }\)
Un élément basé sur cette formulation est convergent lorsqu’il est évalué exactement. Cependant, l’évaluation de cet opérateur \(B\) en chacun des points d’intégration rend cet élément trop couteux pour les applications pratiques, et une forme simplifiée de cet élément s’impose.
Formulation variationnelle du problème#
L’extension de la forme faible du principe variationnel de Hu-Washizuau cas de la mécanique des solides non linéaires s’écrit pour un simple élément \({\Omega}_{e}\)
\(\text{δπ}(v,\stackrel{ˉ}{\dot{\epsilon}},\stackrel{ˉ}{\sigma})={\int}_{{\Omega}_{e}}{\delta \stackrel{ˉ}{\dot{\epsilon}}}^{T}.\sigma d\Omega +\delta {\int}_{{\Omega}_{e}}{\stackrel{ˉ}{\sigma}}^{T}.(\nabla {}_{s}\text{}v-\stackrel{ˉ}{\dot{\epsilon}})d\Omega -\delta {\dot{d}}^{T}.{f}^{\text{ext}}=0\) éq 1.2-1
où \(\delta\) représente une variation, \(v\) le champ de vitesse, \(\stackrel{ˉ}{\dot{\varepsilon}}\) le taux de déformation postulée, \(\stackrel{ˉ}{\sigma}\) la contrainte postulée, \(\sigma\) la contrainte évaluée par la loi constitutive, \(\dot{d}\) les vitesses nodales, \({f}^{\text{ext}}\) les forces nodales externes et \({\nabla}_{s}v\) la partie symétrique du gradient du champ de vitesse.
La formulation «Assumed strain» retenue dans la suite pour construire l’élément est basée sur une forme simplifiée du principe variationnel de Hu-Washizu qui 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 postulée. Ainsi le second terme de l’équation 1.2-1 s’élimine et on obtient:
\(\delta \pi (\stackrel{ˉ}{\dot{\varepsilon}})={\int}_{{\Omega}_{e}}{\delta \stackrel{ˉ}{\dot{\varepsilon}}}^{T}.\sigma d\Omega -\delta {\dot{d}}^{T}.{f}^{\text{ext}}=0\) éq 1.2-2
Sous cette forme, le principe variationnel est indépendant de l’interpolation de la contrainte, puisque la contrainte postulée n’intervient plus et n’a donc pas besoin d’être définie. Les équations discrétisées nécessitent donc la seule interpolation de la vitesse \(v\) et du taux de déformation postulée \(\stackrel{ˉ}{\dot{\varepsilon}}\) dans l’élément. Si \({N}_{I}\) désigne les fonctions de forme isoparamétriques de l’élément, on a:
\(\nu (x,t)=\sum_{I=1}^{m}{N}_{I}(x){\dot{d}}_{I}(t)\) où \(m\) désigne le nombre de nœuds. On en déduit:
\({\nabla}_{s}v(x,t)=B(x).\dot{d}(t)\)
Le taux de déformation postulée \(\stackrel{ˉ}{\dot{\varepsilon}}\) est défini quant à lui par: \(\stackrel{ˉ}{\dot{\varepsilon}}(x,t)=\stackrel{ˉ}{B}(x).\dot{d}(t)\) éq 1.2-3
En remplaçant l’équation 1.2-3 dans le principe variationnel 1.2-2 , on obtient:
\(\delta {\dot{d}}^{T}.{\int}_{{\Omega}_{e}}{\stackrel{ˉ}{B}}^{T}.\sigma d\Omega -\delta {\dot{d}}^{T}.{f}^{\text{ext}}=0\)
Comme \(\delta {\dot{d}}^{T}\) peut être choisi arbitrairement, l’équation précédente conduit à:
\({f}^{int}={f}^{\text{ext}}\)
avec \({f}^{int}={\int}_{{\Omega}_{e}}{\stackrel{ˉ}{B}}^{T}.\sigma (\stackrel{ˉ}{\dot{\varepsilon}})d\Omega\)
Dans l’équation ci-dessus, il est bien précisé que la contrainte \(\sigma\) est calculée par la loi constitutive à partir du taux de déformation postulée \(\stackrel{ˉ}{\dot{\varepsilon}}\) . Pour les problèmes non linéaires, \(\sigma\) peut aussi être une fonction de l’intégrale du taux de déformation postulée et d’autres variables internes. La formulation ainsi obtenue est valable pour des problèmes incluant les deux types de non-linéarités: géométriques et matérielles.
Dans le cas de problèmes linéaires, on a:
\(\sigma =C.\stackrel{ˉ}{\varepsilon}=C.\stackrel{ˉ}{B}.d\)
Les forces internes de l’élément s’écrivent:
\({f}^{int}={K}_{e}.d\) avec \({K}_{e}=\underset{{\Omega}_{e}}{\int}{\stackrel{ˉ}{B}}^{T}.C.\stackrel{ˉ}{B}d\Omega \stackrel{}{}\)
Dans une approche standard en déplacement, le taux de déformation postulée s’identifie à la partie symétrique du gradient de vitesse, ce qui revient à remplacer \(\stackrel{}{\stackrel{ˉ}{B}}\) par \(\stackrel{}{B}\) dans les expressions précédentes. On obtient:
\({K}_{e}=\underset{{\Omega}_{e}}{\int}{\stackrel{}{B}}^{T}.C.\stackrel{}{B}d\Omega \stackrel{}{}\)
Modes de « hourglass » associés à une énergie nulle#
L’écriture matricielle éq 1.1-9 de l’opérateur gradient discrétisé va permettre de comprendre l’origine des modes de hourglass, ou modes de sablier. Comme on va le voir, ces modes cinématiques sont dus à la sous-intégration et sont associés à une énergie nulle alors qu’ils induisent une déformation non nulle. Cette anomalie s’explique par la différence qu’induit la sous-intégration, entre le noyau de l’opérateur de rigidité discrétisé et le noyau de l’opérateur de rigidité continu.
Remarquons d’abord que l’opérateur gradient discrétisé sous-intégré (i.e. associé à un seul point d’intégration situé au centre de l’élément) se réduit à:
\(B=\left[\begin{array}{ccc}{b}_{x}^{T}& 0& 0\\ 0& {b}_{y}^{T}& 0\\ 0& 0& {b}_{z}^{T}\\ {b}_{y}^{T}& {b}_{x}^{T}& 0\\ {b}_{z}^{T}& 0& {b}_{x}^{T}\\ 0& {b}_{z}^{T}& {b}_{y}^{T}\end{array}\right]\) éq 1.3-1
En effet, les termes \({h}_{\alpha ,i}\) de l’équation éq 1.1-9 s’annulent au point d’intégration \(\mid \xi =\eta =\zeta =0\) .
Analysons maintenant le noyau de la matrice de rigidité obtenue par intégration. Dans le cas linéaire, cette matrice élémentaire s’écrit: \({K}_{e}=V{B}^{T}.C.B\) où \(V\) désigne le volume de l’élément.
L’examen du noyau de la rigidité sous-intégrée revient à l’étude du rang de la matrice \(B\) . Il faut donc chercher les modes de vitesse \(\dot{d}\) à déformation nulle, c’est-à-dire vérifiant:
\({\nabla}_{s}v=B.\dot{d}=0\) éq 1.3-2
On doit retrouver dans le noyau de \({K}_{e}\) les modes associés aux mouvements de corps rigide, soit en 3D, 3 translations et 3 rotations.
Le noyau de l’opérateur continu est donc de dimension 6 et se réduit aux 6 vecteurs:
\((\begin{array}{ccc}s& 0& 0\\ 0& s& 0\\ 0& 0& s\end{array})(\begin{array}{ccc}y& z& 0\\ -x& 0& z\\ 0& -x& -y\end{array})\)
Les 3 premiers vecteurs colonne vérifient 1.3-2 du fait de \({b}_{i}^{T}.s=0\) , les 3 derniers du fait de \({b}_{i}^{T}.{x}_{j}={\delta}_{ij}\) .
Mais outre les modes rigides précédents, l’opérateur gradient discrétisé donné en 1.3-1 annule les 12 autres vecteurs suivants:
\((\begin{array}{ccc}{h}_{1}& 0& 0\\ 0& {h}_{1}& 0\\ 0& 0& {h}_{1}\end{array})(\begin{array}{ccc}{h}_{2}& 0& 0\\ 0& {h}_{2}& 0\\ 0& 0& {h}_{2}\end{array})\) \((\begin{array}{ccc}{h}_{3}& 0& 0\\ 0& {h}_{3}& 0\\ 0& 0& {h}_{3}\end{array})(\begin{array}{ccc}{h}_{4}& 0& 0\\ 0& {h}_{4}& 0\\ 0& 0& {h}_{4}\end{array})\)
car \({b}_{i}^{T}.{h}_{\alpha}=0\)
Les vecteurs colonnes ci-dessus sont appelés modes de sablier ou modes de hourglass.
Les 4 modes selon \(\mathit{Ox}\) sont illustrés ci-dessous:
\({h}_{1}\text{warping}\)
\({h}_{1}\text{warping}\)
\({h}_{2}\text{bending}\)
\({h}_{1}\text{warping}\)
z
y
x
\({h}_{3}\text{bending}\)
\({h}_{1}\text{warping}\)
\({h}_{4}\text{nonphysique}\)
\({h}_{1}\text{warping}\)
Modes de sablier selon \(\mathit{Ox}\) d’après Belytschko et Bindeman
L’opérateur gradient discrétisé donné en 1.3-1 et celui intégré exactement donné en 1.1-9 calculent correctement les gradients des champs linéaires et donnent des résultats identiques lorsqu’ils sont appliqués aux 6 modes de corps rigide. Par contre, la formule exacte de \(B\) calcule proprement les autres modes \({h}_{\alpha}\) de déformation, alors que la forme sous-intégrée donne des résultats erronés.
La déficience dans le rang de l’opérateur \(B\) sous-intégré (rang 6 seulement au lieu de 18 pour \(B\) intégré exactement) se traduit par des oscillations irréalistes du maillage conduisant à des solutions qui explosent pour certaines conditions aux limites.
Ces éléments sous-intégrés nécessitent donc des techniques dites de stabilisation.
Stabilisation de type « Assumed Strain Method » (ASM)#
La démarche suivie pour stabiliser l’élément hexaédrique sous-intégré est celle de Belytschko et Bindeman [2]. La méthode développée pour contrôler les modes de hourglass dans un élément sous–intégré est dite Assumed Strain Method. Elle implique en outre la construction d’une forme appropriée de la matrice B permettant d’éviter les blocages numériques.
Formules de Hallquist et de Flanagan-Belytschko#
Jusqu’ici, pour établir les formules 1.1-9 de l’opérateur B, on a utilisé les vecteurs \({b}_{i}\) basés sur la dérivation des fonctions de forme à l’origine (équation 1.1-6 ). Cette écriture de \(B\) est dite forme de Hallquist. Une seconde forme est considérée maintenant, basée sur les valeurs moyennes des dérivées des fonctions de forme de l’élément. Cette forme est due à Flanagan et Belytschko et s’écrit:
\({\stackrel{ˆ}{b}}_{i}^{T}=\frac{1}{V}\underset{{\Omega}_{e}}{\int}{N}_{,i}(\xi ,\eta ,\zeta )d\Omega\) \(i=1,2,3\)
Ces vecteurs peuvent se calculer explicitement. Ils peuvent aussi être intégrés numériquement exactement.
Les conditions d’orthogonalité du §1.2 restent vérifiées exceptée la premièresur les éléments non parallélipipédiques:
\({\stackrel{ˆ}{b}}_{i}^{T}.{h}_{j}\lbrace \begin{array}{c}\text{=0}\text{sur les parallélipipèdes}\\ \mathrm{\ne }0\text{sinon}\end{array}\)
\({\stackrel{ˆ}{b}}_{i}^{T}.{h}_{4}=0\) sur tous les éléments.
En supposant encore vraie cette relation d’orthogonalité, le gradient de déplacement s’écrit comme précédemment:
\({v}_{i,j}=({\stackrel{ˆ}{b}}_{j}^{T}+{h}_{\alpha ,j}{\stackrel{ˆ}{\gamma}}_{\alpha}^{T}).{\dot{d}}_{i}\) avec \({\stackrel{ˆ}{\gamma}}_{\alpha}=\frac{1}{8}\left[{h}_{\alpha}-\sum_{j=1}^{3}({h}_{\alpha}^{T}.{x}_{j}){\stackrel{ˆ}{b}}_{j}\right]\)
On peut continuer à écrire: \({\nabla}_{s}v=\stackrel{ˆ}{B}\cdot \dot{d}\)
où \(\stackrel{ˆ}{B}\) apparait comme la somme d’un terme constant \(\stackrel{ˆ}{{B}_{c}}\) et d’un terme non constant \(\stackrel{ˆ}{{B}_{n}}\) définis par:
\(\stackrel{ˆ}{B}={\stackrel{ˆ}{B}}_{c}+{\stackrel{ˆ}{B}}_{n}\)
\({\stackrel{ˆ}{B}}_{c}=\left[\begin{array}{ccc}{\stackrel{ˆ}{b}}_{x}^{T}& 0& 0\\ 0& {\stackrel{ˆ}{b}}_{y}^{T}& 0\\ 0& 0& {\stackrel{ˆ}{b}}_{z}^{T}\\ {\stackrel{ˆ}{b}}_{y}^{T}& {\stackrel{ˆ}{b}}_{x}^{T}& 0\\ {\stackrel{ˆ}{b}}_{z}^{T}& 0& {\stackrel{ˆ}{b}}_{x}^{T}\\ 0& {\stackrel{ˆ}{b}}_{z}^{T}& {\stackrel{ˆ}{b}}_{y}^{T}\end{array}\right]\) \({\stackrel{ˆ}{B}}_{n}=\left[\begin{array}{ccc}{\stackrel{ˆ}{X}}_{1234}^{T}& 0& 0\\ 0& {\stackrel{ˆ}{Y}}_{1234}^{T}& 0\\ 0& 0& {\stackrel{ˆ}{Z}}_{1234}^{T}\\ {\stackrel{ˆ}{Y}}_{1234}^{T}& {X}_{1234}^{T}& 0\\ {\stackrel{ˆ}{Z}}_{1234}^{T}& 0& {X}_{1234}^{T}\\ 0& {\stackrel{ˆ}{Z}}_{1234}^{T}& {\stackrel{ˆ}{Y}}_{1234}^{T}\end{array}\right]\)
avec \({\stackrel{ˆ}{X}}_{1234}^{T}=\sum_{\alpha =1}^{4}{h}_{\alpha ,x}{\stackrel{ˆ}{{\gamma}_{\alpha}}}^{T}\) \({\stackrel{ˆ}{Y}}_{1234}^{T}=\sum_{\alpha =1}^{4}{h}_{\alpha ,y}{\stackrel{ˆ}{{\gamma}_{\alpha}}}^{T}\) \({\stackrel{ˆ}{Z}}_{1234}^{T}=\sum_{\alpha =1}^{4}{h}_{\alpha ,z}{\stackrel{ˆ}{{\gamma}_{\alpha}}}^{T}\)
Cette formulation, même si elle est moins rigoureuse que celle de Hallquist (puisqu’elle suppose vraies toutes les relations d’orthogonalité) donne de meilleurs résultats en termes de précision et de convergence.
Projection sur un champ de déformation \(\stackrel{ˉ}{B}\)#
On applique une projection à l’opérateur gradient discrétisé \(\stackrel{ˆ}{B}\) pour en déduire un opérateur \(\stackrel{ˉ}{B}\) possédant certaines bonnes propriétés. L’objectif de la projection est double:
Elle permet, d’une part, d’éliminer le blocage volumique de l’élément fini dans le cas incompressible
Elle évite, d’autre part, le blocage dû aux termes de cisaillement transverse excessifs dans les problèmes à flexion dominante.
L’opérateur \(\stackrel{ˆ}{B}\) est remplacé par un opérateur \(\stackrel{ˉ}{B}\) tel que:
\(\lbrace \begin{array}{c}\stackrel{ˆ}{B}={\stackrel{ˆ}{B}}_{c}+{\stackrel{ˆ}{B}}_{n}\\ \stackrel{ˉ}{B}={\stackrel{ˆ}{B}}_{c}+{\stackrel{ˉ}{B}}_{n}\end{array}\)
Seule la partie non constante \({\stackrel{ˆ}{B}}_{n}\) est projetée, la partie constante \({\stackrel{ˆ}{B}}_{c}\) reste inchangée.
Deux projections différentes conduisent aux 2 éléments finis suivants:
Elément ASQBI (Assumed Strain Quintessential Bending Incompressible) !
\(\stackrel{ˉ}{{B}_{n}}=\left[\begin{array}{ccc}{\stackrel{ˆ}{X}}_{1234}^{T}& -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{Y}}_{3}^{T}-\nu {\stackrel{ˆ}{Y}}_{24}^{T}& -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{Z}}_{2}^{T}-\nu {\stackrel{ˆ}{Z}}_{34}^{T}\\ -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{X}}_{3}^{T}-\nu {\stackrel{ˆ}{X}}_{14}^{T}& {\stackrel{ˆ}{Y}}_{1234}^{T}& -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{Z}}_{1}^{T}-\nu {\stackrel{ˆ}{Z}}_{34}^{T}\\ -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{X}}_{2}^{T}-\nu {\stackrel{ˆ}{X}}_{14}^{T}& -\stackrel{ˉ}{\nu}{\stackrel{ˆ}{Y}}_{1}^{T}-\nu {\stackrel{ˆ}{Y}}_{24}^{T}& {\stackrel{ˆ}{Z}}_{1234}^{T}\\ {\stackrel{ˆ}{Y}}_{12}^{T}& {\stackrel{ˆ}{X}}_{12}^{T}& 0\\ {\stackrel{ˆ}{Z}}_{13}^{T}& 0& {\stackrel{ˆ}{X}}_{13}^{T}\\ 0& {\stackrel{ˆ}{Z}}_{23}^{T}& {\stackrel{ˆ}{Y}}_{23}^{T}\end{array}\right]\) éq 2.2-1
Où: \(\stackrel{}{\stackrel{ˉ}{\nu}=\frac{\nu}{1-\nu }}\) , \({\stackrel{ˆ}{X}}_{14}^{T}={h}_{1,x}\stackrel{ˆ}{{\gamma}_{{1}^{T}}}+{h}_{4,x}\stackrel{ˆ}{{\gamma}_{{4}^{T}}}\) \({\stackrel{ˆ}{Z}}_{2}^{T}={h}_{2,z}\stackrel{ˆ}{{\gamma}_{{2}^{T}}}\)
Elément ADS (Assumed Deviatoric Strain) :
\(\stackrel{ˉ}{{B}_{n}}=\left[\begin{array}{c}\frac{2}{3}{\stackrel{ˆ}{X}}_{1234}^{T}-\frac{1}{3}{\stackrel{ˆ}{Y}}_{1234}^{T}-\frac{1}{3}{\stackrel{ˆ}{Z}}_{1234}^{T}\\ -\frac{1}{3}{\stackrel{ˆ}{X}}_{1234}^{T}\text{}\frac{2}{3}{\stackrel{ˆ}{Y}}_{1234}^{T}-\frac{1}{3}{\stackrel{ˆ}{Z}}_{1234}^{T}\text{}\\ -\frac{1}{3}{\stackrel{ˆ}{X}}_{1234}^{T}-\frac{1}{3}{\stackrel{ˆ}{Y}}_{1234}^{T}\text{}\frac{2}{3}{\stackrel{ˆ}{Z}}_{1234}^{T}\\ \text{}{\stackrel{ˆ}{Y}}_{12}^{T}\text{}{\stackrel{ˆ}{X}}_{12}^{T}0\\ \text{}{\stackrel{ˆ}{Z}}_{13}^{T}0{\stackrel{ˆ}{X}}_{13}^{T}\\ 0{\stackrel{ˆ}{Z}}_{23}^{T}\text{}{\stackrel{ˆ}{Y}}_{23}^{T}\end{array}\right]\) éq 2.2-2
Les schémas de projection sont détaillés dans [2]. On indiquera néanmoins les grandes lignes de l’obtention de leurs opérateurs \(\stackrel{ˉ}{{B}_{n}}\) .
Pour l’élément ADS, cela revient dans un premier temps à décomposer l’opérateur \(\stackrel{ˆ}{{B}_{n}}\) en la somme d’un terme sphérique et d’un terme déviatorique. Ensuite, la partie sphérique est sous-intégrée (i.e. évaluée au point \(\xi =\eta =\zeta =0\) , ce qui revient à l’annuler). Cette procédure permet ainsi de traiter les termes diagonaux de la déformation, donc le blocage volumique, et on peut vérifier dans éq 2.2-2 que la somme des trois premiers termes dans chacun des vecteurs colonne est nulle. Pour éviter le blocage en cisaillement transverse, il faut agir cette fois sur les termes non diagonaux de la déformation, et donc annuler ceux responsables d’un cisaillement excessif. Le résultat de cette étape se traduit par la suppression de certains termes dans les trois dernières lignes de l’opérateur \(\stackrel{ˆ}{{B}_{n}}\) (comparer les formules donnant \(\stackrel{ˆ}{{B}_{n}}\) et \(\stackrel{ˉ}{{B}_{n}}\) ).
Pour l’élément ASQBI, le traitement du blocage en cisaillement est le même que pour l’élément ADS. Il en résulte que les trois dernières lignes de l’opérateur \(\stackrel{ˉ}{{B}_{n}}\) sont identiques dans les deux éléments. En revanche, le traitement du blocage volumique est un peu différent pour l’élément ASQBI (voir éq 2.2-1 ). On peut vérifier cependant que, lorsque \(\nu \to \frac{1}{2}\) , donc \(\stackrel{ˉ}{\nu}\to 1\) , la somme des trois premiers termes dans chaque vecteur colonne de \(\stackrel{ˉ}{{B}_{n}}\) est nulle.
Choix des éléments finis#
On a comparé les deux éléments finis ADS et ASQBI sur un certain nombre de tests. Il ressort que ASQBI donne de meilleurs résultats que ADS en élasticité alors que ADS donne de meilleurs résultats en plasticité, son comportement étant un peu trop souple dans les cas élastiques. Le cas test de la poutre console en flexion [V6.04.196] en est une bonne illustration. Ces observations sont en bonne conformité avec celles de Belytschko et Bindeman [2].
On a donc choisi d’implanter dans Aster l’élément ASQBI en élasticité et l’élément ASD en plasticité pour bénéficier dans chaque cas du meilleur élément. De plus, on évite à l’utilisateur d’avoir à choisir l’élément fini. On rappelle que la seule différence entre les deux éléments est la matrice \(\stackrel{ˉ}{{B}_{n}}\) .
Intégration de l’élément dans Code_Aster#
Description et utilisation#
Cet élément s’appuie sur les mailles 3D volumiques HEXA8.
Modélisation#
On affecte la modélisation 3D_SI aux mailles HEXA8 désignées. Les éléments de face habituels de la modélisation 3D sont affectés sur les mailles QUAD4.
Matériau#
Tous les coefficients matériau relatifs aux lois de comportement valides en petites déformations pour les modélisations 3D sont utilisables.
Conditions aux limites et chargement#
Tous les chargements et conditions aux limites disponibles sur les éléments 3D et les éléments de face sont disponibles.
Options de post-traitement#
Toutes les options de post-traitement habituellement disponibles pour les modélisations 3D sont utilisables.
Calcul en flambement linéaire#
L’option RIGI_MECA_GE étant activée dans le catalogue de l’élément, il est possible d’effectuer un calcul de flambement classique après assemblage des matrices de rigidité élastique et géométrique.
Calculs non linéaires#
Tous les comportements disponibles pour les modélisation 3D habituelles sont utilisable, en petites déformations. En ce qui concerne les grandes déformations, la seule option utilisable est l””approximation PETIT_REAC.
L’intégration numérique est réalisée avec un point de Gauss, tout comme en non linéaire matériel.
Implantation#
Les options sont activées dans le catalogue meca_hexs8.cata.
Validation#
Les tests validant cet élément sont, dans la version 9 de Code_Aster :
SSNV196: Poutre 3D en flexion en élasticité et en plasticité
HSNV125G, PERFE01A, SDNV103, SSND105 utilisent un seul élément pour valider des comportements non linéaires, permettant ainsi des économies de temps calcul (1 seul point de Gauss)