This page shows how to implement a behaviour describing the load induced thermal strain (LITS) phenomena in concrete following the work of [1].

This implementation is available here.

Description

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{LITS}}\):

\[ \underline{\varepsilon}^{\mathrm{to}}=\underline{\varepsilon}^{\mathrm{el}}+\underline{\varepsilon}^{\mathrm{LITS}} \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)}\]

Load induced thermal strain

The evolution of the load induced thermal strain \(\underline{\varepsilon}^{\mathrm{LITS}}\) is given by:

\[ \underline{\dot{\varepsilon}}^{\mathrm{LITS}}= {{\displaystyle \frac{\displaystyle \eta{\left(\underline{\sigma}\right)}\,\beta{\left(T\right)}}{\displaystyle \sigma_{u0}}}} \left[ -\nu_{\mathrm{LITS}}{\mathrm{tr}{\left({\left<\underline{\sigma}\right>_{-}}\right)}}\,\underline{I} +{\left(1-\nu_{\mathrm{LITS}}\right)}\,{\left<\underline{\sigma}\right>_{-}} \right]\,\dot{T}_{e} \qquad{(3)}\]

where:

• \({\left<\sigma\right>_{-}}\) is the negative part of the stress tensor \(\underline{\sigma}\). This negative part is computed as follows: \[ {\left<\sigma\right>_{-}} = \sum_{I=1}^{3}{\left<\sigma_{I}\right>_{-}}\,\underline{n}_{I} \] where \(\sigma_{I}\) is the Ith eigenvalue of \(\underline{\sigma}\) and \(\underline{n}_{I}\) the associated eigen tensor.
• \(\eta{\left(\underline{\sigma}\right)} = 1+C_{m}{\left(\underline{\sigma}\right)}\,\gamma\). \(\eta\) describes the effect of the stress triaxiality.
• \(C_{m}{\left(\underline{\sigma}\right)}={{\displaystyle \frac{\displaystyle -{\mathrm{tr}{\left({\left<\underline{\sigma}\right>_{-}}\right)}}}{\displaystyle \sqrt{{\left<\underline{\sigma}\right>_{-}}\,\colon\,{\left<\underline{\sigma}\right>_{-}}}}}}\). \(C_{m}\) is the triaxiality index linked to the principal stresses.
• \(\sigma_{u0}\), \(\nu_{\mathrm{LITS}}\) and \(\gamma\) are material parameters. More precisely \(\sigma_{u0}\) is the compressive strength of the material and \(\gamma\) is a material parameter calibrated to experimental data for the appropriate temperature range.
• \(\beta{\left(T\right)}\) is a four order polynomial of the temperature describing uniaxial temperature-LITS behaviour observed in the literature: \[ \beta{\left(T\right)}=\beta_{0}+\beta_{1}\,T+\beta_{2}\,T^{2}+\beta_{3}\,T^{3}+\beta_{4}\,T^{4} \]

\(\dot{T}_{e}\) is the effective temperature rate defined as follows:

\[ \dot{T}_{e}= \left\{ \begin{aligned} \dot{T} &\quad\text{if}\quad T=T_{\mathrm{max}}\\ 0 &\quad\text{if}\quad T<T_{\mathrm{max}}\\ \end{aligned} \right. \]

where \(T_{\mathrm{max}}\) is the maximal temperature seen by the material over time.

The negative part of the stress \({\left<\sigma\right>_{-}}\) is an isotropic function of the stress. TFEL/Math, the mathematical library behind MFront, offers dedicated functions to compute its value and its derivative \({\displaystyle \frac{\displaystyle \partial {\left<\sigma\right>_{-}}}{\displaystyle \partial \underline{\sigma}}}\).

Implicit scheme

Equation (3) can be discretized using a standard implicit scheme as follows:

\[ \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}=C\,\eta{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{s} \qquad{(4)}\]

where:

• \(\theta\) is a numerical parameter (\(\theta \in \left[0:1\right]\).
• \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\) is evaluated by Equation (2): \[ {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}={\left.\lambda\right|_{t+\theta\,\Delta\,t}}\,{\mathrm{tr}{\left({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}\right)}}+ 2\,{\left.\mu\right|_{t+\theta\,\Delta\,t}}\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}+ \]
• \({\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\theta\,\Delta\,t}}={\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t}}+\theta\,\Delta\,\underline{\varepsilon}^{\mathrm{el}}\)
• \({\left.T\right|_{t+\theta\,\Delta\,t}}={\left.T\right|_{t}}+\theta\,\Delta\,T\)
• \(C={{\displaystyle \frac{\displaystyle \beta{\left({\left.T\right|_{t+\theta\,\Delta\,t}}\right)}\,\Delta\,T_{e}}{\displaystyle \sigma_{u0}}}}\)
• \(\underline{s}= -\nu_{\mathrm{LITS}}{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}\,\underline{I}+{\left(1-\nu_{\mathrm{LITS}}\right)}\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\)

In the previous expressions, \(\Delta\,T_{e}\) is the effective temperature increment.

Equation (4) shows that \(\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}\) is known if \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) is known.

Using Equation (1) in its increment form thus leads to the following non linear equation whose only unknows is \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\):

\[ f_{\underline{\varepsilon}^{\mathrm{el}}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\right)}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}=\underline{0} \qquad{(5)}\]

This non linear equation will be solved by a standard Newton-Raphson algorithm, which requires the jacobian \({\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\):

\({\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is computed as follows:

\[ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{I}+{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{I}+ C\,\left[ \underline{s}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}+ \eta{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}\,{\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \right] \]

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &=\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ &=\theta\,\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ &=\theta\,\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \qquad{(6)}\]

where \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\) is the elastic stiffness tensor at the middle of the time step:

\[ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}={\left.\lambda\right|_{t+\theta\,\Delta\,t}}\,\underline{I}\,\otimes\,\underline{I}+2\,\mu\,\underline{\underline{\mathbf{I}}} \]

Note

The explicit expression of \({\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\) is not given here as it will be computed by an build-in function of the TFEL/Math library.

\({\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\) is computed as follows: \[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}} &={\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}{\left({{\displaystyle \frac{\displaystyle -{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\right)}\\ &={{\displaystyle \frac{\displaystyle -\underline{I}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}+{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}\,{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}{\left({{\displaystyle \frac{\displaystyle -1}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\right)}\\ &={{\displaystyle \frac{\displaystyle -\underline{I}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\,\underline{I}+{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}\,{{\displaystyle \frac{\displaystyle {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle {\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}^{3/2}}}}\, \end{aligned} \]

Finally, \({\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is given by:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&= {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ &\theta\,{\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,\underline{\underline{\mathbf{D}}}\\ \end{aligned} \]

where:

\[ {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}= -\nu_{\mathrm{LITS}}\,\underline{I}\,\otimes\,\underline{I}+{\left(1-\nu_{\mathrm{LITS}}\right)}\,\underline{\underline{\mathbf{I}}}\\ \]

Implementation

Choice of the domain specific language

The domain specific language suitable for implementing an implicit scheme is called Implicit. The @DSL keyword allows specifying the domain specific language to be used:

@DSL Implicit;

Choice of the behaviour name

The @Behaviour allows giving a name the behaviour:

@Behaviour LoadInducedThermalStrain_Torelli2018;

Defining the name of behaviour is required.

Metadata

The @Author, @Date and @Description keywords allows adding metadata to the behaviour (respectively the name of the author of the implementation, the date at which the implementation was written, a short description), as follows:

@Author Giacomo Torelli;
@Date 25 / 03 / 2018;
@Description {
  "A confinement-dependent load-induced thermal strain "
  "constitutive model for concrete subjected "
  "to temperatures up to 500 °C"

  "Torelli, Giacomo and Mandal, Parthasarathi and Gillie, "
  "Martin and Tran, Van-Xuan"

  "International Journal of Mechanical Sciences"

  "http://www.sciencedirect.com/science/article/pii/S0020740317337372"
}

Those metadata, though recommended, are optional.

Use of the StandardElasticity brick

The behaviour treated in this document allows using the StandardElasticity brick. This brick, described in this page, 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 (The axial strain is defined as an additional state variable and the associated equation in the implicit system is added to enforce the plane stess condition).
• Automatic addition of the standard terms 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;

As the StandardElasticity brick provides support for all modelling hypotheses, we can declare that your behaviour indeed supports all modelling hypotheses:

@ModellingHypotheses{".+"};

Choice of the numerical parameter

The implicit scheme requires three numerical parameters:

• the \(\theta\) parameter (default value is \(1/2\)).
• the maximum number of iterations (default value is \(100\)).
• the stopping criteria (default value is \(10^{-8}\)). This may be sufficient in most cases, but not for the computation of the consistent tangent operator by the StandardElasticity brick.

@Theta 0.5;
@IterMax 100;
@Epsilon 1.e-16;

Note

Those values are automatically generates parameters that can be changed at runtime: the names of those parameters are respectively theta, iterMax and epsilon.

Thermoelastic properties, computation of the thermal strain

The @ComputeStiffnessTensor keyword allows computing the elastic stiffness tensor in a variable called D. In this case, the UnAltered version of the stiffness tensor is required (“unaltered” means that the modelling hypothesis was not taken into account).

@ComputeStiffnessTensor<UnAltered>{47000, 0.25};

The @ComputeThermalExpansion keyword allows computing the thermal strain. It is followed by the mean thermal expansion coefficient.

@ComputeThermalExpansion 10e-6;

The thermal strain is automatically removed to the total strain and the total strain increment.

Parameters

The following code declares the various parameters used by the constitutive equations:

@Parameter gamma = 1.5;
@Parameter sigmultimate = 50.;
@Parameter nulits = 0.50;
@Parameter tcrit = 0.;
@Parameter b[5] = {2.7031065533E-05, -1.0209170592E-06, 6.1200423753E-9, //
                   -1.2632648735E-11, 6.9158539621E-15};

Auxiliary state variables

By default, the Implicit DSL declares the elastic strain as a state variable, i.e. a variable which is part of the implicit system and whose value is saved from one time step to the other.

We choose to declare two auxiliary state variables, the load induced thermal strain \(\elits\) and the maximum value of the temperature seen by the material \(T_{\mathrm{max}}\). Those variables are not part of the implicit system but their values are saved from one time step to the other.

Saving maximum value of the temperature seen by the material is required, but its evolution is independent of the mechanics.

Saving the load induced thermal strain is only useful for post-processing. It could be removed for saving some space in memory.

Auxiliary state variables are introduced by the @AuxiliaryStateVariable keyword, as follows:

@AuxiliaryStateVariable StrainStensor elits;
elits.setEntryName("LoadInducedThermaStrain");
@AuxiliaryStateVariable temperature Tmax;
Tmax.setEntryName("MaximalValueOfTheTemperature");

Local Variables

Local variables are available in every code blocks. Here we will introduce two local variables:

• the increment of the load induced thermal strain, that will be used to update the load induced thermal strain at the end of the time step
• the creep coefficient which allows computing the thermal dependency before the implicit resolution.

//! LITS increment
@LocalVariable StrainStensor delits;
//! Creep coefficient
@LocalVariable real C;

Initialization of the local variables

The @InitLocalVariables code block is called before the implicit resolution. It is used here to compute the C variable and update the Tmax auxiliary variable.

@InitLocalVariables {
  Tmax = max(max(tcrit, Tmax), T);
  const auto T_ = T + theta * dT;
  const auto beta = b[0] + T_ * (b[1] + T_ * (b[2] + T_ * (b[3] + T_ * b[4])));
  const auto dTe = max(T + dT - max(tcrit, Tmax), temperature(0));
  C = (beta / (-sigmultimate)) * dTe;
}

Implicit system

The implicit system is defined in the code block introduced by the @Integrator code block.

This code block is called at each iteration of the implicit system, after the stress \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\) at the middle of the time step has been computed and stored in the sig variable (this is handled by the StandardElasticity brick).

@Integrator {

We first define some useful constants.

  constexpr const stress eeps = 1.e-12;
  const stress seps = eeps * young;
  constexpr const auto id = Stensor::Id();
  constexpr const auto esolver = StressStensor::FSESJACOBIEIGENSOLVER;

Then the negative part of the stress \({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\) and its derivative with respect to the stress \({\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\) are computed using the computeIsotropicFunctionAndDerivative method of the stensor class.

Notes

A scalar function \(f\) defines an isotropic function of a symmetric tensor \(\underline{a}\) as follows:

\[ f{\underline{a}} = \sum_{I=1}^{3}{\left<a_{I}\right>_{-}}\,\underline{n}_{I} \]

where \(a_{I}\) are the eigenvalues of \(\underline{a}\) and \(\underline{n}_{I}\) its eigentensors.

This function requires two scalar functions defining the isotropic function and its derivative. Those functions are called np and dnp here:

  // negative part
  const auto np = [](const stress x) -> stress { return x <= 0 ? x : 0; };
  const auto dnp = [&seps](const stress x) -> real {
    return std::abs(x) < seps ? 0.5 : ((x < 0) ? 1 : 0);
  };

Then we call the computeIsotropicFunctionAndDerivative method:

  auto sn = StressStensor{};
  auto dsn = Stensor4{};
  std::tie(sn, dsn) =
      sig.template computeIsotropicFunctionAndDerivative<esolver>(np, dnp,
                                                                  seps * 0.1);

The previous expression can be simplified in C++-17:

  const auto [sn,dsn] = 
      sig.template computeIsotropicFunctionAndDerivative<esolver>(np, dnp,
                                                                  seps * 0.1);

About the choice of the eigensolver

Here we used a the Jacobi algorithm as implemented in [2]. This algorithm is lower but more accurate than the default one (see here for an indepth comparison of available algorithms.

The next step is to compute \({\left.\eta\right|_{t+\theta\,\Delta\,t}}\) and its derivative \({\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.sig\right|_{t+\theta\,\Delta\,t}}}}\) (See Equation (6)):

  //----- Multiaxial correction coefficient
  const auto sn_eq = sqrt(sn | sn);
  const auto isn_eq = 1 / max(seps, sn_eq);
  const auto cm = -trace(sn) * isn_eq;
  const auto dcm_dsig = eval((-isn_eq * id + trace(sn) * power<3>(isn_eq) * sn) | dsn);
  const auto eta = 1 + (cm - 1) * gamma;
  const auto deta_dsig = gamma * dcm_dsig;

We are now able to compute the increment \(\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}\) and its derivative:

  const auto se = eval((1 + nulits) * sn - nulits * trace(sn) * id);
  delits = C * eta * se;
  const auto dse_dsig = (1 + nulits) * dsn - nulits * ((id ^ id) * dsn);
  const auto ddelits_dsig = C * (se ^ deta_dsig) + C * eta * dse_dsig;

Finally, we form the implicit equation associated with the elastic strain and its derivative:

  feel += delits;
  dfeel_ddeel += theta * ddelits_dsig * D;
} // end of @Integrator

Note that we used the fact that feel is implicitly initialized to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\).

Update of the auxiliary state variables

The @UpdateAuxiliaryStateVariables code block is called after the implicit resolution, after the computation of the stress at the end of the time step and after the update of the state variables (here the elastic strain).

We use the last estimate of the its increment to update the elits variable:

@UpdateAuxiliaryStateVariables {
  e_lits += delits;
}

References