1 Description

Let \(\vec{Y}\) be a vector holding all the integration variables¹. The evolution of \(\vec{Y}\) is assumed to satisfy the following differential equation:

\[ \vec{\dot{Y}}=G{\left(\vec{Y},t\right)} \]

where \(t\) is not meant to describe an explicit dependency to the time but rather a placeholder for variables whose evolutions are known (those variables are called external state variables in MFront).

The following notations are used:

• \({\left.\vec{Y}\right|_{t}}\) denotes the value of the integration variables at the beginning of the time step.
• \(\Delta\,\vec{Y}\) denotes the increment of the integration variables.
• \({\left.\vec{Y}\right|_{t+\Delta\,t}}\) denotes the value of the integration variables at the end of the time step: \[ {\left.\vec{Y}\right|_{t+\Delta\,t}}={\left.\vec{Y}\right|_{t}}+\Delta\,\vec{Y} \]

An implicit scheme turns this ordinary differential equation into a non-linear system of equations whose unknown is the increment \(\Delta\,Y\) over a time step \(\Delta\,t\):

\[ \Delta\,\vec{Y}-G{\left({\left.\vec{Y}\right|_{t}}+\theta\,\Delta\,\vec{Y},t+\theta\,\Delta\,t\right)}\,\Delta\,t=0 \]

The increment \(\Delta\,\vec{Y}\) of the state variables satisfies the following implicit system:

\[ F{\left(\Delta\,\vec{Y},\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}=0 \qquad{(1)}\]

Equation (1) implicitly defines \(\Delta\,\vec{Y}\) as an implicit function of the increment of the strain tensor \(\Delta\underline{\varepsilon}^{\mathrm{to}}\) and may be rewritten as:

\[ \vec{F}{\left(\Delta\,\vec{Y}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)},\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}=0 \qquad{(2)}\]

2 Available algorithms

The following algorithms are available:

• NewtonRaphson
• NewtonRaphson_NumericalJacobian
• PowellDogLeg_NewtonRaphson
• PowellDogLeg_NewtonRaphson_NumericalJacobian
• Broyden
• PowellDogLeg_Broyden
• Broyden2
• LevenbergMarquardt
• LevenbergMarquardt_NumericalJacobian

Those algorithms are described in the documentation of the TFEL/Math library. The main steps of those algorithms and the associated code blocks are depicted in Figure 1.

2.1 Notes about updating auxiliary state variable or local variables in the Integrator code blocks when the numerical evaluation of the jacobian is requested

In some cases, it is convenient to update auxiliary state variables in the @Integrator code block, rather than computing them in the @UpdateAuxiliaryStateVariables code block which is only called once the convergence is reached.

However, if the jacobian matrix is computed numerically (at least partially), such updates could be wrong, because they can be based on the perturbed values of the unknowns.

Since TFEL 3.1, this can be circumvented by testing the value of the perturbatedSystemEvaluation boolean value as follows:

// let av be an auxiliary state variable
@AuxiliaryStateVariable StrainStensor av;

@Integrator{
  // put updated value of av in a temporary variable
  const auto av_ = eval(...);  
  ...
  definition of the implicit system
  ...
  if(!perturbatedSystemEvaluation){
    // update auxiliary state variables
    av = av_;
  }
} // end of @Integrator

In many cases, rather than updating auxiliary variables during the Newton iterations, it can be more practical to compute its increment, defined by a local variable and to update the auxiliary variable in the @UpdateAuxiliaryStateVariables code block. The previous trick can be used in this case in a straightforward manner.

3 Computation of the consistent tangent operator

3.1 Computation of the consistent tangent operator for small strain behaviours

In this section, a small strain behaviour is considered.

The increment \(\Delta\,\vec{Y}\) of the state variables satisfies the following implicit system:

\[ F{\left(\Delta\,\vec{Y},\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}=0 \qquad{(3)}\]

Equation (3) implicitly defines \(\Delta\,\vec{Y}\) as an implicit function of the increment of the increment of the strain tensor \(\Delta\underline{\varepsilon}^{\mathrm{to}}\) and may be rewritten as:

\[ \vec{F}{\left(\Delta\,\vec{Y}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)},\Delta\,\underline{\varepsilon}^{\mathrm{to}}\right)}=0 \qquad{(4)}\]

The stress tensor \(\underline{\sigma}\) is assumed to be an explicit function of the integration variables \({\left.\vec{Y}\right|_{t+\Delta\,t}}\) at the end of the time step and the total strain \({\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}}\) at the end of time step:

\[ \underline{\sigma}{\left({\left.\vec{Y}\right|_{t+\Delta\,t}}, {\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}}\right)} \]

3.1.1 Formal derivation of the consistent tangent operator

The consistent tangent operator is thus given by:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\underline{\sigma}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+ {\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\vec{Y}}}\,{\displaystyle \frac{\displaystyle \partial \Delta\,\vec{Y}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

By differentiation:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\vec{F}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\vec{Y}}}\,{\displaystyle \frac{\displaystyle \mathrm{d}\Delta\,\vec{Y}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}+ {\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=0 \]

\({\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\vec{Y}}}\) is the jacobian \(J\) of the implicit system.

Hence,

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\Delta\,\vec{Y}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}} =-J^{-1}\,{\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

Finally, the consistent tangent operator is formally given by:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\underline{\sigma}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}= {\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}- {\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\vec{Y}}}\,J^{-1}\,{\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \qquad{(5)}\]

3.1.2 Discussion

Equation (5) may be extented almost directly to generalised behaviours. However, direct usage of Equation (5) has several practical disadvantages:

• The matrix \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\vec{Y}}}\) is usually sparse. For example, most small strain behaviours are based on the Hooke law, which means that the stress \(\underline{\sigma}\) only depends on the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\). In many behaviours, the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) is only a small subset of the integration variables \(\vec{Y}\).
• It requires to invert the jacobian matrix.
• The matrix \({\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) is also usually sparse. For example, the \(\Delta\,\underline{\varepsilon}^{\mathrm{to}}\) commonly only appears in the implicit equation describing the split of the strain.

Hence, a computation of the consistent tangent based on Equation (5) may imply a significant performance hit which can be avoided in most common cases, as discussed in the following section.

3.1.3 Simplication of the Equation (5) in common cases

Equation (5) can be simplified if the following assumptions are made:

• The integration variable \(\vec{Y}\) is decomposed as: \[ \vec{Y}= \begin{pmatrix} \Delta\,\underline{\varepsilon}^{\mathrm{el}}\\ \Delta\,\vec{z} \end{pmatrix} \] where \(\vec{z}\) is a sub-vector of \(\vec{Y}\) containing the other state variables. The implicit system is also decomposed in a similar manner: \[ \vec{F}= \begin{pmatrix} f_{\underline{\varepsilon}^{\mathrm{el}}}\\ f_{\vec{z}} \end{pmatrix} \]
• The implicit equation \(f_{\underline{\varepsilon}^{\mathrm{el}}}\) associated with the elastic strain holds the split of the strain, which means that is has the form: \[ f_{\underline{\varepsilon}^{\mathrm{el}}}=\Delta\underline{\varepsilon}^{\mathrm{el}}-\Delta\underline{\varepsilon}^{\mathrm{to}}+... \]
• The strain increment \(\Delta\underline{\varepsilon}^{\mathrm{to}}\) only appears in \(f_{\underline{\varepsilon}^{\mathrm{el}}}\). Thus, \({\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) has the form: \[ {\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}=- \begin{pmatrix} 1& 0 & 0 & 0& 0 &0 \\ 0& 1& 0 & 0& 0 &0 \\ 0& 0 & 1 & 0& 0 &0 \\ 0& 0 & 0 & 1& 0 &0 \\ 0& 0 & 0 & 0& 1 &0 \\ 0& 0 & 0 & 0& 0 &1\\ 0& 0 & 0 & 0& 0 &0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0& 0 & 0 & 0& 0 &0 \\ \end{pmatrix} \qquad{(6)}\]
• The stress tensor only depends on the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) through the Hooke law: \[ {\left.\underline{\sigma}\right|_{t+\Delta\,t}}={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}} \qquad{(7)}\] The expression of the consistent tangent operator is then: \[ {\displaystyle \frac{\displaystyle \mathrm{d}\underline{\sigma}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\displaystyle \frac{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \qquad{(8)}\]

Those assumptions are more or the less the basis upon which various bricks in MFront are built (StandardElasticity, StandardElastoViscoplasticity, DDIF2).

Thanks to Equation (6), one sees that the product \(-J^{-1}\,{\displaystyle \frac{\displaystyle \partial \vec{F}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) only contains the \(6\) first columns of the \(J^{-1}\). This allows identifying \({\displaystyle \frac{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) with the \(6\times 6\) upper-left submatrx of \(J^{-1}\). Let \(J_{\underline{\varepsilon}^{\mathrm{el}}}^{-1}\) this submatrix:

\[ J_{\underline{\varepsilon}^{\mathrm{el}}}^{-1}={\displaystyle \frac{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{el}}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}} \]

Thanks to Equation (8), the consistent tangent operator finally reads:

\[ {\displaystyle \frac{\displaystyle \mathrm{d}\underline{\sigma}}{\displaystyle \mathrm{d}\Delta\,\underline{\varepsilon}^{\mathrm{to}}}}={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\cdot\,J_{\underline{\varepsilon}^{\mathrm{el}}}^{-1} \qquad{(9)}\]

3.1.3.1 The getPartialJacobianInvert method

The getPartialJacobianInvert method allows the computation of the \(J_{\underline{\varepsilon}^{\mathrm{el}}}^{-1}\) tensors as follows:

Stensor4 iJe;
getPartialJacobianInvert(iJe);

The getPartialJacobianInvert method assumes that the jacobian matrix is decomposed (using the LU algorithm with partial pivoting in the current version of TFEL), which is guaranteed in the @TangentOperatorBlock. So \(J_{\underline{\varepsilon}^{\mathrm{el}}}^{-1}\) is computed by solving \(6\) linear systems. In particular, the jacobian matrix is not inverted.

A typical computation of the consistent tangent operator is as follows:

    StiffnessTensor De;
    Stensor4 iJe;
    computeElasticStiffness<N,Type>::exe(De,lambda,mu);
    getPartialJacobianInvert(iJe);
    Dt = De * Je;

This code is loosely the one generated by the StandardElastoViscoplasticy brick.

3.1.3.2 Extensions to more general class of behaviours

The getPartialJacobianInvert may take several arguments which will give the derivatives of the other integration variables.

For example, the following code can be used to compute the consistent tangent operator of a behaviour in which the stress computed using a standard damaging law of the form \({\left.\underline{\sigma}\right|_{t+\Delta\,t}}={\left(1-{\left.d\right|_{t+\Delta\,t}}\right)}\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\left.\underline{\varepsilon}^{\mathrm{el}}\right|_{t+\Delta\,t}}\)²\(\mbox{}^{,}\)³:

    StiffnessTensor De;
    Stensor4 iJe;
    Stensor iJd;
    computeUnalteredElasticStiffness<N,Type>::exe(De, lambda, mu);
    getPartialJacobianInvert(iJe, iJd);
    Dt = (1 - d) * De * iJe - ((De * eel) ^ iJd);

3.2 Extension to generalised behaviours

The idea is to extend this method for generalised behaviours.

3.2.1 Discussion

In Section 3.1.2, we discussed various drawbacks of a direct use of Equation (5).

The equivalent of the consistent tangent operator for generalised behaviours are so called tangent operator blocks which can be:

• the derivatives of one of the thermodynamic force with respect to the increment of a gradient.
• the derivative of one of the thermodynamic force with respect to the increment of an external state variable gradient.
• the derivative of a state variable with respect to the increment of a gradient.
• the derivative of a state variable with respect to the increment of an external state variable.

In the spirit of the Section 3.1, the computation of those derivatives boils down to computing the derivatives of the increment of the state variables \(\Delta\,\vec{Y}\) with respect to the increment of the considered variable \(X\) (\(X\) can thus be a gradient or an external state variable).

As the derivative \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,X}}\) can not be assumed to be sparse, two solutions can be considered:

• extend the getPartialJacobianInvert method.
• have access to the blocks of the invert of the jacobian matrix.

In both cases, \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\) matrix must be computed. The next paragraph is devoted to this computation.

3.2.2 Computation of the \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\) matrix.

The \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\) matrix can be decomposed by blocks:

\[ {\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}} = \begin{pmatrix} {\displaystyle \frac{\displaystyle \partial f_{y_{1}}}{\displaystyle \partial \Delta\, X}} \\ \vdots \\\ {\displaystyle \frac{\displaystyle \partial f_{y_{n}}}{\displaystyle \partial \Delta\, X}} \\ \end{pmatrix} \qquad{(10)}\]

where \({\left(y_{1},\cdots\,y_{n}\right)}\) denotes the set of integration variables, i.e.:

\[ \vec{Y}=\begin{pmatrix} y_{1}\\ \vdots\\ y_{n} \end{pmatrix} \quad\text{and}\quad \vec{F}=\begin{pmatrix} f_{y_{1}}\\ \vdots\\ f_{y_{n}} \end{pmatrix} \]

The advantage of the Decomposition (10) is that each blocks \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\, X}}\) is a tensorial object that can be computed using the facilities offered by TFEL/Math library

If one of the block \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\, X}}\) appears in a @TangentOperator code block, MFront automatically defines the matrix \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\) as an hidden variable and initializes it to zero. This does mean that only non zero derivatives must be defined by the user

For a standard integration variable, a view called \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\, X}}\) is also automatically defined. For array of integration variables, a local function object \({\displaystyle \frac{\displaystyle \partial f_{y_{i}}}{\displaystyle \partial \Delta\, X}}\) is defined. This mechanism is similar to the one used to access the jacobian by block.

3.2.3 The getIntegrationVariablesDerivatives_X local function objects

When required in a @TangentOperator code block (i.e. when used), the getIntegrationVariablesDerivatives_X is defined as a local function object which may the blocks of \(-J^{-1}\,{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\)

3.2.4 Block decomposition of the invert of the jacobian

The invert of the jacobian \(J\) can be decomposed as:

\[ J^{-1}= \begin{pmatrix} J^{-1}_{{\left(f_{y_{1}}, y_{1}\right)}}&\cdots&J^{-1}_{{\left(f_{y_{1}}, y_{n}\right)}}\\ \vdots & \ddots &\vdots \\ J^{-1}_{{\left(f_{y_{n}}, y_{1}\right)}}&\cdots&J^{-1}_{{\left(f_{y_{n}}, y_{n}\right)}}\\ \end{pmatrix} \]

When required in the @TangentOperator code block (i.e. when used), for each pair of state variables \(y\), \(z\), a variable called iJ_x_z is defined:

• If \(y\) and \(z\) are not arrays, iJ_x_z is directly a view to the corresponding block of the invert of the jacobian.
• If \(y\) and/or \(z\) are arrays, iJ_x_z is a function object taking indices as arguments and returning a view to the corresponding block of the invert of the jacobian.

For example:

@TangentOperator{
 // iJ_eel_p is the block corresponding to:
 // 1. the elastic strain (row)
 // 2. the equivalent viscoplastic strain
}

3.3 Applications

3.3.1 A first example of a viscoplastic behaviour coupled with heat transfer

3.3.1.1 Description of the behaviour

Let us consider the following simple isotropic viscoplastic behaviour.

The total strain \(\underline{\varepsilon}^{\mathrm{to}}\) is split into an elastic part \(\underline{\varepsilon}^{\mathrm{el}}\), a thermal strain \(\underline{\varepsilon}^{\mathrm{th}}\) and a viscoplastic part \(\underline{\varepsilon}^{\mathrm{vp}}\):

\[ \Delta\,\underline{\varepsilon}^{\mathrm{to}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{th}}+\Delta\,\underline{\varepsilon}^{\mathrm{vp}} \]

The thermal strain is given by:

\[ \underline{\varepsilon}^{\mathrm{th}}=\alpha{\left(T\right)}\,{\left(T-T_{\mathrm{ref}}\right)} \]

where:

• \(T\) is the temperature, assumed to be an external state variable.
• \(T_{\mathrm{ref}}\) is a reference temperature.

The stress tensor \(\underline{\sigma}\) is related to the elastic strain by the Hooke law:

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

The viscoplastic strain rate is given by the famous Norton-Hoff law:

\[ \underline{\dot{\varepsilon}}^{\mathrm{vp}}= \dot{p}\,\underline{n} = A\,\exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle T}}}\right)}\,\sigma_{\mathrm{eq}}^{E}\,\underline{n} \]

where:

• \(\sigma_{\mathrm{eq}}\) is the equivalent von Mises stress.
• \(\underline{n}\) is the von Mises normal.
• \(A\), \(T_{a}\) and \(E\) are material parameters.
• \(p\) is the equivalent plastic strain.

The integration variables are the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) and the equivalent plastic strain \(p\).

The tangent operator blocks to be computed are \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) and \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,T}}\).

The computation is of \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{to}}}}\) readily follows Section 3.1.3.1 and will not be described here.

3.3.1.2 Implicit scheme

\[ \left\{ \begin{aligned} f_{\underline{\varepsilon}^{\mathrm{el}}}&=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{th}}+\Delta\,p\,{\left.\underline{n}\right|_{t+\theta\,\Delta\,t}}-\Delta\underline{\varepsilon}^{\mathrm{to}}\\ f_{p}&=\Delta\,p-A\,\exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle {\left.T\right|_{t+\theta\,\Delta\,t}}}}}\right)}\,{\left({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\right)}^{E} \end{aligned} \right. \]

3.3.2 Computation of \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta T}}\)

\[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta T}}&={\displaystyle \frac{\displaystyle \partial \Delta\underline{\varepsilon}^{\mathrm{th}}}{\displaystyle \partial \Delta T}}={\displaystyle \frac{\displaystyle \partial {\left.\underline{\varepsilon}^{\mathrm{th}}\right|_{t+\Delta\,t}}}{\displaystyle \partial \Delta T}}= {\left(\alpha{\left({\left.T\right|_{t+\Delta\,t}}\right)}+\left.{\displaystyle \frac{\displaystyle \mathrm{d}\alpha}{\displaystyle \mathrm{d}T}}\right|_{{\left.T\right|_{t+\Delta\,t}}}\,{\left({\left.T\right|_{t+\Delta\,t}}-T_{\mathrm{ref}}\right)}\right)}\,\underline{I}\\ {\displaystyle \frac{\displaystyle \partial f_{p}}{\displaystyle \partial \Delta T}}&=-A\,\theta\,{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle {\left.T\right|_{t+\theta\,\Delta\,t}}^{2}}}}\,\exp{\left(-{{\displaystyle \frac{\displaystyle T_{a}}{\displaystyle {\left.T\right|_{t+\theta\,\Delta\,t}}}}}\right)}\,{\left({\left.\sigma_{\mathrm{eq}}\right|_{t+\theta\,\Delta\,t}}\right)}^{E} \end{aligned} \right. \]

∂fεᵉˡ∕∂ΔT = ...;
∂fp∕∂ΔT = ...;

3.3.3 Computing \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,T}}\) using getIntegrationVariablesDerivatives_T

@TangentOperatorBlock{
  ∂fεᵉˡ∕∂ΔT = ...;
  ∂fp∕∂ΔT = ...;
  auto ∂Δεᵉˡ∕∂ΔT = Stensor{};
  getIntegrationVariablesDerivatives_T(∂Δεᵉˡ∕∂ΔT);
  const auto λ =...;
  const auto μ =...;
  const auto ∂λ∕∂T =...;
  const auto ∂μ∕∂T =...;
  const auto De = λ  ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ I₄;  
  const auto ∂De∕∂T = ∂λ∕∂T  ⋅ (I₂ ⊗ I₂) + 2 ⋅ ∂μ∕∂T ⋅ I₄;  
  ∂σ∕∂ΔT = De ⋅ ∂Δεᵉˡ∕∂ΔT + ∂De∕∂T ⋅ εᵉˡ;
}

3.3.4 Computing \({\displaystyle \frac{\displaystyle \partial \underline{\sigma}}{\displaystyle \partial \Delta\,T}}\) using the block decomposition of the invert of jacobian

@TangentOperatorBlock{
  const auto λ =...;
  const auto μ =...;
  const auto De = λ  ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ I₄;  
  // computation of ∂σ∕∂ΔT
  ∂fεᵉˡ∕∂ΔT = ...;
  ∂fp∕∂ΔT = ...;
  const auto ∂Δεᵉˡ∕∂ΔT = -(iJ_fεᵉˡ_εᵉˡ ⋅ ∂fεᵉˡ∕∂ΔT + iJ_fεᵉˡ_p ⋅ ∂fp∕∂ΔT);
  const auto ∂λ∕∂T =...;
  const auto ∂μ∕∂T =...;
  const auto ∂De∕∂T = ∂λ∕∂T  ⋅ (I₂ ⊗ I₂) + 2 ⋅ ∂μ∕∂T ⋅ I₄;  
  ∂σ∕∂ΔT = De ⋅ ∂Δεᵉˡ∕∂ΔT + ∂De∕∂T ⋅ εᵉˡ;
}