This article shows how to implement hyperviscoelastic behaviours in MFront.

The hyperelastic behaviour is left unspecified. We only assume that one can compute \(\underline{S}^{\infty}\) and its derivative \({\displaystyle \frac{\displaystyle \partial \underline{S}^{\infty}}{\displaystyle \partial \underline{C}}}\).

The reader can refer to the following pages describing the implementation of various hyperelastic behaviours:

• Signorini
• Ogden

The formalism used by most hyperelastic behaviours is described in this page.

An example of implementation based on the Signorini hyperelastic behaviour is available here: SignoriniHyperViscoelasticity.mfront

Compatibility with previous versions of TFEL

This implementation presented here is based on the version \(3.0\) of TFEL. However, only a few adaptations are required for used with previous versions:

• the first line of the file must be replace by: @Parser DefaultFiniteStrainParser;
• the power function is not available. It can be replaced by explicit multiplications.
• the computeDeterminantSecondDerivative is not available. It can be replaced by the following code: Stensor4::dsquare(C)-(C^id)-I1*Stensor4::Id()+(id^dI2_dC)
• the use of the auto keyword is only permitted by the C++11 standard, which is default for TFEL-3.0. The user is left with two choices:
  • replace auto with explicit types: real for scalars, Stensor for symmetric tensors, Stensor4 for linear forms on symmetric tensors (“fourth order tensors”). However, in this case, the implementation can be less efficient as it leads to an explicit evaluation of intermediate results.
  • specify the appropriate flag to the C++ compiler. Assuming that the current shell is bash and that the compiler is g++, this can be done by invoking mfront as follows:

CXXFLAGS="`tfel-config --oflags` -std=c++11" mfront --obuild --interface=aster SignoriniHyperViscoelasticity-tfel-2.0.3.mfront`

An implementation where those changes were made is available here

Computation of the Cauchy stress

The hyperelastic behaviour is supposed to decouple the volumetric stress \(\underline{S}^{v\,\infty}\) and the isochoric stress \(S^{i\,\infty}\): \[ \underline{S}^{\infty}=\underline{S}^{v\,\infty}+\underline{S}^{i\,\infty} \]

The second Piola-Kirchhoff stress is given by:

\[ \underline{S}= \underline{S}^{\infty}+\sum_{i=1}^{N}\underline{H}_{i} \]

where \(\underline{H}_{i}\) are the overstresses.

Evolution of the overstresses

The algorithm used to update \(\underline{H}_{i}\) is the following:

\[ {\left.\underline{H}_{i}\right|_{t+\Delta\,t}}=\exp{\left(-{{\displaystyle \frac{\displaystyle dt}{\displaystyle \tau_{i}}}}\right)}{\left.\underline{H}_{i}\right|_{t}}+g_{i}\,\tau_{i}{\left(1-\exp{\left(-{{\displaystyle \frac{\displaystyle dt}{\displaystyle \tau_{i}}}}\right)}\right)}{{\displaystyle \frac{\displaystyle {\left({\left.S^{i\,\infty}\right|_{t+\Delta\,t}}-{\left.S^{i\,\infty}\right|_{t}}\right)}}{\displaystyle dt}}} \]

Computation of second Piola-Kirchhoff stress

// updating the viscoelastic stresses and computing the Second
// Piola-Kirchhoff stress
StressStensor S = Sv+Se;
if(dt>0){
  for(unsigned short i=0;i!=Nv;++i){
    const auto dtr = dt/tau[i];
    e[i] = exp(-dtr);
    H[i] = e[i]*H[i]+g[i]*(1-e[i])/dtr*(Se-Se_1);
    S += H[i];
  }
}

Computation of the Cauchy stress

The function convertSecondPiolaKirchhoffStressToCauchyStress converts the second Piola Kirchhoff stress to the Cauchy stress using the deformation gradient. It is used as follows:

sig = convertSecondPiolaKirchhoffStressToCauchyStress(S,F1);

Computation of the consistent tangent operator

Assuming that \({\displaystyle \frac{\displaystyle \partial \underline{S}^{v\,\infty}}{\displaystyle \partial \underline{C}}}\) and \({\displaystyle \frac{\displaystyle \partial \underline{S}^{i\,\infty}}{\displaystyle \partial \underline{C}}}\) are known, the expression of the consistent tangent operator is:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial \underline{S}^{v\,\infty}}{\displaystyle \partial \underline{C}}} &= {\displaystyle \frac{\displaystyle \partial \underline{S}^{v\,\infty}}{\displaystyle \partial \underline{C}}}+{\displaystyle \frac{\displaystyle \partial \underline{S}^{i\,\infty}}{\displaystyle \partial \underline{C}}}+ \sum_{i=1}^{N}{\displaystyle \frac{\displaystyle \partial {\left.\underline{H}_{i}\right|_{t+\Delta\,t}}}{\displaystyle \partial \underline{C}}}\\ &= {\displaystyle \frac{\displaystyle \partial \underline{S}^{v\,\infty}}{\displaystyle \partial \underline{C}}}+ {\displaystyle \frac{\displaystyle \partial \underline{S}^{i\,\infty}}{\displaystyle \partial \underline{C}}}\,\underbrace{{\left(1+\sum_{i=1}^{N}g_{i}\,\tau_{i}\,{{\displaystyle \frac{\displaystyle {\left(1-\exp{\left(-{{\displaystyle \frac{\displaystyle dt}{\displaystyle \tau_{i}}}}\right)}\right)}}{\displaystyle dt}}}\right)}}_{c} \end{aligned} \]

This expression can be implemented in a straightforward manner:

auto c = real(1);
if(dt>0){
  for(unsigned short i=0;i!=Nv;++i){
    const auto dtr = dt/tau[i];
    c += g[i]*(1-e[i])/dtr;
  }
}
dS_dC = dSv_dC+c*dSi_dC;