This document describes how to implement the Fichant-La Borderie damage behaviour. See [1, 2] for a detailed description.

The Fichant-La Borderie (FLB) model is an extension of the Mazars model, which:

simplifies the number of input parameters defining the tensile and compressive asymmetry.
handles unilateral effects (i.e., crack closure).
features an energy based regularisation method, thus allowing a fracture finite energy in a finite element implementation

The Fichant La Borderie damage behaviour under uniaxial tension and compression

The implementation is available here

Description of the Fichant-La Borderie damage law

Equivalent strain

The equivalent strain \(\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\) is defined as:

\[ \varepsilon^{\mathrm{to}}_{\mathrm{eq}}=\sqrt{{\left<\underline{\varepsilon}^{\mathrm{to}}\right>_{+}}\,\colon\,{\left<\underline{\varepsilon}^{\mathrm{to}}\right>_{+}}}= \sqrt{\sum_{I=1}^{3}{\left<\varepsilon^{\mathrm{to}}_{I}\right>_{+}}^2} \]

where \({\left<\underline{\varepsilon}^{\mathrm{to}}\right>_{+}}\) is the positive part of the total strain and \(\left(\varepsilon^{\mathrm{to}}_{I}\right)_{i \in [1:3]}\) are its eigenvalues.

Derivative of the equivalent strain

The following expression of the derivative of the equivalent strain with respect to the total strain will be useful:

\[ {\displaystyle \frac{\displaystyle \mathrm{d} \varepsilon^{\mathrm{to}}_{\mathrm{eq}}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}}={{\displaystyle \frac{\displaystyle 1}{\displaystyle \varepsilon^{\mathrm{to}}_{\mathrm{eq}}}}}\,{\left<\underline{\varepsilon}^{\mathrm{to}}\right>_{+}}\,\colon\,{\displaystyle \frac{\displaystyle \mathrm{d} {\left<\underline{\varepsilon}^{\mathrm{to}}\right>_{+}}}{\displaystyle \mathrm{d} \underline{\varepsilon}^{\mathrm{to}}}} \qquad{(1)}\]

Damage evolution

The damage evolution is a function of the equivalent strain \(\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\):

Under a given threshold \(\epsilon_{0}\), no damage occurs.
Once this threshold is reached, the damage \(d\) is defined as follows: \[ d=1-{{\displaystyle \frac{\displaystyle \epsilon_{0}}{\displaystyle \varepsilon^{\mathrm{to}}_{\mathrm{eq}}}}}\,\exp{\left(B_{t}\,{\left(\epsilon_{0}-\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right)}\right)} \qquad{(2)}\]

However, the latter expression must be modified to take into account the irreversibility of the damage evolution. Let \({\left.d\right|_{t}}\) be the value of the damage at the beginning of the time step and \({\left.d\right|_{t+\Delta\,t}}\) its value at the end of the time step, \({\left.d\right|_{t+\Delta\,t}}\) is determined as follows:

\[ {\left.d\right|_{t+\Delta\,t}}=\max{\left({\left.d\right|_{t}},1-{{\displaystyle \frac{\displaystyle \epsilon_{0}}{\displaystyle {\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}}\,\exp{\left(B_{t}\,{\left(\epsilon_{0}-{\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}\right)}\right)}\right)} \]

where \({\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}\) is the value of the equivalent strain at the end of the time step.

Hillerborg regularisation (see [3])

Local damage models suffers from spurious mesh dependency. In particular, the dissipated energy tends to \(0\) as the mesh size decreases.

The Hillerborg regularisation consists in introducing the mesh size in the material parameter to get a constant dissipated energy.

In the FLB case, \(B_{t}\) can be deduced from the fracture energy \(G_{f}\) and the mesh size \(h\) as follows:

\[ B_{t} ={{\displaystyle \frac{\displaystyle h\,E\,\epsilon_{0}}{\displaystyle G_{f} −{{\displaystyle \frac{\displaystyle 1}{\displaystyle 2}}}\,E\,\epsilon_{0}^{2}\,h}}} \]

where \(E\) is the Young modulus

Note that \(B_{t}\) must be positive, which, for a given fracture energy \(G_{f}\), imposes a maximal value for the mesh size \(h\).

\[ h < {{\displaystyle \frac{\displaystyle 2\,G_{f}}{\displaystyle E\,\epsilon_{0}^{2}}}} \]

A condition which states that the elastic energy stored inside an element at the onset of damage must be less than dissipated energy by damage.

Derivatives of the damage

In case of a growing damage, the derivative of the damage with respect to the equivalent plastic strain is:

\[ {\displaystyle \frac{\displaystyle \mathrm{d} {\left.d\right|_{t+\Delta\,t}}}{\displaystyle \mathrm{d} {\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}={{\displaystyle \frac{\displaystyle \epsilon_{0}}{\displaystyle {\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}}\exp{\left(B_{t}\,{\left(\epsilon_{0}-{\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}\right)}\right)}\left(B_{t}+{{\displaystyle \frac{\displaystyle 1}{\displaystyle {\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}}\right) \qquad{(3)}\]

Combined with Equation (1), Equation (3) allows computing the derivative of the damage with respect to the total strain.

Computation of the stress

The effective stress tensor \(\underline{\sigma}^{\mathrm{eff}}\) is computed using the standard Hooke law: \[ {\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}={\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\,\colon\,{\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}} \]

The behaviour is assumed isotropic, so that the stiffness tensor \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\) can be related to the first and second Lamé coefficients, denoted respectively \(\lambda\) and \(\mu\), as follows:

\[ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}={\left.\lambda\right|_{t+\Delta\,t}}\,\underline{I}\,\otimes\,\underline{I}+2\,{\left.\mu\right|_{t+\Delta\,t}}\,\underline{\underline{\mathbf{I}}} \]

Note

The restriction to isotropy has an important pratical consequence: \(\underline{\sigma}^{\mathrm{eff}}\) and \(\underline{\varepsilon}^{\mathrm{to}}\) have the same egein basis.

The stress tensor \(\underline{\sigma}\) is then computed as follows: \[ {\left.\underline{\sigma}\right|_{t+\Delta\,t}}={\left(1-{\left.d\right|_{t+\Delta\,t}}\right)}\,{\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{+}}+{\left(1-{\left.d\right|_{t+\Delta\,t}}^{a}\right)}\,{\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{-}} \]

Derivative of the stress with respect to the total strain

In case of damage growth, the derivative of the stress with respect to the total strain is given by: \[ \begin{aligned} {\displaystyle \frac{\displaystyle \mathrm{d} {\left.\underline{\sigma}\right|_{t+\Delta\,t}}}{\displaystyle \mathrm{d} {\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}}}}&= \left[{\left(1-{\left.d\right|_{t+\Delta\,t}}\right)}\, {\displaystyle \frac{\displaystyle \mathrm{d} {\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{+}}}{\displaystyle \mathrm{d} {\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}}}+ {\left(1-{\left.d\right|_{t+\Delta\,t}}^{a}\right)}\,{\displaystyle \frac{\displaystyle \mathrm{d} {\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{-}}}{\displaystyle \mathrm{d} {\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}}}\right]\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\Delta\,t}}\\ &- {\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{+}}\,\otimes\,{\displaystyle \frac{\displaystyle \mathrm{d} {\left.d\right|_{t+\Delta\,t}}}{\displaystyle \mathrm{d} {\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}}}} - a\,{\left.d\right|_{t+\Delta\,t}}^{a-1}\,{\left<{\left.\underline{\sigma}^{\mathrm{eff}}\right|_{t+\Delta\,t}}\right>_{-}}\,\otimes\,{\displaystyle \frac{\displaystyle \mathrm{d} {\left.d\right|_{t+\Delta\,t}}}{\displaystyle \mathrm{d} {\left.\underline{\varepsilon}^{\mathrm{to}}\right|_{t+\Delta\,t}}}} \end{aligned} \qquad{(4)}\]

This expression gives the so-called consistent tangent operator.

Implementation

Choice of the domain specific language

The Fichant-La Borderie damage behaviour does not require any specific integration scheme as the damage evolution is an explicit function of the total strain. The Default domain specific language is suitable for treating that kind of behaviours.

@DSL Default;

Name of the behaviour

The @Behaviour keyword defines the name of the behaviour:

@Behaviour FichantLaBorderieDamageBehaviour;

Metadata

The implementation starts by defining somme metadata associated with the behaviour:

the date.
the names of the authors.
a short description.

@Date 7/12/2019;
@Author A. Gangnant, T. Helfer;
@Description{
  "Implementation of the Fichant-La Borderie damage behaviour"
}

Material properties

Elastic material properties

Most solvers uses the Young modulus and the Poisson ratio to characterize the elasticity of an isotropic material. Here we adopt this convention by defining them as two material properties named young and nu and associate those variables with the appropriate glossary names:

@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");

About glossary names

Associating a glossary name is a way to define the so-called external name of this variable, i.e. the name that will appear on the calling solver side. It also allows MFront to make assumptions about the real meaning of those variables.

The case of the Cast3M solver is a noteworthy example. This solver requires \(4\) elastic properties to be defined for its own use in the case of an isotropic material properties: the Young modulus, the Poisson ratio, the density and the thermal expansion coefficient. Normally, the MFront material properties are appended to those four properties. However, by associating the YoungModulus and the PoissonRatio glossary names to the young and nu material properties respectively, MFront will identify them with the ones required by Cast3M.

For details about the glossary names, the reader may refer to the following page.

Material properties associated with damage

The damage evolution given by Equation (2) requires three material properties \(Bt\), \(e_{0}\) and \(a\).

@MaterialProperty real Bt;
@MaterialProperty strain e0;
@MaterialProperty real a;

State variable

The only state variable of this behaviour is the damage variable \(d\) which we declare as follows:

@StateVariable real d;

For post-processing purposes, we associate this variable with the Damage glossary name.

d.setGlossaryName("Damage");

Local variables

In MFront, a local variable are variables in each code blocks. They are usually evaluated once for all in the @InitLocalVariables code block which is called before the behaviour integration (as defined by the @Integrator code block) or before the computation of the prediction operator (as defined in the @PredictionOperator code block).

In the following, we will use two local variables lambda and mu associated respectively with the first and second Lamé coefficients.

@LocalVariable stress lambda, mu;

Initialization of the local variables

The initialization of the local variables is performed in the @InitLocalVariables code block, as follows:

@InitLocalVariables {
  lambda = computeLambda(young, nu);
  mu = computeMu(young, nu);
}

Here, we use two functions provided by the TFEL/Material library: computeLambda and computeMu.

Prediction operator

The prediction operator is used by the Code_Aster finite element solver at the beginning of each new time step. It is also used by the Abaqus/Explicit solver to retrieve the elastic properties at the packaging step.

@PredictionOperator {
  static_cast<void>(smt);
  Dt = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
}

Here, we only provide the elastic prediction operator. Thus, we discarded the smt variable by casting it to void to avoid a compiler warning.

Behaviour integration

The behaviour integration is performed in the @Integrator code block which is meant to:

update the damage state variable.
compute the stress at the end of the time step.
compute the consistent tangent operator if required.

Note

The computation of the consistent tangent operator is usually performed in the @TangentOperator code block. However, in the case of the Fichant-La Borderie behaviour, this would lead to an unneccessary complexity as a lot of local variables would have been required to store the results of complex computations required in both the @Integrator and the @TangentOperator code blocks.

To state that the consistent tangent operator is indeed computed in the @Integrator code block, one shall use the @ProvidesTangentOperator keyword. Note that this keyword explicitly states that the consistent tangent operator is not symmetric. See the @ProvidesSymmetricTangentOperator keyword to state that the consistent tangent operator is symmetric.

@ProvidesTangentOperator;
@Integrator{
  // a simple alias for the identity
  constexpr const auto id = Stensor::Id();

At the beginning of this block, we define two constant used in the following to decide if two eigenvalues are equal:

The first constant eeps is used when dealing a tensor of the strain type.
The second constant seps is used when dealing a tensor of the stress type.

  constexpr const strain eeps = 1.e-12;
  const auto seps = eeps * young;

Then we define two local functions called pp and square_pp:

pp computes the positive part of a scalar.
square_pp computes the square of the positive part of a scalar.

  // positive part
  const auto pp = [](const real x) { return x > 0 ? x : 0; };
  // square of the posititve part
  auto square_pp = [](const strain v) { return v > 0 ? v * v : 0; };

Note

In C++, those local functions are called lambda functions.

We now compute the total strain and stores it in a variable called e:

  // total strain
  const auto e = eval(eto + deto);

where eto and deto are respectively the total strain at the beginning of the time step and its increment over the time step.

Tensorial operations in the TFEL/Math library are lazy, which means that their evaluation is delayed until it is really needed. This technique is called expression template in C++ and is based on the definition of intermediate objects which hold the operation to be performed.

The eval function is used to force the evaluation of this operation. This is required for the computation of the eigenvalues and the eigenvectors which is done by calling the computeEigenVectors tensor of the stensor class.

In C++-17, this can be called as follows:

  // eigen values and eigen tensors of the total strain
  const auto [e_vp, m] = 
    e.template computeEigenVectors<Stensor::FSESJACOBIEIGENSOLVER>();

In previous C++ version, one must use a more verbose code:

  // eigen values and eigen tensors of the total strain
  auto e_vp = tvector<3u, strain>{};
  auto m = tmatrix<3u, 3u, strain>{};
  e.template computeEigenVectors<Stensor::FSESJACOBIEIGENSOLVER>(e_vp, m);

Here, we must make some important comments:

We use the Jacobi eigen solver which is very accurate but is more expensive than the default eigen solver.
The computation of the equivalent strain only requires the eigenvalues of the total strain. We also compute its eigen vectors because they also are the eigen vectors of the effective stress and its positive part.

The computation of the equivalent strain is then straightforward:

  const auto e_eq = sqrt(square_pp(e_vp[0]) + //
                         square_pp(e_vp[1]) + //
                         square_pp(e_vp[2]));

We are now able to compute the damage at the end of the time step as follows:

  // update the damage, taking irreversibility into account
  const auto Cd = (e0 / e_eq) * exp(Bt * (e0 - e_eq));
  const auto d_p = (e_eq > e0) ? 1 - Cd : 0;
  const auto bp = d_p > d;
  if (bp) {
    d = d_p;
  }

The previous lines also defines two additional variables:

Cd which is used to compute the expression \({{\displaystyle \frac{\displaystyle \epsilon_{0}}{\displaystyle {\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}}\exp{\left(B_{t}\,{\left(\epsilon_{0}-{\left.\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}\right)}\right)}\) only once.
bp which is a boolean stating if the damage increases over the time step. This boolean will be used in the computation of the consistent tangent operator.

The next two lines computes the effective stress:

  // effective stress at the end of the time step
  const auto l_tr_e = lambda * trace(e);
  const auto s = eval(l_tr_e * id + 2 * mu * e);

The next lines computes the positive part of the effective stress:

  // positive part of the effective stress
  const auto s_vp = tvector<3u, real>{l_tr_e + 2 * mu * e_vp[0],  //
                                      l_tr_e + 2 * mu * e_vp[1],  //
                                      l_tr_e + 2 * mu * e_vp[2]};
  const auto sp = StressStensor::computeIsotropicFunction(pp, s_vp, m);

Those lines uses the computeIsotropicFunction which takes the function computing the positive part, the eigenvalues of the effective stress and its eigen vectors.

Finally, the final stress is computed:

  // function of the damage to simplify expressions
  const auto d_a = pow(d, a);
  const auto fpd = (1 - d);
  const auto fpn = (1 - d_a);
  // final stress
  sig = (fpd - fpn) * sp + fpn * s;

Computation of the tangent operator

If one computes the tangent operator in the @Integrator code block, one shall check if this operator was requested using the computeTangentOperator_ boolean value.

  if (computeTangentOperator_) {

Next, we define a local function called dpp which computes the derivative of the positive part:

    // derivative of the positive part
    const auto dpp = [&seps](const stress x) {
      return std::abs(x) < seps ? 0.5 : ((x < 0) ? 0 : 1);
    };

Various kind of tangent operator can be computed. To see which one was requested by the calling solver, one must check the value of the smt variable (smt stands “stiffness matrix type”):

The first case corresponds to the ELASTIC value:

    if (smt == ELASTIC) {
      Dt = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();

The second case corresponds to the SECANTOPERATOR value. This case also corresponds to the consistent tangent operator was no damage increase occurs (i.e. the bp variable’ value is false):

    } else if ((smt == SECANTOPERATOR) || (!bp)) {
      // elastic stiffness
      const auto De = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
      // derivative of the positive part of the stress
      const auto dsp =
          StressStensor::computeIsotropicFunctionDerivative(pp, dpp, s_vp, m, seps * 0.1);
      Dt = ((fpd * dsp + fpn * (Stensor4::Id() - dsp)) + eeps * Stensor4::Id()) * De;

This code uses the computeIsotropicFunctionDerivative which computes \({\displaystyle \frac{\displaystyle \partial {\left<\underline{\sigma}\right>_{+}}}{\displaystyle \partial \underline{\sigma}}}\). We also added the eeps * Stensor4::Id() to avoid singular stiffness matrices and add some coercivity to the global problem.

The last case corresponds to the CONSISTENTANGENTOPERATOR value in case of damage increase

    } else {
      // elastic stiffness
      const auto De = lambda * Stensor4::IxI() + 2 * mu * Stensor4::Id();
      // derivative of the positive part of the stress
      const auto dsp =
          StressStensor::computeIsotropicFunctionDerivative(pp, dpp, s_vp, m, seps * 0.1);
      // derivative with respect to the damage
      const auto dfpd_dd = -1;
      const auto id = 1 / max(eeps, d);
      const auto dfpn_dd = -a * d_a * id;
      // positive part of the total strain
      const auto ep = StrainStensor::computeIsotropicFunction(pp, e_vp, m);
      // derivative of the damage
      const auto dd_deq = Cd * (Bt + 1 / e_eq);
      const auto dep_de =
          StrainStensor::computeIsotropicFunctionDerivative(pp, dpp, e_vp, m, eeps * 0.1);
      const auto deq_dep = ep / e_eq;
      const auto dd_de = dd_deq * deq_dep * dep_de;
      Dt = ((fpd * dsp + fpn * (Stensor4::Id() - dsp)) + eeps * Stensor4::Id()) * De  //
           + (((dfpd_dd - dfpn_dd) * sp + dfpn_dd * s) ^ dd_de);
    }
  }
}

References

Fichant, Stéphanie. Endommagement et anisotropie induite du béton de structures : Modélisations approchées. thesis. Cachan, Ecole normale supérieure, 1996. Available from: http://www.theses.fr/1996DENS0015

Gangnant, Alexandre. Étude de la rupture quasi-fragile d’un béton à l’échelle mésoscopique : Aspects expérimentaux et modélisation. thesis. Bordeaux, 2016. Available from: http://www.theses.fr/2016BORD0401

Hillerborg, A., Modéer, M. and Perterson, P.-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research. 1976. Vol. 6, p. 779–782.

How to implement the Fichant-La Borderie damage behaviour

Thomas Helfer, Alexandre Gagnant

12/2019