This article shows how to implement an isotropic viscoplastic behaviour combining isotropic hardening and multiple kinematic hardenings following an Armstrong-Frederic evolution of the back stress.

Such an example illustrates:

• The usage of StandardElasticity brick (see this page).

The whole implementation is available here.

Description of the behaviour

The behaviour is described by a standard split of the strain \(\underline{\varepsilon}^{\mathrm{to}}\) in an elastic and a plastic parts, respectively denoted \(\underline{\varepsilon}^{\mathrm{el}}\) and \(\underline{\varepsilon}^{\mathrm{vis}}\):

\[ \underline{\varepsilon}^{\mathrm{to}}=\underline{\varepsilon}^{\mathrm{el}}+\underline{\varepsilon}^{\mathrm{vis}} \]

Elastic behaviour

The stress \(\underline{\sigma}\) is related to the the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) by a the standard Hooke behaviour:

\[ \underline{\sigma}= \lambda\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,\underline{\varepsilon}^{\mathrm{el}} \]

Viscoplastic behaviour

The viscoplastic behaviour follows a standard viscoplastic behaviour: \[ \underline{\dot{\varepsilon}}^{\mathrm{vis}}=\left\langle{{\displaystyle \frac{\displaystyle F}{\displaystyle K}}}\right\rangle^{m}\,\underline{n}=\dot{p}\,\underline{n} \]

where \(F\) is the yield surface defined below, \(<.>\) is Macaulay brackets, \(\underline{n}\) is the normal to \(F\) with respect to the stress and \(p\) is the equivalent plastic strain.

The yield surface is defined by: \[ F{\left(\underline{\sigma},\underline{X}_{i},p\right)}={\left(\underline{\sigma}-\sum_{i=1}^{N}\underline{X}_{i}\right)}_{\mathrm{eq}}-R{\left(p\right)}=\underline{\sigma}^{e}_{\mathrm{eq}}-R{\left(p\right)} \]

where:

• \(R{\left(p\right)}\) describes the isotropic hardening as a function of the equivalent viscoplastic strain \(p\).
• the \(N\) tensors \(\underline{X}_{i}\) (i) are backstresses describing the kinematic hardening.
• \({\left(.\right)}_{\mathrm{eq}}\) is the Von Mises norm.

We have introduced an effective stress \(\underline{\sigma}^{e}\) defined by: \[ \underline{\sigma}^{e}=\underline{\sigma}-\sum_{i=1}^{N}\underline{X}_{i} \]

The normal is then given by: \[ \underline{n}={\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \underline{\sigma}}}={{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,{{\displaystyle \frac{\displaystyle \underline{s}^{e}}{\displaystyle \sigma^{e}_{\mathrm{eq}}}}} \] where \(\underline{s}^{e}\) is the deviator of \(\underline{\sigma}^{e}\).

Evolution of the isotropic hardening

The isotropic hardening is defined by: \[ R{\left(p\right)}=R_{\infty} + {\left(R_{0}-R_{\infty}\right)}\,\exp{\left(-b\,p\right)} \]

Evolution of the kinematic hardenings

\[ \underline{X}_{i}={{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,C_{i}\,\underline{a}_{i} \]

The evolution of the kinematic variables \(\underline{a}_{i}\) follows the Armstrong-Frederic rule:

\[ \underline{\dot{a}}_{i}=\underline{\dot{\varepsilon}}^{\mathrm{vis}}-g[i]\,\underline{a}_{i}\,\dot{p}=\dot{p}\,{\left(\underline{n}-g[i]\,\underline{a}_{i}\right)} \]

\(\theta\)-scheme

The previous consitutive equations will be integrated using a \(\theta\)-scheme. The increments of the unknowns satisfy:

\[ \left\{ \begin{aligned} \Delta\,\underline{\varepsilon}^{\mathrm{el}}&=\Delta\,\underline{\varepsilon}^{\mathrm{to}}+\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ \Delta\,p&=\left\langle{{\displaystyle \frac{\displaystyle {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle K}}}\right\rangle^{m}\,\Delta\,t\\ \Delta\,\underline{a}_{i}&=\Delta\,p\,{\left(\underline{n}-g[i]\,{\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}\right)}\\ \end{aligned} \right. \]

where the following notations have been used:

• \({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}={\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{el}}\)
• \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}} = \lambda\,{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}\right)}}\,\underline{I}+2\,\mu\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}\)
• \({\left.p\right|_{t+\theta\,\Delta\,t}} = {\left.p\right|_{t}}+\theta\Delta\,p\)
• \({\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}} = {\left.\underline{a}_{i}\right|_{t}}+\theta\Delta\,\underline{a}\)
• \({\left.\underline{X}_{i}\right|_{t+\theta\,\Delta\,t}} = C_{i}\,{\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}\)
• \({\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}={\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}-\sum_{i=1}^{N}{\left.\underline{X}_{i}\right|_{t+\theta\,\Delta\,t}}\)
• \({\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}={\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}_{\mathrm{eq}}\)
• \({\left.F\right|_{t+\theta\,\Delta\,t}}={\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}-R{\left({\left.p\right|_{t+\theta\,\Delta\,t}}\right)}\)

Computation of the jacobian

Terms related to \(f_{\underline{\varepsilon}^{\mathrm{el}}}\)

Derivative of the normal

In the following, we will make use of the “classical” relationship giving the derivative of the normal: \[ {\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}}={{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)} \]

Here, \(\underline{\underline{\mathbf{M}}}\) is a tensor space defined by: \[ \sigma_{\mathrm{eq}}=\sqrt{\underline{\sigma}\,\colon\,\underline{\underline{\mathbf{M}}}\,\colon\,\underline{\sigma}} \]

\(\underline{\underline{\mathbf{M}}}\) is available in MFront as the result of a static member of the st2tost2 template class. In MFront, the alias Stensor4 is usually used to refer to the st2tost2 class for the current numeric type and space dimension.

Here are the expressions of the term related to \(f_{\underline{\varepsilon}^{\mathrm{el}}}\):

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&=\underline{\underline{\mathbf{I}}}+\Delta\,p{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{\underline{\mathbf{I}}}+\Delta\,p\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}}_{{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}}_{\lambda\,\underline{I}\otimes\underline{I}+2\,\mu\,\underline{\underline{\mathbf{I}}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}}_{\theta\,\underline{\underline{\mathbf{I}}}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}}&={\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{a}_{i}}}&=\Delta\,p\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}}}_{{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}}}}_{-{{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,C_{i}\,\underline{\underline{\mathbf{I}}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,a_{i}}}}_{\theta\,\underline{\underline{\mathbf{I}}}}\\ \end{aligned} \right. \]

Finally,

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &=\underline{\underline{\mathbf{I}}}+{{\displaystyle \frac{\displaystyle 2\,\,\mu\,\theta\,\Delta\,p}{\displaystyle {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,p}} &={\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{a}_{i}}} &=-{{\displaystyle \frac{\displaystyle 2\,C_{i}\,\theta\,\Delta\,p}{\displaystyle 3\,{\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)}\\ \end{aligned} \right. \]

Terms related to \(f_{p}\)

Derivative of the equivalent stress

In the following, we will make use of another “classical” relationship giving the derivative of the equivalent stress: \[ {\displaystyle \frac{\displaystyle \partial {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}}={\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \]

To compute the terms of the jacobian associated with \(f_{p}\), we need the derivatives of \({\left.F\right|_{t+\theta\,\Delta\,t}}\) with respect to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\), \(\Delta\,p\) and \(\Delta\,\underline{a}_{i}\). Assuming that \({\left.F\right|_{t+\theta\,\Delta\,t}}\) is positive, we have:

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &=\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}}_{{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}}}_{\lambda\,\underline{I}\otimes\underline{I}+2\,\mu\,\underline{\underline{\mathbf{I}}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}}_{\theta\,\underline{\underline{\mathbf{I}}}}=2\,\mu\,\theta\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,p}} &=-\theta\,{\displaystyle \frac{\displaystyle \partial R}{\displaystyle \partial p}} \\ {\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{a}_{i}}} &=-=\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}}}_{{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}^{e}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}}}}_{-{{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,C_{i}\,\underline{\underline{\mathbf{I}}}}\,\cdot\,\underbrace{{\displaystyle \frac{\displaystyle \partial {\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{a}_{i}}}}_{\theta\,\underline{\underline{\mathbf{I}}}}=-{{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,C_{i}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \right. \]

The blocks of the jacobian associated with \(f_{p}\) are then: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&={\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}}\,\cdot\,{\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,p}}&=1+{\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}}\,\cdot\,{\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,p}}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta\,\underline{a}_{i}}}&={\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}}\,\cdot\,{\displaystyle \frac{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{a}_{i}}}\\ \end{aligned} \right. \]

where \({\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial {\left.F\right|_{t+\theta\,\Delta\,t}}}}=-{{\displaystyle \frac{\displaystyle m\,\Delta\,t}{\displaystyle K}}}\left\langle{{\displaystyle \frac{\displaystyle {\left.F\right|_{t+\theta\,\Delta\,t}}}{\displaystyle K}}}\right\rangle^{m-1}\,\Delta\,t\)

Terms related to \(f_{a_{i}}\)

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{a}_{i}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}& =-{{\displaystyle \frac{\displaystyle 2\,\mu\,\theta\,\Delta\,p}{\displaystyle {\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)} \\ {\displaystyle \frac{\displaystyle \partial f_{\underline{a}_{i}}}{\displaystyle \partial \Delta\,p}}& =-{\left(\underline{n}-g[i]\,{\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}\right)} \\ {\displaystyle \frac{\displaystyle \partial f_{\underline{a}_{i}}}{\displaystyle \partial \Delta\,\underline{a}_{j}}}& =\delta_{ij}\,\underline{\underline{\mathbf{I}}}+{{\displaystyle \frac{\displaystyle 2\,C_{j}\,\theta\,\Delta\,p}{\displaystyle 3\,{\left.\sigma^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}}\,{\left(\underline{\underline{\mathbf{M}}}-{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\otimes{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)} \\ \end{aligned} \right. \]

Implementation

The Implicit domain specific language is used to integrate a behaviour using a \(\theta\)-scheme:

Choice of the DSL

@DSL Implicit;

Behaviour name

We then give the name of the behaviour:

@Behaviour IsotropicViscoplasticityAmstrongFredericKinematicHardening;

The StandardElasticity brick

To implement this behaviour, we will use the StandardElasticity brick which provides:

• Automatic computation of the stress tensor at various stages of the behaviour integration.
• Automatic computation of the consistent tangent operator.
• Automatic support for plane stress and generalized plane stress modelling hypotheses (definitions of the axial strain as an additional state variables and the associated equation enforcing the plane stess condition).
• Automatic addition of the “standard” terms (see below) associated with the elastic strain state variable.

This behaviour brick is fully described here.

The usage of the StandardElasticity is declared as follows:

@Brick StandardElasticity;

Default numerical parameters

A fully implicit integration scheme is choosen as the default:

@Theta 1;

This value can be dynamically changed at runtime by modifying the theta parameter.

The stopping criterion is chosen low, to ensure the quality of the consistent tangent operator (the default value, \(10^{-8}\) is enough to ensure a good estimation of the state variable evolution, but is not enough to have a proper estimation of the consistent tangent operator):

@Epsilon 1e-13;

This value can be dynamically changed at runtime by modifying the epsilon parameter.

Definition of the state variables

We first define the equivalent viscoplastic strain:

@StateVariable strain    p;
p.setGlossaryName("EquivalentViscoplasticStrain");

We then define the kinematic variables as an array of two state variables:

@StateVariable StrainStensor a[2];
a.setEntryName("KinematicVariables");

Definition of the material properties

@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");

@MaterialProperty stress Rinf;
@MaterialProperty stress R0;
@MaterialProperty real b;
@MaterialProperty stress C[2];
@MaterialProperty real   g[2];
@MaterialProperty real m;
@MaterialProperty real K;

Definition of the local variable

@LocalVariable stress mu;

Initialisation of the local variable

/* Initialize Lame coefficients */
@InitializeLocalVariables{
  mu = computeMu(young,nu);
}

Definition of the implicit system

The definition of the implicit system is done in the @Integrator code block:

@Integrator{

To understand this implementation of this code block, one shall bear in mind the following conventions:

• The current estimation of the stress state \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\) at \(t+\theta\,dt\) is computed by the Standardlasticity brick prior the evaluation of the implicit system and is stored in the sig variable.
• For a given integration variable \(X\), the variable X denotes \({\left.X\right|_{t}}\) and dX the current estimation of the increment of \(X\). Thus, the computation of \({\left.X\right|_{t+\theta\,\Delta\,t}}\) is written X+theta*dX. In the following, we declare a new variable called X_ to hold the current estimation of \({\left.X\right|_{t+\theta\,\Delta\,t}}\).
• For two integration variables \(X\) and \(Y\), the variable fX denotes the equation of the implicit system associated with \(X\) and fY the equation associated with \(Y\). The variable dfX_ddY holds the derivative \({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \Delta\,Y}}\).

First, a few constant values are defined:

  const real two_third     = 2/real(3);
  constexpr const real eps = real(1.e-12);

According the previous conventions, we compute the current estimation of \({\left.p\right|_{t+\theta\,\Delta\,t}}\):

  const auto p_ = p +theta*dp ;

The computation of the current effective stress and the current estimation of the kinematic variables is straightforward:

  auto se     = sig ;
  StrainStensor a_[2];
  for(unsigned short i=0;i!=2;++i){
     a_[i]  = a[i]+theta*da[i];
     se  -= C[i]*a_[i]*two_third;
  }

The current estimation of the equivalent stress is then:

  const auto seq = sigmaeq(se);

The current estimation of \({\left.R\right|_{t+\theta\,\Delta\,t}}\) is then evaluated:

  const auto Rp  = Rinf + (R0-Rinf)*exp(-b*p_) ;

\({\left.F\right|_{t+\theta\,\Delta\,t}}\) is then evaluated:

  const auto F   = seq - Rp ;

If the current estimation of \(F\) is positive, we define the implicit system. Otherwise, the default initialization of the implicit equations

  if(F > 0){
    const auto Fexp = pow(F/K,m-1)/K ;
    const auto iseq = 1/(max(seq,eps*young));
    const auto n    = 3*iseq*deviator(se)/2;
    const auto vp   = Fexp*F;
    const auto Jmn = eval(Stensor4::M()-(n^n));
    feel += dp*n;
    fp   -= vp*dt;
    for(unsigned short i=0;i!=2;++i){
      fa[i] -= dp*(n-g[i]*a_[i]);
    }
    // computation of the jacobian matrix
    dfeel_ddeel += 2*mu*theta*dp*Jmn*iseq ;
    dfeel_ddp    = n;
    dfp_ddeel    = - Fexp*m*dt * 2 * mu* theta * n;
    dfp_ddp     += theta* Fexp *m *dt*b*(Rinf-Rp);
    for(unsigned short i=0;i!=2;++i){
      dfeel_dda(i)   = -C[i]*dp*theta*iseq*two_third*Jmn ;
      dfp_dda(i)     =  Fexp*m*dt*C[i]*theta*two_third*n;
      dfa_ddeel(i)   = -2*mu*theta*dp*Jmn *iseq;
      dfa_ddp(i)     = -n + g[i]*a_[i];
      dfa_dda(i,i)   =  (1+dp*g[i]*theta)*Stensor4::Id()+C[i]*dp*theta*iseq*two_third*Jmn;
    }
    dfa_dda(0,1)   =  C[1]*dp*theta*iseq*two_third*Jmn;
    dfa_dda(1,0)   =  C[0]*dp*theta*iseq*two_third*Jmn;
  }
}

An easy optimization

The last equation can be easily eliminated as \(\Delta\,\underline{a}_{i}\) can be expressed as a simple function of \(\Delta\,p\) and \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\):

\[ \left\{ \begin{aligned} \Delta\,\underline{a}_{i}={{\displaystyle \frac{\displaystyle \Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}-g_{i}{\left.\underline{a}_{i}\right|_{t}}}{\displaystyle 1+g_{i}\,\theta\,\Delta\,p}}} \end{aligned} \right. \]

Such an optimization is strongly encouraged as the reduction in the system size is significant. This optimization can be introducted in a straightforward manner:

• The kinematic variable \(\underline{a}_{i}\) must be defined as an auxiliary state variable (using the @AuxiliaryStateVariable keyword) and updated after the integration using the @UpdateAuxiliaryStateVariables code block.
• The updated values of \({\left.\underline{a}_{i}\right|_{t+\theta\,\Delta\,t}}\) must be computed with the previous update formulae.
• The derivatives with respect to the kinematic variables \(\underline{a}_{i}\) must be added to the derivatives with respect to \(\Delta\,p\) and \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) using the chain rule.

Reduction to a scalar equation

Using the isotropy assumption, this system of equations can be further optimised by a reduction to a scalar equation with \(\Delta\,p\) as the only unknown (see [1] for details). However, with this operation, the StandardElasticity brick can’t be used anymore. Thus plane stress an generalised plane stress hypotheses would require a specific treatment and the expression of the consistent tangent operator, which becomes quite complex, has to be written by hand.

References