3. Local integration and digital implementation#

3.1. System of equations to be solved#

3.1.1. Behavior type MONOCRISTAL#

The local behavior of the monocrystal is defined, at a given moment of discretization in time and at the level of an integration point of a finite element, by the data:

the macroscopic stress tensor at the previous moment \(\Sigma ({t}_{i\mathrm{-}1})\mathrm{=}{\Sigma }^{\mathrm{-}}\),
internal variables at the previous moment, for each sliding system: \({\alpha }_{s}({t}_{i\mathrm{-}1}),{\gamma }_{s}({t}_{i\mathrm{-}1}),{p}_{s}({t}_{i\mathrm{-}1})\),
and the total deformation increase tensor provided by iteration n of the global resolution algorithm \(\Delta E=\Delta {E}_{i}^{n}\) (with the notations of [R5.03.01]).

Integration consists in finding:

the macroscopic stress tensor \(\Sigma =\Sigma ({t}_{i})\)
and the internal variables \({\alpha }_{s}\mathrm{=}{\alpha }_{s}({t}_{i})\), \({\gamma }_{s}\mathrm{=}{\gamma }_{s}({t}_{i})\), \({p}_{s}\mathrm{=}{p}_{s}({t}_{i})\)

verifying the equations of behavior in each sliding system (which are one-dimensional relationships), and the relationships of passage between the macroscopic tensors and all the sliding directions. Notation: equations are written in discretized form in such a way as:

explicit, if the quantities noted \({A}^{\text{+/ -}}\) are evaluated at the moment \({t}_{i\mathrm{-}1}\): \({A}^{\text{+/ -}}\mathrm{=}{A}^{\mathrm{-}}\mathrm{=}A({t}_{i\mathrm{-}1})\)
implicit, if they are evaluated at time \({t}_{i}\): \({A}^{\text{+/ -}}\mathrm{=}{A}^{+}\mathrm{=}A({t}_{i})\)

The equations to be integrated can be put in the following general form:

Given, at a Gauss point, the tensors:

\(\Delta E\): deformation variation at time \({t}_{i}\),

\(E({t}_{i-1})={E}^{-}\): deformation at time \({t}_{i-1}\),

\(\Sigma ({t}_{i\mathrm{-}1})\mathrm{=}{\Sigma }^{\mathrm{-}}\): macroscopic constraint at time \({t}_{i-1}\),

\({\alpha }_{s}({t}_{i\mathrm{-}1}),{\gamma }_{s}({t}_{i\mathrm{-}1}),{p}_{s}({t}_{i\mathrm{-}1})\): internal variables for each sliding system at \({t}_{i-1}\),

You have to find:

\(\Sigma \mathrm{=}\Sigma ({t}_{i})\): macroscopic constraint at time \({t}_{i}\), in the coordinate system corresponding to the global orientation

\({\alpha }_{s}\mathrm{=}{\alpha }_{s}({t}_{i})\)

\({\gamma }_{s}\mathrm{=}{\gamma }_{s}({t}_{i})\)

\({p}_{s}\mathrm{=}{p}_{s}({t}_{i})\)

verifying:

\({D}^{-1}\Sigma =({D}_{-}^{-1}){\Sigma }^{-}+(\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\), where \(\Lambda\) may depend on temperature, and may correspond to orthotropic elasticity.

\(\Delta {E}^{\text{vp}}=\sum _{s}{\mu }_{s}\Delta {\gamma }_{s}\)

for each sliding system (of all the families of systems):

\({n}_{s}\) flow relationships: either in viscoplasticity \(\Delta {\gamma }_{s}\mathrm{=}g({\tau }_{s}^{\text{+/ -}},{\alpha }_{s}^{\text{+/ -}},{\gamma }_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

with \(\Delta {p}_{s}\mathrm{=}f({\tau }_{s}^{\text{+/ -}},{\alpha }_{s}^{\text{+/ -}},{\gamma }_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

or in plasticity \(F({\tau }_{s}^{\text{+/ -}},{\alpha }_{s}^{\text{+/ -}},{\gamma }_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\mathrm{\le }0\), \(F\Delta {p}_{s}\mathrm{=}0\),

with \(\Delta {p}_{s}\mathrm{=}\mathrm{\mid }\Delta {\gamma }_{s}\mathrm{\mid }\)

\({n}_{s}\) evolution equations for kinematic work hardening: \(\Delta {\alpha }_{s}\mathrm{=}h({\tau }_{s}^{\text{+/ -}},{\alpha }_{s}^{\text{+/ -}},{\gamma }_{s}^{\text{+/ -}},{p}_{s}^{\text{+/ -}})\)

\({n}_{s}\) evolution equations for isotropic work hardening: \({R}_{s}({p}_{s}^{\text{+/ -}})\)

This is resolved either explicitly (Runge_Kutta) or implicitly (Newton).

3.1.2. Behavior type POLYCRISTAL#

Discretized behavioral relationships are:

Given (at a Gauss point) the global tensors:

increase in total deformation \(\Delta E\),
total deformation at the previous moment, \(E({t}_{i-1})={E}^{-}\)
constrained at the previous moment: \(\Sigma (t{}_{i-1}\text{})={\Sigma }^{-},\)
all the internal variables \({\alpha }_{s}^{\mathrm{-}},{\gamma }_{s}^{\mathrm{-}},{p}_{s}^{\mathrm{-}}\) relating to the sliding systems of each phase,
the behavior parameters of each phase, the orientations and volume fractions of each phase, and the type of location method.

You need to find \(\Sigma =\Sigma ({t}_{i}),{\alpha }_{s}=\alpha {}_{s}\text{}({t}_{i}),{\gamma }_{s}={\gamma }_{s}({t}_{i}),{p}_{s}={p}_{s}({t}_{i})\) verifying:

at the level of the Gauss point: \(\Sigma =\Lambda ({\Lambda }_{-}^{-1}){\Sigma }^{-}+\Lambda (\Delta E-\Delta {E}^{\text{th}}-\Delta {E}^{\text{vp}})\), in the global coordinate system,
for each phase (or « grain »), defined by an orientation and a proportion \({f}_{g}\), a relationship for the location of the constraints, of the general form (to be expressed in the local coordinate system of each phase)

\({\sigma }_{g}=L(\Sigma ,{E}^{\text{vp}},{\varepsilon }_{g}^{\text{vp}},{\beta }_{g})\)

and for each sliding system in each phase:
\(\Delta {\varepsilon }_{{g}^{\text{vp}}}\mathrm{=}\mathrm{\sum }_{s}{m}_{s}\Delta {\gamma }_{s}\)
\({n}_{s}\) equations: \({\tau }^{s}={\Sigma }_{\text{ij}}:{\mu }_{\text{ij}}^{s}\)
\({n}_{s}\) flow relationships: \(\Delta {\gamma }_{s}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\) with \(\Delta {p}_{s}=\mid \Delta {\gamma }_{s}\mid\)
\({n}_{s}\) evolutions of work hardening: \(\Delta {\alpha }_{s}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)
\(F({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\le \mathrm{0,}F\cdot \Delta {p}_{s}=0\), (in plasticity independent of time)

The homogenization equations are then left: \(\Delta {E}^{\text{vp}}=\sum _{g}{f}_{g}\Delta {\varepsilon }_{g}^{\text{vp}}\)

The viscoplastic behaviors relating to each sliding system are identical to the case of the microstructure.

In the current version of Code_Aster, these behavioral relationships are only explicitly integrated.

3.2. Implicit behavior resolution MONOCRISTAL#

3.2.1. Implicit resolution of phenomenological laws#

A system of the following general form must be solved:

\(R(Y)=\{\begin{array}{}s(\Sigma ,\Delta {E}^{\text{vp}},{\alpha }_{s}^{\text{+}},{\gamma }_{s}^{\text{+}},{p}_{s}^{\text{+}})\\ e(\Sigma ,\Delta {E}^{\text{vp}},{\alpha }_{s}^{\text{+}},{\gamma }_{s}^{\text{+}},{p}_{s}^{\text{+}})\\ {n}_{s}\left\{\begin{array}{}a(\Sigma ,\Delta {E}^{\text{vp}},{\alpha }_{s}^{\text{+}},{\gamma }_{s}^{\text{+}},{p}_{s}^{\text{+}})\\ g(\Sigma ,\Delta {E}^{\text{vp}},{\alpha }_{s}^{\text{+}},{\gamma }_{s}^{\text{+}},{p}_{s}^{\text{+}})\\ p(\Sigma ,\Delta {E}^{\text{vp}},{\alpha }_{s}^{\text{+}},{\gamma }_{s}^{\text{+}},{p}_{s}^{\text{+}})\end{array}\right\}\end{array}\) \(=\{\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

In a more contracted manner, we pose:

\(R(\Delta y)\mathrm{=}0\mathrm{=}\left[\begin{array}{c}s(\Delta y)\\ e(\Delta y)\\ a(\Delta y)\\ g(\Delta y)\\ p(\Delta y)\end{array}\right]\) with \(\Delta y\mathrm{=}(\begin{array}{c}\Delta \Sigma \\ \Delta {E}^{\text{vp}}\\ \Delta {\alpha }_{s}\\ \Delta {\gamma }_{s}\\ \Delta {p}_{s}\end{array})\)

To solve this system of \(6+6+3{n}_{s}\) nonlinear equations (in 3D), we use a Newton method: we build a series of solution vectors in the following way:

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

It is therefore necessary to define the initial values \({Y}_{0}\), and to calculate the Jacobian matrix of the system: \(\frac{\text{dR}}{{\text{dY}}_{k}}\) (this is detailed in the appendix for the viscoplastic behaviors described above). It has the following form:

\(J=\left[\begin{array}{ccccc}\frac{\partial s}{\partial \text{ΔΣ}}& \frac{\partial s}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial s}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial s}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial e}{\partial \text{ΔΣ}}& \frac{\partial e}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial e}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial e}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial a}{\partial \text{ΔΣ}}& \frac{\partial a}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial a}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial a}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial g}{\partial \text{ΔΣ}}& \frac{\partial g}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial g}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial g}{\partial {\mathrm{\Delta p}}_{s}}\\ \frac{\partial p}{\partial \text{ΔΣ}}& \frac{\partial p}{\partial {\mathrm{\Delta E}}^{\text{vp}}}& \frac{\partial p}{\partial {\mathrm{\Delta \alpha }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta \gamma }}_{s}}& \frac{\partial p}{\partial {\mathrm{\Delta p}}_{s}}\end{array}\right]\)

However, this system of equations can be simplified by reducing its size in order to improve integration performance.

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 an equation can be eliminated from the global system to be solved.

Like \(\Delta p=\mid \Delta \gamma \mid\), the equation in \(\Delta p\) can be eliminated.

In the case of MONO_CINE1, \(\Delta \alpha\) is expressed directly as a function of \(\Delta \gamma\). The size of the system to be solved is therefore reduced to \(6+{n}_{s}\) equations in \(\Delta \gamma\), with \({n}_{s}\) the number of sliding systems.

In the case of MONO_CINE2, \(\Delta \alpha\) is calculated numerically as a function of \(\Delta \gamma\), by a secant method.

The equation relating to the stress components is reduced to:

\({R}_{1}^{i}({\Sigma }_{i},{\tau }_{s},\Delta {\gamma }_{s},\Delta {\alpha }_{s},{R}_{s}({p}_{r}))={\Lambda }^{-1}{\Sigma }_{i}-({\Lambda }_{-}^{-1}{\Sigma }_{{i}_{-}})-\Delta \varepsilon +\Delta {\varepsilon }^{\text{th}}+\sum _{s}{\mu }_{s}{g}_{s}({\tau }_{s},\Delta {\gamma }_{s},\Delta {\alpha }_{s},{R}_{s}({p}_{r}))=0\)

The equation relating to the work hardening parameter \(\Delta \gamma\) is written in the form:

\({R}_{2}^{s}({\Sigma }_{i},{\tau }_{s},\Delta {\gamma }_{s},\Delta {\alpha }_{s},{R}_{s}({p}_{r}))=\Delta {\gamma }_{s}-{g}_{s}({\tau }_{s},\Delta {\gamma }_{s},\Delta {\alpha }_{s},{R}_{s}({p}_{r}))=0\)

with \(i\mathrm{=}1à6\), \(s=1\) to \({n}_{s}\) and \(r=1\) to \({n}_{s}\)

In the case of MONO_VISC1:

\({g}_{s}=\Delta t{\langle \frac{\mid {\tau }_{s}-c{\alpha }_{s}\mid -{R}_{s}({p}_{r})}{k}\rangle }^{n}\frac{{\tau }_{s}-{\mathrm{c\alpha }}_{s}}{\mid {\tau }_{s}-c{\alpha }_{s}\mid }\) with \({\tau }_{s}=\Sigma :{\mu }_{s}\)

In the case of MONO_VISC2:

\({g}_{s}=\Delta t{\langle \frac{\mid {\tau }_{s}-c{\alpha }_{s}-a{\gamma }_{s}\mid -{R}_{s}({p}_{r})+\frac{d}{\mathrm{2c}}{(c{\alpha }_{s})}^{2}}{k}\rangle }^{n}\frac{{\tau }_{s}-c{\alpha }_{s}-a{\gamma }_{s}}{\mid {\tau }_{s}-c{\alpha }_{s}-a{\gamma }_{s}\mid }\)

The Jacobian matrix of the reduced system is written in the form:

\(J=\left[\begin{array}{cc}\frac{\partial {R}_{1}^{i}}{\partial {\Sigma }_{j}}& \frac{\partial {R}_{1}^{i}}{\partial \Delta {\gamma }_{s}}\\ \frac{\partial {R}_{2}^{s}}{\partial {\Sigma }_{i}}& \frac{\partial {R}_{2}^{s}}{\partial {\mathrm{\Delta \gamma }}_{r}}\end{array}\right]\) of dimension \((6+{n}_{s})\times (6+{n}_{s})\)

The different terms of the Jacobian matrix are:

\(\frac{\partial {R}_{1}^{i}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\frac{\partial {\tau }_{s}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}{({\mu }_{s})}_{j}\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\)

\(\frac{\partial {R}_{2}^{s}}{\partial {\Sigma }_{i}}=\frac{\partial (-{g}_{s}({\tau }_{s},{\mathrm{\Delta \alpha }}_{s},{\mathrm{\Delta \gamma }}_{s},{R}_{s}({p}_{r})))}{\partial {\Sigma }_{i}}=\frac{\partial (-{g}_{s})}{\partial {\tau }_{s}}\text{.}\frac{\partial {\tau }_{s}}{\partial {\Sigma }_{i}}=-{\mu }_{s}^{i}\text{.}\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\)

\(\frac{\partial {R}_{1}^{i}}{\partial \Delta {\gamma }_{s}}=\sum _{s}{\mu }_{s}^{i}\text{.}\frac{\partial {g}_{s}}{\partial \Delta {\gamma }_{s}}\)

\(\frac{\partial {R}_{2}^{s}}{\partial \Delta {\gamma }_{\mathrm{rs}}}={\delta }_{\mathrm{sr}}-\frac{\partial {g}_{s}}{\partial \Delta {\gamma }_{r}}\) with \(\frac{\partial {g}_{s}}{\partial \Delta {\gamma }_{r}}=\frac{\partial {g}_{s}}{\partial \Delta {\alpha }_{s}}\mathrm{.}\frac{\partial \Delta {\alpha }_{s}}{\partial \Delta {\gamma }_{s}}{\delta }_{\mathit{sr}}+\frac{\partial {g}_{s}}{\partial {R}_{s}}\mathrm{.}\frac{\partial {R}_{s}}{\partial \Delta {\gamma }_{r}}\)

In the expressions above, the terms \(\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\text{,}\frac{\partial {g}_{s}}{\partial \Delta {\alpha }_{s}}\text{,}\frac{\partial \Delta {\alpha }_{s}}{\partial \Delta {\gamma }_{s}}\text{,}\frac{\partial {g}_{s}}{\partial {R}_{s}}\text{,}\frac{\partial {R}_{s}}{\partial \Delta {p}_{r}}\) depend on the laws of flow and work hardening:

In the case of* MONO_VISC1 **:

\(\frac{\partial {g}_{s}}{\partial {\tau }_{s}}=\frac{n\Delta t}{{k}^{n}}{\langle \mid {\tau }_{s}-c{\alpha }_{s}\mid -{R}_{s}\left({p}_{r}\right)\rangle }^{n-1}={H}_{s}\left({\tau }_{s},\Delta {\alpha }_{s},\Delta {\gamma }_{s},{R}_{s}\left({p}_{r}\right)\right)\)

\(\frac{\partial {g}_{s}}{\partial \Delta {\alpha }_{s}}=-c\text{.}{H}_{s}({\tau }_{s},\Delta {\gamma }_{s},{R}_{s})\)

\(\frac{\partial {g}_{s}}{\partial {R}_{s}}=-{H}_{s}\left({\tau }_{s},\Delta {\gamma }_{s},{R}_{s}\right)\text{sgn}s\) with \(\text{sgn}s=\frac{{\tau }_{s}-c{\alpha }_{s}}{\mid {\tau }_{s}-c{\alpha }_{s}\mid }\)

In the case of* MONO_VISC2 **:

\(\frac{\partial {g}_{s}}{\partial {\tau }_{s}}={H}_{s}\left({\tau }_{s},\Delta {\alpha }_{s},\Delta {\gamma }_{s},{R}_{s}\left({p}_{r}\right)\right)=\frac{n}{{k}^{n}}{\langle \mid {\tau }_{s}-c{\alpha }_{s}-a{\gamma }_{s}\mid -{R}_{s}\left({p}_{r}\right)+\frac{d}{\mathrm{2c}}{\left(c{\alpha }_{s}\right)}^{2}\rangle }^{n-1}\)

\(\frac{\partial {g}_{s}}{\partial \Delta {\alpha }_{s}}={H}_{s}({\tau }_{s},\Delta {\gamma }_{s},{R}_{s})(-c\text{sgn}s+d{\alpha }_{s}\Delta {\alpha }_{s})\)

In the case of* MONO_CINE1: **

\(\frac{\partial \Delta {\alpha }_{s}}{\partial \Delta {\gamma }_{r}}=\frac{1-d{\alpha }_{s}^{-}\text{sgn}s}{{\left(1+d\mid \Delta {\gamma }_{s}\mid \right)}^{2}}\)

In the case of* MONO_CINE2: **

\(\frac{\partial {\alpha }_{s}}{\partial \Delta {\gamma }_{r}}\) ; This derivation requires the resolution of a non-linear equation by Newton’s method

In the case of* MONO_ISOT1: **

\(\frac{\partial {R}_{s}}{\partial \Delta {\gamma }_{r}}=\frac{\partial \left[Q\sum _{k}{h}_{\text{sk}}(1-{e}^{-{\text{bp}}_{k}})\right]}{\partial \Delta {\gamma }_{r}}=\frac{\partial }{\partial {p}_{r}}\left[Q\sum _{k}{h}_{\text{sk}}(1-{e}^{-{\text{bp}}_{k}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma }_{r}}={\text{bQh}}_{\text{sr}}{e}^{-{\text{bp}}_{r}}\text{.}\text{sgn}r\)

with \(\text{sgn}r=\frac{\partial {p}_{r}}{\partial \Delta {\gamma }_{r}}=\frac{\partial \mid \Delta {\gamma }_{r}\mid }{\partial \Delta {\gamma }_{r}}=\frac{\Delta {\gamma }_{r}}{\mid \Delta {\gamma }_{r}\mid }\)

In the case of* MONO_ISOT2: **

\(\begin{array}{}\frac{\partial {R}_{s}}{\partial \Delta {\gamma }_{r}}=\frac{\partial \left[{Q}_{1}\sum _{k}{h}_{\text{sk}}(1-{e}^{-{b}_{1}{p}_{k}})+{Q}_{2}(1-{e}^{-{b}_{2}{p}_{r}})\right]}{\partial \Delta {\gamma }_{r}}=\\ \frac{\partial }{\partial {p}_{r}}\left[{Q}_{1}\sum _{k}{h}_{\text{sk}}(1-{e}^{-{b}_{1}{p}_{k}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma }_{r}}+\frac{\partial }{\partial {p}_{r}}\left[{Q}_{2}(1-{e}^{-{b}_{2}{p}_{r}})\right]\frac{\partial {p}_{r}}{\partial \Delta {\gamma }_{r}}=\\ ({b}_{1}{Q}_{1}{h}_{\text{sr}}{e}^{-{b}_{1}{p}_{r}}+{b}_{2}{Q}_{2}{e}^{-{b}_{2}{p}_{r}}{\delta }_{\text{sr}})\text{sgn}r\end{array}\)

3.2.2. Implicit resolution — Model MONO_DD_KR integration#

The system to be solved is reduced to the equations relating to the components of the stresses.

\({R}_{1}^{i}({\Sigma }_{i},{\tau }_{s},{\alpha }_{s})={({\Lambda }^{-1}\Sigma )}_{i}-{({\Lambda }_{\text{-}}^{\text{-1}}{\Sigma }_{\text{-}})}_{i}-\Delta {\varepsilon }_{i}+\Delta {\varepsilon }_{i}^{\text{th}}+\sum _{s}{\mu }_{i}^{s}{g}_{s}({\tau }_{s},{\alpha }_{s})=0\)

and those relating to the work hardening parameter:

\({R}_{2}^{s}({\sum }_{i},{\tau }_{s},{\alpha }_{s})=\Delta {\alpha }_{s}-\mid {g}_{s}({\tau }_{s},{\alpha }_{s})\mid \text{.}h({\alpha }_{s})=0\)

with \(i\) = 1 to 6, \(s\) = 1 to \(\mathrm{ns}\) and \(r\) =1 to \(\mathrm{ns}\)

This makes it possible to obtain a system of equations formulated in an entirely implicit way, where \({R}_{1}^{i}\) and \({R}_{2}^{s}\) only depend on \({\Sigma }_{i},{\alpha }_{s}\) (because \({\tau }_{s}\) is calculated explicitly from \({\Sigma }_{i}\)).

\({g}_{s}({\tau }_{s},{\alpha }_{s})={\dot{\gamma }}_{0}\text{exp}(\frac{-\Delta G({\tau }_{\text{eff}}^{s})}{{k}_{B}T})\text{.}\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid }\text{.}\Delta t\)

\(\Delta G({\tau }_{\text{eff}}^{s})=\Delta {G}_{0}{(1-{(\frac{\langle {\tau }_{\text{eff}}^{s}\rangle }{{\tau }_{R}})}^{p})}^{q}\)

\(h({\alpha }_{s})=(\frac{b}{d}+\frac{\sqrt{\sum _{u\ne s}{\alpha }^{u}}}{K}-\frac{{g}_{c}{\alpha }^{s}}{b})\)

\({\tau }_{\text{eff}}^{s}=\langle \mid {\tau }^{s}\mid -{\tau }_{0}-{\tau }_{\mu }^{s}\rangle\) \(\langle \rangle\): positive part

\({\tau }_{\mu }^{s}=\frac{{(\mu )}^{2}\sum _{u}{a}^{\text{su}}{\alpha }^{u}}{\mid {\tau }^{s}\mid -{\tau }_{0}}\) with \({\alpha }^{u}\mathrm{=}{\rho }^{u}{b}^{2}\)

Note: for each sliding system \(s\) where conditions \({\tau }_{\text{eff}}^{s}>0\mid {\tau }^{s}\mid >{\tau }_{0}\) are not met, we immediately get: \({g}_{s}({\tau }_{s},{\alpha }_{s})=\mathrm{0,}\Delta {\alpha }_{s}=0\)

The Jacobian matrix of the system is written in the form:

\(J=\left[\begin{array}{cc}\frac{\partial {R}_{1}^{i}}{\partial {\Sigma }_{j}}& \frac{\partial {R}_{1}^{i}}{\partial \Delta {\alpha }_{s}}\\ \frac{\partial {R}_{2}^{s}}{\partial {\Sigma }_{i}}& \frac{\partial {R}_{2}^{s}}{\partial \Delta {\alpha }_{r}}\end{array}\right]\) of dimension \((6+{n}_{s})\times (6+{n}_{s})\)

\(\frac{\partial {R}_{1}^{i}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\frac{\partial {\tau }_{s}}{\partial {\Sigma }_{j}}={({\Lambda }^{-1})}_{\mathrm{ij}}+\sum _{s}{({\mu }_{s})}_{i}{({\mu }_{s})}_{j}\frac{\partial {g}_{s}}{\partial {\tau }_{s}}\)

=> You have to calculate \(\frac{\partial {g}^{s}}{\partial {\tau }_{s}}=\text{sgn}({\tau }_{s})\text{.}\frac{\partial \mid {g}^{s}\mid }{\partial {\tau }_{s}}\)

\(\frac{\partial {R}_{2}^{s}}{\partial {\Sigma }_{i}}=\frac{\partial (-\mid {g}_{s}\mid \text{.}{h}_{s})}{\partial {\tau }_{s}}\text{.}\frac{\partial {\tau }_{s}}{\partial {\Sigma }_{i}}=-{\mu }_{s}^{i}\text{.}(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{s}}\text{.}{h}_{s}+\mid {g}_{s}\mid \text{.}\frac{\partial {h}_{s}}{\partial {\tau }_{s}})\)

=> Since the term \(\frac{\partial {h}_{s}}{\partial {\tau }_{s}}\) is zero, it remains to calculate \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{s}}\).

\(\frac{\partial {R}_{1}^{i}}{\partial \Delta {\alpha }_{s}}=\sum _{s}{\mu }_{s}^{i}\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid }\text{.}\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha }_{s}}\)
\(\frac{\partial {R}_{2}^{s}}{\partial \Delta {\alpha }_{r}}={\delta }_{\text{sr}}-\frac{\partial \mid {g}_{s}\mid }{\partial \Delta {\gamma }_{r}}={\delta }_{\text{sr}}-(\frac{\partial \mid {g}_{s}\mid }{\partial \Delta {\alpha }_{r}}{h}_{s}+\mid {g}_{s}\mid \frac{\partial {h}_{s}}{\partial \Delta {\alpha }_{r}})\)

=> You must also calculate \(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha }_{r}}\) and \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{s}}\);

\({h}_{s}=\frac{\Delta {\alpha }_{s}}{\mid \Delta {\gamma }_{s}\mid }\) and \(\frac{\partial {\tau }_{\mu }^{s}}{\partial {\tau }_{s}}=\frac{-{\tau }_{\mu }^{s}}{\mid {\tau }^{s}\mid -{\tau }_{0}}\text{.}\text{sgn}s\)

Calculation of \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{r}}\):

\(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{s}}=-\frac{\Delta t\text{.}{\dot{\gamma }}_{0}}{{k}_{B}T}\text{exp}(-\frac{\Delta G({\tau }_{\text{eff}}^{s})}{{k}_{B}T})\text{.}\frac{\Delta G({\tau }_{\text{eff}}^{s})}{\Delta {\tau }_{s}}=\frac{-\mid {g}_{s}({\tau }_{s},{\alpha }_{s})\mid }{{k}_{B}T}\text{.}\frac{\partial \Delta G({\tau }_{\text{eff}}^{s})}{\partial {\tau }_{s}}\)

where

\(\frac{\partial \Delta G({\tau }_{\text{eff}}^{s})}{\partial {\tau }_{s}}=-\Delta {G}_{0}{(1-{(\frac{{\tau }_{\text{eff}}^{s}}{{\tau }_{R}})}^{p})}^{q-1}\text{.}\frac{q\text{.}p}{{\tau }_{R}}\text{.}{(\frac{{\tau }_{\text{eff}}^{s}}{{\tau }_{R}})}^{p-1}\text{.}\frac{\partial {\tau }_{\text{eff}}^{s}}{\partial {\tau }_{s}}\)

and

\(\frac{\partial {\tau }_{\text{eff}}^{s}}{\partial {\tau }_{s}}=\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid }-\frac{\partial {\tau }_{\mu }^{s}}{\partial {\tau }_{s}}=\text{sgn}({\tau }_{s})+{\mu }^{2}(\sum {a}^{\text{su}}{\alpha }^{u})\text{.}\frac{\text{sgn}({\tau }_{s})}{{(\mid {\tau }_{s}\mid -{\tau }_{0})}^{2}}=\text{sgn}({\tau }_{s})+{\tau }_{\mu }^{s}\text{.}\frac{\text{sgn}({\tau }_{s})}{\mid {\tau }_{s}\mid -{\tau }_{0}}\)

Calculation of \(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha }_{r}}\):

\(\frac{\partial \mid {g}_{s}\mid }{\partial {\alpha }_{r}}=-\frac{\Delta t\text{.}{\dot{\gamma }}_{0}}{{k}_{B}T}\text{exp}(-\frac{\Delta G({\gamma }_{\text{eff}}^{s})}{{k}_{B}T})\text{.}\frac{\Delta G({\tau }_{\text{eff}}^{s})}{\Delta {\alpha }_{r}}=\frac{-\mid {g}_{s}({\tau }_{s},{\alpha }_{s})\mid }{{k}_{B}T}\text{.}\frac{\partial \Delta G({\tau }_{\text{eff}}^{s})}{\partial {\alpha }_{r}}\)

where

\(\frac{\partial \Delta G({\tau }_{\text{eff}}^{s})}{\partial {\alpha }_{r}}=-\frac{\text{qp}}{{\tau }_{R}}\Delta {G}_{0}{(1-{(\frac{{\tau }_{\text{eff}}}{{\tau }_{R}})}^{p})}^{q-1}\text{.}{(\frac{{\tau }_{\text{eff}}}{{\tau }_{R}})}^{p-1}\text{.}\frac{\partial {\tau }_{\text{eff}}^{s}}{\partial {\alpha }_{r}}\)

and

\(\frac{\partial {\tau }_{\text{eff}}^{s}}{\partial {\alpha }_{r}}=-\frac{\partial {\tau }_{\mu }^{s}}{\partial {\alpha }_{r}}=-\frac{{\mu }^{2}}{\mid {\tau }_{s}\mid -{\tau }_{0}}\frac{\partial (\sum _{u}{a}^{\text{su}}{\alpha }^{u})}{\partial {\alpha }_{r}}=-\frac{{\mu }^{2}}{\mid {\tau }_{s}\mid -{\tau }_{0}}{a}^{\text{sr}}\)

Calculation of \(\frac{\partial \mid {g}_{s}\mid }{\partial {\tau }_{s}}\):

\(\frac{\partial {h}^{s}}{\partial {\alpha }_{r}}=\frac{1}{K}\frac{\partial }{\partial {\alpha }_{r}}\sqrt{\sum _{u\ne s}{\alpha }^{u}}-\frac{{g}_{c}}{b}{\delta }_{\text{sr}}=\frac{1}{\mathrm{2K}}\frac{1-{\delta }_{\text{rs}}}{\sqrt{\sum _{u\ne s}{\alpha }^{u}}}-\frac{{g}_{c}}{b}{\delta }_{\text{sr}}\)

3.2.3. Implicit resolution — models MONO_DD_CFC (IRRA)/MONO_DD_FAT#

Expression of the residual Knowing the solution at time \({t}_{i-1}\), noted \({\sigma }^{\text{-}}=\sigma ({t}_{i-1}),{\rho }_{s}^{\text{-}}=\rho {}_{s}\text{}({t}_{i-1}),{\gamma }_{s}^{\text{-}}={\gamma }_{s}({t}_{i-1}),{p}_{s}^{\text{-}}={p}_{s}({t}_{i-1})\),

It is a question of finding, at moment \({t}_{i}\), the quantities \(\sigma =\sigma ({t}_{i}),{\rho }_{s}=\rho {}_{s}\text{}({t}_{i}),{\gamma }_{s}={\gamma }_{s}({t}_{i}),{p}_{s}={p}_{s}({t}_{i})\) verifying:

\({({\Lambda }^{-1}\sigma )}_{i}-{({({\Lambda }^{\text{-}})}^{\text{-1}}{\sigma }^{\text{-}})}_{i}=\Delta {\varepsilon }_{i}-\Delta {\varepsilon }_{i}^{\text{th}}-\sum _{s}({({\mu }_{s})}_{i}\Delta {\gamma }_{s}(\sigma ,\omega ))=0\)

The cission resolved by sliding system is: \({\tau }_{s}=\sigma :{\mu }_{s}=({\sigma }^{\text{-}}+\Delta \sigma ):{\mu }_{s}\)

The flow law after implicit discretization becomes:

\({\mathrm{\Delta \gamma }}_{s}={\dot{\gamma }}_{0}\mathrm{\Delta t}({(\frac{\mid {\tau }_{s}\mid }{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}})}^{n}-1)\text{.}\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid }=\Delta {p}_{s}\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid }\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma }_{0}}\Delta t({(\frac{\mid {\tau }_{s}\mid }{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) if \(∣{\tau }_{s}∣\ge {\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}\); otherwise \(\Delta {p}_{s}=0\)

The law of work hardening is: \({\tau }_{s}^{\mathrm{forest}}(\omega )={\tau }_{s}^{\mathrm{forest}}(\langle {\omega }^{\text{-}}+\Delta \omega \rangle )=\mu C(\omega )\sqrt{{\sum }_{j=\mathrm{1,12}}\langle {a}_{\mathrm{sj}}{\omega }_{j}\rangle }\)

In the presence of irradiation, the law of work hardening becomes:

\({\tau }_{s}^{\mathit{forest}}=\mu \sqrt{\sum _{j=\mathrm{1,12}}C{(\omega )}^{2}{a}_{\mathit{sj}}{\omega }_{j}+{\alpha }^{\mathit{loops}}{\varphi }^{\mathit{loops}}{\omega }_{s}^{\mathit{loops}}+{\alpha }^{\mathit{voids}}{\varphi }_{s}^{\mathit{voids}}{\omega }^{\mathit{voids}}}\)

with \(C(\omega )=0.2+0.8\frac{\mathrm{ln}\left(\alpha \sqrt{\sum _{i=\mathrm{1,12}}\langle {\omega }_{i}\rangle }\right)}{\mathrm{ln}\left(\alpha b\sqrt{{\rho }_{\mathit{ref}}}\right)}\) or, in the case DD_FAT, \(C(\omega )=1.\)

The evolution of the dislocation density is given (always in implicit form) by:

\(\Delta {\omega }_{s}=\Delta {p}_{s}{h}_{s}(\omega )=\Delta {p}_{s}{h}_{s}({\omega }^{\text{-}}+\Delta \omega )\) with:

\({h}_{s}(\omega )=(A\frac{\sum _{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}\langle {\omega }_{j}\rangle }{\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}\langle {\omega }_{j}\rangle }}+BC(\omega )\sum _{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}\langle {\omega }_{j}\rangle }-\frac{y}{b}{\omega }_{s})\)

or, in case DD_FAT, \({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum _{u\ne s}\langle {\omega }_{u}\rangle }}{K}-{g}_{\mathrm{c0}}\frac{\langle {\omega }_{s}\rangle }{b})\)

The evolution of dislocation densities associated with irradiation is of the following form:

\({\omega }_{s}^{\mathit{loops}}(t+\Delta t)={\omega }^{\mathit{sat}}+\left({\omega }_{s}^{\mathit{loops}}(t)-{\omega }^{\mathit{sat}}\right)\mathrm{exp}\left(-\xi \sum _{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

\({\varphi }_{s}^{\mathit{voids}}(t+\Delta t)={\varphi }^{\mathit{sat}}+\left({\varphi }_{s}^{\mathit{voids}}(t)-{\varphi }^{\mathit{sat}}\right)\mathrm{exp}\left(-\zeta \sum _{j\in \mathit{copla}(s)}∣{\Delta \gamma }_{j}∣\right)\)

Using the Kelvin vector notation [chapter 5 of 15) for symmetric tensors (\(\sigma\), \({\mu }_{s}\)) for symmetric tensors (,), and taking into account the possible variations of the elasticity coefficients with temperature, we end up with the system to be solved, comprising 18 equations; \(i=\mathrm{1,}6\), \(s=\mathrm{1,12}\), comprising 18 equations;,, the 18 unknowns being the 6 components of \(\sigma\) and the 12 values \(\omega\).

System to be solved in \(\sigma\) (6 components) and \(\Delta \omega\) (12 components):

\({R}_{i}^{(1)}(\sigma ,\omega )={({\Lambda }^{-1}\sigma )}_{i}-{({({\Lambda }^{\text{-}})}^{\text{-1}}{\sigma }^{\text{-}})}_{i}-\Delta {\varepsilon }_{i}+\Delta {\varepsilon }_{i}^{\text{th}}+\sum _{s}({({\mu }_{s})}_{i}\Delta {p}_{s}(\sigma ,\omega )\text{.}\frac{{\tau }_{s}}{\mid {\tau }_{s}\mid })=0\)

\({R}_{s}^{(2)}(\sigma ,\omega )=\Delta {\omega }_{s}-\Delta {p}_{s}{h}_{s}(\omega )=0\)

with: \({\tau }_{s}=\sigma :{\mu }_{s}=({\sigma }^{\text{-}}+\Delta \sigma ):{\mu }_{s}\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma }_{0}}\Delta t({(\frac{\mid {\tau }_{s}\mid }{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) if \(∣{\tau }_{s}∣\ge {\tau }_{0}+{\tau }_{s}^{f}\); otherwise \(\Delta {p}_{s}=0\)

and if \(\Delta {p}_{s}>0\):

\({\tau }_{s}^{\mathrm{forest}}(\omega )={\tau }_{s}^{\mathrm{forest}}({\omega }^{\text{-}}+\Delta \omega )=\mu C(\omega )\sqrt{{\sum }_{j=\mathrm{1,12}}{a}_{\mathrm{sj}}\langle {\omega }_{j}\rangle }\)

or \({\tau }_{s}^{\mathit{forest}}=\mu \sqrt{\sum _{j=\mathrm{1,12}}C{(\omega )}^{2}{a}_{\mathit{sj}}{\omega }_{j}+{\alpha }^{\mathit{loops}}{\varphi }^{\mathit{loops}}{\omega }_{s}^{\mathit{loops}}+{\alpha }^{\mathit{voids}}{\varphi }_{s}^{\mathit{voids}}{\omega }^{\mathit{voids}}}\)

\(C(\omega )=0.2+0.8\frac{\mathrm{ln}(\alpha \sqrt{\sum _{i=\mathrm{1,12}}\langle {\omega }_{i}\rangle })}{\mathrm{ln}(\alpha b\sqrt{{\rho }_{\mathrm{ref}}})}\)

In case DD_FAT, this system becomes:

\({R}_{s}^{(2)}(\sigma ,\omega )=\Delta {\omega }_{s}-\Delta {p}_{s}{h}_{s}(\omega )=0\)

with: \({\tau }_{s}=\sigma :{\mu }_{s}=({\sigma }^{\text{-}}+\Delta \sigma ):{\mu }_{s}\)

\(\Delta {p}_{s}(\sigma ,\omega )=\dot{{\gamma }_{0}}\Delta t({(\frac{\mid {\tau }_{s}\mid }{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )})}^{n}-1)\) if \(∣{\tau }_{s}∣\ge {\tau }_{0}+{\tau }_{s}^{f}\); otherwise \(\Delta {p}_{s}=0\)

and if \(\Delta {p}_{s}>0\):

\({\tau }_{s}^{\mathrm{forest}}(\omega )={\tau }_{s}^{\mathrm{forest}}({\omega }^{\text{-}}+\Delta \omega )=\mu \sqrt{{\sum }_{j=\mathrm{1,12}}{a}_{\mathrm{sj}}\langle {\omega }_{j}\rangle }\)

\({h}_{s}(\omega )=(\frac{b}{d}+\frac{\sqrt{\sum _{u\ne s}\langle {\omega }_{u}\rangle }}{K}-{g}_{\mathrm{c0}}\frac{\langle {\omega }_{s}\rangle }{b})\)

The Jacobian matrix of the system can then be calculated for integration by Newton’s method (R5.03.14).

Jacobian matrix

This paragraph only concerns laws DD_CFC and DD_CC: for law DD_FAT, the Jacobian matrix has not been programmed and the matrix calculated by disturbance is used. The Jacobian matrix of the system is written in the form:

\(J=\left[\begin{array}{cc}{(\frac{\partial {R}_{i}^{(1)}}{\partial {\sigma }_{j}})}_{i=\mathrm{1,6};j=\mathrm{1,6}}& {(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega }_{s}})}_{i=\mathrm{1,6};s=\mathrm{1,12}}\\ {(\frac{\partial {R}_{s}^{(2)}}{\partial {\sigma }_{j}})}_{s=\mathrm{1,12};j=\mathrm{1,6}}& {(\frac{\partial {R}_{s}^{(2)}}{\partial \Delta {\omega }_{r}})}_{s=\mathrm{1,12};r=\mathrm{1,12}}\end{array}\right]\) of dimension \(18\times 18\)

\(\begin{array}{c}\frac{\partial {R}_{i}^{(1)}}{\partial {\sigma }_{j}}={\left({\Lambda }^{-1}\right)}_{\mathit{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{\partial \left(\Delta {p}_{s}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\right)}{\partial {\sigma }_{j}}={\left({\Lambda }^{-1}\right)}_{\mathit{ij}}+\sum _{s}{({\mu }_{s})}_{i}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}\frac{\partial {\tau }_{s}}{\partial {\sigma }_{j}}=\\ {\left({\Lambda }^{-1}\right)}_{\mathit{ij}}+\sum _{s}{({\mu }_{s})}_{i}{({\mu }_{s})}_{j}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}\end{array}\)

with \(\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}=n\frac{\left(\Delta {p}_{s}+\dot{{\gamma }_{0}}\Delta t\right)}{{\tau }_{s}}\)

indeed: \(\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}=\frac{n\dot{{\gamma }_{0}}\Delta t}{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )}{(\frac{∣{\tau }_{s}∣}{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )})}^{n-1}\cdot \frac{{\tau }_{s}}{∣{\tau }_{s}∣}=n\frac{(\Delta {p}_{s}+\dot{{\gamma }_{0}}\Delta t)}{∣{\tau }_{s}∣}\cdot \frac{{\tau }_{s}}{∣{\tau }_{s}∣}\)

because \({\tau }_{s}^{\mathrm{forest}}(\omega )\) doesn’t depend on \({\tau }_{s}\) and \(\Delta {p}_{s}+\dot{{\gamma }_{0}}\Delta t=\dot{{\gamma }_{0}}\Delta t{(\frac{\mid {\tau }_{s}\mid }{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}(\omega )})}^{n}\)

\(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega }_{s}}=\sum _{u=\mathrm{1,12}}{({\mu }_{u})}_{i}\frac{{\tau }_{u}}{∣{\tau }_{u}∣}\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega }_{s}}\)

with \(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega }_{s}}=-n\dot{{\gamma }_{0}}\Delta t∣{\tau }_{u}∣{(\frac{∣{\tau }_{u}∣}{{\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega )})}^{n-1}\frac{\frac{\partial {\tau }_{u}^{\mathrm{forest}}}{\partial \Delta {\omega }_{s}}}{{({\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega ))}^{2}}\) or

\(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega }_{s}}=-n\dot{{\gamma }_{0}}\Delta t{(\frac{∣{\tau }_{u}∣}{{\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega )})}^{n}\frac{\frac{\partial {\tau }_{u}^{\mathrm{forest}}}{\partial \Delta {\omega }_{s}}}{({\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega ))}\)

\(\frac{\partial \Delta {p}_{u}}{\partial \Delta {\omega }_{s}}=-n\frac{(\Delta {p}_{u}+\dot{{\gamma }_{0}}\Delta t)}{({\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega ))}\frac{\partial {\tau }_{u}^{\mathrm{forest}}}{\partial \Delta {\omega }_{s}}\)

It remains to be expressed: \(\frac{\partial {\tau }_{u}^{\mathit{forest}}}{\partial \Delta {\omega }_{s}}\). Without irradiation, we obtain:

\(\frac{\partial {\tau }_{u}^{\mathrm{forest}}}{\partial \Delta {\omega }_{s}}=\mu \frac{\partial C(\omega )}{\partial \Delta {\omega }_{s}}\cdot \sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathrm{uj}}{\omega }_{j}}+\mu C(\omega )\frac{{a}_{\mathrm{us}}}{2\sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathrm{uj}}{\omega }_{j}}}\frac{\langle {\omega }_{s}\rangle }{{\omega }_{s}}\)

and with irradiation: \(\frac{\partial ({\tau }_{u}^{\mathit{forest}})}{\partial (\Delta {\omega }_{s})}=\frac{{\mu }^{2}C(\omega )}{2{\tau }_{u}^{\mathit{forest}}}\left(2\frac{\partial C(\omega )}{\Delta {\omega }_{s}}\sum _{j=\mathrm{1,12}}{a}_{\mathit{uj}}{\omega }_{j}+C(\omega ){a}_{\mathit{us}}\right)\)

with \(\frac{\partial C(\omega )}{\partial \Delta {\omega }_{s}}=\frac{0.8}{2\mathrm{ln}(\alpha b\sqrt{{\rho }_{\mathrm{ref}}})}\frac{1}{\sum _{k=\mathrm{1,12}}{\omega }_{k}}\frac{\langle {\omega }_{s}\rangle }{{\omega }_{s}}\)

Which provides the derivative: \(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega }_{s}}=-\sum _{u=\mathrm{1,12}}{({\mu }_{u})}_{i}\frac{{\tau }_{u}}{∣{\tau }_{u}∣}n\frac{(\Delta {p}_{u}+\dot{{\gamma }_{0}}\Delta t)}{({\tau }_{f}+{\tau }_{u}^{\mathrm{forest}}(\omega ))}\frac{\partial {\tau }_{u}^{\mathrm{forest}}}{\partial \Delta {\omega }_{s}}\)

\(\frac{\partial {R}_{s}^{2}}{\partial {\sigma }_{i}}=-{({\mu }_{s})}_{i}\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}{h}_{s}(\omega )\)

with \(\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}=n\frac{\left(\Delta {p}_{s}+\dot{{\gamma }_{0}}\Delta t\right)}{{\tau }_{s}}\)

\(\frac{\partial {R}_{s}^{(2)}}{\partial {\omega }_{r}}={\delta }_{\mathrm{sr}}-\Delta {p}_{s}\frac{\partial {h}_{s}(\omega )}{\partial {\omega }_{r}}-{h}_{s}(\omega )\frac{\partial \Delta {p}_{s}}{\partial {\omega }_{r}}\)
We’ve already calculated \(\frac{\partial \Delta {p}_{s}}{\partial \Delta {\omega }_{r}}\). It remains to be drifted: \(\frac{\partial {h}_{s}(\omega )}{\partial \Delta {\omega }_{r}}\):
- \(\frac{\partial }{\partial \Delta {\omega }_{r}}(A\frac{\sum _{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\omega }_{j}}{\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}{\omega }_{j}}}+BC(\omega )\sum _{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\omega }_{j}}-\frac{y}{b}{\omega }_{s})=A\frac{\partial {T}_{1}(\omega )}{\partial \Delta {\omega }_{r}}+B\frac{\partial {T}_{2}(\omega )}{\partial \Delta {\omega }_{r}}-\frac{y}{b}{\delta }_{\mathrm{sr}}\frac{\langle {\omega }_{r}\rangle }{{\omega }_{r}}\)

\(\frac{\partial }{\partial \Delta {\omega }_{r}}({T}_{1}(\omega ))=\frac{\sqrt{{a}_{\mathrm{sr}}}}{\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}{\omega }_{j}}}\cdot {I}_{s}^{\mathrm{forest}}(r)\frac{\langle {\omega }_{r}\rangle }{{\omega }_{r}}-\frac{{a}_{\mathrm{sr}}\sum _{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\omega }_{j}}{2\sqrt{{a}_{\mathrm{sr}}{\omega }_{r}}{(\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}{\omega }_{j}})}^{2}}\frac{\langle {\omega }_{r}\rangle }{{\omega }_{r}}\)

\(\frac{\partial }{\partial \Delta {\omega }_{r}}{T}_{2}(\omega )=\frac{\partial C}{\partial \Delta {\omega }_{r}}\sum _{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\omega }_{j}}+C(\omega )\frac{{a}_{\mathrm{sr}}}{2\sqrt{{a}_{\mathrm{sr}}{\omega }_{r}}}{I}_{s}^{\mathrm{copla}}(r)\frac{\langle {\omega }_{r}\rangle }{{\omega }_{r}}\)

\({I}_{s}^{\mathrm{copla}}(r)=1\) if \(r\in \mathrm{copla}(s)\), \(=0\) otherwise and \({I}_{s}^{\mathrm{forest}}(r)=1\) if \(r\in \mathrm{forest}(s)\), \(0\) otherwise

3.2.4. Implicit resolution — models MONO_DD_CC/MONO_DD_CC_IRRA#

We choose to dimension the unknowns by asking \(\omega =\rho {b}^{2}\). The proposed algorithm (cf. 17) consists in eliminating the dependence of \({R}^{s}\) to \({\tau }_{\mathrm{eff}}^{s}\) by approaching \(\Delta G\) by:

\(\Delta {G}_{\mathit{app}}=\mathit{min}\left(\Delta {G}_{0};kT\mathrm{ln}\left[\frac{{\omega }_{m}^{s}H}{\sqrt{{\omega }_{\mathit{tot}}^{s}}{\dot{\epsilon }}_{\mathit{eq}}}\right]\right)\) with \({\omega }_{\mathit{tot}}^{s}(\omega )={\omega }_{F}^{s}(\omega )=\sum _{j\ne s}{\omega }^{j}+{\omega }_{\mathit{irr}}^{s}\)

In this expression, \(\dot{{\varepsilon }_{\mathrm{eq}}}\) is a fixed parameter.

\({R}^{s}(\omega )=\frac{\mu b}{2{\tau }_{0}{\left(1-\frac{\Delta {G}_{\mathit{app}}}{\Delta {G}_{0}}\right)}^{2}}\) by checking that \(\Delta {G}_{\mathit{app}}\le \Delta {G}_{0}\) (otherwise, \({R}^{s}\) is chosen large)
\(\frac{1}{{\lambda }^{s}(\omega )+D}=\mathit{min}\left(\frac{\sqrt{{\omega }_{\mathit{tot}}^{s}(\omega )}}{b};(D+2{R}^{s}(\omega ))\frac{{\omega }_{\mathit{tot}}^{s}(\omega )}{{b}^{2}}\right)\)
\({l}^{s}(\omega )=\mathit{max}\left({\lambda }^{s}(\omega )-2{\alpha }_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega );{l}_{c}\right)\)
\({\alpha }_{\mathit{AT}}^{s}(\omega )=\sqrt{\sum _{j\ne s}{h}^{\mathit{sj}}\frac{{\omega }^{j}}{{\omega }_{\mathit{tot}}^{s}}+{a}_{\mathit{irr}}\frac{{\omega }_{\mathit{irr}}^{s}}{{\omega }_{\mathit{tot}}^{s}}}\) with \({h}_{\mathrm{sj}}\) the current term in the interaction matrix
\({\tau }_{\text{LT}}^{s}(\omega )=\text{max}\left(0;{\alpha }_{\mathit{AT}}^{s}(\omega )\mu b\left(\frac{1}{{\lambda }^{s}(\omega )}-\frac{1}{2{\alpha }_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega )+{l}_{c}(T)}\right)\right)\)
\({\tau }_{\mathit{LR}}^{s}(\omega )=\mu \sqrt{{h}^{\mathit{ss}}{\omega }^{s}}\)
\({\tau }_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})=∣{\tau }^{s}∣-{\tau }_{c}(\omega )\) and \({\tau }_{c}^{s}={\tau }_{F}+\text{}\sqrt{({\tau }_{\mathit{LR}}^{s}{(\omega )}^{2}+{\tau }_{\text{LT}}^{s}{(\omega )}^{2})}\)
\({\Delta \gamma }_{\mathit{nuc}}^{S}(\omega ,{\tau }^{s})=\frac{\Delta t{\omega }_{m}^{s}H}{b}\cdot {l}^{s}(\omega )\mathrm{exp}\left(\frac{-\Delta G(\omega ,{\tau }^{s})}{{k}_{B}T}\right)\mathit{sgn}({\tau }_{s})\) with \(\Delta G(\omega ,{\tau }^{s})=\Delta {G}_{0}\left(1-\sqrt{\frac{\text{<}{\tau }_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})\text{>}}{{\tau }_{0}}}\right)\) if \({\tau }_{\mathit{eff}}^{s}\le {\tau }_{0}\), otherwise \(\Delta G=\Delta {G}_{0}\)
\({\Delta \gamma }_{\mathit{prob}}^{s}(\omega ,{\tau }^{s})=\Delta t\dot{{\gamma }_{0}}{\left(\frac{∣{\tau }^{s}∣}{{\tau }_{c}^{s}(\omega )}\right)}^{n}\mathit{sgn}({\tau }^{s})\)
\(\frac{1}{{\Delta \gamma }^{s}(\omega ,{\tau }^{s})}=\frac{1}{{\Delta \gamma }_{\mathit{nuc}}^{s}(\omega ,{\tau }^{s})}+\frac{1}{{\Delta \gamma }_{\mathit{prob}}^{s}(\omega ,{\tau }^{s})}\)
\(\Delta {\omega }^{s}(\omega ,{\tau }^{s})=∣\Delta {\gamma }^{s}(\omega ,{\tau }^{s})∣\left(\frac{b}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})\sqrt{{h}^{\mathit{ss}}{\omega }^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha }_{\mathit{AT}}^{s}(\omega ){c}_{\mathit{eff}}^{s}(\omega ,{\tau }^{s}){\lambda }_{s}(\omega ){\omega }_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-{y}^{s}(\omega ,{\tau }^{s})\frac{{\omega }^{s}}{b}\right)\) with

\({c}_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})=\left(1-\frac{\text{<}{\tau }_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})\text{>}}{{\tau }_{0}}\right)\) \(\frac{1}{{y}^{s}}(\omega ,{\tau }^{s})=\frac{1}{{y}_{\mathit{AT}}^{s}}+\frac{2\pi {\tau }_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})}{\mu b}\)

In the case DD_CC_IRRA, \(\dot{{\omega }_{\mathit{irr}}^{s}}=-\xi {\omega }_{\mathit{irr}}^{s}∣\dot{{\gamma }^{s}}(\omega )∣\), that is: \({\omega }_{\mathit{irr}}^{s}(t+\Delta t)={\omega }_{\mathit{irr}}^{s}(t)\text{exp}(-\xi ∣\Delta {\gamma }_{s}∣)\)

So you have to solve a system of 12 equations in \(\Delta {\omega }^{s}\)

Note: If you want to completely solve the above equations without approximation \(\Delta {G}_{\mathit{app}}\), you must then add the 12 unknowns \({\tau }_{\mathit{eff}}\). In fact, the dependence on \({\tau }_{\mathit{eff}}\) is implicit in the expression for \({R}_{s}\) if we replace \(\Delta {G}_{\mathit{app}}\) by \(\Delta G\).

The Jacobian matrix of the system can then be calculated for integration by Newton’s method.

Jacobian matrix

With the approximation chosen, the system of equations to be solved is:

\({R}_{i}^{(1)}\left(\Delta {\epsilon }^{\mathrm{el}},\omega \right)=\Delta {\epsilon }_{i}^{\mathit{el}}-\Delta {\epsilon }_{i}+\Delta {\epsilon }_{i}^{\text{th}}+\sum _{s}\left({\mu }_{i}^{s}\Delta {\gamma }^{s}\left({\tau }^{s},\omega \right)\right)=0\) with \({\tau }^{s}=\sigma \mathrm{:}{\mu }^{s}=\Lambda {\epsilon }^{\mathit{el}}\mathrm{:}{\mu }^{s}\)

\({R}_{s}^{(2)}\left({\tau }^{s},\omega \right)=\Delta {\omega }_{s}-\Delta {p}_{s}\left(\omega ,{\tau }^{s}\right){H}_{s}\left(\omega ,{\tau }^{s}\right)=0\) with \(\Delta {p}_{s}(\omega ,{\tau }^{s})=∣(\Delta {\gamma }_{s}(\omega ,{\tau }^{s}))∣\) and

\({H}^{s}(\omega ,{\tau }^{s})=\left(\frac{b}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}^{s}(\omega ,{\tau }^{s})\sqrt{{h}^{\mathit{ss}}{\omega }^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha }_{\mathit{AT}}^{s}(\omega ){c}_{\mathit{eff}}^{s}(\omega ,{\tau }^{s}){\lambda }_{s}(\omega ){\omega }_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-{y}^{s}(\omega ,{\tau }^{s})\frac{{\omega }^{s}}{b}\right)\)

The Jacobian matrix of the system is written in the form:

\(J=\left[\begin{array}{cc}{\left(\frac{\partial {R}_{i}^{(1)}}{\partial {\epsilon }_{j}^{\mathit{el}}}\right)}_{i=\mathrm{1,6};j=\mathrm{1,6}}& {\left(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega }_{s}}\right)}_{i=\mathrm{1,6};s=\mathrm{1,12}}\\ {\left(\frac{\partial {R}_{s}^{(2)}}{\partial {\epsilon }_{j}^{\mathit{el}}}\right)}_{s=\mathrm{1,12};j=\mathrm{1,6}}& {\left(\frac{\partial {R}_{r}^{(2)}}{\partial \Delta {\omega }_{s}}\right)}_{r=\mathrm{1,12};s=\mathrm{1,12}}\end{array}\right]\) of dimension \(18\mathrm{\times }18\)

\(\frac{\partial {R}_{i}^{(1)}}{\partial {\epsilon }_{j}^{\mathit{el}}}={\delta }_{\mathit{ij}}+\sum _{s}{\mu }_{i}^{s}\frac{\partial {\tau }^{s}}{\partial {\sigma }_{j}}={\delta }_{\mathit{ij}}+\sum _{s}\frac{\partial \Delta {\gamma }^{s}}{\partial {\tau }^{s}}{\mu }_{i}^{s}{\mu }_{k}^{s}{\Lambda }_{\mathit{kj}}\)

with \(\frac{\partial \Delta {\gamma }^{s}}{\partial {\tau }^{s}}={(\Delta {\gamma }^{s})}^{2}\left[\frac{1}{{(\Delta {\gamma }_{\mathit{nuc}}^{s})}^{2}}\frac{\partial \Delta {\gamma }_{\mathit{nuc}}^{s}}{\partial {\tau }_{s}}+\frac{1}{{(\Delta {\gamma }_{\mathit{prob}}^{s})}^{2}}\frac{\partial \Delta {\gamma }_{\mathit{prob}}^{s}}{\partial {\tau }_{s}}\right]\)

\(\frac{\partial \Delta {\gamma }_{\mathit{nuc}}^{s}}{\partial {\tau }_{s}}=\frac{\Delta {\gamma }_{\mathit{nuc}}^{s}\Delta {G}_{0}}{2{k}_{b}T\sqrt{{\tau }_{0}{\tau }_{\mathit{eff}}}}\mathit{sgn}({\tau }^{s})\) \(\mathit{si}{\tau }_{\mathit{eff}}>0\), 0 otherwise

\(\frac{\partial \Delta {\gamma }_{\mathit{prob}}^{s}}{\partial {\tau }_{s}}=\frac{n\Delta {\gamma }_{\mathit{prob}}^{s}}{∣{\tau }_{s}∣}\)

\(\frac{\partial {{R}^{s}}^{(2)}}{\partial {\epsilon }_{i}^{\mathit{el}}}=-\left[\frac{\partial \Delta {p}^{s}}{\partial {\tau }^{s}}{H}^{s}+\Delta {p}^{s}\frac{\partial {H}^{s}}{\partial {\tau }^{s}}\right]{\mu }_{k}^{s}{\Lambda }_{\mathit{kj}}\)

with \(\frac{\partial \Delta {p}_{s}}{\partial {\tau }_{s}}=\frac{\partial \Delta {p}^{s}}{\partial \Delta {\gamma }^{s}}\frac{\partial \Delta {\gamma }^{s}}{\partial {\tau }_{s}}=\frac{\Delta {\gamma }^{s}}{∣\Delta {\gamma }^{s}∣}\frac{\partial \Delta {\gamma }^{s}}{\partial {\tau }_{s}}\) \(\frac{\partial \Delta {\gamma }^{s}}{\partial {\tau }_{s}}\) already calculated previously.

and \(\frac{\partial {H}^{s}}{\partial {\tau }^{s}}=\left(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau }^{s}}\frac{\sqrt{{h}^{\mathit{ss}}{\omega }^{s}}}{{K}_{\mathit{self}}}+\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau }^{s}}\frac{{\alpha }_{\mathit{AT}}^{s}(\omega ){\lambda }_{s}(\omega ){\omega }_{\mathit{tot}}^{s}(\omega )}{b{K}_{f}}-\frac{\partial {y}^{s}}{\partial {\tau }^{s}}\frac{{\omega }^{s}}{b}\right)\)

\(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial {\tau }^{s}}=\frac{-\partial {\tau }_{\mathit{eff}}^{s}}{\partial {\tau }^{s}}\frac{1}{{\tau }_{0}}=\frac{-{\tau }^{s}}{∣{\tau }^{s}∣}\frac{1}{{\tau }_{0}}=\frac{-\mathit{sgn}({\tau }^{s})}{{\tau }_{0}}\) if \({\tau }_{\mathit{eff}}^{s}>0\), 0 otherwise

\(\frac{\partial {y}^{s}}{\partial {\tau }^{s}}=-{({y}^{s})}^{2}\frac{2\pi }{\mu b}\frac{\partial {\tau }_{\mathit{eff}}^{s}}{\partial {\tau }^{s}}=-{({y}^{s})}^{2}\frac{2\pi }{\mu b}\mathit{sgn}({\tau }^{s})\)

\(\frac{\partial {R}_{i}^{(1)}}{\partial \Delta {\omega }^{s}}=\sum _{r=\mathrm{1,12}}{\mu }_{i}^{r}\frac{\partial \Delta {\gamma }^{r}}{\partial \Delta {\omega }^{s}}\)

with \(\frac{\partial \Delta {\gamma }^{r}}{\partial \Delta {\omega }^{s}}={(\Delta {\gamma }^{r})}^{2}\left[\frac{1}{{(\Delta {\gamma }_{\mathit{nuc}}^{r})}^{2}}\frac{\partial \Delta {\gamma }_{\mathit{nuc}}^{r}}{\partial \Delta {\omega }^{s}}+\frac{1}{{(\Delta {\gamma }_{\mathit{prob}}^{r})}^{2}}\frac{\partial \Delta {\gamma }_{\mathit{prob}}^{r}}{\partial \Delta {\omega }^{s}}\right]\)

The various terms used above are:

\(\frac{\partial \Delta {\gamma }_{\mathit{nuc}}^{r}}{\partial \Delta {\omega }^{s}}=\Delta {\gamma }_{\mathit{nuc}}^{r}\left[\frac{1}{{l}^{r}}\frac{\partial {l}^{r}}{\partial \Delta {\omega }^{s}}-\frac{1}{{k}_{b}T}\frac{\partial \Delta G}{\partial \Delta {\omega }^{s}}\right]\)

with

\(\frac{\partial {l}^{r}}{\partial \Delta {\omega }^{s}}=\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }^{s}}-2\frac{\partial {\alpha }_{\mathit{AT}}^{r}}{\partial \Delta {\omega }^{s}}{R}^{r}-2{\alpha }_{\mathit{AT}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega }^{s}}\mathit{si}{\lambda }^{s}(\omega )-2{\alpha }_{\mathit{AT}}^{s}(\omega ){R}^{s}(\omega )>{l}_{c};0\mathit{sinon}\)

The calculation of \(\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }^{s}}\) is divided into two possibilities:

If \(\mathit{min}\left(\frac{\sqrt{{\omega }_{\mathit{tot}}^{r}(\omega )}}{b};(D+2{R}^{r}(\omega ))\frac{{\omega }_{\mathit{tot}}^{r}(\omega )}{{b}^{2}}\right)=\frac{\sqrt{{\omega }_{\mathit{tot}}^{r}(\omega )}}{b}\) then \(\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }^{s}}=\frac{-b}{2{({\omega }^{r})}^{3/2}}(1-{\delta }_{\mathit{rs}})\)

If not \(\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }^{s}}=\frac{-{b}^{2}}{{(D+2{R}^{r})}^{2}{({\omega }_{\mathit{tot}}^{r})}^{2}}\left[2{\omega }_{\mathit{tot}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega }^{s}}+(D+2{R}^{r})(1-{\delta }_{\mathit{rs}})\right]\)

\(\frac{\partial {\alpha }_{\mathit{AT}}^{r}}{\partial \Delta {\omega }_{s}}=\frac{1}{2{\alpha }_{\mathit{AT}}^{r}}\frac{1-{\delta }_{\mathit{rs}}}{{\omega }_{\mathit{tot}}^{r}}\left[{h}^{\mathit{rs}}-{({\alpha }_{\mathit{AT}}^{r})}^{2}\right]\)

\(\frac{\partial {R}^{r}}{\partial \Delta {\omega }_{s}}=\frac{2{R}^{r}}{\Delta {G}_{0}-\Delta {G}_{\mathit{app}}}\frac{\partial \Delta {G}_{\mathit{app}}}{\partial \Delta {\omega }_{s}}\)

\(\frac{\partial \Delta {G}_{\mathit{app}}}{\partial \Delta {\omega }_{s}}=\frac{-kT(1-{\delta }_{\mathit{rs}})}{2{\omega }_{\mathit{tot}}^{r}}\)

\(\frac{\partial \Delta G}{\partial \Delta {\omega }^{s}}=\frac{\Delta {G}_{0}}{2\sqrt{{\tau }_{0}{\tau }_{\mathit{eff}}}}\frac{\partial {\tau }_{c}^{r}}{\partial \Delta {\omega }^{s}}\)

\(\frac{\partial \Delta {\gamma }_{\mathit{prob}}^{r}}{\partial {\tau }^{s}}=-\Delta {\gamma }_{\mathit{prob}}^{r}\frac{n}{∣{\tau }_{c}^{r}∣}\frac{\partial {\tau }_{c}^{r}}{\partial \Delta {\omega }^{s}}\)

with \(\frac{\partial {\tau }_{c}^{r}}{\partial \Delta {\omega }^{s}}=\frac{{\tau }_{\mathit{LR}}^{r}\frac{\partial {\tau }_{\mathit{LR}}^{r}}{\partial \Delta {\omega }^{s}}+{\tau }_{\text{LT}}^{r}\frac{\partial {\tau }_{\text{LT}}^{r}}{\partial \Delta {\omega }^{s}}}{{\tau }_{c}^{r}-{\tau }_{F}}\)

\(\frac{\partial {\tau }_{\mathit{LR}}^{r}}{\partial \Delta {\omega }^{s}}={\delta }_{\mathit{rs}}\frac{\mu }{2}\sqrt{\frac{{h}^{\mathit{rr}}}{{\omega }^{r}}}=\frac{{\delta }_{\mathit{rs}}}{2}\frac{{\tau }_{\mathit{LR}}^{r}}{{\omega }^{r}}\)

\(\frac{\partial {\tau }_{\text{LT}}^{r}}{\partial \Delta {\omega }^{s}}=\frac{\partial {\alpha }_{\mathit{AT}}^{r}}{\partial \Delta {\omega }_{s}}\frac{{\tau }_{\text{LT}}^{r}}{{\alpha }_{\mathit{AT}}^{r}}-\mu b\frac{{\alpha }_{\mathit{AT}}^{r}}{{({\lambda }^{r})}^{2}}\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }^{s}}+2\mu b\frac{{\alpha }_{\mathit{AT}}^{r}}{{(2{\alpha }_{\mathit{AT}}^{r}{R}^{r}+{L}_{c})}^{2}}\left(\frac{\partial {\alpha }_{\mathit{AT}}^{r}}{\partial \Delta {\omega }_{s}}{R}^{r}+{\alpha }_{\mathit{AT}}^{r}\frac{\partial {R}^{r}}{\partial \Delta {\omega }_{s}}\right)\) if \({\tau }_{\text{LT}}^{s}>0\), 0

If not

\(\frac{\partial {R}_{r}^{(2)}}{\partial \Delta {\omega }^{s}}={\delta }_{\mathit{rs}}-\Delta {p}_{r}\frac{\partial {H}_{r}}{\partial \Delta {\omega }_{s}}-{H}_{r}\frac{\partial \Delta {p}_{r}}{\partial \Delta {\omega }_{s}}\)

\(\frac{\partial \Delta {p}_{r}}{\partial \Delta {\omega }_{s}}=\frac{{\gamma }^{r}}{∣{\gamma }^{r}∣}\frac{\partial \Delta {\gamma }^{r}}{\partial \Delta {\omega }_{s}}\) \(\frac{\partial \Delta {\gamma }^{s}}{\partial \Delta {\omega }_{s}}\) has already been calculated previously.

It remains to calculate:

\(\begin{array}{c}\frac{\partial {H}_{r}}{\partial \Delta {\omega }_{s}}=\frac{\partial {c}_{\mathit{eff}}^{r}}{\partial \Delta {\omega }_{s}}\frac{\sqrt{{h}^{\mathit{rr}}{\omega }^{r}}}{{K}_{\mathit{self}}}+\frac{{c}_{\mathit{eff}}^{r}\sqrt{{h}^{\mathit{rr}}}{\delta }_{\mathit{rs}}}{2\sqrt{{\omega }^{r}}{K}_{\mathit{self}}}-\frac{\partial {y}^{r}}{\partial \Delta {\omega }_{s}}\frac{{\omega }^{s}}{b}-\frac{{y}^{r}}{b}{\delta }_{\mathit{rs}}\\ \frac{+1}{b{K}_{f}}\left(\frac{\partial {c}_{\mathit{eff}}^{r}}{\partial \Delta {\omega }_{s}}{\alpha }_{\mathit{AT}}^{r}{\lambda }^{r}{\omega }_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}\frac{\partial {\alpha }_{\mathit{AT}}^{r}}{\partial \Delta {\omega }_{s}}{\lambda }^{r}{\omega }_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}{\alpha }_{\mathit{AT}}^{r}\frac{\partial {\lambda }^{r}}{\partial \Delta {\omega }_{s}}{\omega }_{\mathit{tot}}^{r}+{c}_{\mathit{eff}}^{r}{\alpha }_{\mathit{AT}}^{r}{\lambda }^{r}(1-{\delta }_{\mathit{rs}})\right)\end{array}\)

with:

\(\frac{\partial {c}_{\mathit{eff}}^{s}}{\partial \Delta {\omega }^{s}}=\frac{1}{{\tau }_{0}}\frac{\partial {\tau }_{c}^{r}}{\partial \Delta {\omega }^{s}}\)

\(\frac{\partial {y}^{r}}{\partial \Delta {\omega }^{s}}={({y}^{r})}^{2}\frac{2\pi }{\mu b}\frac{\partial {\tau }_{c}^{r}}{\partial \Delta {\omega }^{s}}\)

3.2.5. Implicit integration algorithm in large deformations#

From \({U}_{n}\) and \(\Delta U\), we can calculate \({F}_{n}\), \(\Delta F\) and \({F}_{n+1}=\Delta F{F}_{n}\). In addition, we know the results of the time increment \(n\): Cauchy constraints \({\sigma }_{n}\) internal variables.

It is a question of determining \({\sigma }_{n+1}\) and the internal variables \({({\gamma }_{s})}_{n+1}\) (or \({({\alpha }_{s})}_{n+1}\)) solutions of the system of \(6+{n}_{s}\) equations:

\({S}_{n+1}=\Lambda \mathrm{.}\frac{1}{2}({{F}_{n+1}^{e}}^{T}{F}_{n+1}^{e}–{I}_{d})\)

and for each sliding system:

\({(\Delta {\gamma }_{s})}_{n+1}=\Delta tg({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\) or \(\Delta {\alpha }_{s}=\Delta th({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)

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

Indeed, according to 7 and 13, Piola-Kirchhoff stresses measured from the intermediate configuration are considered here, so only the elastic part of the transformation gradient is considered.

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

Several choices are possible for the calculation of \({(\Delta {F}^{p})}^{\text{-1}}\):

\({(\Delta {F}^{p})}^{\text{-1}}\mathrm{=}\mathrm{exp}(\mathrm{-}\mathrm{\sum }\Delta {\gamma }_{s}(S,{\omega }_{s}){m}_{s}\mathrm{\otimes }{n}_{s})\) 12 which has the advantage of keeping the zero trace of plastic deformations, but which can present an additional cost;
\({(\Delta {F}^{p})}^{\text{-1}}={({I}_{d}+\sum \Delta {\gamma }_{s}(S,{\omega }_{s}){m}_{s}\otimes {n}_{s})}^{\text{-1}}{\left(\mathit{det}({I}_{d}+\sum \Delta {\gamma }_{s}(S,{\omega }_{s}){m}_{s}\otimes {n}_{s})\right)}^{-\frac{1}{3}}\) 13
\({(\Delta {F}^{p})}^{\text{-1}}=({I}_{d}-\sum \Delta {\gamma }_{s}(S,{\omega }_{s}){m}_{s}\otimes {n}_{s}){\left(\mathit{det}({I}_{d}-\sum \Delta {\gamma }_{s}(S,{\omega }_{s}){m}_{s}\otimes {n}_{s})\right)}^{-\frac{1}{3}}\)

\({F}^{\text{e}}\) gradients are stored as internal variables (by removing the identity tensor)

In post-processing we calculate the Cauchy constraints: \({\sigma }_{n+1}=\frac{1}{\mathrm{det}{F}_{n+1}^{e}}{F}_{n+1}^{e}{(S)}_{n+1}{{F}_{n+1}^{e}}^{\text{T}}\)

or Kirchhoff’s constraints: \({\tau }_{n+1}=\mathrm{det}{F}_{n+1}{\sigma }_{n+1}={F}_{n+1}^{e}{(S)}_{n+1}{{F}_{n+1}^{e}}^{\text{T}}\)

Car \(\mathrm{det}{F}_{n+1}=\mathrm{det}{F}_{n+1}^{e}\mathrm{det}{F}_{n+1}^{p}=\mathrm{det}{F}_{n+1}^{e}\)

The purpose of the problem to be solved at each integration point is therefore to know the internal variables. \({({\gamma }_{s})}_{n}\) \({({\alpha }_{s})}_{n}\) \({({p}_{s})}_{n}\) \({F}_{n}^{e}\) to determine \({\sigma }_{n+1}\) and the internal variables \({({\gamma }_{s})}_{n+1}\) (\({({\alpha }_{s})}_{n+1}\), \({F}_{n+1}^{e}\)) solutions of the system to \(6+{n}_{s}\) equations:

\({R}_{1}={(\Lambda )}^{\text{-1}}{S}_{n+1}-\frac{1}{2}({{F}_{n+1}^{e}}^{T}{F}_{n+1}^{e}–{I}_{d})=0\)

and for each sliding system:

\({R}_{2}={(\Delta {\gamma }_{s})}_{n+1}-\Delta tg({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})=0\) or \({R}_{2}=\Delta {\alpha }_{s}-\Delta t\cdot h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})=0\)

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

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

The resolution of this system of nonlinear equations can be done classically by Newton’s algorithm [R5.03.14].

The unknowns vector is composed of: \({S}_{n+1}\) (6 components) and \(\Delta {\alpha }_{s}\) (\({n}_{s}\) components).

The calculation of the Jacobian matrix is described in the appendix.

Initialization,

The algorithm is initialized as follows:

We ask \({F}_{n+1}^{\text{e-trial}}={F}_{n+1}{({F}_{n}^{p})}^{\text{-1}}\) which is the same as assuming that \({F}_{n+1}^{p}=\Delta {F}^{p}{F}_{n}^{p}={I}_{d}{F}_{n}^{p}\)

We calculate \({R}_{1}\) with \({F}_{n+1}^{\text{e}}={F}_{n+1}^{\text{e-trial}}\) and then all the equations \({R}_{2}\) If no threshold is reached, then \({F}_{n+1}^{\text{e}}={F}_{n+1}^{\text{e-trial}}\). Otherwise you have to solve system \({R}_{1}\) \({R}_{2}\).

We can choose to keep the tensor \({F}_{n+1}^{e}\) as an internal variable to more easily calculate \({F}_{n+1}^{e}=\Delta F{F}_{n}^{e}{(\Delta {F}^{p})}^{\text{-1}}\)

In practice we store \({F}^{e}-{I}_{d}\) in the internal variables \(V(6+3\ast {n}_{s}+1)\) to \(V(6+3\ast {n}_{s}+9)\)

We also keep \({F}^{p}-{I}_{d}\) in the internal variables \(V(6+3\ast {n}_{s}+10)\) to \(V(6+3\ast {n}_{s}+19)\)

In post-processing, tensors \({m}_{s}\) and \({n}_{s}\) can be updated using the rotation resulting from polar decomposition \({F}^{e}={R}^{e}{U}^{e}\): \({m}_{s}={R}^{e}{{m}_{s}}_{0}\) \({n}_{s}={R}^{e}{{n}_{s}}_{0}\)

For the calculation of polar decomposition, an exact algorithm is proposed in 12.

Next 13, they can also be carried in the current configuration by:

\({m}_{s}={F}^{e}{{m}_{s}}_{0}\) and \({n}_{s}={F}^{e-T}{{n}_{s}}_{0}\)

3.2.6. Convergence criteria used for implicit resolutions#

The criterion for stopping iterations relates to the relative nullity of the residue: \(\frac{{∣∣R(Y)∣∣}_{h}}{{∣∣{R}_{\mathrm{ref}}∣∣}_{h}}<{\varepsilon }_{c}\). where \({\varepsilon }_{c}\) represents the tolerance given by RESI_INTE_RELA. The problem is choosing \({R}_{\mathrm{ref}}\) well as the \({∣∣\mathrm{.}∣∣}_{h}\) standard

For consistency with Newton’s global algorithm used in STAT_NON_LINE [R5.03.02], and to avoid unnecessary iterations when the initial residue is very small, we choose:

\(\frac{\underset{i=\mathrm{1,6}}{\mathrm{max}}∣{R}_{i}(Y)∣}{\underset{i=\mathrm{1,6}}{\mathrm{max}}∣{E}_{i}^{\mathrm{tr}}∣}<{\varepsilon }_{c}\) with \({E}^{\mathrm{tr}}={\Lambda }^{\text{-1}}{\Sigma }^{\text{-}}+\Delta E-\Delta {E}^{\mathrm{th}}\)

\(\frac{\underset{i=\mathrm{7,}n}{\mathrm{max}}∣{R}_{i}(Y)∣}{\underset{i=\mathrm{7,}n}{\mathrm{max}}∣{Y}_{i}^{\text{-}}+\Delta {Y}_{i}∣}<{\varepsilon }_{c}\)

The two previous criteria must be verified in order to obtain convergence.

If convergence is not reached after the maximum number of iterations, the stationarity of the solution is also tested:

\(\mid {Y}_{k+1}-{Y}_{k}\mid <\varepsilon\)

The method used allows a local re-division of the time step, either systematic or in case of non-convergence (keyword ITER_INTE_PAS): this allows easier integration, but then losing the benefit of a coherent tangent operator at the end of the integration of the behavior.

It is often preferable, in case of local non-convergence after ITER_INTE_MAXI iterations, to proceed with a global redistribution of the time step (using DEFI_LIST_INST), in order to maintain a coherent tangent matrix.

3.3. Explicit resolution#

Another solution method, which is very simple to implement to solve the differential equations of monocrystalline behavior, is explicit resolution. For it to be digitally effective, it is essential to associate it with automatic step control. As in [R5.03.14], we use the Runge and Kutta method. The calculation of internal variables at the moment

is only a function of the values of their derivatives \(\frac{\text{dY}}{\text{dt}}=F(Y,t)\):

\(\Delta Y=\{\begin{array}{}\Delta {\alpha }_{s}\\ \Delta {\gamma }_{s}\\ \Delta {p}_{s}\\ \Delta {E}^{\mathrm{vp}}\end{array}=\{\begin{array}{}h({\tau }_{s}^{\text{-}},{\alpha }_{s}^{\text{-}},{\gamma }_{s}^{\text{-}},{p}_{s}^{\text{-}})\\ g({\tau }_{s}^{\text{-}},{\alpha }_{s}^{\text{-}},{\gamma }_{s}^{\text{-}},{p}_{s}^{\text{-}})\\ ∣\Delta {\gamma }_{s}∣=f({\tau }_{s}^{\text{-}},{\alpha }_{s}^{\text{-}},{\gamma }_{s}^{\text{-}},{p}_{s}^{\text{-}})\\ \sum _{s}{m}_{s}\Delta {\gamma }_{s}\end{array}\)

with \({\tau }_{s}=\Sigma :{\mu }_{s}=({\Sigma }^{\text{-}}+\Lambda (\Delta E-\Delta {E}^{\mathrm{th}}-\Delta {E}^{\mathrm{vp}})):{\mu }_{s}\)

We integrate according to the following diagram:

\({Y}_{t+h}={Y}^{\text{(2)}}\) if the accuracy criterion is met

\({Y}^{(2)}\mathrm{=}Y+\frac{h}{2}\left[F(Y,t)+F({Y}^{(1)},t+h)\right]\) with \({Y}^{(1)}=Y+hF(Y,t)\)

The difference between \({Y}^{(2)}\) (2nd order diagram) and \({Y}^{(1)}\) (1st order diagram, Euler) provides an estimate of the integration error and allows you to control the size of the time step

which is set to \(\Delta {t}_{i}\) for the first attempt. The step control strategy is defined based on a standard of the gap between the two integration methods: \(∥{Y}^{(2)}-{Y}^{(1)}∥\) and the precision required by the user \(\eta\) (keyword: RESI_INTE_RELA). The criterion used is as follows, where we note

\(\delta Y(t)=\underset{j=\mathrm{1,}N}{\text{sup}}\left\{\frac{∣{y}_{i}^{(2)}-{y}_{i}^{(1)}∣}{\mathrm{max}\left[\varepsilon ,∣{y}_{j}^{(t)}∣\right]}\right\}<\eta\)

Parameter \(\varepsilon\) is set to 0.001. The desired integration precision

should be consistent with the level of precision required for the overall stage.

If the criterion is not verified, the time step is re-divided according to a heuristic method (number of sub-steps defined by the user via the ITER_INTER_PAS keyword). When the time step becomes too short (h < 1.10—20), the calculation is stopped with an error message.

Explicit integration algorithm in large deformations

At each moment and at each integration point, we can calculate \({F}_{n}\), \(\Delta F\) and \({F}_{n+1}=\Delta F{F}_{n}\). We know the Cauchy constraints \({\sigma }_{n}\) and the internal variables at time \(n\)

It’s about determining \(S\) and the internal variables \(({\gamma }_{s})\), \(({\alpha }_{s})\), \({F}^{p}\) at time \(n+1\). To do this, we can use a Runge-Kutta method [cf. R5.03.14] on the system of differential equations:

\({\dot{F}}^{p}=(\sum _{s}\dot{{\gamma }_{s}}{m}_{s}\otimes {n}_{s}){F}^{p}\)

\(\dot{{\gamma }_{s}}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\) for each sliding system

\(\dot{{\alpha }_{s}}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\) for each sliding system

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

\((S)=\Lambda \mathrm{.}\frac{1}{2}({{F}^{e}}^{T}{F}^{e}–{I}_{d})\) and \({F}^{e}=F{({F}^{p})}^{\text{-1}}\)