This page describe how to implement the Signorini hyperelastic behaviour (see [2]).

The whole implementation is available here: Signorini.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 Signorini.mfront`

An implementation where those changes were made is available here

Description

The Signorini behaviour is described by the following decoupled potential: \[ W{\left(\underline{C}\right)}=W^{v}{\left(J\right)}+W^{i}{\left(\bar{I}_{1},\bar{I}_{2}\right)} \] where

• \(\underline{C}\) is 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}}\))
• \(\bar{I}_{1}\) and \(\bar{I}_{2}\) are the first and second invariants of the modified right Cauchy tensor \(\underline{\bar{C}}\) defined by: \(\underline{\bar{C}}=J^{-2/3}\,\underline{C}= I_{3}^{-1/3}\,\underline{C}\)

A general derivation of the stress and 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 Signorini behaviour mostly refers to the following choice of \(W^{i}\): \[ W^{i}= C_{10}\,{\left(\bar{I}_{1}-3\right)}+2\,C_{20}\,{\left(\bar{I}_{1}-3\right)}^{2}+C_{01}\,{\left(\bar{I}_{2}-3\right)} \]

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;

We use parameters for defining the material constants:

@Parameter K   = 2.939e9;
@Parameter C10 = 2.668e6;
@Parameter C20 = 0.446e6;
@Parameter C01 = 0.271e6;

Finally, we use a local variable to store the consistent tangent operator:

@LocalVariable StiffnessTensor dS_dC;

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  = C2-I1*C+I2*id;

In TFEL 3.1, one may use the computeDeterminantDerivative function to compute \({\displaystyle \frac{\displaystyle \partial I_{3}}{\displaystyle \partial \underline{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

We first compute the invariants of the modified right Cauchy tensor \(\underline{\bar{C}}\) and their derivatives with respect to \(\underline{C}\).

\(\bar{I}_{1}\) and \(\bar{I}_{2}\) are related to the first and second invariants of the right Cauchy tensor by:

\[ \begin{aligned} \bar{I}_{1} &= I_{3}^{-1/3}\,I_{1} = \bar{c}_{I_{3}}\,I_{1}\\ \bar{I}_{2} &= I_{3}^{-2/3}\,I_{2} = \bar{c}_{I_{3}}^{2}\,I_{2}\\ \end{aligned} \]

The expression of the derivatives of \(\bar{I}_{1}\) and \(\bar{I}_{2}\) is detailled here.

The computation of all those terms is implemented as follows:

const auto iJb        =  1/cbrt(I3);
const auto iJb2       =  power<2>(iJb);
const auto iJb4       =  iJb2*iJb2;
const auto diJb_dI3   = -iJb4/3;
const auto diJb_dC    = diJb_dI3*dI3_dC;
const auto I1b        =  I1*iJb;
const auto dI1b_dC    =  iJb*id+I1*diJb_dC;
const auto dI2b_dC    = iJb2*dI2_dC+2*I2*iJb*diJb_dC;

By chain rule, we have: \[ \underline{S}^{i} = 2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \underline{C}}}=2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{1}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{1}}{\displaystyle \partial \underline{C}}}+2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{2}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{2}}{\displaystyle \partial \underline{C}}} \]

In the case of the Signorini behaviour, the derivatives of \(W^{i}\) are: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{1}}} &= C_{10}+2\,C_{20}\,(\bar{I}_{1}-3)\\ {\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{2}}} &= C_{01}\\ \end{aligned} \right. \]

The computation of \(\underline{S}^{i}\) is:

const auto dPi_dI1b    = C10+2*C20*(I1b-3);
const auto dPi_dI2b    = C01;
const StressStensor Si = 2*(dPi_dI1b*dI1b_dC+dPi_dI2b*dI2b_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 second derivatives of the invariants of \(\underline{C}\) and \(\underline{\bar{C}}\)

const auto d2I3_dC2 = computeDeterminantSecondDerivative(C);
const auto d2I2_dC2   = (id^id)-Stensor4::Id();
const auto iJb7       =  iJb4*power<3>(iJb);
const auto d2iJb_dI32 = 4*iJb7/9;
const auto d2iJb_dC2  = d2iJb_dI32*(dI3_dC^dI3_dC)+ diJb_dI3*d2I3_dC2;
const auto d2I1b_dC2  = (id^diJb_dC)+(diJb_dC^id)+I1*d2iJb_dC2;
const auto d2I2b_dC2  = 2*iJb*(dI2_dC^diJb_dC)+ iJb2*d2I2_dC2+
                        2*iJb*(diJb_dC^dI2_dC)+
                        2*I2*(diJb_dC^diJb_dC)+
                        2*I2*iJb*d2iJb_dC2;

Computation of the volumetric part of the consistent tangent operator

\[ \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} \]

For the Signorini behaviour: \[ {\displaystyle \frac{\displaystyle \partial^{2} W^{v}}{\displaystyle \partial J^{2}}}=K \]

const auto d2Pv_dJ2 = K;
const auto dSv_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

\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial^{2} \underline{S}^{i}}{\displaystyle \partial C^{2}}}=& 2\,{\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\displaystyle \partial \bar{I}_{1}^{2}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{1}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{I}_{1}}{\displaystyle \partial \underline{C}}}+ 2\,{\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\partial \bar{I}_{1}\partial \bar{I}_{2}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{1}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{I}_{2}}{\displaystyle \partial \underline{C}}}+ 2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{1}}}{\displaystyle \frac{\displaystyle \partial^{2} \bar{I}_{1}}{\displaystyle \partial \underline{C}^{2}}}+\\ & 2\,{\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\displaystyle \partial \bar{I}_{2}^{2}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{2}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{I}_{2}}{\displaystyle \partial \underline{C}}}+ 2\,{\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\partial \bar{I}_{1}\partial \bar{I}_{2}}}\,{\displaystyle \frac{\displaystyle \partial \bar{I}_{2}}{\displaystyle \partial \underline{C}}}\otimes{\displaystyle \frac{\displaystyle \partial \bar{I}_{1}}{\displaystyle \partial \underline{C}}}+ 2\,{\displaystyle \frac{\displaystyle \partial W^{i}}{\displaystyle \partial \bar{I}_{2}}}{\displaystyle \frac{\displaystyle \partial^{2} \bar{I}_{2}}{\displaystyle \partial \underline{C}^{2}}} \end{aligned} \]

For the Signorini behaviour: \[ \left\{ \begin{aligned} {\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\displaystyle \partial \bar{I}_{1}^{2}}} &= 2\,C_{20} \\ {\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\displaystyle \partial \bar{I}_{2}^{2}}} &= {\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\partial \bar{I}_{1}\partial \bar{I}_{2}}} = {\displaystyle \frac{\displaystyle \partial^{2} W^{i}}{\partial \bar{I}_{2}\partial \bar{I}_{1}}}=0 \end{aligned} \right. \]

const auto d2Pi_dI1b2 = 2*C20;
const auto dSi_dC = 2*(d2Pi_dI1b2*(dI1b_dC^dI1b_dC)+dPi_dI1b*d2I1b_dC2+
                                                   +dPi_dI2b*d2I2b_dC2);

Tests

The comparison of previous derivation of the consistent tangent operator to a numerical approximation is made in the Hyperelasticity.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^{-7}\,C10\).

References