This article shows how to implement an orthotropic plastic behaviour with isotropic linear hardening in MFront.

The implementation is available here.

Description of the behaviour

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

\[ \underline{\varepsilon}^{\mathrm{to}}=\underline{\varepsilon}^{\mathrm{el}}+\underline{\varepsilon}^{\mathrm{p}} \qquad{(1)}\]

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}} \qquad{(2)}\]

Plastic flow

The plastic part of the behaviour is described by the following yield surface: \[ f{\left(\sigma_{\mathrm{eq}},p\right)} = \sigma_{\mathrm{eq}}-\sigma_{0}-R\,p \qquad{(3)}\]

where \(\sigma_{\mathrm{eq}}\) is the vonMises stress defined below, \(p\) is the equivalent plastic strain, \(\sigma_{0}\) is the yield stress and \(R\) is the hardening slope.

The von Mises stress is defined by: \[ \sigma_{\mathrm{eq}}=\sqrt{{{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,\underline{s}\,\colon\,\underline{s}}=\sqrt{3\,J_{2}} \] where:

• \(\underline{s}\) is the deviatoric stress defined as follows: \[ \underline{s}=\underline{\sigma}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,{\mathrm{tr}{\left(\underline{\sigma}\right)}}\,\underline{I} \]
• \(J_{2}\) is the second invariant of \(\underline{s}\).

The previous expression can be rewritten by introducing a fourth order tensor called \(\underline{\underline{\mathbf{M}}}\): \[ \sigma_{\mathrm{eq}}=\sqrt{\sigma\,\colon\,\underline{\underline{\mathbf{M}}}\,\colon\,\underline{\sigma}} \]

The tensor \(\underline{\underline{\mathbf{M}}}\) is given by: \[ \underline{\underline{\mathbf{M}}}={{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,{\left(\underline{\underline{\mathbf{I}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\underline{I}\,\otimes\,\underline{I}\right)} \]

The plastic flow is assumed to be associated, so the flow direction \(\underline{n}\) is given by \({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{\sigma}}}\):

\[ \underline{n} = {\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{\sigma}}} = {{\displaystyle \frac{\displaystyle 3}{\displaystyle 2\,\sigma_{\mathrm{eq}}}}}\,\underline{s} \]

One may readily check that:

\[ \underline{n}\,\colon\underline{n}={{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}} \qquad{(4)}\]

Integration algorithm

The previous constitutive equations will be integrated using a standard implicit scheme which, due to the isotropy of the elastic behaviour, the quadratic nature of the equivalent stress and the linearity of the isotropic hardening rule can be solved analytically.

Expression of the plastic strain increment

Assuming that a plastic flow occurs, the plastic increment \(\Delta\,\underline{\varepsilon}^{\mathrm{p}}\) over the time step is approximated by:

\[ \Delta\,\underline{\varepsilon}^{\mathrm{p}}= \Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \]

The slip of strain given by Equation (1) may be written incrementally:

\[ \Delta\,\underline{\varepsilon}^{\mathrm{to}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{p}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \qquad{(5)}\]

and projected on the deviatoric space as follows:

\[ \Delta\,\underline{s}^{\mathrm{to}} = \Delta\,\underline{s}^{\mathrm{el}}+\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \]

The isotropy of the elastic behaviour projected on the deviatoric space leads to:

\[ \begin{aligned} {\left.\underline{s}\right|_{t+\theta\,\Delta\,t}} &= {{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,{\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ &= 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{s}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{s}^{\mathrm{el}}\right)} \\ &= 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{s}^{\mathrm{el}}\right|_{t}} +\theta\,\Delta\,\underline{s}^{\mathrm{to}} -\theta\,\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\right)} \end{aligned} \]

The previous equations leads to:

\[ {\left({{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,{\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}+\theta\,\Delta\,p\right)}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} = 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{s}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{s}^{\mathrm{to}}\right)} \qquad{(6)}\]

The right hand side is the value of the deviatoric part of the stress tensor that would be obtained in a pure elastic loading. The latter is commonly called the trial stress \(\underline{\sigma}^{\mathrm{tr}}\):

\[ \underline{\sigma}^{\mathrm{tr}}= {\left.\lambda\right|_{t+\theta\,\Delta\,t}}\,{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+ 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)} \]

Its deviatoric part is denoted \({\left.\underline{s}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}\):

\[ {\left.\underline{s}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}=2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left({\left.\underline{s}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{s}^{\mathrm{to}}\right)} \]

Taking the von Mises norm of both side of Equation (6) yields:

\[ {\left({{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,{\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}+\theta\,\Delta\,p\right)}\,\sqrt{{{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\colon{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}} = {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}} \]

where \({\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}\) is the so-called trial equivalent stress.

It is important to note that \(\underline{s}^{\mathrm{tr}}\) and \({\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}\) can be computed a priori, only with the known values of the elastic strain at the beginning of the time step and the total strain increment.

The previous expressions leads to a closed form expression of the normal:

\[ {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} = {{\displaystyle \frac{\displaystyle 3}{\displaystyle 2\,{\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}}}}\underline{s}^{\mathrm{tr}} \qquad{(7)}\]

Derivatives of the equivalent trial stress \({\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}\) and normal \({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\)

The following classical relationships will be useful in the following:

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}}}&={\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ {\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.\underline{\sigma}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}}}}&= {{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\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. \qquad{(8)}\]

Finally, taking into account Equation (4), the von Mises stress satisfies:

\[ {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}} = {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}} - 3\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,\theta\,\Delta\,p \]

Using the fact that, during plastic loading, the stress state must be on the Yield Surface (3), the von Mises stress also satisfies:

\[ {\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}=\sigma_{0}+R\,{\left.p\right|_{t+\theta\,\Delta\,t}}=\sigma_{0}+R\,{\left({\left.p\right|_{t}}+\theta\,\Delta\,p\right)} \]

Combining those equations yields this closed-form expression of the plastic strain increment:

\[ \Delta\,p = {{\displaystyle \frac{\displaystyle {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}} - \sigma_{0} - R\,{\left.p\right|_{t}}}{\displaystyle {\left(3\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\, + R\right)}\,\theta}}} \qquad{(9)}\]

Expression of the elastic strain increment

The elastic strain increment is given in closed form by Equations (5), (7) and (9):

\[ \Delta\underline{\varepsilon}^{\mathrm{el}}= \Delta\,\underline{\varepsilon}^{\mathrm{to}}-\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \qquad{(10)}\]

Expression of the consistent tangent operator

The consistent tangent operator is the derivative \({\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) of the stress at the end of the time step with respect the increment of the strain.

Using Equation (2), this derivative can be expressed as:

\[ {\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

where \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}={\left.\lambda\right|_{t+\Delta\,t}}\,\,\underline{\underline{\mathbf{I}}}+2\,{\left.\mu\right|_{t+\Delta\,t}}\,\underline{I}\,\otimes\,\underline{I}\) denotes the stiffness matrix.

The derivative \({\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) can be deduced from Equation (10):

\[ {\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= \underline{\underline{\mathbf{I}}} -{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} -\Delta\,p\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

Using Equation (8), the derivative \({\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) can be deduced from Equation (9) as follows:

\[ {\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left(3\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\, + R\right)}\,\theta}}}\, {\displaystyle \frac{\displaystyle \partial {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\right|_{t+\theta\,\Delta\,t}} }{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} = {{\displaystyle \frac{\displaystyle 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}}{\displaystyle 3\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\, + R}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \qquad{(11)}\]

Using Equations (8) and (7), the derivative \({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) is also given by:

\[ {\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {{\displaystyle \frac{\displaystyle 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,\theta}{\displaystyle {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\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)} \qquad{(12)}\]

Equations (11) and (12), shows that the derivatives \({\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) and \({\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) are linear combinations of the two fourth order tensors \(\underline{M}\) and \({\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\) which only acts on the deviatoric parts of the symmetric tensors. The following properites then holds:

\[ \left\{ \begin{aligned} {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,\underline{\underline{\mathbf{M}}} &= 2\,{\left.\mu\right|_{t+\Delta\,t}}\,\underline{\underline{\mathbf{M}}}\\ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} &= 2\,{\left.\mu\right|_{t+\Delta\,t}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \right. \]

The consistent tangent operator is then:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} &= {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}} -{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\left( {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} +\Delta\,p\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \right)}\\ &={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}} -2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left( {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \Delta\,p}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} +\Delta\,p\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \right)} \end{aligned} \]

and finaly,

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} ={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}} -4\,{\left({\left.\mu\right|_{t+\theta\,\Delta\,t}}\right)}^{2}\, &\left( {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\, + R}}}\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}} \otimes\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}+\right.\\ &\left.{{\displaystyle \frac{\displaystyle \theta\Delta\,p}{\displaystyle {\left.\sigma_{\mathrm{eq}}^{\mathrm{tr}}\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)}\right)\\ \end{aligned} \qquad{(13)}\]

Fully implicit implementation

We now describe how to make a fully implicit implementation, i.e. we fix the \(\theta\) parameter to \(1\).

This choice is consistent with the description of a plastic flow as it imposes that the material lies on the yield surface at the end of the time step if a plastic loading occurs.

With other choices, the material may be outside the yield surface at the end of the time step and spurious oscillations in the stress in time can be observed, for example during a simple uniaxial tensile test test.

Choice of the domain specific language

The behaviour integration is trivial as all the computated quantities can be obtained thanks in closed form:

• The values of the state variables \(\underline{\varepsilon}^{\mathrm{el}}\), \(p\) and the value of the stress \(\underline{\sigma}\) at the end of the time step.
• The consistent tangent operator \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\).

In this case, an appropriate domain specific language is the Default:

@DSL Default;

Contrary to other DSLs, the Default does not provide a built-in support for a given integration scheme.

Name of the behaviour

The name of the behaviour is introduced by the @Behaviour keyword:

@Behaviour IsotropicLinearHardeningPlasticity;

Metadata

The following lines defines the author of the implementation, the date and provides a small description:

@Author Thomas Helfer;
@Date   14/10/2016;
@Description{
  An implicit implementation of a simple
  isotropic plasticity behaviour with
  isotropic linear hardening.

  The yield surface is defined by:
  "\["
  "  f(\sigmaeq,p) = \sigmaeq-s_{0}-H\,p"
  "\]"
}

Declaration of state parameters

The behaviour used two state variables, the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) and the equivalent plastic strain \(p\).

Those state variables are declared as follows:

@StateVariable StrainStensor eel;
eel.setGlossaryName("ElasticStrain");
@StateVariable strain p;
p.setGlossaryName("EquivalentPlasticStrain");

The state variables have been associated to glossary names which are meant to:

• automatically declare the external name of the variable (i.e. the name as seen from the calling solver if this solver uses the metadata exported by MFront).
• automatically document the variable.

Increments of the state variables

The Default automatically declares the increments of those variables, called respectively deel and dp. The user is free to use those increments or to directly update the values of the state variables.

If the increment of a state variable is defined, it is added to the state variable after the @Integrator code block.

In our implementation, we will directly update the elastic strain and use the plastic strain increment.

Declaration of parameters

In this implementation, we choose to use parameters associated respectively with the Young’ modulus, Poisson’ ratio, hardening slope and yield stress, as follows:

@Parameter stress young =  70.e9;
young.setGlossaryName("YoungModulus");
@Parameter real nu =   0.34;
nu.setGlossaryName("PoissonRatio");
@Parameter stress H =  10.e9;
H.setEntryName("HardeningSlope");
@Parameter stress s0 =  300.e6;
s0.setGlossaryName("YieldStress");

Using material properties

The use of parameters could be replaced by the use of material properties as follows:

@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");
@MaterialProperty stress H;
H.setEntryName("HardeningSlope");
@MaterialProperty stress s0;
s0.setGlossaryName("YieldStress");

Contrary to parameters, material properties does not have default values and must be provided by the calling solvers, which generally uses values provided by the user in its input file.

Computation of the prediction operator

The current implementation only provides the elastic operator:

@PredictionOperator{
  const auto lambda = computeLambda(young, nu);
  const auto mu     = computeMu(young, nu);
  Dt = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
}

We could also provide a tangent operator, but this would mean saving an auxiliary state variable stating if a plastic loading occured at the previous time step.

Behaviour integration

The behaviour integration is implemented in the @Integrator code block.

Computation of the consistent tangent operator

The Default DSL allows to compute the consistent tangent operator either inside the @Integrator code block or in a dedicated code block name @TangentOperator.

The first choice is simplier and more efficient for the considered implementation, but this requires to explicitely state that the consistent tangent operator is defined in the @Integrator code block using the @ProvidesSymmetricTangentOperator keyword, as follows:

@ProvidesSymmetricTangentOperator;

The consistent tangent operator may or may not be requested by the calling solver at runtime.

Computation of the Lamé’ coefficients

The @Integrator code blocks begins by computing the Lamé’ coefficients as follows:

@Integrator{
  const auto lambda = computeLambda(young, nu);
  const auto mu     = computeMu(young, nu);

Elastic loading

The next lines assume a purely elastic loading and:

• Updates the value of the elastic strain with the increment of the total strain.
• Initializes the value of the consistent tangent operator to the elastic stiffness matrix.

  eel += deto;
  if(computeTangentOperator_){
    Dt = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
  }

If required, a plastic correction will be applied to the elastic strain and the consistent tangent operator in the next part of the code.

Computation of the trial stress

To check of if a plastic correction must be applied, we first compute the deviator of the trial stress, the equivalent trial stress and checks if its value is outside the elastic domain at the beginning of the time step.

  const auto se     = 2 * mu * deviator(eel);
  const auto seq_e  = sigmaeq(se);
  const auto plastic_loading = seq_e > s0 + H * p;

The plastic_loading boolean value is thus true if a plastic correction must be applied, false otherwise.

Plastic loading

Computation of the normal

If a plastic correction must be applied, we first compute the normal using Equation (7):

  if(plastic_loading){
    const auto iseq_e = 1 / seq_e;
    const auto n      = eval(3 * se / (2 * seq_e));

The eval function

TFEL/Math uses expression templates to perform lazy evaluations of most tensorial operations. See its documentation for details.

The eval function evaluates an expression to get its results. We could have removed it here, but the normal is used multiple times, so we choose to evaluate it once.

Update of the state variables

The increment of the equivalent plastic strain is computed using Equation (9):

    const auto cste   = 1 / (H + 3 * mu);
    dp   = (seq_e - s0 - H * p) * cste;

The value of the equivalent plastic strain will be updated automatically after the @Integrator code block.

Once the equivalent plastic strain is known, the elastic strain can be corrected using the split of the strain (Equation (1)):

    eel -= dp * n;

Update of the consistent tangent operator

The consistent tangent operator can be corrected using Equation (13):

    if((computeTangentOperator_) && (smt==CONSISTENTTANGENTOPERATOR)){
      Dt -= -4 * mu * mu * (dp * iseq_e * (Stensor4::M()-(n ^ n)) + cste * (n ^ n));
    }
  } // end of if(plastic_loading)

Computation of the stress at the end of the time step

At this stage, the elastic strain have been updated and hols its value at the end of the time step. The last step of the implementation is to use it to compute the stress at the end of the time step:

  sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
} // end of @Integrator