3. Architecture de l’intégration#

3.1. Catalogues#

Les lois de comportement cristallines sont utilisables dans C_COMPORTEMENT.capy via RELATION = “MONOCRISTAL” ou RELATION=”POLYCRISTAL”. Les données relatives aux lois cristallines spécifiques sotn quant à elle définies dans sd_compor, issue de DEFI_COMPOR et fournie sous le mot-clé COMPOR dans c_comportement.capy.

catapy.commun/c_comportement.capy

b_monox = BLOC(condition = « RELATION == “MONOCRISTAL” « ,

fr=tr(« SD issue de DEFI_COMPOR »),

COMPOR =SIMP(statut=”o”,typ=compor_sdaster,max=1),),

b_polyx = BLOC(condition = « RELATION == “POLYCRISTAL” « ,

fr=tr(« SD issue de DEFI_COMPOR »),

COMPOR =SIMP(statut=”o”,typ=compor_sdaster,max=1),),

Les catalogues des lois de comportement sont:

bibpyt/Comportement/monocristal.py

from cata_comportement import LoiComportement

loi = LoiComportement(

nom = “MONOCRISTAL”,

doc = « «  »Ce modèle permet de décrire le comportement d’un monocristal dont les relations de comportement sont fournies via le concept compor, issu de DEFI_COMPOR. Le nombre de variables internes est fonction des choix effectués dans DEFI_COMPOR ; pour plus de précisions voir [R5.03.11]. » » »,

num_lc = 32,

nb_vari = 0,

nom_vari = None,

mc_mater = None,

modelisation = (“3D”,”AXIS”,”D_PLAN”),

deformation = (“PETIT”, “PETIT_REAC”, “SIMO_MIEHE”),

nom_varc = (“TEMP”),

algo_inte = (“NEWTON”,”NEWTON_RELI”,”RUNGE_KUTTA”,”NEWTON_PERT”,),

type_matr_tang = (“PERTURBATION”, “VERIFICATION”),

proprietes = None, )

bibpyt/Comportement/polycristal.py

from cata_comportement import LoiComportement

loi = LoiComportement(

nom = “POLYCRISTAL”,

doc = » » »Comportement poly-cristallin homogénéisé, défini par DEFI_COMPOR » » »,

num_lc = 37,

nb_vari = 0,

nom_vari = None,

mc_mater = None,

modelisation = (“3D”,”AXIS”,”D_PLAN”),

deformation = (“PETIT”, “PETIT_REAC”, “GROT_GDEP”),

nom_varc = (“TEMP”),

algo_inte = (“RUNGE_KUTTA”),

type_matr_tang = (“PERTURBATION”, “VERIFICATION”),

proprietes = None, )

3.2. Routines d’intégration des lois#

L’intégration des lois de comportement cristallines se fait comme pour toutes les lois de comportement au niveau de chaque point d’intégration des éléments finis pour les options non linéaires (FULL_MECA, RAPH_MECA, RIGI_MECA_TANG).

La routine LC0032 (MONOCRISTAL) est appelée pour chaque point l’intégration par les routines suivantes:

  • si DEFORMATION=”PETIT” ou =”PETIT_REAC”:

  • en 3D: TE0139/NMPL3D/NMCOMP/REDECE/LC0000/ LC0032

  • en 2D: TE0100/NMPL2D/NMCOMP/REDECE/LC0000/ LC0032

  • si DEFORMATION=”SIMO_MIEHE” :

  • en 2D ou en 3D, NMPL2D ou NMPL3D sont remplacées par NMGPFI

La routine LC0032 fait appel suivant l’algorithme choisi, soit à PLASTI (intégration implicite par la méthode de Newton ou une variante) soit à NMVPRK (intégration par Runge-Kutta).

La routine LC0037 (POLYCRISTAL) est appelée pour chaque point l’intégration par les routines suivantes:

  • si DEFORMATION=”PETIT” ou =”PETIT_REAC”:

  • en 3D: TE0139/NMPL3D/NMCOMP/REDECE/LC0000/ LC0032

  • en 2D: TE0100/NMPL2D/NMCOMP/REDECE/LC0000/ LC0032

La routine LC0037 fait appel à NMVPRK (intégration par Runge-Kutta).

3.3. Récupération des coefficients matériau#

Quelque soit le type d’intégration choisi (implicite ou explicite) la récupération des caractéristiques matériau et comportement, issues de DEFI_MATERIAU et DEFI_COMPOR, se fait par l’intermédiaire de la routine LCMMAT, appelée par LCMATE, pour le MONOCRISTAL , et par la routine LCMMAP, appelée également par LCMATE, pour le POLYCRISTAL.

Ces routines ont plusieurs fonctions :

  • Récupération des valeurs des mot-clés définissant les paramètres, principalement à l’aide de la routine générale RCVALB, et stockage dans deux tableaux (différents seulement si les coefficients dépendent de la température) : MATERD définissant les paramètres à l’instant précédent, c’est à dire au début du pas de temps, et MATERF à l’instant actuel, donc à la fin du pas de temps). Ces tableaux permettent de passer les paramètres matériau aux routines de résolution ;

  • Lecture de la sd_compor (issue de DEFI_COMPOR) et stockage des informations (familles de systèmes de glissement, matrice d’interaction, …) dans des tableaux utilisés lors de la résolution.

Architecture de LCMMAT :

LCMMJV : lecture de la sd_compor issue de DEFI_COMPOR

pour chaque famille de systèmes :

LCMMSG fournit le nombre de systèmes de glissement

LCMMJS (si il s’agit d’une famille «utilisateur» )

LCMAFL récupère les coefficients matériau relatifs à l’écoulement

LCMHSR+LCMHDD : appel spécifique dans cette routine pour MONO_DD_KR

LCMAEC coefficients matériau relatifs à l’écrouissage cinématique,

LCMAEI coefficients matériau relatifs à l’écrouissage isotrope,

LCMHSR : calcul ou lecture de la matrice d’interaction

DMAT3D, D1MA3D : opérateur d’élasticité et son inverse

CALCMM LCMMSG.: calcul et stockage des tenseurs d’orientation de tous les systèmes, définis dans le repère global, pour optimiser les performances.

Dans le cas d’un nouveau comportement monocristallin, il suffit a priori d’intervenir dans les routines LCMAFL, LCMAEI et éventuellement LCMAEC, en ajoutant dans chacune de ces routines le bloc d’instructions nécessaires à la récupération des coefficients matériau de ce comportement.

Il convient également de lui attribuer un numéro : en effet, pour optimiser les performances, il est préférable de lire et comparer des entiers plutôt que ces chaînes de caractères ; dans les routines d’intégration, plutôt que de tester :

IF (NECOUL.EQ.”MONO_VISC1”) THEN..

on testera :

IF (NUCOUL.EQ.1) THEN….

La nomenclature des numéros de lois d’écoulement est définie dans LCMAFL:

Nom de la loi d’écoulement

Numéro associé

MONO_VISC1

1

MONO_VISC2

2

MONO_DD_KR

4

MONO_DD_CFC

5

MONO_DD_CFC_IRRA

8

MONO_DD_FAT

6

MONO_DD_CC

7

MONO_DD_CC_IRRA

7

Tableau 3.3-1

Les numéros de lois d’écrouissage isotrope sont définis dans LCMAEI:

Nom de la loi d’écoulement

Numéro associé

MONO_ISOT1

1

MONO_ISOT2

2

MONO_DD_CFC

3

MONO_DD_CFC_IRRA

8

MONO_DD_FAT

4

MONO_DD_CC

7

MONO_DD_CC_IRRA

7

Tableau 3.3-2

Les numéros de lois d’écrouissage isotrope sont définis dans LCMAEC:

Nom de la loi d’écoulement

Numéro associé

MONO_CINE1

1

MONO_CINE2

2

Tableau 3.3-3

Architecture de LCMMAP :

Lecture de la sd_compor de type polycristal

Pour chaque comportement monocristal (5 au maximum) utilisé par l’ensemble des grains

pour chaque famille de systèmes : lecture des caractéristiques comme LCMMAT

LCMMSG fournit le nombre de systèmes de glissement

LCMAFL récupère les coefficients matériau relatifs à l’écoulement

LCMAEC coefficients matériau relatifs à l’écrouissage cinématique,

LCMAEI coefficients matériau relatifs à l’écrouissage isotrope, et matrice d’interaction

DMAT3D, D1MA3D : opérateur d’élasticité et son inverse.

Stockage des informations relatives aux monocristaux utilisés pour chaque phase dans des tableaux spécifiques (nbcomm, coeft/materf,cpmono, décrits en fin de ce document).

Juste avant l’appel à la résolution explicite par Runge-Kutta (routine GERPAS), appel à la routine spécifique CALCMS permettant de stocker dans un unique tableau les systèmes de glissement relatifs à tous les grains, définis en repère global, pour optimiser les performances.

3.4. Intégration explicite – schéma de RUNGE - KUTTA#

3.4.1. Cas du monocristal :#

C’est la façon la plus rapide (mais moins optimale que l’intégration implicite, pour un système de quelques dizaines d’équations) d’introduire un nouveau comportement monocristallin en petites déformations : il suffit de d’écrire les dérivées des variables internes dans la routine LCMMON, appelée par RDIF01.

La routine LCMMONa pour but de calculer les dérivées des variables internes. On doit résoudre un système de \(6+3\times {n}_{s}\) équations différentielles, du type (cf. [R5.03.11]): \(\frac{\text{dY}}{\text{dt}}=F(Y,t)\) , où \(Y\) représente l’ensemble des variables internes: \(Y=\lbrace \begin{array}{c}{\alpha}_{s}\\ {\gamma}_{s}\\ {p}_{s}\\ {E}^{\mathit{vp}}\end{array}\)

Le système d’équations différentielles à résoudre est:

  • \(\dot{{\epsilon}^{\text{vp}}}=\sum_{s}{\mu}_{s}\dot{{\gamma}_{s}}\) 6 équations

  • pour chaque système de glissement (sur l’ensemble des familles de systèmes) 3 relations :

      • \(\dot{{\gamma}_{s}}=\dot{{p}_{s}}\left({\tau}_{s}(\sigma ),{\alpha}_{s},{\gamma}_{s},{R}_{s}(p)\right)\eta ({\tau}_{s},{\alpha}_{s})\) où \(\eta ({\tau}_{s},{\alpha}_{s})=\pm 1\)

      • \(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

      • \({R}_{s}(p)\)

où \({\tau}_{s}={\mu}_{s}:\sigma\) où \(\sigma\) est déduite de la relation : \(\sigma =\Lambda (\varepsilon –{\varepsilon}^{\mathrm{vp}})\)

En pratique, la résolution s'effectue de la façon suivante :


La routine LCMMONcalcule les contraintes par la relation d'élasticité (isotrope ou orthotrope)


CALL CALSIG(...) ce qui fournit le tenseur SIG = :math:`\sigma`


Puis elle calcule les :math:`6+3\times {n}_{s}`dérivées issues des équations différentielles ci-dessus, et les stocke dans un tableau DVIN :

  • Boucle sur les familles de systèmes de glissement :

DO IFA=1,NBFSYS

…

Récupération du nombre de systèmes de glissement CALL LCMMSG(NOMFAM,NBSYS,0,PGL,MS,NG,LG,0,Q)

Boucle sur les systèmes de glissement de la famille IFA :

DO IS=1,NBSYS

CALL LCMMSG(..) relecture du tenseur d’orientation MuS

CALCUL DE LA CISSION REDUITE TAUS= SIG(I)*MuS(I) pour i=1,6 \({\tau}_{s}={\mu}_{s}:\sigma\)

CALL LCMMFI( => RP) routine de calcul de l’écrouissage isotrope:math:{R}_{s}(p)

CALL LCMMFE(=> DGAMMA, DP) routine de calcul de l’écoulement:math:dot{{gamma}_{s}},dot{{p}_{s}}

CALL LCMMEC( => DALPHA)routine de calcul de l’écrouissage cinématique \(\dot{{\alpha}_{s}}\)

Calcul de la déformation viscoplastique globale DO ITENS=1,6 DEVI(ITENS)=DEVI(ITENS)+MuS(ITENS)*DGAMMA \(\dot{{\varepsilon}^{\text{vp}}}=\sum_{s}{\mu}_{s}\dot{{\gamma}_{s}}\) ENDDO

stockage des dérivées des variables internes pour le système de glissement IS DVIN(NUVI-2)=DALPHA DVIN(NUVI-1)=DGAMMA DVIN(NUVI )=DP

ENDDO

ENDDO

stockage du tenseur dérivée de la déformation viscoplastique.

DO ITENS=1,6

DVIN(ITENS)= DEVI(ITENS)

ENDDO

L’algorithme de Runge-Kutta gère alors l’intégration de ces équations différentielles, en contrôlant l’erreur, et en raffinant le pas de temps jusqu’à obtenir une erreur inférieure à la précision demandée (RESI_INTE_RELA) [R5.03.14].

A priori, la structure de la routine LCMMON, ne doit pas évoluer lors de l’ajout d’un nouveau comportement monocristallin ; seules les routinesLCMMFI, LCMMFE et LCMMEC sont à modifier.

Elles sont construites de la façon suivante :

LCMMFI

C——————————————————————–

C POUR UN NOUVEAU TYPE D’ECROUISSAGE ISOTROPE, AJOUTER UN BLOC IF

C——————————————————————–

C IF (NECRIS.EQ.”MONO_ISOT1”) THEN

IF (NUEISO.EQ.1) THEN

……

RP=…

C ELSEIF (NECRIS.EQ.”MONO_ISOT2”) THEN

ELSEIF (NUEISO.EQ.2) THEN

…………

RP=…

C ELSEIF (NECRIS.EQ.”MONO_DD_CFC”) THEN

ELSEIF (NUEISO.EQ.3) THEN

……….

RP=MU*SQRT(RP)*CEFF

ENDIF

Remarque : on utilise ici les numéros associés à chaque type de comportement plutôt que les noms, pour optimiser les performances.

De même la structure de la routine LCMMFC est :

C———————————————————————-

C POUR UN NOUVEAU TYPE D’ECROUISSAGE CINEMATIQUE, AJOUTER UN BLOC IF

C———————————————————————-

C IF (NECRCI.EQ.”MONO_CINE1”) THEN

IF (NUECIN.EQ.1) THEN

DALPHA=…

C ELSEIF (NECRCI.EQ.”MONO_CINE2”) THEN

ELSEIF (NUECIN.EQ.2) THEN

ENDIF

Et de façon similaire, la structure de la routine LCMMFC est :

C————————————————————-

C POUR UN NOUVEAU TYPE D’ECOULEMENT, CREER UN BLOC IF

C————————————————————

C IF (NECOUL.EQ.”MONO_VISC1”) THEN

IF (NUECOU.EQ.1) THEN

DP=…

DGAMMA= …

ELSEIF (NUECOU.EQ.2) THEN

ENDIF

Remarque : Les routines LCMMFE, LCMMFI, LCMMFC sont également utilisées par l’intégration explicite du polycristal et par l’intégration implicite. Ceci simplifie la mise en œuvre d’un nouveau comportement cristallin.

3.4.2. Cas du polycristal :#

L’intégration du polycristal s’appuie largement sur celle du monocristal : on calcule là encore les dérivées des variables internes pour chaque monocristal de chaque grain \(g\) , dans la routine LCMMOP, appelée par RDIF01.

On doit résoudre un système de \(6+{n}_{g}(6+3\times {n}_{s}(g))\) équations différentielles, du type :

  • pour chaque grain défini par une orientation et une proportion \({f}_{g}\) , une relation de localisation des contraintes, de la forme générale :

\({\sigma}_{g}=L(\Sigma ,{E}^{\text{vp}},{\varepsilon}_{g}^{\text{vp}},{\beta}_{g})\) avec, \(\Sigma =\Lambda ({\Lambda}_{-}^{-1}){\Sigma}^{-}+\Lambda (\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\)

  • pour chacun des \({n}_{s}(g)\) systèmes de glissement de chaque grain \(g\) , 3 relations :

      • \(\dot{{\gamma}_{s}}=\dot{{p}_{s}}({\tau}_{s}(\sigma ),{\alpha}_{s},{\gamma}_{s},{R}_{s}(p))\eta ({\tau}_{s},{\alpha}_{s})\) où \(\eta ({\tau}_{s},{\alpha}_{s})=\pm 1\)

      • \(\dot{{\alpha}_{s}}=h({\tau}_{s},{\alpha}_{s},{\gamma}_{s},{p}_{s})\)

      • \({R}_{s}(p)\)

  • à l’échelle du grain, lcalcul de la déformation plastique : \(\dot{{\varepsilon}_{g}^{\text{vp}}}=\sum_{s}{\mu}_{s}\dot{{\gamma}_{s}}\)

  • calcul de la déformation plastique macroscopique : \(\dot{{E}^{\text{vp}}}=\sum_{g}{f}_{g}\dot{{\epsilon}_{g}^{\text{vp}}}\)

En pratique, la résolution s'effectue de la façon suivante :


La routine LCMMOPcalcule les contraintes par la relation d'élasticité (isotrope ou orthotrope)

CALL CALSIG(...) ce qui fournit le tenseur SIG (:math:`\Sigma`)


Puis elle calcule les :math:`6+{n}_{g}(6+3\times {n}_{s}(g))`dérivées issues des équations différentielles ci-dessus, et les stocke dans un tableau DVIN :

  • Boucle sur les grains :

DO IGRAIN=1,NGRAIN

CALL LCLOCA() relation de localisation permettant de calculer SIGG :math:`({\sigma}_{g})`

  • Boucle sur les familles de systèmes de glissement :

DO IFA=1,NBFSYS

…

Récupération du nombre de systèmes de glissement CALL LCMMSG(…)

Boucle sur les systèmes de glissement de la famille IFA : DO IS=1,NBSYS

CALL LCMMSG(..) relecture du tenseur d’orientation MuS(IGRAIN, IS)

CALCUL DE LA CISSION REDUITE TAUS= SIGG(I)*MuS(I) pour i=1,6

CALL LCMMFI( => RP) routine de calcul de l’écrouissage isotrope

CALL LCMMFE(=> DGAMMA, DP) routine de calcul de l’écoulement

CALL LCMMEC( => DALPHA )routine de calcul de l’écrouissage cinématique

Calcul de la déformation viscoplastique globale DO ITENS=1,6

c DEVG(ITENS)=DEVG(ITENS)+MuS(ITENS)*DGAMMA \((\text{DEVG}={\varepsilon}_{g}^{\mathit{vp}})\)

ENDDO

stockage des dérivées des variables internes pour le système de glissement IS DVIN(NUVI-2)=DALPHA DVIN(NUVI-1)=DGAMMA DVIN(NUVI )=DP

ENDDO

ENDDO

homogénéisation des déformations viscoplastiques

DO I=1,6

DEVI(I)=DEVI(I)+FV*DEVG(I)

ENDDO

stockage du tenseur dérivée de la déformation viscoplastique.

DO ITENS=1,6

DVIN(ITENS)= DEVI(ITENS)

ENDDO

La septième variable interne contient la déformation viscoplastique équivalente cumulée

DVIN(7)= DVINEQ

L’algorithme de Runge-Kutta gère alors l’intégration de ces équations différentielles, en contrôlant l’erreur, et en raffinant le pas de temps jusqu’à obtenir une erreur inférieure à la précision demandée (RESI_INTE_RELA) [R5.03.14].

La structure de la routine LCMMOP, ne doit pas évoluer lors de l’ajout d’un nouveau comportement monocristallin ; seules les routinesLCMMFI, LCMMFE et LCMMEC sont à modifier (ce qui est déjà fait en principe pour l’intégration du monocristal). Une modification supplémentaire est à effectuer dans la routine LCLOCA lors de l’ajout d’une nouvelle règle de localisation. Enfin, pour certains post-traitements spécifiques au niveau du point d’intégration, il faut intervenir dans la routine LCDPEQ.

3.5. Intégration implicite du monocristal par Newton dans PLASTI#

On intègre cette fois le comportement monocristallin par une méthode de Newton. Cette méthode est programmée dans PLASTI [R5.03.14]. Il faut donc fournir à cet algorithme écrire le système d’équations à résoudre sous forme purement implicite, de la façon suivante : \(R(Y)=0\)

\(R(Y)=\lbrace \begin{array}{}{\Lambda}^{\text{-1}}\Sigma -({\Lambda}_{-}^{\text{- 1}}){\Sigma}^{-}-(\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\\ \Delta {E}^{\text{vp}}-{\sum}_{s}{\mu}_{s}\Delta {\gamma}_{s}\\ {n}_{s}\left\lbrace \begin{array}{}\Delta {\alpha}_{s}-h({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ \Delta {\gamma}_{s}-g({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\\ {\mathrm{\Delta p}}_{s}-f({\tau}_{s}^{+},{\alpha}_{s}^{+},{\gamma}_{s}^{+},{p}_{s}^{+})\end{array}\right\rbrace \end{array}=0\)

C’est un système de \(6+6+3{n}_{s}\) équations non linéaires (de même taille que celui qui est intégré par la méthode de Runge-Kutta en explicite). Mais ce n’est pas le système effectivement résolu: afin d’optimiser les performances, on résout en fait un système d’équations réduit, de taille \(6+{n}_{s}\) , construit de la façon suivante [R5.03.11] :

  • Dans l’expression des 6 composantes du tenseur des contraintes, \(\Delta {E}^{\text{vp}}\) peut être exprimée en fonction de \(\sum_{s}{\mu}_{s}\Delta {\gamma}_{s}\) donc 6 équations peuvent être éliminées du système global à résoudre.

  • Comme \(\Delta {p}_{s}=\mid \Delta {\gamma}_{s}\mid\) , l’équation en \(\Delta {p}_{s}\) peut être éliminée

  • pour tous les comportements cristallins considérés actuellement, soit \(\Delta {\alpha}_{s}\) s’exprime directement en fonction de \(\Delta {\gamma}_{s}\) , dans le cas de l’écrouissage cinématique, soit \(\Delta {\gamma}_{s}\) s’exprime en fonction de \(\Delta {\alpha}_{s}\) qui représente alors la densité de dislocation (à un facteur près) pour les comportements de type MONO_DD_*. Il y a donc une seule inconnue par système de glissement: \(\Delta {\beta}_{s}\) , qui représente soit \(\Delta {\gamma}_{s}\) , soit \(\Delta {\alpha}_{s}\)

On obtient donc le système suivant :

  • Petites déformations:

\(\begin{array}{}{R}_{1}(\sigma ,\Delta \beta )={\Lambda}^{\text{-1}}.\Delta \sigma -\Delta \varepsilon +\Delta {\varepsilon}^{\mathrm{th}}+\sum_{s}\Delta {\gamma}_{s}{\mu}_{s}=0\\ {R}_{2}(\sigma ,\Delta \beta )=\Delta {\beta}_{s}-{k}_{s}({\tau}_{s}(\sigma ),\Delta \beta )=0\text{}\end{array}\) où l’inconnue est : \(Y=\left[\begin{array}{c}\sigma \\ \Delta \beta \end{array}\right]\)

avec \({\tau}_{s}(\sigma )=\sigma :{\mu}_{s}\)

  • Grandes déformations

\(\begin{array}{}{R}_{1}(S,\Delta \beta )={\Lambda}^{\text{-1}}.S-\frac{1}{2}({F}^{\text{eT}}{F}^{\text{e}}-{I}_{d})=0\\ {R}_{2}(S,\Delta \beta )=\Delta {\beta}_{s}-{k}_{s}({\tau}_{s}(S),\Delta \beta )=0\text{}\end{array}\) où l’inconnue est : \(Y=\left[\begin{array}{c}S\\ \Delta \beta \end{array}\right]\)

avec \({\tau}_{s}(S)=\left[(2{\Lambda}^{\text{-1}}S+{I}_{d})S\right]:{m}_{s}\otimes {n}_{s}\) et \({F}_{n+1}^{e}=\Delta F{F}_{n}^{e}{(\Delta {F}^{p}(\Delta {\gamma}_{s}))}^{\text{-1}}\)

Et, suivant le comportement considéré,

  • \(\Delta {\gamma}_{s}=\Delta {p}_{s}({\tau}_{s},\Delta {\beta}_{s}){\xi}_{s}\) et \({\xi}_{s}=\frac{{\tau}_{s}}{∣{\tau}_{s}∣}\mathrm{ou}\frac{{\tau}_{s}-f(\alpha )}{∣{\tau}_{s}-f(\alpha )∣}\) correspond au signe de l’écoulement

  • \(\Delta {\beta}_{s}\) représente soit l’incrément de glissement plastique \(\Delta {\gamma}_{s}\) , pour les lois MONO_VISC* , soit la variation de densité de dislocations \(\Delta {\omega}_{s}\) pour les lois MONO_DD_*

Remarque : l’extraction des inconnues du système à partir des variables internes (qui restent au nombre de 3 par système de glissement) se fait dans la routine LCAFYD. Inversement, après la résolution par NEWTON, le calcul des 3 variables internes par système de glissement en fonction de la variable principale utilisée dans la résolution se fait dans la routine LCPLNF.

La forme générale de l’algorithme résolu par Newton est

\({Y}_{k+1}={Y}_{k}-(\frac{\text{dR}}{{\text{dY}}_{k}}{)}^{-1}R({Y}_{k})\)

Il faut donc définir les valeurs initiales \(\Delta {Y}_{0}\) (0 par défaut, dans la routine LCMMIN), et calculer le résidu \({Y}_{k}\) ,ainsi que lamatrice jacobienne du système: \(\frac{\text{dR}}{{\text{dY}}_{k}}\)

3.5.1. Architecture générale de PLASTI:#

L’algorithme de Newton utilisé dans PLASTI est décrit en [R5.03.14].

Lecture des coefficients matériau et stockage dans les objets NBCOMM, MATERF/MATERD, CPMONO

CALL LCMATE => LCMMAT, identique à RUNGE_KUTTA

Prédiction élastique

CALL LCELAS

Calcul du SEUIL

CALL LCCNVX => routine LCMMVX évaluation du seuil pour MONOCRISTAL.


Calcul de la solution élasto-visco-plastique par la méthode de Newton :

IF ( SEUIL .GE. 0.D0 ) THEN

CALL LCPLAS / CALL LCPLNL

ENDIF


Calcul de l'opérateur tangent :

IF ( OPT .EQ. 'RIGI_MECA_TANG' .OR. OPT .EQ. 'FULL_MECA' ) THEN

CALL LCJPLC => CALL LCMMJP

ENDIF

Remarque : L’opérateur tangent est calculé automatiquement en fonction de la matrice jacobienne du système d’équations local [R5.03.11]. Ceci est réalisé en petites et en grandes déformations dans la routine LCMMJP. Il n’y a donc a priori rien à modifier pour ce calcul lors de l’ajout d’un nouveau comportement cristallin.

La routine LCCNVX permet de détecter si le seuil est franchi pour au moins un système de glissement. Elle appelle dans le cas du monocristal la routine LCMMVX . Sa structure est la suivante :

SEUIL=0.D0

DO IFA=1,NBFSYS

DO IS=1,NBSYS

CALL LCMMFI

C ECOULEMENT VISCOPLASTIQUE

CALL LCMMFE => DP calculé à partir de la prédiction élastique

IF (DP.GT.0.D0) SEUIL=1.D0

ENDDO

ENDDO

On utilise donc les mêmes routines LCMMFE et LCMMFI que pour l’intégration explicite et implicite.

La routine LCPLNL réalise la boucle de Newton . Sa structure (générique à l’ensemble des lois de comportement sous PLASTI) est la suivante :

LCPLNL

    • LCAFYD(extraction des variables internes utiles au système réduit)

    • LCINIT=> LCMMIN: initialisation de \(\Delta {Y}_{0}\) , 0 par défaut

    • LCRESI=> LCMMRE: calcul du résidu

    • Soit LCJACB=> LCMMJA : calcul de la matrice jacobienne

      • Soit (si ALGO_INTE=NEWTON_PERT) LCJACPcalcul de la matrice jacobienne par perturbation, qui appelle LCRESI

    • MGAUSSrésolution

    • LCRELIrecherche linéaire (si ALGO_INTE=NEWTON_RELI), qui appelleLCRESI…

    • LCCONV=> LCMMCVcritère de convergence

    • LCPLNF=> LCDPECcalcul des toutes les variables internes.

3.5.2. Calcul du résidu#

La routine LCMMRE calcule le résidu. Sa structure est la suivante :

DO IFA=1,NBFSYS

DO IS=1,NBSYS

CALTAUcalcul de \({\tau}_{s}\) (petites ou grandes déformations)

LCMMLCcalcule des quantités relatives au système de glissement :

CALL LCMMFI( => RP)        routine de calcul de l'écrouissage isotrope

CALL LCMMFE( => DGAMMA, DP)        routine de calcul de l'écoulement

CALL LCMMEC( => DALPHA )routine de calcul de l'écrouissage cinématique

calcul de :math:`{k}_{s}`et stockage dans :math:`{R}_{2}(\sigma ,\Delta \beta )=\Delta {\beta}_{s}-{k}_{s}({\tau}_{s}(\sigma ),\Delta \beta )`


En petites déformations: déformation viscoplastique globale

DO ITENS=1,6

DEVI(ITENS)=DEVI(ITENS)+MuS(ITENS)*DGAMMA

ENDDO

En grandes déformations, calcul des termes nécessaires à :math:`\Delta {F}^{p}(\Delta {\gamma}_{s})`

ENDDO

ENDDO

  • en petites déformations, calcul de \({R}_{1}(\sigma ,\Delta \beta )={\Lambda}^{\text{-1}}.\Delta \sigma -\Delta \varepsilon +\Delta {\varepsilon}^{\mathrm{th}}+\sum_{s}\Delta {\gamma}_{s}{\mu}_{s}=0\)

  • en grandes déformations

CALCFEcalcul de \({F}_{n+1}^{e}=\Delta F{F}_{n}^{e}{(\Delta {F}^{p}(\Delta {\gamma}_{s}))}^{\text{-1}}\)

LCGRLAcalcul de \({E}_{\mathrm{GL}}^{e}=\frac{1}{2}({{F}^{e}}^{T}{F}^{e}–{I}_{d})\) \(S=\Lambda :{E}_{\mathrm{GL}}^{e}\)

et \({R}_{1}(S,\Delta \beta )={\Lambda}^{\text{-1}}.S-\frac{1}{2}({F}^{\text{eT}}{F}^{\text{e}}-{I}_{d})\)

On constate donc que l’on utilise là encore les même routines que pour l’intégration explicite : LCMMFI, LCMMFE, LCMMEC sont a priori les seules routines à modifier lors de l’ajout d’un comportement, que ce soit en petites ou en grandes déformations, pour le calcul du résidu, lors de l’intégration implicite.

3.5.3. Calcul de la matrice jacobienne#

La routine LCMMJAcalcule la matrice jacobienne (sauf si ALGO_INTE=”NEWTON_PERT”).

Sa structure est la suivante :

DO IFA=1,NBFSYS

DO IS=1,NBSYS

LCMMJB: calcul des termes dérivés

LCMMJ2: calcul des termes dérivés pour MONO_DD_KR

LCMMJD: calcul des termes dérivés pour MONO_DD_CFC, MONO_DD_CC(plus _IRRA)

LCMMJ1: calcul des termes dérivés pour MONO_VISC1, MONO_VISC2

ENDDO

ENDDO

Les dérivées de ces équations pour le calcul de la matrice jacobienne peuvent être écrites de façon générale :

HPP

GDEF

\({J}_{11}=\frac{\partial {R}_{1}{(\sigma ,\Delta \beta )}_{i}}{\partial {\sigma}_{j}}\)

\({J}_{12}=\frac{\partial {R}_{1}{(\sigma ,\Delta \beta )}_{i}}{\partial {\beta}_{s}}\)

\({J}_{11}=\frac{\partial {R}_{1}{(S,\Delta \beta )}_{i}}{\partial {S}_{j}}\)

\({J}_{12}=\frac{\partial {R}_{1}{(S,\Delta \beta )}_{i}}{\partial {\beta}_{s}}\)

\({J}_{21}=\frac{\partial {R}_{2}{(\sigma ,\Delta \beta )}_{s}}{\partial {\sigma}_{i}}\)

\({J}_{22}=\frac{\partial {R}_{2}{(\sigma ,\Delta \beta )}_{s}}{\partial {\beta}_{r}}\)

\({J}_{21}=\frac{\partial {R}_{2}{(S,\Delta \beta )}_{s}}{\partial {S}_{i}}\)

\({J}_{22}=\frac{\partial {R}_{2}{(S,\Delta \beta )}_{s}}{\partial {\beta}_{r}}\)

Tableau 3.5.3-1

Les termes intervenant dans chaque sous-matrice ont en commun des termes spécifiques à chaque comportement, qui sont calculés dans les routines LCMMJ* [R5,03,11 annexe 5] :

  • \(\frac{\partial \Delta {\gamma}_{s}}{\partial {\tau}_{s}}=\frac{\partial \Delta {p}_{s}}{\partial {\tau}_{s}}{\xi}_{s}\) ,

  • \(\frac{\partial \Delta {\gamma}_{r}}{\partial \Delta {\beta}_{s}}=\frac{\partial \Delta {p}_{r}}{\partial \Delta {\beta}_{s}}{\xi}_{r}\) ,

  • \(\frac{\partial {k}_{s}}{\partial {\tau}_{s}}\)

  • \(\frac{\partial {k}_{r}}{\partial \Delta {\beta}_{s}}\)

Dans le cas d’un nouveau comportement, il faut donc soit ajouter le calcul de ces termes dérivés dans une routine LCMMJ* existante, soit en ajouter une nouvelle.

Remarque : pour un premier test, on peut se passer du calcul de la matrice jacobienne, en utilisant la construction automatique de la matrice jacobienne (par perturbation, ALGO_INTE=”NEWTON_PERT”). De ce fait, il est très rapide d’introduire un nouveau comportement cristallin dans l’environnement PLASTI: il suffit de calculer le résidu dans la routine LCMMRE, appelée par LCRESI. Par contre, le temps calcul ne sera optimisé qu’avec une matrice jacobienne programmée.

3.5.4. Critère de convergence et post-traitements#

Le critère de convergence est générique a priori, et ne dépend pas du comportement, mais peut être éventuellement modifié :

    • LCCONV=> LCMMCVcritère de convergence

La dernière routine à modifier (éventuellement, si on veut calculer et ajouter aux variables internes en sortie des valeurs utiles au post-traitement) est :

    • LCPLNF=> LCDPECqui recalcule toutes les variables internes à partir de la solution du système réduit. Elle fait appel encore une fois à la routine de comportement LCMMLC.