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 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)}}=s^{e}_{\mathrm{eq}}-R{{\left(p\right)}} \]

where:

We have introduced an effective deviatoric stress \({\underline{s}}^{e}\) defined by: \[ {\underline{s}}^{e}={\underline{s}}-\sum_{i=1}^{N}{\underline{X}}_{i} \] where \({\underline{s}}\) is the deviatoric part of the stress.

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 s^{e}_{\mathrm{eq}}}}}} \]

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:

Implicit system formulation and implementation

The previous equations can be translated in an implicit system in a straightforward manner:

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

Computation of the jacobian

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{s}}\right|_{t+\theta\,\Delta\,t}}}}}}={{{\displaystyle \frac{\displaystyle 1}{\displaystyle {{\left.s^{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.s^{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{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{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{s}}\right|_{t+\theta\,\Delta\,t}}}}}}}_{{{{\displaystyle \frac{\displaystyle 1}{\displaystyle {{\left.s^{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{s}}\right|_{t+\theta\,\Delta\,t}}}}{\displaystyle \partial {{\left.{\underline{a}}_{i}\right|_{t+\theta\,\Delta\,t}}}}}}}_{-{{{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}}\,C_{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.s^{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.s^{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. \]

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.s^{e}_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}}}{\displaystyle \partial {{\left.{\underline{s}}\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{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{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 {\underline{s}}}}}}_{{{\left.{\underline{n}}\right|_{t+\theta\,\Delta\,t}}}}\,\cdot\,\underbrace{{{\displaystyle \frac{\displaystyle \partial {{\left.{\underline{s}}\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{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\)

\[ \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.s^{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.s^{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:

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 `(???) code block:

@Integrator{

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

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| Stensor4::M() )*two_third;
    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}}}}{\displaystyle 1+g_{i}\,{{\left.{\underline{a}}_{i}\right|_{t}}}\,\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:

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 Chaboche and Cailletaud (1996) 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

Chaboche, J.L., and G. Cailletaud. 1996. “Integration Methods for Complex Plastic Constitutive Equations.” Computer Method in Applied Mechanics and Engineering 133: 125–55.