This document describes how to build a behaviour of concrete describing the damage evolution following Fichant-La Borderie and the load induced thermal strain (LITS) following Torelli et al (See [1–4] for details).

We strongly advice readers to refer to the following documents before continuing:

• The implementation of the Fichant-La Borderie damage behaviour is detailled here.
• The implementation of the the load induced thermal strain is detailled here.

Description of the coupled behaviour

The behaviour is described by a standard split of the strain \(\underline{\varepsilon}^{\mathrm{to}}\) in an elastic and a viscoplastic parts, respectively denoted \(\underline{\varepsilon}^{\mathrm{el}}\) and \(\underline{\varepsilon}^{\mathrm{LITS}}\):

\[ \underline{\varepsilon}^{\mathrm{to}}=\underline{\varepsilon}^{\mathrm{el}}+\underline{\varepsilon}^{\mathrm{LITS}} \qquad{(1)}\]

Damage behaviour

The stress \(\underline{\sigma}\) is related to the the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) by the Fichant-La Borderie damage behaviour:

\[ \underline{\sigma}= f_{FLB}{\left(d,\underline{\varepsilon}^{\mathrm{el}}\right)} \qquad{(2)}\]

where \(d\) is a scalar damage variable. See this page for details.

Damage evolution

The damage evolution is a function of the equivalent elastic strain \(\varepsilon^{\mathrm{el}}_{\mathrm{eq}}\), which will be defined hereafter:

• 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{el}}_{\mathrm{eq}}}}}\,\exp{\left(B_{t}\,{\left(\epsilon_{0}-\varepsilon^{\mathrm{el}}_{\mathrm{eq}}\right)}\right)} \qquad{(3)}\]

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{el}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}}}}\,\exp{\left(B_{t}\,{\left(\epsilon_{0}-{\left.\varepsilon^{\mathrm{el}}_{\mathrm{eq}}\right|_{t+\Delta\,t}}\right)}\right)}\right)} \]

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

Equivalent strain

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

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

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

Computation of the stress

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>_{-}} \]

where 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{el}}\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}}} \]

Load induced thermal strain

The evolution of the load induced thermal strain \(\underline{\varepsilon}^{\mathrm{LITS}}\) is given by:

\[ \underline{\dot{\varepsilon}}^{\mathrm{LITS}}= {{\displaystyle \frac{\displaystyle \eta{\left(\underline{\sigma}\right)}\,\beta{\left(T\right)}}{\displaystyle \sigma_{u0}}}} \left[ -\nu_{\mathrm{LITS}}{\mathrm{tr}{\left({\left<\underline{\sigma}\right>_{-}}\right)}}\,\underline{I} +{\left(1-\nu_{\mathrm{LITS}}\right)}\,{\left<\underline{\sigma}\right>_{-}} \right]\,\dot{T}_{e} \qquad{(4)}\]

where:

• \({\left<\sigma\right>_{-}}\) is the negative part of the stress tensor \(\underline{\sigma}\). This negative part is computed as follows: \[ {\left<\sigma\right>_{-}} = \sum_{I=1}^{3}{\left<\sigma_{I}\right>_{-}}\,\underline{n}_{I} \] where \(\sigma_{I}\) is the Ith eigenvalue of \(\underline{\sigma}\) and \(\underline{n}_{I}\) the associated eigen tensor.
• \(\eta{\left(\underline{\sigma}\right)} = 1+C_{m}{\left(\underline{\sigma}\right)}\,\gamma\). \(\eta\) describes the effect of the stress triaxiality.
• \(C_{m}{\left(\underline{\sigma}\right)}={{\displaystyle \frac{\displaystyle -{\mathrm{tr}{\left({\left<\underline{\sigma}\right>_{-}}\right)}}}{\displaystyle \sqrt{{\left<\underline{\sigma}\right>_{-}}\,\colon\,{\left<\underline{\sigma}\right>_{-}}}}}}\). \(C_{m}\) is the triaxiality index linked to the principal stresses.
• \(\sigma_{u0}\), \(\nu_{\mathrm{LITS}}\) and \(\gamma\) are material parameters. More precisely \(\sigma_{u0}\) is the compressive strength of the material and \(\gamma\) is a material parameter calibrated to experimental data for the appropriate temperature range.
• \(\beta{\left(T\right)}\) is a four order polynomial of the temperature describing uniaxial temperature-LITS behaviour observed in the literature: \[ \beta{\left(T\right)}=\beta_{0}+\beta_{1}\,T+\beta_{2}\,T^{2}+\beta_{3}\,T^{3}+\beta_{4}\,T^{4} \]

\(\dot{T}_{e}\) is the effective temperature rate defined as follows:

\[ \dot{T}_{e}= \left\{ \begin{aligned} \dot{T} &\quad\text{if}\quad T=T_{\mathrm{max}}\\ 0 &\quad\text{if}\quad T<T_{\mathrm{max}}\\ \end{aligned} \right. \]

where \(T_{\mathrm{max}}\) is the maximal temperature seen by the material over time.

Implicit algorithm

Following the implementation of Torelli’ LITS behaviour (see here for details), the implicit scheme can be reduced to one non linear equation whose only unknown is the elastic strain increment:

\[ f_{\underline{\varepsilon}^{\mathrm{el}}}=\Delta\,\underline{\varepsilon}^{\mathrm{el}}+\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}{\left(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\right)}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}=\underline{0} \qquad{(5)}\]

Computing the jacobian \({\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is straightforward:

\[ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}=\underline{\underline{\mathbf{I}}}+{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \]

The derivative \({\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\) has been computed on this page while the derivative \({\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is computed on this page.

Implementation

Choice of the domain specific language

The domain specific language suitable for implementing an implicit scheme is called Implicit. The @DSL keyword allows specifying the domain specific language to be used:

@DSL Implicit;

Choice of the behaviour name

The @Behaviour allows giving a name the behaviour:

@Behaviour FichantLaBorderieDamageAndTorelliLoadInducedThermalStrainBehaviour;

Defining the name of behaviour is required.

Metadata

The @Author, @Date and @Description keywords allows adding metadata to the behaviour (respectively the name of the author of the implementation, the date at which the implementation was written, a short description), as follows:

@Author Thomas Helfer;
@Date 15 / 12 / 2019;
@Description {
  "Implementation of a behaviour combining Fichant-La Borderie' "
  "description of damage and Torelli'LITS"
}

Numerical parameters

We specify two numerical parameters impacting the implicit scheme:

• the \(\theta\) parameter (default value is \(1/2\)). In pratice, we set to \(1\) to indicate that the resolution is purely implicit, but we won’t use this value in the implementation.
• the stopping criteria (default value is \(10^{-8}\)). This may be sufficient in most cases, but not for the computation of the consistent tangent operator.

@Epsilon 1.e-14;
@Theta 1;

Thermoelastic properties, computation of the thermal strain

The @ElasticMaterialProperties keyword allows defining the material properties. When those are defined, the Implicit DSL automatically defines the first and second Lamé coefficients in local variables called:

• lambda and mu at time \(t+\theta\,\Delta\,t\)
• lambda_tdt and mu_tdt at time \(t+\Delta\,t\)

@ElasticMaterialProperties{47000, 0.25};

The @ComputeThermalExpansion keyword allows computing the thermal strain. It is followed by the mean thermal expansion coefficient.

@ComputeThermalExpansion 10e-6;

// parameters of the Fichant-La Borderie part
@Parameter Bt = 3690.070983;
@Parameter e0 = 1.03e-04;
@Parameter a = 3;

// parameters of the LITS part
@Parameter gamma = 1.5;
@Parameter sigmultimate = 50.;
@Parameter nulits = 0.50;
@Parameter tcrit = 0.;
@Parameter b[5] = {2.7031065533E-05, -1.0209170592E-06, 6.1200423753E-9,  //
                   -1.2632648735E-11, 6.9158539621E-15};

@AuxiliaryStateVariable real d;
d.setGlossaryName("Damage");
@AuxiliaryStateVariable StrainStensor elits;
elits.setEntryName("LoadInducedThermaStrain");
@AuxiliaryStateVariable temperature Tmax;
Tmax.setEntryName("MaximalValueOfTheTemperature");

//! LITS increment
@LocalVariable StrainStensor delits;
//! Creep coefficient
@LocalVariable real C;

@LocalVariable StiffnessTensor dsig_ddeel;
@LocalVariable real d_p;
@LocalVariable StiffnessTensor De;

@InitLocalVariables {
  Tmax = max(max(tcrit, Tmax), T);
  const auto T_ = T + theta * dT;
  const auto beta = b[0] + T_ * (b[1] + T_ * (b[2] + T_ * (b[3] + T_ * b[4])));
  const auto dTe = max(T + dT - max(tcrit, Tmax), temperature(0));
  C = (beta / (-sigmultimate)) * dTe;
  De = lambda_tdt * Stensor4::IxI() + 2 * mu_tdt * Stensor4::Id();
} // end of @InitLocalVariables

@Integrator {
  constexpr const auto id = Stensor::Id();
  constexpr const auto esolver = StressStensor::FSESJACOBIEIGENSOLVER;
  constexpr const stress eeps = 1.e-12;
  const stress seps = 1.e-12 * young;

  // positive part
  const auto pp = [](const real x) { return x > 0 ? x : 0; };
  // derivative of the positive part
  const auto dpp = [&eeps](const real x) { return std::abs(x) < eeps ? 0.5 : ((x < 0) ? 0 : 1); };
  // square of the posititve part
  auto square_ppos = [](const strain& v) { return v > 0 ? v * v : 0; };
  // elastic strain at the midle of the time step
  const auto e = eval(eel + deel);
  // eigen values and eigen tensors of the elastic strain
  auto e_vp = tvector<3u, strain>{};
  auto m = tmatrix<3u, 3u, strain>{};
  e.template computeEigenVectors<esolver>(e_vp, m);
  // update the damage
  const auto e_eq = sqrt(square_ppos(e_vp[0]) + square_ppos(e_vp[1]) + square_ppos(e_vp[2]));
  // effective stress at t+dt
  const auto Cd = (e0 / e_eq) * exp(Bt * (e0 - e_eq));
  d_p = (e_eq > e0) ? 1 - Cd : 0;
  const auto bp = d_p > d;
  const auto de = bp ? d_p : d;
  // derivative with respect to the damage
  const auto dde_ddeel = [&]() -> Stensor {
    if (!bp) {
      return Stensor(real(0));
    }
    // positive part of the total strain
    const auto ep = StrainStensor::computeIsotropicFunction(pp, e_vp, m);
    // derivative of the damage
    const auto dde_deq = Cd * (Bt + 1 / e_eq);
    const auto dep_ddeel =
        StrainStensor::computeIsotropicFunctionDerivative(pp, dpp, e_vp, m, eeps * 0.1);
    const auto deq_dep = ep / e_eq;
    return dde_deq * deq_dep * dep_ddeel;
  }();
  // function of the damage to simplify expressions
  const auto de_a = pow(de, a);
  const auto fpd = (1 - de);
  const auto fpn = (1 - de_a);
  // effective stress at the end of the time step
  const auto l_tr_e = lambda_tdt * trace(e);
  const auto s = eval(l_tr_e * id + 2 * mu_tdt * e);
  const auto s_vp = tvector<3u>{l_tr_e + 2 * mu_tdt * e_vp[0],  //
                                l_tr_e + 2 * mu_tdt * e_vp[1],  //
                                l_tr_e + 2 * mu_tdt * e_vp[2]};
  const auto sp = StressStensor::computeIsotropicFunction(pp, s_vp, m);
  const auto dsp = StressStensor::computeIsotropicFunctionDerivative(pp, dpp, s_vp, m, seps * 0.1);
  const auto sn = eval(s - sp);
  const auto dsn = eval(Stensor4::Id() - dsp);
  // final stress
  sig = fpd * sp + fpn * sn;
  // derivative of the stress
  dsig_ddeel = ((fpd - fpn) * dsp + fpn * Stensor4::Id()) * De;
  if (bp) {
    const auto ide = 1 / max(eeps, de);
    const auto dfpd_dd = -1;
    const auto dfpn_dd = -a * de_a * ide;
    dsig_ddeel += dfpd_dd * (sp ^ dde_ddeel) + dfpn_dd * (sn ^ dde_ddeel);
  }

  // LITS part
  const auto sn_eq = sqrt(sn | sn);
  const auto isn_eq = 1 / max(seps, sn_eq);
  const auto cm = -trace(sn) * isn_eq;
  const auto dcm_dsig = eval((-isn_eq * id + trace(sn) * power<3>(isn_eq) * sn) | dsn);
  const auto eta = 1 + (cm - 1) * gamma;
  const auto se = eval((1 + nulits) * sn - nulits * trace(sn) * id);
  delits = C * eta * se;
  const auto deta_dsig = gamma * dcm_dsig;
  const auto dse_dsig = (1 + nulits) * dsn - nulits * ((id ^ id) * dsn);
  const auto ddelits_dsig = C * (se ^ deta_dsig) + C * eta * dse_dsig;

  // elasticity
  feel += delits;
  dfeel_ddeel += ddelits_dsig * dsig_ddeel;
} // end of @Integrator

@UpdateAuxiliaryStateVariables {
  d = max(d, d_p);
  elits += delits;
}

@TangentOperator {
  if (smt == ELASTIC) {
    Dt = De;
  } else if (smt==CONSISTENTTANGENTOPERATOR){
    Stensor4 ddeel_ddeto;
    getPartialJacobianInvert(ddeel_ddeto);
    Dt = dsig_ddeel * ddeel_ddeto;
  }
}

References