This page describe how to implement the Ogden hyperelastic behaviour.

The whole implementation is available here: Ogden.mfront

This implementation is compatible with TFEL version 3.1 which has introduced several overloaded versions of the computeIsotropicFunction and computeIsotropicFunctionDerivative methods which reduce the number of evaluation of the pow function.

An implementation compatible with TFEL version 3.0 is available here: Ogden-tfel-3.0.0.mfront

Description

Let \(\underline{C}\) be the right Cauchy tensor. The three invariants of \(\underline{C}\) are defined by: \[ \left\{ \begin{aligned} I_{1} &= {\mathrm{tr}{\left(C\right)}} \\ I_{2} &= {{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\left({\left({\mathrm{tr}{\left(\underline{C}\right)}}\right)}^{2}-{\mathrm{tr}{\left(\underline{C}^{2}\right)}}\right) \\ I_{3} &= \det{\left(\underline{C}\right)} \end{aligned} \right. \]

\(J\) is the determinant of deformation gradient (\(J=\sqrt{I_{3}}\))

The Ogden hyperelastic behaviour is described by a decoupled potential \(W\) decomposed in an volumetric part \(W^{\text{v}}\) and an isochoric part \(W^{\text{i}}\):

\[ W{\left(\underline{C}\right)}=W^{v}{\left(J\right)}+\bar{W}^{i}{\left(\bar{\lambda}_{1},\bar{\lambda}_{2},\bar{\lambda}_{3}\right)} \]

where:

• \({\left(\bar{\lambda}_{i}\right)}_{i\in \{1,2,3\}}\) are the eigenvalues of the isochoric right Cauchy tensor \(\underline{\bar{C}}\). \({\left(\bar{\lambda}_{i}\right)}_{i\in \{1,2,3\}}\) are defined by: \(\bar{\lambda}_{i}=\bar{c}_{I_{3}}\,\lambda_{i}\)
• \(\lambda_{i}\) are the right Cauchy stress tensor \(\underline{C}\)
• \(\bar{c}_{I_{3}}=I_{3}^{-1/3}\)
• \(\underline{\bar{C}}=J^{-2/3}\,\underline{C}= I_{3}^{-1/3}\,\underline{C}= \bar{c}_{I_{3}}\,\underline{C}\)

A general derivation of the stress and of the consistent tangent operator for this kind of hyperelastic behaviours is described here. The reader shall refer to this page for intermediate computations that won’t be detailled here.

For our implementation, we choose \(W^{v}\) as follows: \[ W^{v}{\left(J\right)}={{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,K\,(J-1)^{2} \] where \(K\) is the bulk modulus.

The Ogden behaviour mostly refers to the following choice of \(\bar{W}^{i}\): \[ \bar{W}^{i}{\left(\bar{\lambda}_{1},\bar{\lambda}_{2},\bar{\lambda}_{3}\right)}=\sum_{p=1}^{N}\bar{W}^{i}_{p} =\sum_{p=1}^{N}{{\displaystyle \frac{\displaystyle \mu_{p}}{\displaystyle \alpha_{p}}}}{\left(\bar{\lambda}_{1}^{\alpha_{p}/2}+\bar{\lambda}_{2}^{\alpha_{p}/2}+\bar{\lambda}_{3}^{\alpha_{p}/2}-3\right)} =\sum_{p=1}^{N}{{\displaystyle \frac{\displaystyle \mu_{p}}{\displaystyle \alpha_{p}}}}\bar{c}_{I_{3}}^{\alpha_{p}/2}\,f_{p}{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}-{{\displaystyle \frac{\displaystyle 3\,\mu_{p}}{\displaystyle \alpha_{p}}}} \]

with \(f_{p}{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}=\lambda_{1}^{\alpha_{p}/2}+\lambda_{2}^{\alpha_{p}/2}+\lambda_{3}^{\alpha_{p}/2}\)

In the following, we will make the simple choice (\(p=1\)): \[ \bar{W}^{i}{\left(\bar{\lambda}_{1},\bar{\lambda}_{2},\bar{\lambda}_{3}\right)}=\sum_{p=1}^{N}\bar{W}^{i}_{p} ={{\displaystyle \frac{\displaystyle \mu}{\displaystyle \alpha}}}\bar{c}_{I_{3}}^{\alpha/2}\,f{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}-{{\displaystyle \frac{\displaystyle 3\,\mu}{\displaystyle \alpha}}} \]

The general expression of the second Piola-Kirchhoff stress is:

\[ \underline{S}=2\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial \underline{C}}}+2{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \underline{C}}}=\underline{S}^{v}+\underline{S}^{i} \]

where \(\underline{S}^{v}=2\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial \underline{C}}}\) is the volumetric part of the the second Piola-Kirchhoff stress and \(\underline{S}^{i}=2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \underline{C}}}\) is the isochoric part.

Implementation

The implementation of this behaviour is decomposed in two parts:

• The computation of the Cauchy stress.
• The computation of the consistent tangent operator.

We use the DefaultFiniteStrain implementation as this behaviour does not require an integration algorithm.

@DSL DefaultFiniteStrain;

Then, we define the material parameters:

@Parameter alpha =     28.8;
@Parameter mu    =    27778;
@Parameter K     = 69444444;

In the following, we will mostly use \(\alpha/2\):

const auto a = alpha/2;

Computation of the Cauchy stress

Computation of the right Cauchy tensor \(\underline{C}\), the invariants and the invariants derivatives

The identity tensor is stored in the id variable for a shorter and cleaner code:

const auto id = Stensor::Id();

Then, we compute the determinant of the deformation gradient at the end of the time step, called F1:

const auto J  = det(F1);

Then we compute the right Cauchy tensor \(\underline{C}\):

const auto C  = computeRightCauchyGreenTensor(F1);

For the computation of the invariants of \(\underline{C}\) and theirs derivatives, we need to compute \(\underline{C}^{2}\). Since the multiplication of two symmetric tensors leads to a non symmetric tensors, we use the special function square which is a more efficient equivalent of syme(C*C):

const auto C2 = square(C);

The first invariant \(I_{1}\) is the trace of \(\underline{C}\):

const auto I1 = trace(C);

The second invariant \(I_{2}\) is computed as follows:

const auto I2 = (I1*I1-trace(C2))/2;

The derivative of \({\displaystyle \frac{\displaystyle \partial I_{2}}{\displaystyle \partial \underline{C}}}\) is equal to \({\mathrm{tr}{\left(\underline{C}\right)}}\,\underline{I}-\underline{C}\), which is readily written as:

const auto dI2_dC = I1*id-C;

At this stage, \({\displaystyle \frac{\displaystyle \partial I_{2}}{\displaystyle \partial \underline{C}}}\) is not evaluated: the associated variable dI2_dC is an evaluation tree that stores the computation that will be performed latter. This is an essential optimization techniques used in MFront to achieve optimal performances. If we wanted or needed to explicitly evaluate, we could either use the eval function or explicitly precise the type of dI2_dC.

The derivative \({\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}}\) of the third invariant is given by: \[ {\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}} = \underline{C}^{2}-I_{1}\,\underline{C}+I_{2}\,\underline{I} \]

The computation of the third invariant and its derivatives is straightforward:

const auto I3      = J*J;
const auto dI3_dC  = computeDeterminantDerivative(C);

Computation of the volumetric part of the second Piola Kirchhoff stress

The volumetric part of the second Piola Kirchhoff stress \(\underline{S}^{v}\) is given by: \[ \underline{S}^{v} = {{\displaystyle \frac{\displaystyle 1}{\displaystyle J}}}\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}\,{\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}} \]

Here, \[ {\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}=K\,(J-1) \]

We can implement it as follows:

const auto dPv_dJ   = K*(J-1);
const StressStensor Sv = dPv_dJ/J*dI3_dC;

Here, we have explicitly given the type of Sv to force its evaluation.

Computation of the isochoric part of the second Piola Kirchhoff stress

By chain rule, one have: \[ \begin{aligned} \underline{S}^{i}=2\,{\displaystyle \frac{\displaystyle \partial \bar{W}^{i}}{\displaystyle \partial \underline{C}}}= \mu\,\bar{c}_{I_{3}}^{\alpha/2-1}\,f{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}{\displaystyle \frac{\displaystyle \partial \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}}}+ {{\displaystyle \frac{\displaystyle 2\,\mu\,\bar{c}_{I_{3}}^{\alpha/2}}{\displaystyle \alpha}}}\,{\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\\ \end{aligned} \]

\({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\) is an isotropic function of \(\underline{C}\), thus we can use the various methods provided by the TFEL library to evaluate \({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\) and its derivative \({\displaystyle \frac{\displaystyle \partial^{2} f}{\displaystyle \partial \underline{C}^{2}}}\) using the computeIsotropicFunction and computeIsotropicFunctionDerivative static methods of the Stensor class.

First, we need to compute \(\bar{c}_{I_{3}}\) and its derivatives. To minimize the number of call to the pow function, we also compute \(\bar{c}_{I_{3}}^{a-2}\):

const auto iJb        =  1/cbrt(I3);
const auto iJb2       =  power<2>(iJb);
const auto iJb4       =  iJb2*iJb2;
const auto iJb7       =  iJb4*power<3>(iJb);
const auto c          = pow(iJb,a-2);
// derivatives
const auto diJb_dI3   = -iJb4/3;
const auto d2iJb_dI32 = 4*iJb7/9;
const auto diJb_dC    = diJb_dI3*dI3_dC;

Then, we can compute the eigenvalues and the eigentensors of \(\underline{C}\):

tvector<3u,real> vp;
tmatrix<3u,3u,real> m;
std::tie(vp,m) = C.computeEigenVectors();

In C++-17, the previous code can be replaced by:

const auto [vp,m] = C.computeEigenVectors();

We can now evaluate \(f\) and \({\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\).:

const tvector<3u,real> pwv  = {pow(vp(0),a-2),pow(vp(1),a-2),pow(vp(2),a-2)};
const tvector<3u,real> d2fv = {a*(a-1)*pwv(0),a*(a-1)*pwv(1),a*(a-1)*pwv(2)};
const tvector<3u,real> dfv  = {a*vp(0)*pwv(0),a*vp(1)*pwv(1),a*vp(2)*pwv(2)};
const auto fv    = vp(0)*vp(0)*pwv(0)+vp(1)*vp(1)*pwv(1)+vp(2)*vp(2)*pwv(2);
const auto df_dC = Stensor::computeIsotropicFunction(dfv,m);

Finally, we can compute \(\underline{S}^{i}\):

const StressStensor Si = mu*c*iJb*((fv*diJb_dC+(iJb/a)*df_dC));

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(Sv+Si,F1);

Computation of the consistent tangent operator

The most direct expression the consistent tangent operator is given by \({\displaystyle \frac{\displaystyle \partial \underline{S}}{\displaystyle \partial \underline{C}}}\). We let MFront make the appropriate conversion to the consistent tangent operator expected by the solver.

Here, \[ {\displaystyle \frac{\displaystyle \partial \underline{S}}{\displaystyle \partial \underline{C}}}={\displaystyle \frac{\displaystyle \partial \underline{S}^{v}}{\displaystyle \partial \underline{C}}}+{\displaystyle \frac{\displaystyle \partial \underline{S}^{i}}{\displaystyle \partial \underline{C}}} \]

Computation of the \({\displaystyle \frac{\displaystyle \partial^{2} I_{3}}{\displaystyle \partial \underline{C}^{2}}}\)

To compute \({\displaystyle \frac{\displaystyle \partial^{2} I_{3}}{\displaystyle \partial \underline{C}^{2}}}\), one may use the computeDeterminantSecondDerivative function:

const auto d2I3_dC2 = computeDeterminantSecondDerivative(C);

Computation of the volumetric part of the consistent tangent operator

\({\displaystyle \frac{\displaystyle \partial \underline{S}^{v}}{\displaystyle \partial \underline{C}}}\) is given by:

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial \underline{S}^{v}}{\displaystyle \partial \underline{C}}} &={\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial J}}{\left({{\displaystyle \frac{\displaystyle 1}{\displaystyle J}}}\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}\right)}\,{\displaystyle \frac{\displaystyle \partial J}{\displaystyle \partial I_{3}}}\,{\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}}+{{\displaystyle \frac{\displaystyle 1}{\displaystyle J}}}\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}\,{\displaystyle \frac{\displaystyle \partial^{2} I_{3}}{\displaystyle \partial \underline{C}^{2}}}\\ &={\left({\displaystyle \frac{\displaystyle \partial^{2} W^{v}}{\displaystyle \partial J^{2}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle J}}}\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}\right)}\,\frac{1}{2\,I_{3}}\,{\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{C}}}+{{\displaystyle \frac{\displaystyle 1}{\displaystyle J}}}\,{\displaystyle \frac{\displaystyle \partial W^{v}}{\displaystyle \partial J}}\,{\displaystyle \frac{\displaystyle \partial^{2} I_{3}}{\displaystyle \partial \underline{C}^{2}}}\\ \end{aligned} \]

Here: \[ {\displaystyle \frac{\displaystyle \partial^{2} W^{v}}{\displaystyle \partial J^{2}}}=K \]

Those expression are straightforward to implement:

const auto d2Pv_dJ2 = K;
dS_dC = ((d2Pv_dJ2-dPv_dJ/J)/(2*I3)*(dI3_dC^dI3_dC)
        +dPv_dJ/J*d2I3_dC2);

Computation of the isochoric part of the consistent tangent operator

As detailled here, the expression of \({\displaystyle \frac{\displaystyle \partial S^{i}}{\displaystyle \partial \underline{C}}}\) is given by:

\[ \begin{aligned} {{\displaystyle \frac{\displaystyle 1}{\displaystyle \mu}}}\,{\displaystyle \frac{\displaystyle \partial S^{i}}{\displaystyle \partial \underline{C}}} &=\frac{\alpha-2}{2}\,\bar{c}_{I_{3}}^{\alpha/2-2}\,f{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}{\displaystyle \frac{\displaystyle \partial \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}}}\\ &+\bar{c}_{I_{3}}^{\alpha/2-1}\,f{\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)}\,{\displaystyle \frac{\displaystyle \partial^{2} \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}^{2}}}+\bar{c}_{I_{3}}^{\alpha/2-1}\,{\displaystyle \frac{\displaystyle \partial \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\\ &+\bar{c}_{I_{3}}^{\alpha/2-1}\,{\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{c}_{I_{3}}}{\displaystyle \partial \underline{C}}}+ 2\,{{\displaystyle \frac{\displaystyle \bar{c}_{I_{3}}^{\alpha/2}}{\displaystyle \alpha}}}\,{\displaystyle \frac{\displaystyle \partial^{2} f}{\displaystyle \partial \underline{C}^{2}}} \end{aligned} \]

This expression requires the computation of \({\displaystyle \frac{\displaystyle \partial^{2} f}{\displaystyle \partial \underline{C}^{2}}}\), which is a direct application of the static method computeIsotropicFunctionDerivative from the Stensor class:

const tvector<3u,real> d2fv = {a*(a-1)*pwv(0),a*(a-1)*pwv(1),a*(a-1)*pwv(2)};
const auto d2f_dC2    = Stensor::computeIsotropicFunctionDerivative(dfv,d2fv,vp,m,1.e-12);

With this result, the computation of \({\displaystyle \frac{\displaystyle \partial^{2} f}{\displaystyle \partial \underline{C}^{2}}}\) is a direct translation from the previous equation:

const auto d2iJb_dI32 = 4*iJb7/9;
const auto d2iJb_dC2  = d2iJb_dI32*(dI3_dC^dI3_dC)+ diJb_dI3*d2I3_dC2;
dS_dC += mu*c*((a-1)*fv*(diJb_dC^diJb_dC)+iJb*(fv*d2iJb_dC2+
               ((diJb_dC^df_dC)+(df_dC^diJb_dC))+iJb/a*d2f_dC2));

Tests

The comparison of previous derivation of the consistent tangent operator to a numerical approximation is made in the Ogden.cxx test that can be found in the tests/Material directory of the TFEL sources.

The components of the consistent tangent operator match their numerical approximation with an accuracy lower than \(10^{-9}\max_{i,j}{\displaystyle \frac{\displaystyle \partial D}{\displaystyle \partial \underline{C}}}(i,j)\).