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:

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:

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;