r5.03.14 Intégration implicite et explicite des relations de comportements non linéaires#

Résumé

Ce document décrit deux méthodes d’intégration (une implicite et une explicite) pour la résolution de problèmes avec comportements non linéaires dans les opérateurs STAT_NON_LINE [R5.03.01] et DYNA_NON_LINE [R5.05.05].

La méthode numérique implicite présentée est celle de Newton, avec ou sans recherche linéaire.

La méthode numérique explicite est celle de Runge Kutta d’ordre 2 à pas adaptatif.

Intégration explicite d’une relation de comportement incrémentale#

Pour résoudre le système non linéaire relatif à la loi de comportement,on peut utiliser une méthode d’intégration explicite de type Runge Kutta d’ordre 2 (option ‘RUNGE_KUTTA’ de l’opérande ALGO_INTE du mot-clé facteur COMPORTEMENT) pour les problèmes incrémentaux de comportement non linéaire.

Les relations concernées actuellement sont:

VISCOCHAB	Comportement élasto-viscoplastique de Chaboche, à deux centres, avec effet de mémoire et de restauration
VENDOCHAB	Comportement viscoplastique avec endommagement de Chaboche
POLYCRISTAL	Comportement élasto-viscoplastique polycristallin homogénéisé
MONOCRISTAL	Comportements viscoplastiques monocristallins des métaux
HAYHURST	Comportement viscoplastique avec endommagement de Hayhurst

La relation POLYCRISTAL n’est accessible que par la méthode explicite présentée ci-après.

Ce type d’intégration permet d’implanter très facilement un nouveau modèle de comportement [bib2]. On décrit le calcul du champ de contraintes à partir d’un incrément de déformation donné en suivant l’évolution des variables internes.

Utilisation d’une méthode explicite#

On suppose que le système à résoudre peut s’exprimer sous la forme :

\(\frac{dY}{dt}=F(Y,t;\Delta {\varepsilon}_{i}^{n},{\sigma}_{i-1})\)

où \(Y\) représente l’ensemble des variables internes du modèle, y compris la variation de déformation plastique \(\Delta {\varepsilon}^{p}\) .

Leurs dérivées en temps ne sont supposées dépendre que des variables internes et des contraintes à l’instant précédent, et de l’accroissement de déformation totale fourni par l’algorithme global de Newton \(\Delta {\varepsilon}_{i}^{n}\) .

Cette hypothèse repose en particulier sur la relation linéaire (loi de Hooke) entre le tenseur des contrainte à l’instant actuel et la partie élastique du tenseur des déformations : \(\Delta \sigma =A(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{p})\) , ce qui permet d’exprimer \(\sigma\) en fonction de \(\Delta {\varepsilon}^{p}\) dans les équations différentielles \(\frac{dY}{dt}=F(Y,t;\sigma (\Delta {\varepsilon}_{i}^{n},{\sigma}_{i-1},Y))\)

L’une des techniques les plus simples à mettre en œuvre pour résoudre ces équations différentielles est l’utilisation de méthodes explicites. Pour qu’elles restent efficaces numériquement, il est indispensable de leur associer un contrôle de pas automatique pour conserver un bon compromis coût précision.

Les méthodes de Runge et Kutta explicites et emboîtées [bib1], [bib2] sont sans doute les schémas d’intégration les plus simples respectant ces critères. Leur principe est d’associer deux schémas d’intégration d’ordre différent pour contrôler le pas de temps en fonction d’une précision requise. En fonction de l’ordre d’intégration choisi, plusieurs algorithmes sont disponibles et le schéma le plus simple est une méthode d’ordre deux.

On intègre selon le schéma suivant :

\({Y}_{t+h}={Y}^{(2)}\) si le critère de précision est satisfait ;
\({Y}^{(2)}=Y+\frac{h}{2}[F(Y,t)+F({Y}^{(1)},t+h)]\) avec \({Y}^{(1)}=Y+hF(Y,t)\)

La différence entre \({Y}^{(2)}\) (schéma d’ordre 2) et \({Y}^{(1)}\) (schéma d’ordre 1, Euler) fournit une estimation de l’erreur d’intégration et permet de contrôler la taille du pas de temps \(h\) qui est initialisé à \(\Delta {t}_{i}\) pour la première tentative. Ainsi, la méthode reste efficace si le comportement reste élastique au cours de l’incrément et on dispose naturellement de sous pas locaux indépendants de l’incrément global pour intégrer, avec une meilleure précision, l’évolution des variables internes aux points de Gauss où la non linéarité du comportement est la plus significative.

Contrôle du pas de temps#

La stratégie du contrôle du pas est définie à partir d’une norme de l’écart entre les deux méthodes d’intégration: \(\parallel {Y}^{(2)}-{Y}^{(1)}\parallel\) et de la précision requise par l’utilisateur \(\eta\) (mot-clé: RESI_INTE_RELA). Le critère retenu est le suivant, où l’on note \(Y=({y}_{1,}{y}_{2,}\mathrm{...}{y}_{N})\) :

\(\deltay (t)=\underset{j=\mathrm{1...N}}{\sup}\left\lbrace \frac{\mid {y}_{j}^{(2)}-{y}_{j}^{(1)}\mid }{\max[\varepsilon ,∣{y}_{j}(t)∣]}\right\rbrace <\eta\)

Le paramètre \(\eta\) est fixé à 0,001. La précision d’intégration souhaitée \(\eta\) doit être cohérente avec le niveau de précision requis pour l’étape globale. La valeur par défaut \(\eta ={10}^{-6}\) s’avère souvent pénalisante en termes de temps calcul. Elle peut être souvent être augmentée d’un ordre ou deux, en particulier pour les modèles élasto-visco-plastiques, sans réelle perte de précision (à voir au cas par cas).

Si le critère n’est pas vérifié, le pas de temps est re-découpé selon une méthode heuristique. Lorsque le pas de temps devient trop faible (\(h<\text{}\) précision machine), un code retour est émis, et le re-découpage global du pas de temps est activé, si demandé par l’utilisateur.

Influence sur l’étape globale#

Dans sa version actuelle, la méthode ne fournit pas de matrice tangente ; son évaluation est toutefois possible (mais coûteuse) par une technique de perturbations (activée sous COMPORTEMENT le mot-cléTYPE_MATR_TANG=”PERTURBATION”) . Si on utilise la matrice élastique pour l’équilibre global, le nombre d’itérations de Newton pour converger peut devenir très important, ce qui amène souvent à nombre d’incréments plus important pour que la résolution globale soit facilitée.

Pour améliorer la convergence, il est possible d’initialiser les itérations de Newton à partir de la solution précédemment calculée : au lieu de considérer l’incrément de déplacement pour la première itération à l’instant courant (obtenu via la matrice élastique), on initialise l’incrément de déplacement recherché pour la première itération \(\Delta {U}_{1}\) à partir de la valeur \(\Delta {U}_{n}^{-}\) , solution à l’incrément précédent (cf. figure page suivante). Le mot-clé PREDICTION=‘EXTRAPOLE’ (sous NEWTON) [U4.51.03] permet d’initialiser avec l’incrément convergé du pas précédent (pondéré par le rapport des pas de temps). Cette estimation est projetée sur le champ des déplacements cinématiquement admissibles pour que la solution finale vérifie bien les conditions aux limites de Dirichlet.

../../../../_images/1000468A000069BB000045C30B87E6F3E3C607B7.svg

Intégration implicite d’une relation de comportement incrémentale.#

Pour résoudre les systèmes d’équations non linéaires relatifs aux problèmes incrémentaux de comportement non linéaire (opérande RELATION du mot clé facteur COMPORTEMENT). il est possible d’utiliser une méthode d’intégration implicite de type Newton-Raphson, à laquelle on peut adjoindre une étape de recherche linéaire (respectivement option ‘NEWTON’ et option “NEWTON_RELI” de l’opérande ALGO_INTE du mot-clé facteur CONVERGENCE).

Pour plus de détails consulter le document [U4.51.03] du manuel utilisateur.

Les relations concernées sont par exemple :

VISCOCHAB	Comportement élasto-viscoplastique de J.L. Chaboche, à deux centres, avec effet de mémoire et de restauration
MONOCRISTAL	Comportements viscoplastiques monocristallins des métaux
IRRAD3M	Comportement élasto-plastique sous irradiation des aciers inoxydables
…

Les relations VISCOCHAB et MONOCRISTAL peuvent également être intégrées par la méthode explicite exposée au § 2 .

Intégration implicite par la méthode de Newton#

L’environnement PLASTI permet dans Code_Aster d’intégrer par une méthode de Newton, de façon systématique, des modèles de comportement non linéaire. La démarche pour résoudre les équations du modèle peut se schématiser de la manière suivante: connaissant toutes les variables (contraintes, variables internes) à l’instant précédent \({t}_{i-1}\) , et l’accroissement de déformation totale fourni au point considéré par l’algorithme global de Newton \(\varepsilon (\Delta {u}_{i}^{n})\) , noté aussi \(\Delta {\varepsilon}_{i}^{n}\) , il faut trouver la solution du système :

\(f(\Delta Y)=0=\left[\begin{array}{c}g(\Delta y)\\ l(\Delta y)\\ i(\Delta y)\\ j(\Delta y)\\ k(\Delta y)\end{array}\right]\) avec \(\Delta y=(\begin{array}{c}\Delta \sigma \\ \Delta {\varepsilon}^{p}\\ \Delta {X}_{1}\\ \Delta {X}_{2}\\ \Delta p\end{array})={y}_{i}-{y}_{i-1}\) .

La première équation représente par exemple le relation contrainte-déformation élastique \(g(\Delta y)=\Delta \sigma -A(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{p})\) , la deuxième exprimant l’écoulement plastique, et les suivantes correspondent aux lois d’évolution des différentes variables internes, vectorielles ou scalaires, notées ici \(\Delta {X}_{1}\) , \(\Delta {X}_{2}\) et \(\Delta p\) .

Le nombre d’équations à résoudre varie évidemment en fonction du comportement (voir par exemple les comportements MONOCRISTAL [R5.03.11], VISCOCHAB [R5.03.12], IRRAD3M [R5.03.23]).

On résout ce système par la méthode de Newton proposée dans l’environnement PLASTI, soit:

\(\lbrace \begin{array}{c}\frac{\partial {f}^{l}}{\partial \Delta {y}_{k}}d(\Delta {y}_{k})\text{= -}{f}^{l}(\Delta {y}_{k})\\ \Delta {y}_{k+1}=\Delta {y}_{k}+d(\Delta {y}_{k})\end{array}\) en itérant en \(k\) jusqu’à convergence.

Il reste à calculer la matrice jacobienne du système: \(\frac{\text{df}}{{\text{dy}}_{k}}\) . Elle a la forme suivante:

\(J=\left[\begin{array}{ccccc}\frac{\partial g}{\partial \Delta \sigma }& \frac{\partial g}{\partial \Delta {\varepsilon}^{p}}& \frac{\partial g}{\partial \Delta {X}_{1}}& \frac{\partial g}{\partial \Delta {X}_{2}}& \frac{\partial g}{\partial \Delta p}\\ \frac{\partial l}{\partial \Delta \sigma }& \frac{\partial l}{\partial \Delta {\varepsilon}^{p}}& \frac{\partial l}{\partial \Delta {X}_{1}}& \frac{\partial l}{\partial \Delta {X}_{2}}& \frac{\partial l}{\partial \Delta p}\\ \frac{\partial i}{\partial \Delta \sigma }& \frac{\partial i}{\partial \Delta {\varepsilon}^{p}}& \frac{\partial i}{\partial \Delta {X}_{1}}& \frac{\partial i}{\partial \Delta {X}_{2}}& \frac{\partial i}{\partial \Delta p}\\ \frac{\partial j}{\partial \Delta \sigma }& \frac{\partial j}{\partial \Delta {\varepsilon}^{p}}& \frac{\partial j}{\partial \Delta {X}_{1}}& \frac{\partial j}{\partial \Delta {X}_{2}}& \frac{\partial j}{\partial \Delta p}\\ \frac{\partial k}{\partial \Delta \sigma }& \frac{\partial k}{\partial \Delta {\varepsilon}^{p}}& \frac{\partial k}{\partial \Delta {X}_{1}}& \frac{\partial k}{\partial \Delta {X}_{2}}& \frac{\partial k}{\partial \Delta p}\end{array}\right]\)

Le critère d’arrêt des itérations porte, suivant les comportements, sur la nullité du résidu: \(\frac{f({y}_{k})}{f({y}_{0})}<\epsilon\) ou la stationnarité de la solution: \(\mid {y}_{k+1}-{y}_{k}\mid <\epsilon\) .

Les valeurs initiales \({y}_{0}\) sont choisies nulles par défaut, mais peuvent prendre des valeurs particulières suivant les comportements.

Plus précisément, les tests d’arrêt et les valeurs initiales sont :

Comportements	Critère d’arrêt	Initialisation
VISCOCHAB	\(\mid {y}_{k+1}-{y}_{k}\mid <\epsilon \mid {y}_{k+1}\mid\)	Solution explicite et modification de l’initialisation en cas de mauvaise convergence
MONOCRISTAL	\(∣f({y}_{k})∣<\epsilon ∣f({y}_{0})∣\)	Solution nulle
IRRAD3M	\(∣f({y}_{k})∣<\epsilon ∣f({y}_{0})∣\)	Solution élastique

La méthode utilisée permet un re-découpage local du pas de temps, soit systématique, soit en cas de non convergence.

Opérateur tangent#

L’opérateur tangent peut être obtenu directement à partir du système précédent. En effet, le système formé des équations du modèle à convergence, \((f(\Delta {y}_{\mathrm{CV}})=0)\) est vérifié en fin d’incrément. Pour une petite variation de \(f\) , en considérant cette fois l’accroissement de déformation totale \(\Delta {\varepsilon}_{i}^{n}\) comme variable et non comme paramètre, le système reste à l’équilibre et on vérifie \(\mathrm{df}=0\) , c’est-à-dire:

\(\frac{\partial f}{\partial \Delta \sigma }\delta \Delta \sigma +\frac{\partial f}{\partial \Delta {\varepsilon}_{i}^{n}}\delta \Delta {\varepsilon}_{i}^{n}+\frac{\partial f}{\partial \Delta {\varepsilon}^{p}}\delta \Delta {\varepsilon}^{p}+\frac{\partial f}{\partial \Delta {X}_{1}}\delta \Delta {X}_{1}+\frac{\partial f}{\partial \Delta {X}_{2}}\delta \Delta {X}_{2}+\frac{\partial f}{\partial \Delta p}\delta \Delta p=0\)

Ce système peut encore s’écrire:

\(\frac{\partial f}{\partial y}\delta (y)=X,\text{avec}Y=\left[\begin{array}{c}\Delta \sigma \\ \Delta {Z}_{s}\end{array}\right]\text{et}X=\left[\begin{array}{c}\delta \Delta {\varepsilon}_{i}^{n}\\ 0\end{array}\right]\)

où \(\Delta Z\) représente les inconnues complémentaires de \(\Delta \sigma\) dans \(y\) .

En écrivant la matrice jacobienne sous la forme:

\(J.\mathrm{\delta Y}=\left[\begin{array}{cc}{Y}_{0}& {Y}_{1}\\ {Y}_{2}& {Y}_{3}\end{array}\right]\left[\begin{array}{c}\Delta \sigma \\ \Delta Z\end{array}\right]\)

puis en opérant par éliminations et substitutions successives, l’opérateur tangent recherchépeut s’écriredirectement:

\({(\frac{\partial \Delta \sigma }{\partial \Delta \varepsilon })}_{i}^{n}={({Y}_{0}-{Y}_{1}{Y}_{3}^{-1}{Y}_{2})}^{-1}\)

Les équations précédentes montrent que l’on est conduit à réutiliser la même matrice jacobienne \(J\) que précédemment pour évaluer l’opérateur tangent.

Intégration implicite par la méthode de Newton avec recherche linéaire#

Afin d’améliorer la robustesse de la méthode de Newton, on offre la possibilité d’activer la recherche linéaire.

Présentation de l’algorithme général de recherche linéaire#

Soit un système d’équations non-linéaires à résoudre: \(f(x)=0\)

On associe à \(f\) une fonctionnelle \(F\) .

La méthode de Newton avec recherche linéaire consiste à trouver à chaque itération un pas d’avancement suivant la direction de descente qui minimise (au sens fort ou au sens faible) la fonctionnelle. Cela se ramène en un enchainement des étapes suivantes:

calcul de la direction de descente: \({d}_{k}=-{[\nabla f({x}_{k})]}^{-1}.f({x}_{k})\)
calcul du pas d’avancement: \(\rho\) tel que \(F({x}_{k}+\rho .{d}_{k})<F({x}_{k})\) , le meilleur étant le \(\rho\) qui minimise la fonctionnelle \(\rho =\mathrm{argmin}[F({x}_{k}+\rho .{d}_{k})]\)
actualisation des inconnues: \({x}_{k+1}={x}_{k}+\rho {d}_{k}\)

On associe à \(f\) la fonctionnelle \(F(\rho )=\frac{1}{2}{\parallel f({x}_{k}+\rho .{d}_{k})\parallel }^{2}\) .

Contrairement à la recherche linéaire de Newton global où l’on détermine \(\rho\) par la résolution de \(F'(\rho )=0\) , on cherche ici un \(\rho \in [0,1],\rho \ne 0\) qui diminue suffisamment la fonctionnelle \(F\) . On ne va pas chercher un pas d’avancement supérieur au pas de Newton (\(\rho =1\) ). Le principe est décrit dans [1].

Intervalle de confiance:

Pour savoir si un \(\rho\) diminue suffisamment la fonctionnelle \(F\) , on utilise la règle d’Armijo. Celle-ci permet de caractériser un intervalle de confiance en disant que \(\rho\) convient si \(F(\rho )\le F(0)+\omega F'(0)\rho\) où \(\omega\) est un réel fixé inférieur à 1. Cette expression peut se représenter graphiquement par la : la courbe en pointillés bleu est le graphe de la fonctionnelle \(h(\rho )={(7\rho -4)}^{4}+{(7\rho -5)}^{3}\) , celle est violet est la tangente à l’origine, de pente \(F'(0)\) et celle en pointillés rouge est la droite de pente \(\omega F'(0)\) avec \(\omega =0,1\) . Dans cet exemple, \(\rho\) diminue suffisamment la fonctionnelle si \(\rho \in [0;0,871]\) .

Principe du calcul de \(\rho\) :

On commence d’abord par essayer le pas de Newton (\(\rho =1\) ). Si ce pas satisfaisait la règle d’Armijo, il est retenu et on continue les itérations de Newton. Sinon, on fait une interpolation quadratique connaissant \(F(0)\) , \(F'(0)\) et \(F(1)\) . On détermine alors un nouveau \(\rho ={\rho}_{q}\) . Si ce pas satisfaisait la règle d’Armijo, il est retenu et on continue les itérations de Newton. Sinon, on réalise une interpolation cubique connaissant \(F(0)\) , \(F'(0)\) , \(F(1)\) et \(F({\rho}_{q})\) . On détermine alors un nouveau \(\rho ={\rho}_{c}\) . Si ce pas satisfaisait la règle d’Armijo, il est retenu et on continue les itérations de Newton. Sinon, on continue les interpolations cubiques entre 0 et les 2 derniers \(\rho\) . Le paragraphe ci-dessous détaille le calcul de \(\rho\) .

Détails du calcul de \(\rho\) :

On commence d’abord par essayer le pas de Newton (\(\rho =1\) ). À ce stade, on connait \(F(0)=\frac{1}{2}{\parallel f({x}_{k})\parallel }^{2}\) et \(F'(0)=-{\parallel f({x}_{k})\parallel }^{2}\) .

On calcule \(F(1)=\frac{1}{2}{\parallel f({x}_{k}+{d}_{k})\parallel }^{2}\) .

Si la règle d’Armijo n’est pas satisfaite, c’est-à-dire si \(F(1)>F(0)+\omega F'(0)\) , alors on on modélise \(F(\rho )\) par un polynôme d’ordre 2 (interpolation quadratique):

\(q(\rho )=[F(1)-F(0)-F'(0)]{\rho}^{2}+F'(0)\rho +F(0)\)

Le minimum de ce polynôme est obtenu pour \({\rho}_{q}=\frac{-F'(0)}{2[F(1)-F(0)-F'(0)]}\) . Notons que si \(F(1)\) est beaucoup plus grand que \(F(0)\) , alors le minimum du polynôme est atteint en 0. Or on ne souhaite pas avoir un pas d’avancement trop petit, sous peine de devoir faire un très grand nombre d’itérations de Newton. Pour cela, on impose un \({\rho}_{\min}>0\) (par exemple \({\rho}_{\min}=0,1\) ). On peut aussi imposer une borne supérieure \({\rho}_{\max}\) (par exemple \({\rho}_{\max}=0,6\) ). Ainsi, si \({\rho}_{q}<{\rho}_{\min}\) alors \({\rho}_{q}={\rho}_{\min}\) et si \({\rho}_{q}>{\rho}_{\max}\) alors \({\rho}_{q}={\rho}_{\max}\) .

Ensuite on calcule \(F({\rho}_{q})=\frac{1}{2}{\parallel f({x}_{k}+{\rho}_{q}.{d}_{k})\parallel }^{2}\) .

Si la règle d’Armijo n’est pas satisfaite, c’est-à-dire si \(F({\rho}_{q})>F(0)+\omega F'(0){\rho}_{q}\) , alors on on modélise \(F(\rho )\) par un polynôme d’ordre 3 (interpolation cubique):

\(c(\rho )=a{\rho}^{3}+b{\rho}^{2}+F'(0)\rho +F(0)\)

avec

\(\left[\begin{array}{}a\end{array}\right]=\frac{1}{{\rho}_{1}-{\rho}_{2}}\left[\begin{array}{cc}\frac{1}{{\rho}_{1}^{2}}& \frac{-1}{{\rho}_{2}^{2}}\\ \frac{-{\rho}_{2}}{{\rho}_{1}^{2}}& \frac{{\rho}_{1}}{{\rho}_{2}^{2}}\end{array}\right]\left[\begin{array}{}F({\rho}_{1})-F(0)-F'(0){\rho}_{1}\\ F({\rho}_{2})-F(0)-F'(0){\rho}_{2}\end{array}\right]\)

où \({\rho}_{1}\) et \({\rho}_{2}\) sont les deux dernières valeurs échouées de \(\rho\) .

Le point minimal de \(c(\rho )\) est:

\({\rho}_{c}=\frac{-b+\sqrt{{b}^{2}-\mathrm{3a}F'(0)}}{\mathrm{3a}}\) .

On continue les interpolations cubiques tant que la règle d’Armijo n’est pas satisfaite. Dans la pratique, on limite à 2 le nombre d’interpolations cubiques.

Figure 3.3.1-a : illustration de la règle d’Armijo

Mise en œuvre de la recherche linéaire dans PLASTI#

Le système à résoudre est:

\(f(\Delta y)=0\)

On fait les itérations de Newton suivantes:

mise à jour des variables à l’instant «+»
calcul du résidu \(f(\Delta {y}_{i})=0\)
calcul de la Jacobienne \(\frac{\partial f}{\partial y}(\Delta {y}_{i})\)
résolution du système linéaire \(\frac{\partial f}{\partial y}(\Delta {y}_{i}).\delta \Delta {y}_{i}=-f(\Delta {y}_{i})\)
recherche linéaire: calcul de \(\rho\)
actualisation de la solution: \(\Delta {y}_{i+1}=\Delta {y}_{i}+\rho \delta \Delta {y}_{i}\)

Pour l’étape de recherche linéaire:

On introduit la fonctionnelle \(F(\rho )=\frac{1}{2}{\parallel f({y}_{i}+\rho .\delta \Delta {y}_{i})\parallel }^{2}\) .

Les paramètres de la recherche linéaire sont fixés en dur:

Paramètre de la règle d’Armijo	\(\omega =0,1\)
Borne min pour le rabattement	\({\rho}_{\min}=0,1\)
Borne max pour le rabattement	\({\rho}_{\max}=0,5\)
Nombre d’itérations d’interpolations cubiques	\({\mathrm{it}}_{\max}=2\)

Remarque:

Certains auteurs prennent des valeurs encore plus petites pour \(\omega\) (par exemple \({10}^{-4}\) ).

Implantation d’un nouveau modèle de comportement#

Méthodes possibles#

En utilisant les architectures précédentes, un nouveau modèle de comportement peut être développé relativement facilement :

la méthode la plus simple est l’intégration explicite, par la méthode de RUNGE_KUTTA (voir § 2 ) ;
l’intégration implicite complète par le méthode de NEWTON dans l’environnement PLASTI (voir § 3 ) demande un peu plus de développements, en particulier le calcul de la matrice jacobienne ;
il est toutefois possible d’utiliser directement le développement réalisé pour l’intégration explicite (2 routines) , sans routine supplémentaire, pour effectuer des calculs avec intégration implicite, la matrice jacobienne étant alors calculée par perturbation.

Implantation d’un nouveau modèle de comportement par la méthode explicite#

Pour implanter un modèle de comportement XXX, utilisable avec ALGO_INTE_=”RUNGE_KUTTA”, il suffit d’écrire 2 routines :

une routine récupération des données matériau, appelée par la routine d’aiguillage LCMATE
une routine calculant les dérivés temporelles des variables internes, appelée par la routine d’aiguillage LCDVIN.

La mise à jour des catalogues concerne uniquement DEFI_MATERIAU et C_COMPORTEMENT, ainsi que le catalogue de la loi de comportement (cf. [D5.04.01]).

Implantation complète d’un nouveau comportement dans PLASTI#

Sans faire une description informatique précise des développements (voir à ce propos [D5.04.01]), on peut toutefois décrire l’architecture de PLASTI et les lieux où intervenir pour intégrer un nouveau comportement.

Architecture générale de PLASTI :

- Appel aux routines récupérant les paramètres matériau (à l’instant \({t}_{i}\) car dépendant éventuellement des variables de commandes) : LCMATE => routines spécifiques de récupération du matériau, communes avec RUNGE_KUTTA (§3)
Intégration élastique sur \(\Delta t\) : LCELAS=> par défaut, élasticité linéaire isotrope ou anisotrope ; si nécessaire, appel à des routines spécifiques pour le calcul élastique XXXELA.
Évaluation du seuil pour l’intégration élastique : LCCNVX => appel aux routines d’évaluation du seuil XXXCVX
Si \(\mathrm{seuil}>0\) , alors résolution par Newton : LCPLASqui appelle LCPLNL, qui réalise la boucle de Newton. Sa structure est la suivante :
- Initialisation : LCINIT => par défaut initialisation à zéro : \({y}^{0}=0\) ; si nécessaire, appel aux routines spécifiques d’initialisation XXXINI
- Itérations de Newton :
- Incrémentation de la solution \({y}^{k}={y}_{i-1}+\Delta {y}^{k}\)
- calcul du résidu \(-f(\Delta {y}^{k})\) LCRESI => appel aux routines calculant le résidu XXXRES
- calcul de la matrice jacobienne \(\frac{\text{df}}{{\text{dy}}_{k}}\) : LCJACB => appel aux routines spécifiques XXXJAC
- résolution du système linéaire \(\frac{\partial {f}^{l}}{\partial \Delta {y}_{k}}d(\Delta {y}_{k})\text{= -}{f}^{l}(\Delta {y}_{k})\) par la méthode de Gauss
- incrémentation de \(\Delta {y}_{k+1}=\Delta {y}_{k}+d(\Delta {y}_{k})\)
- test de convergence LCCONV => par défaut, test sur l’erreur relative sur la solution \(\mid {y}_{k+1}-{y}_{k}\mid <\epsilon \mid {y}_{k+1}\mid\) ; si nécessaire, appel aux routines spécifiques XXXCVG
- si non convergence, poursuite des itérations, sauf si le nombre maximum est atteint, auquel cas on sort de PLASTI avec un code retour, pour activer le re-découpage local ou global du pas de temps (suivant les données utilisateur : ITER_INTE_PAS pour le re-découpage local et DEFI_LIST_INST pour le re-découpage global).
- Si convergence, incrémentation \({y}^{k}={y}_{i-1}+\Delta {y}^{k}\) et sortie.

Pour implanter un modèle de comportement XXX, il faut donc écrire a minima les routines :

XXXMAT: récupération des données matériau, appelée par LCMATE
XXXCVX pour l’évaluation du seuil initial
XXXRES: calcul du résidu, appelée par LCRESI
XXXJAC : calcul de la matrice jacobienne, appelée par LCJACB.
On peut aussi modifier les routines suivantes, si besoin, pour définir de façon particulière :
l’élasticité XXXELA appelée par LCELAS
l’initialisation XXXINI appelée par LCINIT
le test de convergence XXXCVG, appelé par LCCONV.

La mise à jour des catalogues concerne uniquement DEFI_MATERIAU et C_COMPORTEMENT, ainsi que le catalogue de la loi de comportement (voir [D5.04.01]).

Implantation implicite facile d’un nouveau modèle de comportement dans l’environnement PLASTI#

Il est possible d’utiliser directement les 2 routines (récupération des données matériau, dérivées des variables internes) utilisées avec ALGO_INTE_=”RUNGE_KUTTA” pour réaliser un intégration implicite. En effet, le système d’équation différentielles résolu par RUNGE_KUTTA peut s’écrire :

\(\sigma ={\sigma}_{i-1}+\Delta \sigma ={\sigma}_{i-1}+A(\Delta {\varepsilon}_{i}^{n}-\Delta {\varepsilon}^{\text{th}}-\Delta {\varepsilon}^{p}(Y))=G(Y)\)

\(\frac{dY}{dt}=F(Y,t;\sigma )\)

où \(Y\) représente l’ensemble des variables internes du modèle. On suppose ici que les premières valeurs de \(Y\) sont les composantes de la variation de déformation plastique \(\Delta {\varepsilon}^{p}\) . La relation entre le tenseur des contraintes et la partie élastique du tenseur des déformations est généralement linéaire, mais peut être évaluée de façon non linéaire par une expression spécifique.

Une fois programmée la routine permettant de calculer \(\frac{dY}{dt}=F(Y,t;\sigma )\) , il est possible de l’utiliser pour une intégration implicite, ce qui consiste à résoudre (cf. 3.1 ) :

\(R(\Delta Z)=0=\left[\begin{array}{c}{R}_{1}(\Delta Z)\\ {R}_{2}(\Delta Z)\end{array}\right]\) avec \(\Delta Z=(\begin{array}{c}\Delta \sigma \\ \Delta Y\end{array})={Z}_{i}-{Z}_{i-1}\) .

Le premier système d’équations représente la relation contrainte - déformation élastique : \({R}_{1}(\Delta Z)=\sigma -G(Y)=0\) ;
Le deuxième exprime les lois d’évolution des différentes variables internes : \({R}_{2}(\Delta Z)=\Delta Y-\Delta t.F(Y,t;\sigma )=0\) , , où \(G\) et \(F\) sont calculés par la routine explicite.

Pour résoudre ce système par la méthode de Newton proposée dans l’environnement PLASTI, soit:

\(\lbrace \begin{array}{c}\frac{\partial {f}^{l}}{\partial \Delta {Z}_{k}}d(\Delta {Z}_{k})\text{= -}{f}^{l}(\Delta {Z}_{k})\\ \Delta {Z}_{k+1}=\Delta {Z}_{k}+d(\Delta {Z}_{k})\end{array}\) , il faut disposer de la matrice jacobienne \(\frac{\partial {f}^{l}}{\partial \Delta {Z}_{k}}\) .

Son calcul peut être effectué automatiquement par perturbation : il suffit d’utiliser ALGO_INTE=”NEWTON_PERT”.

De plus, la matrice tangente est alors calculées directement à partir de la matrice jacobienne suivant la méthode décrite au § 3.2 .

En résumé, ce procédé permet, avec les deux seules routines nécessaires à l’intégration explicite, d’utiliser une intégration implicite, et de bénéficier d’un matrice tangente.

La mise à jour des catalogues concerne uniquement DEFI_MATERIAU et C_COMPORTEMENT, ainsi que le catalogue de la loi de comportement (cf. [D5.04.01]).

La mise en œuvre de cette méthode est illustrée dans le test SSNV225, pour une loi de fluage.

Bibliographie#

CROUZEIXM., MIGNOTA.L.: “Analyse numérique des équations différentielles”, Masson, 1989.
CAILLETAUD G., PILVIN P. : “Utilisation de modèles polycristallins pour le caclul par éléments finis”, Revue européenne des éléments finis, vol.3, n°4, pp. 515-541, 1994.
ZIANI M., «Accélération de la convergence des méthodes de type Newton pour la résolution des systèmes non-linéaires», thèse de doctorat, 2008.

Description des versions du document#

Version Aster	Auteur(s) Organisme(s)	Description des modifications
4	P. GEYER, E. LORENTZ, C. VOGEL (EDF/RNE/EMA, IMA/MMN)	Texte initial (intégration explicite)
10	S. GENIAUT, J.M. PROIX	Ajout de l’intégration implicite, et de la recherche linéaire.
\(11.2\)	J.M. PROIX	Ajout du §:ref:4.4 <RefNumPara__31362832>