2. Formulation of mono and polycrystalline behaviors#

2.1. Principle of the laws of behavior of single crystals#

The behaviors relating to single-crystal sliding systems are (in all the behaviors envisaged) of the elasto-visco-plastic type. They can be built on phenomenological bases or on physical bases (dynamics of dislocations). For each of the sliding directions, the behavior is*mondimensional*. It can be broken down into three types of equations:

Flow relationship: \(\dot{{\gamma }_{s}}\) represents the plastic sliding of the \(s\) system

\(\dot{{\gamma }_{s}}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\), with \(\dot{{p}_{s}}=∣\dot{{\gamma }_{s}}∣\) and:

elastoplastic case, a criterion of the type: \(F({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\le 0\) and \(F\cdot \dot{{p}_{s}}=0\)
elasto-viscoplastic case, \(\dot{{p}_{s}}=f({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)
\(\dot{{\alpha }_{s}}\) represents kinematic work hardening in the case of a phenomenological law and the dislocation density for physically-based behavior.

Its evolution is the form: \(\dot{{\alpha }_{s}}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)

for phenomenological laws, isotropic work hardening is defined by a function: \(R({p}_{s})\).

These relationships become, after discretization in time:

Flow relationship:

\(\Delta {\gamma }_{s}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\), with, \(\Delta {p}_{s}=∣\Delta {\gamma }_{s}∣\)

for elastoplastic behavior a criterion of the type: \(F({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\le 0\) and \(F\cdot \Delta {p}_{s}=0\)

and for elasto-viscoplastic behavior, \(\Delta {p}_{s}=f({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\Delta t\)

Evolution of \(\dot{{\alpha }_{s}}\): \(\Delta {\alpha }_{s}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)
Evolution of isotropic work hardening: \(R({p}_{s})\)

Quantities \(({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\) are evaluated at the current moment for an implicit discretization and at the previous moment for an explicit discretization.

The available monocrystalline laws are described below. The names of these relationships correspond to their names in the DEFI_MATERIAU command [U4.43.01].

2.2. Single crystal behavior laws available#

2.2.1. Phenomenological laws#

They are constructed from a flow law, an isotropic work hardening law and a kinematic work hardening law. The combination of the laws MONO_VISC1, MONO_ISOT1, MONO_CINE1 corresponds to the crystalline law known as Meric-Cailletaud 1.

2.2.1.1. Flow law MONO_VISC1#

\(\begin{array}{}\Delta {\gamma }_{s}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})=\Delta {p}_{s}\frac{{\tau }_{s}-c{\alpha }_{s}}{∣{\tau }_{s}-c{\alpha }_{s}∣}\\ \Delta {p}_{s}=\Delta t{\langle \frac{∣{\tau }_{s}-c{\alpha }_{s}∣-{R}_{s}({p}_{s})}{k}\rangle }^{n}\end{array}\)

The settings are: \((c,k,n)\).

2.2.1.2. Flow law MONO_VISC2#

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

The settings are: \(c,k,n,a,d\).

2.2.1.3. Kinematic work-hardening law MONO_CINE1#

\(\Delta {\alpha }_{s}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})=\Delta {\gamma }_{s}-d{\alpha }_{s}\Delta {p}_{s}\)

The setting is: \(d\).

2.2.1.4. Kinematic work-hardening law MONO_CINE2#

\(\Delta {\alpha }_{s}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})=\Delta {\gamma }_{s}-d{\alpha }_{s}\Delta {p}_{s}-{(\frac{∣c{\alpha }_{s}∣}{M})}^{m}\frac{{\alpha }_{s}}{∣{\alpha }_{s}∣}\)

The parameters are then: \(d,M,m,c\).

2.2.1.5. Isotropic work hardening law MONO_ISOT1#

\({R}_{s}(p)={R}_{0}+Q(\sum _{r=1}^{\mathrm{ns}}{h}_{\mathrm{sr}}(1-{e}^{{\mathrm{bp}}_{r}}))\); the parameters are: \({R}_{\mathrm{0,}}Q,b,{h}_{\mathrm{sr}}\).

2.2.1.6. Isotropic work hardening law MONO_ISOT2#

\({R}_{s}(p)={R}_{0}+{Q}_{1}\sum _{r}{h}_{\mathrm{sr}}(1-{e}^{-{b}_{1}{p}_{r}})+{Q}_{2}(1-{e}^{-{b}_{2}{p}_{s}})\); the parameters are: \({R}_{\mathrm{0,}}{Q}_{\mathrm{1,}}{b}_{\mathrm{1,}}{b}_{\mathrm{2,}}{h}_{\mathrm{sr}},{Q}_{2}\).

2.2.1.7. Interaction matrix#

\({h}_{\mathrm{sr}}\) a symmetric square matrix of interaction between sliding systems (see annex 4).

In the case where only one coefficient \(H\) is provided, the matrix \({h}_{\mathrm{sr}}\) takes a very simple form:

\({h}_{\mathrm{rs}}=1\mathrm{si}r=s\) and \({h}_{\mathrm{rs}}=H\mathrm{si}r\ne s\), with by default \(H=0\) (identity matrix).

2.2.2. Law resulting from the Dynamics of Dislocations: MONO_DD_KR#

Model MONO_DD_KR (Kocks-Rauch) is applicable to \(\mathit{CC}\) type steels (Centered Cubic Crystal Structure). Its main internal variable is the density of dislocations on each sliding system. It also offers the advantage of being valid for a wide range of temperatures [bib7]. The flow is formulated as follows:

If \(∣{\tau }^{s}∣>{\tau }_{0}\)

\({\tau }_{\mathrm{eff}}^{s}=\langle ∣{\tau }^{s}∣-{\tau }_{0}-{\tau }_{\mu }^{s}\rangle\) \(\langle \rangle\) positive part

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

otherwise \({\tau }_{\mathrm{eff}}^{s}=0\). In these expressions, \(\mu\) designates the shear modulus of the material, \(b\) is the Burgers constant, is the Burgers constant, \({\rho }^{u}\) is the dislocation density, \({a}^{\mathrm{su}}\) is the interaction matrix between the 24 sliding systems (its particular shape, depending on 5 coefficients, is provided in the appendix). The diagonal terms represent the self-working of a sliding system, the other terms represent latent work hardening.

If \(∣{\tau }^{s}∣>{\tau }_{0}\) and \({\tau }_{\mathrm{eff}}^{s}>0\)

\(\Delta {\gamma }_{s}=g({\tau }_{s},{\gamma }_{s})=\dot{{\gamma }_{0}}\mathrm{exp}(\frac{-\Delta G({\tau }_{\mathrm{eff}}^{s})}{{k}_{B}T})\cdot \frac{{\tau }_{s}}{∣{\tau }_{s}∣}\cdot \Delta t\)

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

If not

\(\Delta {\gamma }_{s}=0\)

\(\Delta {\alpha }_{s}=0\)

\(\Delta G({\tau }_{\mathrm{eff}}^{s})\) is an enthalpy activation corresponding to a free activation in an athermal domain. \({\gamma }_{\mathrm{0,}}{\tau }_{R},p,q\) is material data, \({k}_{B}\) is the Boltzman constant, and \(T\) is the temperature.

The evolution of the dislocation density is proportional to the square root of the sum of the dislocation densities of all the other sliding systems:

\(\dot{{\rho }_{s}}=\frac{∣\dot{{\gamma }^{s}}∣}{b}(\frac{1}{{\Lambda }^{s}}-{g}_{c}(T){\rho }^{s})\) with \(\frac{1}{{\Lambda }^{s}}=\frac{1}{d}+\frac{\sqrt{\sum _{u\ne s}{\rho }^{u}}}{K(T)}\)

and \({g}_{c}(T)={g}_{\mathrm{c0}}\mathrm{exp}\left[-\frac{{E}_{\mathrm{gc}}}{{k}_{B}T}\right]\) where \({E}_{\mathrm{gc}},{k}_{B},b,d,{g}_{\mathrm{c0}},K(T)\) are material parameters,

From where

\(\Delta {\rho }^{s}=∣\Delta {\gamma }^{s}∣(\frac{1}{\mathrm{bd}}+\frac{\sqrt{\sum _{u\ne s}{\alpha }^{u}}}{{b}^{2}K}-\frac{{g}_{c}{\alpha }^{s}}{{b}^{3}})\) and therefore:

\(\Delta {\alpha }^{s}=∣\Delta {\gamma }^{s}∣(\frac{b}{d}+\frac{\sqrt{\sum _{u\ne s}{\alpha }^{u}}}{K}-\frac{{g}_{c}{\alpha }^{s}}{b})=∣g({\tau }_{s},{\alpha }_{s})∣\cdot h({\alpha }_{s})\)

The settings are: \(\dot{{\gamma }_{\mathrm{0,}}}{k}_{B},T,{\tau }_{R},\Delta {G}_{\mathrm{0,}}p,q,{\tau }_{\mathrm{0,}}\mu ,b,d,K,{g}_{c}\)

\({k}_{B}\) is the Boltzman constant.

Note: Newton’s resolution can lead to negative values for \(\Delta {\alpha }^{u}\), so the terms under the root must be protected. To do this, we choose the following modification: \(\sqrt{\sum _{u\ne s}\langle {\alpha }^{u}\rangle }\) where \(\langle \text{.}\rangle\) represents the positive side: \(\langle x\rangle =\{\begin{array}{}x\text{si}x>0\\ 0\text{sinon}\end{array}\)

Note also that the expression for \(\Delta {\gamma }_{s}\), which is a function of \(\Delta G({\tau }_{\mathrm{eff}}^{s})\), leads to a discontinuity in the vicinity of \({\tau }_{\mathrm{eff}}^{s}=0\): for \({\tau }_{\mathrm{eff}}^{s}<0\), \(\Delta {\gamma }_{s}=0\), and \(\underset{{\tau }_{\mathrm{eff}}^{s}\to 0\text{*}}{\mathrm{lim}}\Delta {\gamma }_{s}=\dot{{\gamma }_{0}}\mathrm{exp}(\frac{-\Delta {G}_{0}}{{k}_{B}T})\cdot \frac{{\tau }_{s}}{∣{\tau }_{s}∣}\cdot \Delta t\)

2.2.3. Laws from the Dynamics of Dislocations: MONO_DD_CFC/MONO_DD_CFC_IRRA/MONO_DD_FAT#

The formulation of laws DD_CFC [8], [9], and DD_FAT [18] is constructed from calculations of dislocation dynamics. It applies to materials with a face-centered cubic crystalline structure (\(\mathit{CFC}\)) such as austenitic steels. In principle, model DD_CFC is not compatible with a change in loading path, in particular the parameters are not adapted to cyclic solicitation (when a dislocation « retraces its steps » its kinetics are different due to a different interaction with obstacles).

The equations below are written for each sliding system \(s\).

The direction of flow is given by:

\(\dot{{\gamma }_{s}}=\dot{{p}_{s}}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}=\dot{{p}_{s}}\text{signe}({\tau }_{s})\) and the intensity of the flow by:

\(\dot{{p}_{s}}=\dot{{\gamma }_{0}}({(\frac{∣{\tau }_{s}∣}{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}})}^{n}-1)\text{si}∣{\tau }_{s}∣\ge {\tau }_{f}+{\tau }_{s}^{\mathrm{forest}},\text{sinon}\dot{{p}_{s}}=0\)

with \(n\) large (to represent an almost elasto-plastic model: for CFC, the time dependence is in fact negligible).

Compared to the initial model [8], the threshold was introduced as well as the term \(–\dot{{\gamma }_{0}}\) to ensure continuity when the threshold is crossed.

In practice, the parameter \(\dot{{\gamma }_{0}}\) is small and the term \({(\frac{∣{\tau }_{s}∣}{{\tau }_{f}+{\tau }_{s}^{\mathrm{forest}}})}^{n}\) very quickly becomes much greater than \(1\), because of the large values of the exponent \(n\).

\({\tau }_{s}^{\mathrm{forest}}\) is a work hardening function whose coefficients are provided by the keyword MONO_DD_ *(*= CFC, FAT or CC), and which is described in the following §.

In DD_CFC, the evolution of the dislocation density is given by:

\(\dot{{\rho }_{s}}=\frac{\dot{{p}_{s}}}{b}(A\frac{\sum _{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}}{\rho }_{j}}{\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho }_{j}}}+B\sum _{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho }_{j}}-y{\rho }_{s})\)

\(j\in \mathrm{forest}(s)\) are the \({n}_{j}\ne {n}_{s}\) distinct normal systems
\(j\in \mathrm{copla}(s)\) are the same normal systems \({n}_{j}={n}_{s}\)

The initial value for the dislocation density is: \({\rho }_{s}(t=0)={\rho }_{s}^{0}\) (field provided by the user. Carefully check whether it is an overall density (for all 12 sliding systems) or a density per system).

Like \({a}_{\mathrm{ij}}^{\mathrm{eff}}=C{(\rho )}^{2}\mathrm{.}{a}_{\mathrm{ij}}\), \({\tau }_{s}^{\mathrm{forest}}(\rho )=\mu bC(\rho )\sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathrm{sj}}{\rho }_{j}}\), and \(\dot{{\rho }_{s}}=\frac{\dot{{p}_{s}}}{b}(A\frac{\sum _{j\in \mathrm{forest}(s)}\sqrt{{a}_{\mathrm{sj}}}{\rho }_{j}}{\sum _{j=\mathrm{1,12}}\sqrt{{a}_{\mathrm{sj}}{\rho }_{j}}}+BC(\rho )\sum _{j\in \mathrm{copla}(s)}\sqrt{{a}_{\mathrm{sj}}{\rho }_{j}}-y{\rho }_{s})\)

To avoid any problem with conditioning the Jacobian matrix, we choose to size the system by replacing the dislocation density by: \({\omega }_{s}={b}^{2}\ast {\rho }_{s}\) and solve in \({\omega }_{s}\). This gives:

\(\dot{{\omega }_{s}}=\dot{{p}_{s}}{h}_{s}(\omega )\text{avec}{h}_{s}(\omega )=(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})\)

NOTE: Physically, the dislocation density should remain positive. Numerically, for all the terms in \(\sqrt{\rho }\) to be legal, it is necessary to test the positivity of \(\rho\) and either recut the time step if this value becomes negative (by « hesitation » of the algorithm during the iterations) or introduce a positive part into the evolution equations. The latter solution is chosen for explicit integration as well as for implicit integration.

The set of differential equations governing viscoplastic flow in model DD_CFC can therefore be written as:

\(\dot{{\gamma }_{s}}=\dot{{p}_{s}}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\) (12 equations)

\(\dot{{\omega }_{s}}=\dot{{p}_{s}}{h}_{s}(\langle \omega \rangle )\) (12 equations) with: \(\langle {\omega }_{i}\rangle =\{\begin{array}{c}{\omega }_{i}\text{si}{\omega }_{i}>0\\ 0\text{sinon}\end{array}\)

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

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

\({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}\langle {\omega }_{s}\rangle )\)

with \({\omega }_{s}={b}^{2}\ast {\rho }_{s}\)

In case MONO_DD_CFC_IRRA, the evolution of the variables related to irradiation is given by:

\({\dot{\rho }}_{s}^{\mathit{loops}}=-\xi \left(\sum _{j\in \mathit{copla}(s)}∣{\dot{\gamma }}_{j}∣\right)\left({\rho }_{s}^{\mathit{loops}}-{\rho }^{\mathit{sat}}\right)\) \({\dot{\varphi }}_{s}^{\mathit{voids}}=-\zeta \left(\sum _{j\in \mathit{copla}(s)}∣{\dot{\gamma }}_{j}∣\right)\left({\varphi }_{s}^{\mathit{voids}}-{\varphi }^{\mathit{sat}}\right)\)

The expression for \({\tau }_{s}^{\mathit{forest}}\) becomes: \({\tau }_{s}^{\mathit{forest}}=\mu \sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathit{sj}}^{\mathit{eff}}{\omega }_{j}+{\alpha }^{\mathit{loops}}{\varphi }^{\mathit{loops}}{\rho }_{s}^{\mathit{loops}}{b}^{2}+{\alpha }^{\mathit{voids}}{\varphi }_{s}^{\mathit{voids}}{\rho }^{\mathit{voids}}{b}^{2}}\)

Which, after multiplying by \({b}^{\mathrm{2 }}\): \({\omega }_{s}^{\mathit{loops}}={b}^{2}\ast {\rho }_{s}^{\mathit{loops}}\) gives:

\({\dot{\omega }}_{s}^{\mathit{loops}}=-\xi \left(\sum _{j\in \mathit{copla}(s)}∣{\dot{\gamma }}_{j}∣\right)\left({\omega }_{s}^{\mathit{loops}}-{\omega }^{\mathit{sat}}\right)\) and \({\tau }_{s}^{\mathit{forest}}=\mu \sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathit{sj}}^{\mathit{eff}}{\omega }_{j}+{\alpha }^{\mathit{loops}}{\varphi }^{\mathit{loops}}{\omega }_{s}^{\mathit{loops}}+{\alpha }^{\mathit{voids}}{\varphi }_{s}^{\mathit{voids}}{\omega }^{\mathit{voids}}}\)

we can notice in passing that the evolution equations of \({\omega }_{s}^{\mathit{loop}}\) and \({\varphi }_{s}^{\mathit{voids}}\) integrate analytically:

\({\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)\)

In case DD_FAT, the evolution of the dislocation density is given by:

\(\dot{{\rho }_{s}}\mathrm{=}\frac{\dot{{p}_{s}}}{b}(\frac{1}{d}+\frac{\sqrt{\mathrm{\sum }_{u\mathrm{\ne }s}{\rho }_{u}}}{K}\mathrm{-}{g}_{\mathit{c0}}{\rho }_{s})\)

By dimensioning as above, that is to say by asking:

\({\omega }_{s}\mathrm{=}{b}^{2}\mathrm{\times }{\rho }_{s}\)

we get:

\(\dot{{\omega }_{s}}\mathrm{=}\dot{{p}_{s}}{h}_{s}(\omega )\) with \({h}_{s}(\omega )\mathrm{=}(\frac{b}{d}+\frac{\sqrt{\mathrm{\sum }_{u\mathrm{\ne }s}{\omega }_{u}}}{K}\mathrm{-}{g}_{\mathit{c0}}\frac{{\omega }_{s}}{b})\)

Hence by introducing the positive parts:

\(\dot{{\gamma }_{s}}=\dot{{p}_{s}}\frac{{\tau }_{s}}{∣{\tau }_{s}∣}\) (12 equations)

\(C(\omega )=1.\)

\({\tau }_{s}^{\mathrm{forest}}(\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})\)

with \({\omega }_{s}={b}^{2}\ast {\rho }_{s}\)

The work-hardening law for DD_CFC models is: \({\tau }_{s}^{\mathrm{forest}}=\mu b\sqrt{\sum _{j=\mathrm{1,}{n}_{s}}{a}_{\mathrm{sj}}^{\mathrm{eff}}{\rho }_{j}}\)

and \({a}_{\mathrm{ij}}^{\mathrm{eff}}=C{(\rho )}^{2}\mathrm{.}{a}_{\mathrm{ij}}\text{avec}C(\rho )=0.2+0.8\frac{\mathrm{ln}(\alpha b\sqrt{\sum _{i=\mathrm{1,12}}{\rho }_{i}})}{\mathrm{ln}(\alpha b\sqrt{{\rho }_{\mathrm{ref}}})}\)

where \({a}_{\mathrm{ij}}\) represents the interaction matrix between sliding systems, whose particular form is provided in appendix 4.

The work-hardening law of model DD_FAT is: \({\tau }_{s}^{\mathrm{forest}}(\rho )=\mu b\sqrt{\sum _{j=\mathrm{1,12}}{a}_{\mathrm{sj}}{\rho }_{j}}\), which corresponds to \(C(\rho )\mathrm{=}1.\) in the previous model

2.2.4. Laws from the Dynamics of Dislocations: MONO_DD_CC/MONO_DD_CC_IRRA#

The formulation of laws DD_CC (with or without irradiation) is constructed from calculations of dislocation dynamics [17]. It applies to materials with a Cubic Centered crystal structure (\(\mathit{CC}\))

The equations below are written for each sliding system \(s\).

Visco-plastic flow is obtained by harmonic averaging of two regimes (low temperature, and high temperature).

Here, we choose to show the unknowns \(\rho\) (representing the set of dislocation densities \({\rho }^{s}\text{et}{\rho }_{\mathit{irr}}^{s}\)) and \({\tau }_{\mathrm{eff}}^{s}\) (which is an additional unknown a priori. We will see that the proposed resolution algorithm allows everything to be expressed as a function of \(\rho\)).

\({\rho }_{\mathit{irr}}^{s}\) refers to the density of dislocations induced by irradiation (case of MONO_DD_CC_IRRA), and is not used for the MONO_DD_CC behavior in order to limit the number of internal variables.

\({\tau }^{s}\) being the split resolved on system \(s\):

plastic sliding is obtained by: \(\frac{1}{{\dot{\gamma }}^{s}}=\frac{1}{{\dot{\gamma }}_{\mathrm{nuc}}^{s}}+\frac{1}{{\dot{\gamma }}_{\mathrm{prob}}^{s}}\)
\({\dot{\gamma }}_{\mathrm{nuc}}^{S}={\rho }_{m}^{s}bH\cdot {l}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})\cdot \mathrm{exp}(\frac{-\Delta G({\tau }_{\mathrm{eff}}^{s})}{{k}_{B}T})\mathrm{sgn}({\tau }_{s})\)
\(\Delta G({\tau }_{\mathrm{eff}}^{s})=\Delta {G}_{0}(1-{(\frac{\langle {\tau }_{\mathrm{eff}}^{s}\rangle }{{\tau }_{0}})}^{0.5})\) (with: \({\tau }_{\mathrm{eff}}^{s}<{\tau }_{0}\) and \(0<\Delta G<\Delta {G}_{0}\))
\({\tau }_{\mathrm{eff}}^{s}({\tau }^{s},\rho )=∣{\tau }^{s}∣-{\tau }_{c}\) \({\tau }_{c}={\tau }_{F}+\sqrt{({\tau }_{\mathit{LR}}^{s}{(\rho )}^{2}+{\tau }_{\text{LT}}^{s}{(\rho ,{\tau }_{\mathit{eff}}^{s})}^{2})}\)
\({\tau }_{\mathrm{LR}}^{s}(\rho )=\mu b\sqrt{{a}_{\mathrm{AT}}^{s}(\rho ){\rho }^{s}}\)
\({\tau }_{\text{LT}}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})=\text{max}(0;{\alpha }_{\mathrm{AT}}^{s}(\rho )\mu b(\frac{1}{{\lambda }^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})}-\frac{1}{2{\alpha }_{\mathrm{AT}}^{s}(\rho ){R}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})+{l}_{c}(T)}))\)
\({\rho }_{\mathrm{tot}}^{s}(\rho )=\sum _{j\ne s}{\rho }^{j}+{\rho }_{\mathrm{irr}}^{s}\)
\({\alpha }_{\mathit{AT}}^{s}(\rho )=\sqrt{\sum _{j\ne s}{a}_{\mathit{AT}}^{\mathit{sj}}\frac{{\rho }^{j}}{{\rho }_{\mathit{tot}}^{s}}+{a}_{\mathit{irr}}\frac{{\rho }_{\mathit{irr}}^{s}}{{\rho }_{\mathit{tot}}^{s}}}\)
\(\frac{1}{{\lambda }^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})+D}=\mathrm{min}(\sqrt{{\rho }_{\mathrm{tot}}^{s}(\rho )};(D+2{R}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})){\rho }_{\mathrm{tot}}^{s}(\rho ))\)
\({R}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})=\frac{\mu b}{2{\tau }_{0}{(1-\frac{\Delta G({\tau }_{\mathrm{eff}}^{s})}{\Delta {G}_{0}})}^{2}}\)
\({l}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})=\mathrm{max}({\lambda }^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s})-2{\alpha }_{\mathrm{AT}}^{s}(\rho ){R}^{s}(\rho ,{\tau }_{\mathrm{eff}}^{s});{l}_{c})\)
\({\dot{\gamma }}_{\mathrm{prob}}^{s}=\dot{{\gamma }_{0}}{(\frac{∣{\tau }^{s}∣}{{\tau }_{c}})}^{n}\mathrm{sgn}({\tau }^{s})\)

The parameters are: \({\rho }_{m}^{s},b,H,\Delta {G}_{\mathrm{0,}}{\tau }_{\mathrm{0,}}{\tau }_{F},\mu ,{a}_{\mathit{irr}},{l}_{c}(T),D,\dot{{\gamma }_{\mathrm{0,}}}n\) and the interaction matrix \({\alpha }_{\mathrm{AT}}^{\mathrm{sr}}\)

The laws of evolution of dislocation densities (\({\rho }^{s}\) and \({\rho }_{\mathrm{irr}}^{s}\)) are: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~

\(\dot{{\rho }^{s}}=\frac{∣\dot{{\gamma }^{s}}∣}{b}\left(\frac{1}{{d}_{\mathit{lath}}}+\frac{{c}_{\mathit{eff}}\sqrt{{a}_{\mathit{AT}}^{\mathit{ss}}{\rho }^{s}}}{{K}_{\mathit{self}}}+\frac{{\alpha }_{\mathit{AT}}^{s}{c}_{\mathit{eff}}{\lambda }^{s}{\rho }^{s}}{{K}_{f}}-{y}^{s}{\rho }^{s}\right)\)
\(\dot{{\rho }_{\mathrm{irr}}^{s}}=-\xi {\rho }_{\mathrm{irr}}^{s}∣\dot{{\gamma }^{s}}∣\)

with:

\(\frac{1}{{y}^{s}}=\frac{1}{{y}_{\mathrm{AT}}^{s}}+\frac{2\pi {\tau }_{\mathrm{eff}}^{s}(\rho )}{\mu b}\) and \({c}_{\mathit{eff}}^{s}=\left(1-\frac{{\tau }_{\mathit{eff}}^{s}}{{\tau }_{0}}\right)\),

with parameters \({d}_{\mathit{lath}}{K}_{\mathit{self}}{K}_{f}\xi {y}_{\mathit{AT}}^{s}\)

Remarks:

Physically, the dislocation density should remain positive. Numerically, for all the terms in \(\sqrt{\rho }\) to be legal, it is necessary to test the positivity of \(\rho\) and either recut the time step if this value becomes negative (by « hesitation » of the algorithm during the iterations) or introduce a positive part into the evolution equations. The latter solution is chosen for explicit integration as well as for implicit integration.
In the case where the elasticity is not isotropic, the value of \(\mu\) used above is independent of the elasticity coefficients, because it is obtained by homogenization.

2.3. Sliding systems and the overall behavior of the single crystal#

A monocrystal is composed of one or more families of sliding systems (cubic, octahedral, basal, prismatic,…), each family comprising a certain number of systems (12 for the octahedral family for example).

Each family of sliding systems is associated with a flow law, a type of kinematic and isotropic work hardening, and values of the parameters for these laws. In other words, it is not expected to vary behavioral relationships or coefficients within the same family of sliding systems. On the other hand, from one family to another, the laws of behavior can change, as well as the values of the parameters.

A sliding system is determined by an orientation tensor (or Schmid tensor) constructed from crystallographic definitions of:

the sliding direction of system \(s\) (defined by the unit vector \({m}_{s}\))
and the normal to the sliding plane (defined by the unit vector \({n}_{s}\)).

2.3.1. Small deformations, initial configuration#

In small deformations, the orientation tensor is defined by:

\({\mu }_{\mathrm{ij}}^{s}=\frac{1}{2}({m}_{i}\otimes {n}_{j}+{n}_{j}\otimes {m}_{i})\)

From the point of view of behavior at the material point of view, this tensor is used to calculate the resolved cission.

\({\tau }^{s}=\Sigma :{\mu }^{s}\)

and that of the global viscoplastic deformation rate \({E}^{\mathrm{vp}}\), defined from the knowledge of sliding speeds \(\dot{{\gamma }^{s}}\) for all sliding systems:

\({\dot{E}}_{\mathrm{ij}}^{\mathrm{vp}}=\sum _{s\in g}\dot{{\gamma }^{s}}{\mu }_{\mathrm{ij}}^{s}\)

In addition, the monocrystal can be oriented with respect to the global axes for defining the coordinates. This orientation is defined for each mesh or group of elements (typically for each grain) by the data of 3 nautical angles or 3 Euler angles.

The components of the orientation tensor \({\mu }^{s}\), defined in the coordinate system linked to the monocrystal, are then expressed in the global coordinate system using these nautical angles. In small deformations, they are fixed in the initial configuration.

In addition, at the level of a Gauss point, we postulate an elasticity relationship on the global tensors:

total macroscopic deformation \(E\)
macroscopic viscoplastic deformation \({E}^{\mathrm{vp}}\)
macroscopic constraint: \(\Sigma\)

\(\Sigma =\Lambda (E-{E}^{\mathrm{th}}-{E}^{\mathrm{vp}})\)

where \(\Lambda\) represents the elasticity operator (isotropic or orthotropic)

2.3.2. Small deformations, rotation of the crystal lattice#

An approximation of the effect of deformation on the behavior of the single crystal consists in taking into account the rotation of the crystal lattice. This rotation is interesting for two reasons:

its effective consideration in the calculations is consistent with deformations that are not infinitely small (in current calculations, it is not uncommon to reach 10%). It is a preliminary to a complete model of large deformations on the single crystal.
This is interesting information for post-treatment, as it can be compared with experimental results.

We describe here (following 9, 10, 11) the method for taking into account the rotation of a crystal lattice in the general framework of small deformations. The main physical hypotheses are:

knowing the vorticity tensor \(\stackrel{ˆ}{\dot{\omega }}\) or rotation rate, which is the antisymmetric part of the tensor \(L=\dot{F}{F}^{\text{-1}}\)), an additive decomposition is postulated into a part called plastic rotation \(\stackrel{ˆ}{\dot{{\omega }^{\mathrm{p}}}}\) (antisymmetric part of the plastic rotation) and an « elastic » part \(\stackrel{ˆ}{\dot{{\omega }^{e}}}\), which will make it possible to calculate the rotation of the crystal lattice;
we assume for the moment that we can calculate the rotation rate tensor, while integrating the behavior with a linearized deformation tensor. This limitation falls naturally with the complete model in large deformations, but it is assumed here that one can reason in a decoupled manner. So, in incremental form, we have:

\(\Delta \stackrel{ˆ}{\omega }=\frac{1}{2}(\Delta F{F}^{\text{-1}}-{(\Delta F{F}^{\text{-1}})}^{T})\) and \(\Delta \stackrel{ˆ}{{\omega }^{e}}=\Delta \stackrel{ˆ}{\omega }-\Delta \stackrel{ˆ}{{\omega }^{p}}\) which are antisymmetric tensors.

The calculation of \(\Delta \stackrel{ˆ}{{\omega }^{p}}\) is performed at the end of integration: it is stored and used at the next iteration (a purely implicit calculation would involve solving the system of equations described in paragraph 3 by adding all the equations on the orientation of the network, which would greatly burden the calculation).

Once the rotation increment \(\Delta \stackrel{ˆ}{{\omega }^{e}}\) is known, we extract the « axial » vector [voir par exemple R5.03.40]:: \(\Delta {\omega }^{e}\) such as \(\Delta {\omega }^{e}\wedge U=\Delta \stackrel{ˆ}{{\omega }^{e}}\cdot U\) and we calculate the angle of rotation and the vector that corresponds to the axis of rotation:

\(\Delta \theta =\parallel \Delta {\omega }^{e}\parallel\) and \({n}_{\omega }=\frac{\Delta {\omega }^{e}}{\Delta \theta }\).

The rotation increment matrix is then obtained by the Euler-Rodrigues formula:

\(\Delta Q={I}_{d}+\text{sin}(\Delta \theta )\stackrel{ˆ}{{n}_{\omega }}+(1-\text{cos}(\Delta \theta ))\stackrel{ˆ}{{n}_{\omega }}\cdot \stackrel{ˆ}{{n}_{\omega }}\) where \(\stackrel{ˆ}{{n}_{\omega }}\) is the antisymmetric matrix from \({\mathrm{n}}_{\omega }\)

This allows you to update the network rotation matrix: \({Q}_{\text{+}}=\Delta Q\cdot {Q}_{\text{-}}\)

It allows you to update the orientation tensor by: \(\begin{array}{}{n}^{s}={Q}_{\text{+}}{n}_{0}^{s}\\ {m}^{s}={Q}_{\text{+}}{m}_{0}^{s}\end{array}\).

The behavior is then integrated with this definition of the crystal lattice.

It then remains to update the plastic rotation increment: \(\Delta \stackrel{ˆ}{{\omega }^{p}}=\frac{1}{2}{\sum }_{s}\Delta {\gamma }^{s}({m}^{s}\otimes {n}^{s}-{n}^{s}\otimes {m}^{s})\)

Additional internal variables are: \(Q,\Delta {\omega }^{p},\Delta {\omega }^{e},\theta\)

In practice, network rotation is taken into account via the ROTA_RESEAU keyword of the DEFI_COMPOR [U4.43.06] command. It is only effective for the implicit integration of MONOCRISTAL behavior.

2.3.3. Large deformations#

Expected properties: the formalism chosen must be objective, and must maintain plastic incompressibility (\(\mathrm{det}({F}^{p})=1\)). Multiplicative decomposition \(F={F}^{e}{F}^{p}\) is used with Mandel’s relaxed isocline intermediate configuration.

The law of (hyper-) elasticity chosen uses the second Piola-Kirchhoff tensor \(S={J}^{e}{F}^{e-1}\sigma {F}^{e-T}\) and the Green-Lagrange deformation associated with the elastic deformation gradient:

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

On each sliding system \(s\), the resolved cission is written with the Mandel tensor \(M\) defined by \(M={J}_{e}{{F}^{e}}^{T}\sigma {{F}^{e}}^{\text{-T}}\) in the following way: \({\tau }_{s}\mathrm{=}M\mathrm{:}{m}_{s}\mathrm{\otimes }{n}_{s}\) (cf. 7, 14).

Indeed, the power of inner efforts can be written as: \(\frac{1}{\rho }\sigma :D=\frac{1}{\rho }\sigma :L=\frac{1}{\rho }\sigma :\dot{F}{F}^{\text{-1}}\)

and \(\dot{F}{F}^{\text{-1}}=(\dot{{F}^{\text{e}}}{F}^{p}+{F}^{\text{e}}\dot{{F}^{p}}){({F}^{p})}^{\text{-1}}{({F}^{e})}^{\text{-1}}={L}^{\text{e}}+{F}^{\text{e}}{L}^{p}{({F}^{e})}^{\text{-1}}\)

So \(\frac{1}{\rho }\sigma :D=\frac{1}{\rho }\left[\sigma :{L}^{\text{e}}+\sigma :({F}^{\text{e}}{L}^{p}{({F}^{\text{e}})}^{\text{-1}})\right]=\frac{1}{{\rho }_{0}}S:\dot{{E}_{\mathrm{GL}}^{\text{e}}}+\frac{1}{{\rho }_{0}}M:{L}^{p}\)

The power dissipated plastically is written: \(M:{L}^{p}\) and also: \(\sum _{s}\dot{{\gamma }_{s}}{\tau }_{s}\), so \(M:{L}^{p}=M:\sum _{s}\dot{{\gamma }_{s}}{\mu }_{s}=\sum _{s}\dot{{\gamma }_{s}}M:{\mu }_{s}\), from where \({\tau }_{s}=M:{\mu }_{s}\)

For each law of crystalline behavior, the equations of flow and the evolution of the internal variables are classically written on each sliding system:

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

On the other hand, the relationship linking plastic deformation to the sliding of active systems is now written: \({\dot{F}}^{p}=(\sum _{s}\dot{{\gamma }_{s}}{m}_{s}\otimes {n}_{s}){F}^{p}\)

2.4. Behaviour of homogenized polycrystal#

In the case of homogenized polycrystal, each mono-crystalline phase must be defined by its orientation, its proportion (volume fraction) and the associated behavior. An overall orientation of the polycrystal can be defined. In addition, a location rule must be defined.

The monocrystalline behavior is constructed as before from the previous elasto-visco-plastic behavior and from the data of families of sliding systems.

In addition to the monocrystalline behavior described above, a modeling scale is added, which represents the assembly of the phases.

At the level of a Gauss point, we postulate an elasticity relationship on the global tensors:

total macroscopic deformation \(E\)
macroscopic viscoplastic deformation \({E}^{\mathrm{VP}}\)
macroscopic constraint: \(\Sigma\)

\(\Sigma =\Lambda (E-{E}^{\mathrm{th}}-{E}^{\mathrm{VP}})\)

\(\Lambda\) represents the elasticity operator (isotropic or orthotropic)

In addition, knowing all the internal variables 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,
for each monocrystalline phase (or « grain »), defined by an orientation and a \({f}_{g}\) proportion, a stress location relationship, of the general form (to be expressed in the local coordinate system of each phase)

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

and for each sliding system and each phase, behavioral relationships relating to each sliding system, similar to the case of the monocrystal:

Flow relationship:

\(\dot{{\gamma }_{s}}=g({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\), with, \(\dot{{p}_{s}}=∣\dot{{\gamma }_{s}}∣\) and \(\dot{{p}_{s}}=f({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)

Evolution of kinematic work hardening or dislocation density: \(\dot{{\alpha }_{s}}=h({\tau }_{s},{\alpha }_{s},{\gamma }_{s},{p}_{s})\)
Evolution of isotropic work hardening defined by a function: \(R({p}_{s})\)
Viscoplastic deformations of the phase: \({\dot{\varepsilon }}_{\mathrm{ij}}^{v{p}_{g}}=\sum _{s\in g}\dot{{\gamma }^{s}}{\mu }_{\mathrm{ij}}^{s}\)

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

Two \({\sigma }_{g}=L(\Sigma ,{E}^{\mathrm{vp}},{\varepsilon }_{g}^{\mathrm{vp}},{\beta }_{g})\) location relationships are available:

The Berveiller-Zaoui relationship [bib5] based on the concept of self-coherence. This relationship is validated under certain conditions, namely: isotropy of the material (recent work uses extensions to cubic elasticity), homogeneous elastic behavior and monotonic loading:

\({\sigma }_{\mathit{ij}}^{g}={\Sigma }_{\mathit{ij}}+\alpha {\mu }^{\mathit{loca}}\left({E}_{\mathit{ij}}^{\mathit{vp}}-{\epsilon }_{\mathit{ij}}^{v{p}_{g}}\right)\) \(\frac{1}{\alpha }=1+\frac{3}{2}{\mu }^{\mathit{loca}}\frac{∥{E}_{\mathit{ij}}^{\mathit{vp}}∥}{{J}_{2}({\Sigma }_{\mathit{ij}})}\)

The second relationship, inspired by the previous one, and developed more particularly for cyclic loadings [bib4] makes it possible to give a good description to schematize the interactions between grains:
\({\sigma }_{\mathit{ij}}^{g}={\Sigma }_{\mathit{ij}}+{\mu }^{\mathit{loca}}\left({B}_{\mathit{ij}}-{\beta }_{\mathit{ij}}^{g}\right)\) \({B}_{\mathrm{ij}}=\sum _{g}{f}_{g}{\beta }_{\mathrm{ij}}^{g}\)
\({\dot{\beta }}_{\mathrm{ij}}^{g}={\dot{\varepsilon }}_{\mathrm{ij}}^{v{p}_{g}}-D({\beta }_{\mathrm{ij}}^{g}-\delta {\varepsilon }_{\mathrm{ij}}^{v{p}_{g}})∥{\dot{\varepsilon }}_{\mathrm{ij}}^{v{p}_{g}}∥\)

where \({\mu }^{\mathit{loca}}\), \(D\), and \(\delta\) are characteristic parameters of the material and of the temperature.