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:
The implementation is available here
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)\]
The damage evolution is a function of the equivalent strain \(\varepsilon^{\mathrm{to}}_{\mathrm{eq}}\):
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.
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.
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;
The @Behaviour
keyword defines the name of the behaviour:
@Behaviour FichantLaBorderieDamageBehaviour;
The implementation starts by defining somme metadata associated with the behaviour:
@Date 7/12/2019;
@Author A. Gangnant, T. Helfer;
@Description{
"Implementation of the Fichant-La Borderie damage behaviour"
}
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;
"YoungModulus");
young.setGlossaryName(@MaterialProperty real nu;
"PoissonRatio"); nu.setGlossaryName(
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, theMFront
material properties are appended to those four properties. However, by associating theYoungModulus
and thePoissonRatio
glossary names to theyoung
andnu
material properties respectively,MFront
will identify them with the ones required byCast3M
.For details about the glossary names, the reader may refer to the following page.
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;
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.
"Damage"); d.setGlossaryName(
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;
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
.
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);
2 * mu * Stensor4::Id();
Dt = lambda * Stensor4::IxI() + }
Here, we only provide the elastic prediction operator. Thus, we discarded the smt
variable by casting it to void
to avoid a compiler warning.
The behaviour integration is performed in the @Integrator
code block which is meant to:
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:
eeps
is used when dealing a tensor of the strain type.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 calledlambda
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] =
template computeEigenVectors<Stensor::FSESJACOBIEIGENSOLVER>(); e.
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>{};
template computeEigenVectors<Stensor::FSESJACOBIEIGENSOLVER>(e_vp, m); e.
Here, we must make some important comments:
The computation of the equivalent strain is then straightforward:
const auto e_eq = sqrt(square_pp(e_vp[0]) + //
1]) + //
square_pp(e_vp[2])); square_pp(e_vp[
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], //
2 * mu * e_vp[1], //
l_tr_e + 2 * mu * e_vp[2]};
l_tr_e + 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;
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) {
2 * mu * Stensor4::Id(); Dt = lambda * Stensor4::IxI() +
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 =
s_vp, m, seps * 0.1);
StressStensor::computeIsotropicFunctionDerivative(pp, dpp, 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 =
s_vp, m, seps * 0.1);
StressStensor::computeIsotropicFunctionDerivative(pp, dpp, // 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 =
0.1);
StrainStensor::computeIsotropicFunctionDerivative(pp, dpp, e_vp, m, eeps * 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);
}
} }