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.
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 @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 Equation (1), the parameter theta is considered
fixed. If this parameter varies, the solution \(\Delta\,\vec{Y}\) of the implicit system
changes, which can be written as follows:
\[ F{\left(\Delta\,\vec{Y}{\left(\theta\right)},\theta\right)}=0 \qquad{(3)}\]
Using the implicit function theorem, one can evaluate the sensibility of the solution with respect to \(\theta\) as follows:
\[ {\displaystyle \frac{\displaystyle \partial \Delta\,\vec{Y}}{\displaystyle \partial \theta}} = -{\left({\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,\vec{Y}}}\right)}^{-1}\,{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \theta}} = -J^{-1}\,{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \theta}} \qquad{(4)}\]
where \(J={\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Delta\,\vec{Y}}}\) is the jacobian of the implicit system.
Equation (4) shows that \(\Delta\,\vec{Y}\) is a smooth function of \(\theta\) as long as the jacobian is invertible which is the requirements of most algorithms available.
Equation (4) could be integrated using an explicit algorithm and theorically find \(\Delta\,\vec{Y}{\left(\theta\right)}\) to the Implicit System (1) if \(\Delta\,\vec{Y}{\left(0\right)}\) is known.
Note
Equation (4) can be integrated as follows:
\[ \Delta\,\vec{Y}{\left(\theta\right)} = \Delta\,\vec{Y}{\left(0\right)}-\int_{0}^{\theta}J^{-1}\,{\displaystyle \frac{\displaystyle \partial F}{\displaystyle \partial \Theta}}\,{{\mathrm{d}}}\,\Theta \qquad{(5)}\]
For time dependent behaviours, Equation (3) is trival to solve for \(\theta=0\), \(\Delta\,\vec{Y}{\left(0\right)}\) corresponds to the fully explicit Euler approximation:
\[ \Delta\,\vec{Y}{\left(0\right)} = G{\left({\left.\vec{Y}\right|_{t}},t\right)}\,\Delta\,t \]
The case of time independent behaviours is more involved. Let us consider for instance, a plastic behaviour with isotropic hardening for which the implicit equation \(f_{p}\) associated with the plastic multiplier \(p\) is:
\[ f_{p} = {\left.\seq\right|_{t+\theta\,\Delta\,t}}-R{\left({\left.p\right|_{t+\theta\,\Delta\,t}}\right)} \]
In this case, the derivatives \(f_{p}\) with respect to the internal state variables are null for \(\theta=0\).
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{(6)}\]
Equation (6) 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{(7)}\]
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{(8)}\]
Equation (8) may be extented almost directly to generalised behaviours. However, direct usage of Equation (8) has several practical disadvantages:
Hence, a computation of the consistent tangent based on Equation (8) may imply a significant performance hit which can be avoided in many most common cases, as discussed in the following section.
Equation (8) 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 (9), 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 (11), 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{(12)}\]
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 (8).
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{(13)}\]
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 (13) 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 ⋅ εᵉˡ;
}