Implicit domain specific languagegetIntegrationVariablesDerivatives_TLet \(\vec{Y}\) be a vector holding all the integration variables1. 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:
An implicit scheme turns this ordinary differential equation into a non linear system of equations whose unknowns are is 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 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)\]
The following algorithms are available:
NewtonRaphsonNewtonRaphson_NumericalJacobianPowellDogLeg_NewtonRaphsonPowellDogLeg_NewtonRaphson_NumericalJacobianBroydenPowellDogLeg_BroydenBroyden2LevenbergMarquardtLevenbergMarquardt_NumericalJacobianThoses algorithms are described in the documention of the TFEL/Math library. The main steps of those algorithms and the associated code blocks are depicted in Figure 1.
Integrator code blocks when the numerical evaluation of the jacobian is requestedIn 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 perturbated values of the unknowns.
In 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 @IntegratorIn many cases, rather than updating auxiliary variables during the Newton iterations, it can be more pratical to compute its increment, defined in by 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.
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}\) are 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)} \]
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)\]
Equation (5) may be extented almost directly to generalised behaviours. However, direct usage of Equation (5) has several practical disadvantages:
Hence, a computation of the consistent tangent based on Equation (5) may imply a significant performance hit which can be avoided in many most common cases, as discussed in the following section.
Equation (5) can be simplified if the following assumptions are made:
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 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)\]
getPartialJacobianInvert methodThe 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.
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}}\)2\(\mbox{}^{,}\)3:
StiffnessTensor De;
Stensor4 iJe;
Stensor iJd;
computeUnalteredElasticStiffness<N,Type>::exe(De, lambda, mu);
getPartialJacobianInvert(iJe, iJd);
Dt = (1 - d) * De * iJe - ((De * eel) ^ iJd);The idea is to extend this method for generalised behaviours.
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:
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:
getPartialJacobianInvert method.In both cases, \({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\, X}}\) matrix must be computed. The next paragraph is devoted to this computation.
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.
getIntegrationVariablesDerivatives_X local function objectsWhen 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}}\)
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:
iJ_x_z is directly a view to the corresponding block of the invert of the jacobian.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
}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:
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:
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.
\[ \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. \]
\[ \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 = ...;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 ⋅ εᵉˡ;
}@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 ⋅ εᵉˡ;
}In MFront, the integration variables refers to variables which are part of the implicit systems, i.e. variables the increments of which are the solutions of the implicit system. A state variable, declared with @StateVariable keyword, is saved from one step to the other and is automatically added to the list of integration variables. An auxiliary state variable, declared with @AuxiliaryStateVariable keyword, is also saved from one step to the other, but is not added to the list of the integration variables. The @IntegrationVariable keyword allows to append a variable to the integration variables but its value wil not be saved from one step to the other.↩︎
This assumes that the damage \(d\) is the second integration variables.↩︎
This code uses the fact the state variables have been updated to their values at the end of the time step before the call to @TangentOperator code block.↩︎