In fracture mechanics, it is classical to replace the plastic behaviour by an equivalent non linear elastic behaviour.

This document describes the implementation of non linear elastic behaviour of the Ramberg-Osgood type using the Default domain specific language of MFront.

The full implementation is available here.

Note

This implementation is compatible with MFront-3.3 and uses unicode characters (see this page for details).

Description of the behaviour

The behaviour is based on the following relationship between the strain \(\underline{\varepsilon}^{\mathrm{to}}\) and the stress \(\underline{\sigma}\):

\[ \underline{\varepsilon}^{\mathrm{to}}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle 3\,K}}}\,{\mathrm{tr}{\left(\underline{\sigma}\right)}}\,\underline{I}+ {{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}}{\displaystyle 3\,\mu}}}\,\underline{n}+ \alpha\,{\left({{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}}{\displaystyle \sigma_{0}}}}\right)}^{n}\,\underline{n} \qquad{(1)}\]

where:

• \(K\) is the bulk modulus.
• \(\mu\) is the shear modulus.
• \(\sigma_{0}\), \(\alpha\) and \(n\) are material properties describing the pseudo-plastic part of the behaviour
• \(\underline{s} = \underline{\sigma}-{{\displaystyle \frac{\displaystyle {\mathrm{tr}{\left(\underline{\sigma}\right)}}}{\displaystyle 3}}}\,\underline{I}\) is the deviatoric part of the stress tensor.
• \(\sigma_{\mathrm{eq}}\) is the von Mises equivalent stress.
• \(\underline{n}={{\displaystyle \frac{\displaystyle 3\,\underline{s}}{\displaystyle 2\,\sigma_{\mathrm{eq}}}}}\) is the von Mises normal.

Algorithm

Computation of the stress

Knowing the strain \(\underline{\varepsilon}^{\mathrm{to}}\), integration of the behaviour, which boils down to the computation of the stress \(\underline{\sigma}\), requires the inversion of Equation (1).

The projection of this equation on the volumetric axis yields:

\[ {\mathrm{tr}{\left(\underline{\sigma}\right)}}=K\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}} \]

Taking the deviatoric part of Equation (1) yields:

\[ \underline{s}^{e}= {\left({{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}}{\displaystyle 3\,\mu}}}+ \alpha\,{\left({{\displaystyle \frac{\displaystyle \sigma_{\mathrm{eq}}}{\displaystyle \sigma_{0}}}}\right)}^{n}\right)}\,\underline{n} = f{\left(\sigma_{\mathrm{eq}}\right)}\,\underline{n} \qquad{(2)}\]

Contracting each side of his equations with itself leads to:

\[ \underline{s}^{e}\,\colon\,\underline{s}^{e} = f^{2}{\left(\sigma_{\mathrm{eq}}\right)}\, \underline{n}\,\colon\,\underline{n} = {{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\,f^{2}{\left(\sigma_{\mathrm{eq}}\right)} \]

where we used the classical result: \(\underline{n}\,\colon\,\underline{n}={{\displaystyle \frac{\displaystyle 3}{\displaystyle 2}}}\).

Finally, \(\sigma_{\mathrm{eq}}\) satisfies:

\[ f{\left(\sigma_{\mathrm{eq}}\right)}=\varepsilon_{\mathrm{eq}} \qquad{(3)}\]

with \(\varepsilon_{\mathrm{eq}}=\sqrt{{{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,\underline{s}^{e}\,\colon\,\underline{s}^{e}}\).

The scalar Equation (3) allows resolving iteratively \(\sigma_{\mathrm{eq}}\) using a standard Newton-Raphson algorithm. Let \(\sigma_{\mathrm{eq}}^{(i)}\) be the current estimate of \(\sigma_{\mathrm{eq}}\), the next estimate is given by:

\[ \sigma_{\mathrm{eq}}^{(i+1)}=\sigma_{\mathrm{eq}}^{(i)}-{{\displaystyle \frac{\displaystyle f{\left(\sigma_{\mathrm{eq}}^{(i)}\right)}}{\displaystyle {\displaystyle \frac{\displaystyle \mathrm{d} f}{\displaystyle \mathrm{d} \sigma_{\mathrm{eq}}}}{\left(\sigma_{\mathrm{eq}}^{(i)}\right)}}}} \]

The initial estimation is choosen as follows:

\[ \sigma_{\mathrm{eq}}^{(0)}=\sigma_{0}\,\sqrt[n]{{{\displaystyle \frac{\displaystyle \varepsilon_{\mathrm{eq}}}{\displaystyle \alpha}}}} \]

Once \(\sigma_{\mathrm{eq}}\) is known, Equations (2) and (3) allow computing the normal \(\underline{n}\):

\[ \underline{n}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle f{\left(\sigma_{\mathrm{eq}}\right)}}}}\,\underline{s}^{e}= {{\displaystyle \frac{\displaystyle 1}{\displaystyle \varepsilon_{\mathrm{eq}}}}}\,\underline{s}^{e} \]

The following relation allows computing the stress:

\[ \underline{\sigma} =K\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+\underline{s} =K\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+{{\displaystyle \frac{\displaystyle 2\,\sigma_{\mathrm{eq}}}{\displaystyle 3}}}\underline{n} =K\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+{{\displaystyle \frac{\displaystyle 2\,\sigma_{\mathrm{eq}}}{\displaystyle 3\,\varepsilon_{\mathrm{eq}}}}}\,\underline{s}^{e} \]

Finally: \[ \underline{\sigma}=K\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{to}}\right)}}\,\underline{I}+\sigma_{\mathrm{eq}}\,\underline{n}^{e} \qquad{(4)}\] with \(\underline{n}^{e}={{\displaystyle \frac{\displaystyle 2}{\displaystyle 3\,\varepsilon_{\mathrm{eq}}}}}\,\underline{s}^{e}\).

Computation of the tangent operator

The following relations will be useful: \[ \begin{aligned} {\displaystyle \frac{\displaystyle \mathrm{d} \varepsilon_{\mathrm{eq}}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}}&=\underline{n}^{e}\\ {\displaystyle \frac{\displaystyle \mathrm{d} \underline{n}^{e}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}}&= {{\displaystyle \frac{\displaystyle 1}{\displaystyle \varepsilon_{\mathrm{eq}}}}}{\left({{\displaystyle \frac{\displaystyle 2}{\displaystyle 3}}}\,\underline{\underline{\mathbf{P}}}-\underline{n}^{e}\,\otimes\,\underline{n}^{e}\right)}\\ \end{aligned} \] with \(\underline{\underline{\mathbf{P}}}=\underline{\underline{\mathbf{I}}}-{{\displaystyle \frac{\displaystyle 1}{\displaystyle 3}}}\,\underline{I}\,\otimes\,\underline{I}\)

Equation (4) allows computing the tangent operator as follows:

\[ {\displaystyle \frac{\displaystyle \mathrm{d} \underline{\sigma}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}}= K\,\underline{I}\,\otimes\,\underline{I}+ {\displaystyle \frac{\displaystyle \mathrm{d} \sigma_{\mathrm{eq}}}{\displaystyle \mathrm{d} \varepsilon_{\mathrm{eq}}}}\,\underline{n}^{e}\,\otimes\,\underline{n}^{e}+ \sigma_{\mathrm{eq}}\,{\displaystyle \frac{\displaystyle \mathrm{d} \underline{n}^{e}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}} \qquad{(5)}\]

\({\displaystyle \frac{\displaystyle \mathrm{d} \sigma_{\mathrm{eq}}}{\displaystyle \mathrm{d} \varepsilon_{\mathrm{eq}}}}\) is obtained by derivation of Equation (3) as follows:

\[ f{\left(\sigma_{\mathrm{eq}}{\left(\varepsilon_{\mathrm{eq}}\right)}\right)}=\varepsilon_{\mathrm{eq}}\Rightarrow {\displaystyle \frac{\displaystyle \mathrm{d} \sigma_{\mathrm{eq}}}{\displaystyle \mathrm{d} \varepsilon_{\mathrm{eq}}}}={{\displaystyle \frac{\displaystyle 1}{\displaystyle {\displaystyle \frac{\displaystyle \partial f}{\displaystyle \partial \sigma_{\mathrm{eq}}}}}}} \qquad{(6)}\]

Implementation

Choice of domain specific language

While not mandatory (the @DSL keyword can be place anywhere in the file), its is convenient to start the implementation by declaring the domain specific language to be used. The Default DSL, which does not set any particular integration scheme, is the most appropriate here:

@DSL Default;

Name of the behaviour

The @Behaviour keyword is used to give the name of the behaviour.

@Behaviour RambergOsgoodNonLinearElasticity;

Metadata

The following instructions give some information about the author, the date and some description of the behaviour.

@Author L. Gelebart/V. Duc Nguyen;
@Description {
  "Ramberg-Osgood model (non-(linear elastic formulation)"
  "Eps = S:Sig +β.σ_eq^n. N"
  "with:"
  " - N = (3/2) Sig_dev / σ_eq"
  " - σ_eq = rac((3/2)σ_dev:σ_dev))"
  " - β = α⋅σ₀/E"
}

Material properties

This behaviour is meant to be generic, which means that the specific characteristics of the material are to be given by the calling solver.

Here, we declare \(5\) material properties:

• Two properties associated with the elastic part of the behaviour: the Young modulus \(E\) and the Poisson ratio \(\nu\). Although this set of elastic properties is not the most convenient for this specific case (the bulk modulus and the shear modulus are more appropriate), it is compatible with the requirement of the Cast3M solver.
• Three properties associated with the plastic part of the behaviour: \(n\), \(\alpha\) and the \(\sigma_{0}\).

@MaterialProperty stress E;
E.setGlossaryName("YoungModulus");
@MaterialProperty real ν;
ν.setGlossaryName("PoissonRatio");
@MaterialProperty real n;
@MaterialProperty real α;
α.setEntryName("alpha");
@MaterialProperty stress σ₀;
σ₀.setGlossaryName("YieldStrength");

Here, we use the subset of UTF-8 characters supported by MFront (See this page for details) to define the variables. Note that it is important to associate so called external names to those variables using one of the setEntryName and setGlossaryName methods: it is highly recommended that the external names only contains ASCII characters.

A numerical parameter

The following numerical parameter will be used as threshold to avoid division by zero.

//! Numerical threshold 
@Parameter real e_ε = 1e-12;
e_ε.setEntryName("NumericalThreshold");

Local variables

The following local variables are used:

//! the shear modulus
@LocalVariable stress μ;
//! the bulk modulus
@LocalVariable stress K;
//! an helper variable
@LocalVariable real β;
//! von Mises equivalent stress
@LocalVariable stress σₑ;
/*!
 * derivative of the von Mises equivalent stress
 * with respect to the equivalent strain
 */
@LocalVariable stress ∂σₑ∕∂εₑ;
//! the equivalent strain
@LocalVariable strain εₑ;
//! inverse of the equivalent strain
@LocalVariable real iεₑ;
//! normal associated with the deviatoric part of the strain
@LocalVariable StrainStensor ne;

Local variables are available in each code blocks.

Initialisation of local variables

The initialisation of local variables is devoted to the @InitLocalVariables.

@InitLocalVariables {
  μ = computeMu(E, ν);
  K = E / (3 ⋅ (1 - 2 ⋅ ν));
  β = α ⋅ σ₀ / E;
  σₑ = stress{};
  ∂σₑ∕∂εₑ = stress{};
}

Note that we use here the symbol ⋅ in place of * as the multiplication sign. Both are equivalent.

The statement σₑ = stress{}; initializes σₑ to a default value (which is 0).

Prediction operator

The prediction operator can used by code_aster to initialize the algorithm used to find the equilibrium. Here we only implement an elastic prediction operator.

@PredictionOperator{
  const auto P = I₄ - (I₂ ⊗ I₂) / 3;
  Dt = K ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ P;
}

where:

• I₄ is the fourth order identity tensor (of type st2tost2, i.e. a fourth order tensor which transforms a symmetric tensor to a a symmetric tensor).
• I₂ is the second order identity symmetric tensor
• ⊗ is the tensorial product. It is equivalent to the ^ characters.

The behaviour integration

The behaviour integration is @Integrator.

@Integrator{

The integrator part is divided in three parts:

• The first part computes some useful variables.
• The second part deals with the computation of the von Mises stress σₑ and its derivative ∂σₑ∕∂εₑ.
• The last part computes the stress tensor.

First part of the computation

The first part computes the strain at the end of the time step, its deviatoric part, the equivalent strain and the normal.

  const auto ε = εᵗᵒ + Δεᵗᵒ;
  const auto se = deviator(ε);
  εₑ = sqrt(2 ⋅ (se | se) / 3);
  iεₑ = 1 / max(εₑ, e_ε);
  ne = 2 ⋅ se ⋅ (iεₑ / 3);

Computation of the von Mises stress and its derivative

The following code computes the von Mises stress and its derivative:

  if (εₑ < e_ε) {
    σₑ = 3 ⋅ μ ⋅ εₑ;
    ∂σₑ∕∂εₑ = 3 ⋅ μ;
  } else {
    const auto e_σ = E ⋅ e_ε;
    const auto i3μ = 1 / (3 ⋅ μ);
    auto f = real{};
    auto fidf = [&]() {
      const auto rσₑⁿ = pow(σₑ / σ₀, n);
      const auto iσₑ = 1 / max(e_σ, σₑ);
      f = σₑ ⋅ i3μ + β ⋅ rσₑⁿ - εₑ;
      ∂σₑ∕∂εₑ = 1 / max(i3μ + n ⋅ β ⋅ rσₑⁿ ⋅ iσₑ, i3μ ⋅ e_ε);
    };
    σₑ = σ₀ ⋅ pow(εₑ / β, 1 / n);
    auto iter = int{};
    fidf();
    while (abs(f) > e_ε) {
      fidf();
      σₑ -= f ⋅ ∂σₑ∕∂εₑ;
      if (++iter > 20) {
        throw(DivergenceException());
      }
    }
  }

If the equivalent strain is below the numerical threshold, a solution corresponding to an elastic behaviour is returned.

In the other case, a local Newton algorithm is used to determine σₑ and ∂σₑ∕∂εₑ, as described by Section 2.1. Except for the usage of a lambda function called ifdf, the code is quite easy to read. A lambda function can be seen as a local function which allows to factorise the code computing the residual and its derivative. The syntax [&] means that the lambda function of has access to all the variables of the surrounding scope by references. Note that we used Equation (6) to identify \({{\displaystyle \frac{\displaystyle 1}{\displaystyle {\displaystyle \frac{\displaystyle \mathrm{d} f}{\displaystyle \mathrm{d} \sigma_{\mathrm{eq}}}}}}}\) to ∂σₑ∕∂εₑ.

Computation of the stress

The computation of the stress is a direct translation of Equation (4):

  σ = K ⋅ trace(ε) ⋅ I₂ + σₑ ⋅ ne;

Finalisation

The @Integrator code block ends with a final curly brace.

} // end of @Integrator

Computation of the tangent operator

The following block compute the tangent operator using Equation (5):

@TangentOperator{
  const auto P = I₄ - (I₂ ⊗ I₂) / 3;
  if ((smt == ELASTIC) || (smt == SECANTOPERATOR)) {
    Dt = K ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ P;
  } else {
    if (εₑ < e_ε) {
      Dt = K ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ P;
    } else {
      Dt = K ⋅ (I₂ ⊗ I₂) + ∂σₑ∕∂εₑ ⋅ (ne ⊗ ne) +
           σₑ ⋅ iεₑ ⋅ (2 ⋅ P / 3 - (ne ⊗ ne));
    }
  }
}

Examples of usage

Uniaxial tensile test