3. Integration architecture#

3.1. Catalogs#

The laws of crystalline behavior can be used in C_ COMPORTEMENT .capy via RELATION = “MONOCRISTAL” or RELATION =” POLYCRISTAL “. Data relating to specific crystal laws are defined in sd_compor, from DEFI_COMPOR and provided under the COMPOR keyword in c_behaviour.capy.

catapy.commun/c_behavior.capy

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

fr=tr (« SD from DEFI_COMPOR « ),

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

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

fr=tr (« SD from DEFI_COMPOR « ),

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

The catalogs of the laws of behavior are:

bibpyt/behavior/monocrystal.py

from cata_behavior import LoiBehavior

law = law/Behaviour (

name = “MONOCRISTAL”,

doc = «  » « This model makes it possible to describe the behavior of a single crystal whose behavior relationships are provided via the compor concept, from DEFI_COMPOR. The number of internal variables depends on the choices made in DEFI_COMPOR; for more details see [R5.03.11] . « «  »,

num_lc = 32,

nb_vari = 0,

vari name = None,

mc_mater = None,

modeling = (“3D”, “AXIS”, “D_ PLAN “),

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

varc_name = (“TEMP”),

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

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

properties = None,)

bibpyt/behavior/polycristal.py

from cata_behavior import LoiBehavior

law = law/Behaviour (

name = “POLYCRISTAL”,

doc = « Homogenized polycrystalline behavior, defined by DEFI_COMPOR « «  »,

num_lc = 37,

nb_vari = 0,

vari name = None,

mc_mater = None,

modeling = (“3D”, “AXIS”, “D_ PLAN “),

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

varc_name = (“TEMP”),

algo_inte = (“RUNGE_KUTTA”),

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

properties = None,)

3.2. Law integration routines#

The integration of crystalline behavior laws is done as for all behavior laws at the level of each finite element integration point for non-linear options (FULL_MECA, RAPH_MECA, RIGI_MECA_TANG).

The routine LC0032 (MONOCRISTAL) is called for each integration point by the following routines:

  • if DEFORMATION =” PETIT “or =” PETIT_REAC”:

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

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

  • if DEFORMATION =” SIMO_MIEHE “:

  • in 2D or 3D, NMPL2D or NMPL3D are replaced by NMGPFI

Depending on the algorithm chosen, routine LC0032 uses either PLASTI (implicit integration by Newton’s method or a variant) or NMVPRK (integration by Runge-Kutta).

The routine LC0037 (POLYCRISTAL) is called for each integration point by the following routines:

  • if DEFORMATION =” PETIT “or =” PETIT_REAC”:

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

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

The LC0037 routine uses NMVPRK (integration by Runge-Kutta).

3.3. Recovery of material coefficients#

Whatever the type of integration chosen (implicit or explicit), the recovery of the material and behavioral characteristics, from DEFI_MATERIAU and DEFI_COMPOR, is done through the routine LCMMAT, called by LCMATE, for the MONOCRISTAL, and by the routine LCMMAP, and by the routine, also called by LCMATE, for the POLYCRISTAL.

These routines have several functions:

  • Retrieving the values of the keywords defining the parameters, mainly using the general routine RCVALB, and storing them in two tables (different only if the coefficients depend on the temperature): MATERD defining the parameters at the previous moment, i.e. at the beginning of the time step, and MATERF at the current moment, so at the end of the time step). These tables make it possible to pass the material parameters to the resolution routines;

  • Reading the sd_compor (from DEFI_COMPOR) and storing information (families of sliding systems, interaction matrix,…) in tables used during the resolution.

Architecture of LCMMAT:

LCMMJV: reading the sd_compor from DEFI_COMPOR

for each family of systems:

LCMMSG provides the number of sliding systems

LCMMJS (if it is a « user » family)

LCMAFL get the material coefficients relating to the flow

LCMHSR + LCMHDD: specific call in this routine for MONO_DD_KR

LCMAEC material coefficients relating to kinematic work hardening,

LCMAEI material coefficients relating to isotropic work hardening,

LCMHSR: calculating or reading the interaction matrix

DMAT3D, D1 MA3D: elasticity operator and its inverse

CALCMM LCMMSG.: calculation and storage of the orientation tensors of all systems, defined in the global frame of reference, to optimize performance.

In the case of a new monocrystalline behavior, it is sufficient a prima facie to intervene in the routines LCMAFL, LCMAEI and possibly LCMAEC, by adding to each of these routines the block of instructions necessary for the recovery of the material coefficients of this behavior.

It should also be given a number: in fact, to optimize performance, it is preferable to read and compare integers rather than these character strings; in integration routines, rather than testing:

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

we will test:

IF (NUCOUL .EQ.1) THEN…

The nomenclature of flow law numbers is defined in LCMAFL:

Flow law name

Associated number

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

Table 3.3-1

The numbers of isotropic work hardening laws are defined in LCMAEI:

Flow law name

Associated number

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

Table 3.3-2

The numbers of isotropic work hardening laws are defined in LCMAEC:

Flow law name

Associated number

MONO_CINE1

1

MONO_CINE2

2

Table 3.3-3

Architecture of LCMMAP:

Reading the polycristal-type sd_compor

For each single-crystal behavior (5 at most) used by all the grains

for each family of systems: reading characteristics like LCMMAT

LCMMSG provides the number of sliding systems

LCMAFL get the material coefficients relating to the flow

LCMAEC material coefficients relating to kinematic work hardening,

LCMAEI material coefficients relating to isotropic work hardening, and interaction matrix

DMAT3D, D1 MA3D: elasticity operator and its inverse.

Storage of information relating to the monocrystals used for each phase in specific tables (nbcomm, coeft/materf, cpmono, described at the end of this document).

Just before the call for explicit resolution by Runge-Kutta (routine GERPAS), call to the specific routine CALCMS allowing the sliding systems relating to all grains, defined in a global coordinate system, to be stored in a single table, in order to optimize performance.

3.4. Explicit integration — diagram of RUNGE - KUTTA#

3.4.1. Single crystal case:#

This is the fastest way (but less optimal than implicit integration, for a system of a few dozen equations) to introduce a new monocrystalline behavior in small deformations: all you have to do is write the derivatives of the internal variables in routine LCMMON, called by RDIF01.

The purpose of routine LCMMONa is to calculate the derivatives of internal variables. We must solve a system of \(6+3\mathrm{\times }{n}_{s}\) differential equations, of the type (cf. [R5.03.11]): \(\frac{\text{dY}}{\text{dt}}=F(Y,t)\), where \(Y\) represents all the internal variables: \(Y=\{\begin{array}{c}{\alpha }_{s}\\ {\gamma }_{s}\\ {p}_{s}\\ {E}^{\mathit{vp}}\end{array}\)

The system of differential equations to be solved is:

  • \(\dot{{\epsilon }^{\text{vp}}}=\sum _{s}{\mu }_{s}\dot{{\gamma }_{s}}\) 6 equations

  • for each sliding system (across all families of systems) 3 relationships:

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

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

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

where \({\tau }_{s}={\mu }_{s}\mathrm{:}\sigma\) where \(\sigma\) is inferred from the relationship: \(\sigma =\Lambda (\varepsilon –{\varepsilon }^{\mathrm{vp}})\)

In practice, the resolution is carried out in the following way:


Routine LCMMONcalcule the constraints through the elasticity relationship (isotropic or orthotropic)


CALL CALSIG (...) which provides the tensor SIG = :math:`\sigma`


Then it calculates the :math:`6+3\mathrm{\times }{n}_{s}` derivatives from the above differential equations, and stores them in a DVIN table:

  • Loop on the families of sliding systems:

DO IFA =1, NBFSYS

...

Retrieving the number of sliding systems CALL LCMMSG (NOMFAM, NBSYS, 0, PGL,, MS, NG, LG,0, Q)

Loop on the sliding systems of the IFA family:

DO IS=1, NBSYS

CALL LCMMSG (..) rereading the orientation tensor MuS

CALCUL OF THE CISSION REDUITE TAUS = SIG (I) *MuS (I) for i=1.6 \({\tau }_{s}={\mu }_{s}:\sigma\)

CALL LCMMFI (=> RP) routine for calculating isotropic work hardening \({R}_{s}(p)\)

CALL LCMMFE (=> DGAMMA, DP) flow calculation routine \(\dot{{\gamma }_{s}},\dot{{p}_{s}}\)

CALL LCMMEC (=> DALPHA) routine for calculating kinematic work hardening \(\dot{{\alpha }_{s}}\)

Calculation of the global viscoplastic deformation DO ITENS =1.6 DEVI (ITENS) = DEVI (ITENS) +MuS (ITENS) * DGAMMA \(\dot{{\varepsilon }^{\text{vp}}}\mathrm{=}\mathrm{\sum }_{s}{\mu }_{s}\dot{{\gamma }_{s}}\) ENDDO

storage of derivatives of internal variables for the IS sliding system DVIN (NUVI -2) = DALPHA DVIN (NUVI -1) = DGAMMA DVIN (NUVI) =DP

ENDDO

ENDDO

storage of the tensor derived from viscoplastic deformation.

DO ITENS =1.6

DVIN (ITENS) = DEVI (ITENS)

ENDDO

The Runge-Kutta algorithm then manages the integration of these differential equations, by controlling the error, and by refining the time step until an error less than the requested precision is obtained (RESI_INTE_RELA) [R5.03.14].

In principle, the structure of routine LCMMON should not change when adding a new monocrystalline behavior; only routines LCMMFI, LCMMFE and LCMMEC need to be modified.

They are built in the following way:

LCMMFI

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

C POUR A NOUVEAU TYPE OF ECROUISSAGE ISOTROPE, AJOUTER A 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

Note: the numbers associated with each type of behavior are used here rather than the names, to optimize performance.

Likewise the structure of routine LCMMFC is:

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

C POUR A NOUVEAU TYPE OF ECROUISSAGE CINEMATIQUE, AJOUTER A 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

And similarly, the structure of routine LCMMFC is:

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

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

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

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

IF (NUECOU .EQ.1) THEN

DP=…

DGAMMA =…

ELSEIF (NUECOU .EQ.2) THEN

ENDIF

Note: The routines LCMMFE, LCMMFI, LCMMFC are also used by the explicit integration of the polycrystal and by the implicit integration. This simplifies the implementation of a new crystalline behavior.

3.4.2. Polycrystal case:#

The integration of the polycrystal is largely based on that of the monocrystal: here again, the derivatives of the internal variables are calculated for each monocrystal of each grain \(g\), in the LCMMOP routine, called by RDIF01.

We need to solve a system of \(6+{n}_{g}(6+3\mathrm{\times }{n}_{s}(g))\) differential equations, such as:

  • for each grain defined by an orientation and a \({f}_{g}\) proportion, a constraint location relationship, of the general form:

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

  • for each of the \({n}_{s}(g)\) sliding systems of each \(g\) grain, 3 relationships:

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

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

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

  • at the grain scale, plastic deformation calculation: \(\dot{{\varepsilon }_{g}^{\text{vp}}}=\sum _{s}{\mu }_{s}\dot{{\gamma }_{s}}\)

  • calculation of macroscopic plastic deformation: \(\dot{{E}^{\text{vp}}}=\sum _{g}{f}_{g}\dot{{\epsilon }_{g}^{\text{vp}}}\)

In practice, the resolution is carried out in the following way:


Routine LCMMOPcalcule the constraints through the elasticity relationship (isotropic or orthotropic)

CALL CALSIG (...) which provides the tensor SIG (:math:`\Sigma`)


Then it calculates the :math:`6+{n}_{g}(6+3\mathrm{\times }{n}_{s}(g))` derivatives from the above differential equations, and stores them in a DVIN table:

  • Loop on the beans:

DO IGRAIN =1, NGRAIN

CALL LCLOCA () location relationship to calculate SIGG :math:`({\sigma }_{g})`

  • Loop on the families of sliding systems:

DO IFA =1, NBFSYS

...

Retrieving the number of sliding systems CALL LCMMSG (…)

Loop on the sliding systems of the IFA family: DO IS=1, NBSYS

CALL LCMMSG (..) rereading the orientation tensor MuS (IGRAIN, IS)

CALCUL OF THE CISSION REDUITE TAUS = SIGG (I) *MuS (I) for i=1.6

CALL LCMMFI (=> RP) routine for calculating isotropic work hardening

CALL LCMMFE (=> DGAMMA, DP) flow calculation routine

CALL LCMMEC (=> DALPHA) routine for calculating kinematic work hardening

Calculation of the global viscoplastic deformation DO ITENS =1.6

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

ENDDO

storage of derivatives of internal variables for the IS sliding system DVIN (NUVI -2) = DALPHA DVIN (NUVI -1) = DGAMMA DVIN (NUVI) =DP

ENDDO

ENDDO

homogenization of viscoplastic deformations

DO I=1.6

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

ENDDO

storage of the tensor derived from viscoplastic deformation.

DO ITENS =1.6

DVIN (ITENS) = DEVI (ITENS)

ENDDO

The seventh internal variable contains the cumulative equivalent viscoplastic deformation

DVIN (7) = DVINEQ

The Runge-Kutta algorithm then manages the integration of these differential equations, by controlling the error, and by refining the time step until an error less than the requested precision is obtained (RESI_INTE_RELA) [R5.03.14].

The structure of routine LCMMOP must not change when a new monocrystalline behavior is added; only routines LCMMFI, LCMMFE and LCMMEC need to be modified (which is already done in principle for the integration of the monocrystal). An additional change must be made in routine LCLOCA when adding a new location rule. Finally, for some specific post-treatments at the level of the integration point, it is necessary to intervene in routine LCDPEQ.

3.5. Implicit integration of the single crystal by Newton in PLASTI#

This time, the monocrystalline behavior is integrated by a Newton method. This method is programmed in PLASTI [R5.03.14]. It is therefore necessary to provide this algorithm with writing the system of equations to be solved in purely implicit form, in the following way: \(R(Y)=0\)

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

It is a system of \(6+6+3{n}_{s}\) nonlinear equations (of the same size as the one that is integrated by the Runge-Kutta method in explicit). But this is not the system actually solved: in order to optimize performance, we in fact solve a reduced system of equations, of size \(6+{n}_{s}\), built in the following way [R5.03.11]:

  • In the expression for the 6 components of the stress tensor, \(\Delta {E}^{\text{vp}}\) can be expressed as a function of \(\sum _{s}{\mu }_{s}\Delta {\gamma }_{s}\) so 6 equations can be eliminated from the global system to be solved.

  • Like \(\Delta {p}_{s}=\mid \Delta {\gamma }_{s}\mid\), the equation in \(\Delta {p}_{s}\) can be eliminated

*for all the crystalline behaviors currently considered, either \(\Delta {\alpha }_{s}\) is expressed directly as a function of \(\Delta {\gamma }_{s}\), in the case of kinematic work hardening, or \(\Delta {\gamma }_{s}\) is expressed as a function of \(\Delta {\alpha }_{s}\), which then represents the dislocation density (with one factor or less) for MONO_DD_* behaviors. So there is only one unknown per sliding system: \(\Delta {\beta }_{s}\), which represents either \(\Delta {\gamma }_{s}\) or \(\Delta {\alpha }_{s}\)

The following system is therefore obtained:

  • Small deformations:

\(\begin{array}{}{R}_{1}(\sigma ,\Delta \beta )={\Lambda }^{\text{-1}}\mathrm{.}\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}\) where the unknown is: \(Y=\left[\begin{array}{c}\sigma \\ \Delta \beta \end{array}\right]\)

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

  • Large deformations

\(\begin{array}{}{R}_{1}(S,\Delta \beta )={\Lambda }^{\text{-1}}\mathrm{.}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}\) where the unknown is: \(Y=\left[\begin{array}{c}S\\ \Delta \beta \end{array}\right]\)

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

And, depending on the behavior considered,

  • \(\Delta {\gamma }_{s}=\Delta {p}_{s}({\tau }_{s},\Delta {\beta }_{s}){\xi }_{s}\) and \({\xi }_{s}=\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\mathrm{ou}\frac{{\tau }_{s}-f(\alpha )}{∣{\tau }_{s}-f(\alpha )∣}\) correspond to the sign of the flow

:math:`Delta {beta }_{s}` represents either the plastic slip increment :math:`Delta {gamma }_{s}`, for laws MONO_VISC, or the variation in dislocation density \(\Delta {\omega }_{s}\) for laws MONO_DD_*

Note: the extraction of system unknowns from the internal variables (which remain at 3 per drag system) is done in routine LCAFYD. Conversely, after the resolution by NEWTON, the calculation of the 3 internal variables by sliding system according to the main variable used in the resolution is done in the LCPLNF routine.

The general form of the algorithm solved by Newton is

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

It is therefore necessary to define the initial values \(\Delta {Y}_{0}\) (0 by default, in routine LCMMIN), and calculate the residue \({Y}_{k}\), as well as the Jacobian matrix of the system: \(\frac{\text{dR}}{{\text{dY}}_{k}}\)

3.5.1. General architecture of PLASTI:#

The Newton algorithm used in PLASTI is described in [R5.03.14].

Reading material and storage coefficients in objects NBCOMM, MATERF/MATERD, CPMONO

CALL LCMATE => LCMMAT, same as RUNGE_KUTTA

Elastic prediction

CALL LCELAS

Calculating SEUIL

CALL LCCNVX => routine LCMMVX threshold assessment for MONOCRISTAL.


Calculation of the elasto-visco-plastic solution by the Newton method:

IF (SEUIL .GE. 0.D0) THEN

CALL LCPLAS/CALL LCPLNL

ENDIF


Calculation of the tangent operator:

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

CALL LCJPLC => CALL LCMMJP

ENDIF

Note: The tangent operator is automatically calculated based on the Jacobian matrix of the local system of equations [R5.03.11]. This is achieved in small and large deformations in routine LCMMJP. Therefore, in principle, there is nothing to modify for this calculation when adding a new crystalline behavior.

Routine LCCNVX makes it possible to detect if the threshold has been crossed for at least one sliding system. In the case of single crystals, it calls routine LCMMVX. Its structure is as follows:

SEUIL =0.D0

DO IFA =1, NBFSYS

DO IS=1, NBSYS

CALL LCMMFI

C ECOULEMENT VISCOPLASTIQUE

CALL LCMMFE => DP calculated from elastic prediction

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

ENDDO

ENDDO

We therefore use the same routines LCMMFE and LCMMFI as for explicit and implicit integration.

Routine LCPLNL realizes the Newton loop. Its structure (generic to the set of laws of behavior under PLASTI) is as follows:

LCPLNL

    • LCAFYD (extraction of internal variables useful to the reduced system)

    • LCINIT => LCMMIN: initialization of \(\Delta {Y}_{0}\), 0 by default

    • LCRESI => LCMMRE: residue calculation

    • Let LCJACB => LCMMJA: calculation of the Jacobian matrix

      • Let (if ALGO_INTE = NEWTON_PERT) LCJACPcalcul from the Jacobian matrix by disturbance, which calls LCRESI

    • MGAUSSrésolution

    • LCRELIrecherche linear (if ALGO_INTE = NEWTON_RELI), which calls LCRESI…

    • LCCONV => LCMMCVcritère of convergence

    • LCPLNF => LCDPECcalcul of all internal variables.

3.5.2. Residue calculation#

Routine LCMMRE calculates the residue. Its structure is as follows:

DO IFA =1, NBFSYS

DO IS=1, NBSYS

CALTAUcalcul of \({\tau }_{s}\) (small or large deformations)

LCMMLCcalcule of the quantities relating to the sliding system:

CALL LCMMFI (=> RP) routine for calculating isotropic work hardening

CALL LCMMFE (=> DGAMMA, DP) flow calculation routine

CALL LCMMEC (=> DALPHA) routine for calculating kinematic work hardening

calculating :math:`{k}_{s}` and storing it in :math:`{R}_{2}(\sigma ,\Delta \beta )=\Delta {\beta }_{s}-{k}_{s}({\tau }_{s}(\sigma ),\Delta \beta )`


In small deformations: global viscoplastic deformation

DO ITENS =1.6

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

ENDDO

In large deformations, calculate the terms needed for :math:`\Delta {F}^{p}(\Delta {\gamma }_{s})`

ENDDO

ENDDO

  • in small deformations, calculation of \({R}_{1}(\sigma ,\Delta \beta )={\Lambda }^{\text{-1}}\mathrm{.}\Delta \sigma -\Delta \varepsilon +\Delta {\varepsilon }^{\mathrm{th}}+\sum _{s}\Delta {\gamma }_{s}{\mu }_{s}=0\)

  • in large deformations

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

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

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

We therefore note that the same routines are used here again as for explicit integration: LCMMFI, LCMMFE, LCMMEC are in principle the only routines to be modified when adding a behavior, whether in small or large deformations, for the calculation of the residue, during the implicit integration.

3.5.3. Calculation of the Jacobian matrix#

Routine LCMMJAcalcule the Jacobian matrix (unless ALGO_INTE =” NEWTON_PERT “).

Its structure is as follows:

DO IFA =1, NBFSYS

DO IS=1, NBSYS

LCMMJB: calculating derived terms

LCMMJ2: calculating derived terms for MONO_DD_KR

LCMMJD: calculating derived terms for MONO_DD_CFC, MONO_DD_CC (plus _ IRRA)

LCMMJ1: calculating derived terms for MONO_VISC1, MONO_VISC2

ENDDO

ENDDO

The derivatives of these equations for calculating the Jacobian matrix can be written in general terms:

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}}\)

Table 3.5.3-1

The terms used in each sub-matrix have in common terms specific to each behavior, which are calculated in the 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}}\)

In the case of a new behavior, you must therefore either add the calculation of these derived terms into an existing LCMMJ * routine, or add a new one.

Note: for a first test, we can do without calculating the Jacobian matrix, using the automatic construction of the Jacobian matrix (by disturbance, ALGO_INTE =” NEWTON_PERT “). As a result, it is very fast to introduce a new crystalline behavior into environment PLASTI: all you have to do is calculate the residue in routine LCMMRE, called by LCRESI. On the other hand, the calculation time will only be optimized with a programmed Jacobian matrix.

3.5.4. Convergence criteria and post-treatments#

The convergence criterion is generic a prima facie, and does not depend on behavior, but can possibly be modified:

    • LCCONV => LCMMCVcritère of convergence

The last routine to be modified (possibly, if we want to calculate and add values useful for post-processing to the internal variables at the output) is:

    • LCPLNF => LCDPECqui recalculates all internal variables based on the reduced system solution. Once again, it uses the LCMMLC behavior routine.