4. Fonctionnement de la macro_commande#

Cette macro_commande a pour but de restreindre au strict nécessaire les données relatives à une simulation sur un point matériel pour un modèle de comportement incrémental.

Le fonctionnement interne allège donc le fichier de commandes de l’utilisateur, en réalisant des opérations répétitives pour ce genre de situations.

4.1. Cas SUPPORT=”ELEMENT”#

Le fonctionnement est :

  1. création d’un maillage d’un seul élément à un seul point de Gauss (un tétraèdre à quatre nœuds en 3D, un triangle à trois nœuds en 2D) (voir par exemple [V6.04.176]).

  2. affectation d’un modèle 3D ou C_PLAN ou D_PLAN

  3. affectation du matériau sur ce maillage ;

  4. affectation des chargements :

  5. en ce qui concerne les déformations imposées, pour chaque composante affectée via un des mots-clés de EPSI_IMPOSE, création d’un chargement unitaire en déformation qui sera multiplié par la fonction du temps fournie pour cette composante par l’utilisateur ;

  6. en ce qui concerne les contraintes imposées, pour chaque composante affectée via un des mots-clés de SIGM_IMPOSE, création d’un chargement unitaire en contraintes qui sera multiplié par la fonction du temps fournie pour cette composante par l’utilisateur ;

  7. Appel à STAT_NON_LINE. Tous les mots-clés ayant des valeurs par défaut sont utilisés, sauf s’ils sont surchargés par l’utiliseur (NEWTON, CONVERGENCE, SUIVI_DDL, ARCHIVAGE, RECH_LINEAIRE) et de l’état initial.

L’ensemble des résultats (six ou quatre composantes de contraintes et de déformations, variables internes) sont stockés dans une table (tabres). Pour chaque composante (colonne de la table) figure l’évolution en fonction du temps.

4.2. Cas SUPPORT=”POINT”#

Dans ce cas, plutôt que d’utiliser un élément fini (même unique) pour réaliser le calcul, SIMU_POIN_MAT fait appel à une commande dédiée, CALC_POINT_MAT, qui englobe en Fortran l’appel direct à la routine 3D d’intégration des comportements, NMCOMP. Ceci n’est disponible que dans le cas des petites déformations.

Rappelons que NMCOMP est la routine générale d’intégration des lois de comportement, appelée par tous les éléments finis 3D et 2D. Elle permet le calcul en un point (ce point étant le point d’intégration pour un élément fini) des contraintes et variables internes à l’instant actuel, connaissant les contraintes et variables internes à l’instant précédent; et l’accroissement de déformation actuel. (cf. [D5.04.01] et [R5.03.14] ). Plus précisément, à l’instant \({t}_{i},\) et à l’itération \(n\) le tenseur des contraintes \({\sigma}_{i}^{n}\) , en un point est calculé à partir de \(({\sigma}_{i-1},{\alpha}_{i-1})\) et de l’incrément de déformation \(\Delta {\varepsilon}_{i}^{n}\) .

Quand toutes les composantes du tenseur des déformations sont fournies, l’algorithme est immédiat : il s’agit d’une simple boucle en temps, contenant pour chaque incrément temporel la donnée de l’état mécanique de l’incrément précédent et le tenseur (symétrique) correpondant à l’accroissement de déformation connu.

Mais dans le cas contraire, soit que l’on fournit seulement \(n\) composantes de l’histoire des déformations, \(n<6\) , soit que l’on fournit \(n\) composantes de l’histoire des contraintes, l’algorithme est le suivant :

  • par défaut, toute composante non spécifiée correspond à une composante de contrainte assujettie à rester nulle (condition de Neumann)

  • les équations à résoudre sont (en utilisant la notation sous forme vectorielle des tenseurs symétriques d’ordre 2) :

  • \({\sigma}_{i}=F(\Delta {\varepsilon}_{i};{\sigma}_{i-1},{\alpha}_{i-1})\) où \(F\) représente le résultat de l’intégration du comportement par NMCOMP

  • pour j variant de 1 à 6 :

  • soit \({({\sigma}_{i})}_{j}={g}_{j}({t}_{i})\)

  • soit \({({\varepsilon}_{i})}_{j}={g}_{j}({t}_{i})\)

  • où \(\stackrel{ˉ}{{\sigma}_{j}}(t)\) et \(\stackrel{ˉ}{{\varepsilon}_{j}}(t)\) sont données par SIGM_IMPOSE / EPSI_IMPOSE.

Ceci peut s’écrire encore :

pour chaque instant \({t}_{i},\) résoudre :

\(R({Y}_{i})=0\) avec \({Y}_{i}=Y({t}_{i})=\left\lbrace \begin{array}{}\left[{\sigma}_{i}\right]\\ \left[{\varepsilon}_{i}\right]\end{array}\right\rbrace\) et \(R({Y}_{i})=\left\lbrace \begin{array}{}\left[{\sigma}_{i}-F(\Delta {\varepsilon}_{i};{\sigma}_{i-1},{\alpha}_{i-1})\right]\\ \left[{C}_{1}\right]{\sigma}_{i}+\left[{C}_{2}\right]{\varepsilon}_{i}-g({t}_{i})\end{array}\right\rbrace\)

qui est un système non linéaire d’ordre 12.

La dernière relation traduit les condtitions de contraintes ou déformations imposées : les matrices \({C}_{1}\) et \({C}_{2}\) ne contiennent que des termes sur la diagonale, valant 1 si la composante correspondante est imposée, sachant que l’on ne peut avoir à la fois contrainte et déformation imposée.

Par exemple, si la déformation \({\varepsilon}_{yy}\) est imposée, la dernière relation s’écrit :

\(\left[{C}_{1}\right]{\sigma}_{i}+\left[{C}_{2}\right]{\varepsilon}_{i}-g({t}_{i})=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right](\begin{array}{}{\sigma}_{xx}\\ {\sigma}_{yy}\\ {\sigma}_{zz}\\ {\sigma}_{xy}\\ {\sigma}_{xz}\\ {\sigma}_{yz}\end{array})+\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right](\begin{array}{}{\varepsilon}_{xx}\\ {\varepsilon}_{yy}\\ {\varepsilon}_{zz}\\ {\varepsilon}_{xy}\\ {\varepsilon}_{xz}\\ {\varepsilon}_{yz}\end{array})=(\begin{array}{}0\\ {\varepsilon}_{yy}({t}_{i})\\ 0\\ 0\\ 0\\ 0\end{array})\)

Dans le cas où l’utilisateur a précisé (via les mots clés MATR_C1 ,MATR_C2 et VECT_IMPO) des relations supplémentaires, celles-ci sont prises en compte directement dans les matrices \({C}_{1}\) et \({C}_{2.}\)

La résolution de ce système d’équations non linéaires s’effectue par une méthode de Newton :

  • initialisation : \({Y}_{i}^{0}={Y}_{i-1}+{[{K}_{i}^{0}]}^{\text{-1}}\left[\begin{array}{}0\\ -\left[{C}_{1}\right]{\sigma}_{i-1}-\left[{C}_{2}\right]{\varepsilon}_{i-1}+g({t}_{i})\end{array}\right]\)

  • a l’itération n :

        • \(\delta {Y}_{i}^{n+1}=-{[{K}_{i}^{n}]}^{\text{-1}}R({Y}_{i}^{n})={[{K}_{i}^{n}]}^{\text{-1}}\left[\begin{array}{c}F(\Delta {\varepsilon}_{i}^{\mathrm{n}};{\sigma}_{i-1},{\alpha}_{i-1})-{\sigma}_{i}^{\mathrm{n}}\\ -\left[{C}_{1}\right]{\sigma}_{i}^{\mathrm{n}}-\left[{C}_{2}\right]{\varepsilon}_{i}^{\mathrm{n}}+\mathrm{g}({t}_{i})\end{array}\right]\) ;

\(\Delta {Y}_{i}^{n}=\delta {Y}_{i}^{n}+\Delta {Y}_{i}^{n-1}\) ; \({Y}_{i}^{n}=\Delta {Y}_{i}^{n}+{Y}_{i-1}\)

avec

\([{K}_{i}^{n}]=\left[\frac{\partial R}{\partial Y}\right]=\left[\begin{array}{cc}1& -{(\frac{\partial \sigma }{\partial \varepsilon })}_{i}^{n}\\ \left[{C}_{1}\right]& \left[{C}_{2}\right]\end{array}\right]\) et \([{K}_{i}^{0}]=\left[\begin{array}{cc}1& -{(\frac{\partial \sigma }{\partial \varepsilon })}_{i}^{0}\\ \left[{C}_{1}\right]& \left[{C}_{2}\right]\end{array}\right]\)

où \({(\frac{\partial \sigma }{\partial \varepsilon })}_{i}^{0}\) représente l’opérateur tangent de prédiction (option RIGI_MECA_TANG, cf. [R5.03.01, R5.03.02]) et \({(\frac{\partial \sigma }{\partial \varepsilon })}_{i}^{n}\) représente l’opérateur tangent cohérent (option FULL_MECA), cf. [R5.03.01, R5.03.02]). Ces opérateurs peuvent être remplacés par l’opérateur d’élasticité suivant les mots-clés PREDICTION , REAC_ITER.

La seule non linéarité du problème provient du comportement : \(F(\Delta {\varepsilon}_{i}^{n};{\sigma}_{i-1},{\alpha}_{i-1})\) .

Dans le cas d’un comportement linéaire, on vérifie que la solution du problème est obtenue à l’issue de phase de prédiction.

La convergence des itérations est vérfiée :

  • soit en valeur relative, (mot clé RESI_GLOB_RELA) :

    • \(\max(\frac{\underset{j=1,6}{\max}∣{({R}_{i}^{n})}_{j}∣}{\underset{j=1,6}{\max}∣{({\sigma}_{i}^{0})}_{j}∣};\frac{\underset{j=7,12}{\max}∣{({R}_{i}^{n})}_{j}∣}{\underset{j=7,12}{\max}∣{({R}_{i}^{0})}_{j}∣})<\text{RESI\_GLOB\_RELA}\)

Dans ce cas, les termes en contraintes et en déformations sont séparés pour l’examen du critère de convergence pour éviter les problèmes dus aux différences d’ordres de gandeur.

  • soit en valeur absolue (mot clé RESI_GLOB_MAXI)ouvavleur du dénominateur proche de zéro dans le critère relatif ci-dessus :

\(\underset{j=1,12}{\max}∣{({R}_{i}^{n})}_{j}∣<\text{RESI\_GLOB\_MAXI}\)

Les options de calcul de la rigidité tangente par perturbation et la gestion automatique du pas de temps sont également activées, comme dans [U4.51.03].

Dans la résolution précédente, les termes en contraintes sont adimensionnalisés, pour éviter un mauvais conditionnement de la matrice jacobienne. On divise donc pour la résolution tous les termes en contraintes par le max des termes diagonaux de l’opérateur d’élasticité ; il faut donc fournir dans DEFI_MATERIAU le mot-clé ELAS ou ELAS_ORTH ou ELAS_ISTR.

Dans le cas SUPPORT=”POINT”, certains mots clés n’ont pas d’utilité :

  • dans le cas de la recherche linéaire, celle-ci n’est pas programmée dans la version actuelle

  • dans le cas de l’archivage, seul le mot-clé LIST_INST est pris en compte

  • dans le cas CONVERGENCE / RESI_REFE_RELA, ce mot-clé sans objet : le résidu par valeur de référence n’a pas de sens pour un point matériel.

  • Dans le cas où la valeur du mot clé DEFORMATIONn’est pas PETIT, on alarme l’utilisateur en précisant que le type de DEFORMATIONchoisi est incompatible avec SUPPORT=POINT, et que l’on utilise donc SUPPORT=ELEMENT, sauf dans le cas où on fournit toutes les composantes du gradient de transformation (mot clé GRAD_IMPOSE ).