StandardElasticity
brickThis page shows how to implement a behaviour describing the load induced thermal strain (LITS) phenomena in concrete following the work of [1].
This implementation is available here.
The behaviour is described by a standard split of the strain \(\underline{\varepsilon}^{\mathrm{to}}\) in an elastic and a plastic 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)\]
The stress \(\underline{\sigma}\) is related to the the elastic strain \(\underline{\varepsilon}^{\mathrm{el}}\) by a the standard Hooke behaviour:
\[ \underline{\sigma}= \lambda\,{\mathrm{tr}{\left(\underline{\varepsilon}^{\mathrm{el}}\right)}}\,\underline{I}+2\,\mu\,\underline{\varepsilon}^{\mathrm{el}} \qquad(2)\]
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(3)\]
where:
\(\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.
The negative part of the stress \({\left<\sigma\right>_{-}}\) is an isotropic function of the stress. TFEL/Math
, the mathematical library behind MFront
, offers dedicated functions to compute its value and its derivative \({\displaystyle \frac{\displaystyle \partial {\left<\sigma\right>_{-}}}{\displaystyle \partial \underline{\sigma}}}\).
Equation (3) can be discretized using a standard implicit scheme as follows:
\[ \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}=C\,\eta{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}\,\underline{s} \qquad(4)\]
where:
In the previous expressions, \(\Delta\,T_{e}\) is the effective temperature increment.
Equation (4) shows that \(\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}\) is known if \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\) is known.
Using Equation (1) in its increment form thus leads to the following non linear equation whose only unknows is \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}\):
\[ 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)\]
This non linear equation will be solved by a standard Newton-Raphson algorithm, which requires the jacobian \({\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\):
\({\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is computed as follows:
\[ {\displaystyle \frac{\displaystyle \partial f_{\underline{\varepsilon}^{\mathrm{el}}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{I}+{\displaystyle \frac{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{LITS}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} =\underline{I}+ C\,\left[ \underline{s}\,\otimes\,{\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}+ \eta{\left({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right)}\,{\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} \right] \]
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial \eta}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}} &=\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,{\displaystyle \frac{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\\ &=\theta\,\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ &=\theta\,\gamma\,{\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,{\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\\ \end{aligned} \qquad(6)\]
where \({\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}\) is the elastic stiffness tensor at the middle of the time step:
\[ {\left.\underline{\underline{\mathbf{D}}}\right|_{t+\theta\,\Delta\,t}}={\left.\lambda\right|_{t+\theta\,\Delta\,t}}\,\underline{I}\,\otimes\,\underline{I}+2\,\mu\,\underline{\underline{\mathbf{I}}} \]
Note
The explicit expression of \({\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\) is not given here as it will be computed by an build-in function of the
TFEL/Math
library.
\({\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\) is computed as follows: \[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial C_{m}}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}} &={\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}{\left({{\displaystyle \frac{\displaystyle -{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\right)}\\ &={{\displaystyle \frac{\displaystyle -\underline{I}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}+{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}\,{\displaystyle \frac{\displaystyle \partial }{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}{\left({{\displaystyle \frac{\displaystyle -1}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\right)}\\ &={{\displaystyle \frac{\displaystyle -\underline{I}}{\displaystyle \sqrt{{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}}}\,\underline{I}+{\mathrm{tr}{\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}}\,{{\displaystyle \frac{\displaystyle {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle {\left({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\,\colon\,{\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\right)}^{3/2}}}}\, \end{aligned} \]
Finally, \({\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}\) is given by:
\[ \begin{aligned} {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\underline{\varepsilon}^{\mathrm{el}}}}&= {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\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}}}}\\ &\theta\,{\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}\,\colon\,{\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\,\colon\,\underline{\underline{\mathbf{D}}}\\ \end{aligned} \]
where:
\[ {\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}}= -\nu_{\mathrm{LITS}}\,\underline{I}\,\otimes\,\underline{I}+{\left(1-\nu_{\mathrm{LITS}}\right)}\,\underline{\underline{\mathbf{I}}}\\ \]
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;
The @Behaviour
allows giving a name the behaviour:
@Behaviour LoadInducedThermalStrain_Torelli2018;
Defining the name of behaviour is required.
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 Giacomo Torelli;
@Date 25 / 03 / 2018;
@Description {
"A confinement-dependent load-induced thermal strain "
"constitutive model for concrete subjected "
"to temperatures up to 500 °C"
"Torelli, Giacomo and Mandal, Parthasarathi and Gillie, "
"Martin and Tran, Van-Xuan"
"International Journal of Mechanical Sciences"
"http://www.sciencedirect.com/science/article/pii/S0020740317337372"
}
Those metadata, though recommended, are optional.
StandardElasticity
brickThe behaviour treated in this document allows using the StandardElasticity
brick. This brick, described in this page, provides:
This behaviour brick is fully described here.
The usage of the StandardElasticity
is declared as follows:
@Brick StandardElasticity;
As the StandardElasticity
brick provides support for all modelling hypotheses, we can declare that your behaviour indeed supports all modelling hypotheses:
@ModellingHypotheses{".+"};
The implicit scheme requires three numerical parameters:
StandardElasticity
brick.@Theta 0.5;
@IterMax 100;
@Epsilon 1.e-16;
Note
Those values are automatically generates parameters that can be changed at runtime: the names of those parameters are respectively
theta
,iterMax
andepsilon
.
The @ComputeStiffnessTensor
keyword allows computing the elastic stiffness tensor in a variable called D
. In this case, the UnAltered
version of the stiffness tensor is required (“unaltered” means that the modelling hypothesis was not taken into account).
@ComputeStiffnessTensor<UnAltered>{47000, 0.25};
The @ComputeThermalExpansion
keyword allows computing the thermal strain. It is followed by the mean thermal expansion coefficient.
@ComputeThermalExpansion 10e-6;
The thermal strain is automatically removed to the total strain and the total strain increment.
The following code declares the various parameters used by the constitutive equations:
@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}; -
By default, the Implicit
DSL
declares the elastic strain as a state variable, i.e. a variable which is part of the implicit system and whose value is saved from one time step to the other.
We choose to declare two auxiliary state variables, the load induced thermal strain \(\elits\) and the maximum value of the temperature seen by the material \(T_{\mathrm{max}}\). Those variables are not part of the implicit system but their values are saved from one time step to the other.
Saving maximum value of the temperature seen by the material is required, but its evolution is independent of the mechanics.
Saving the load induced thermal strain is only useful for post-processing. It could be removed for saving some space in memory.
Auxiliary state variables are introduced by the @AuxiliaryStateVariable
keyword, as follows:
@AuxiliaryStateVariable StrainStensor elits;
"LoadInducedThermaStrain");
elits.setEntryName(@AuxiliaryStateVariable temperature Tmax;
"MaximalValueOfTheTemperature"); Tmax.setEntryName(
Local variables are available in every code blocks. Here we will introduce two local variables:
//! LITS increment
@LocalVariable StrainStensor delits;
//! Creep coefficient
@LocalVariable real C;
The @InitLocalVariables
code block is called before the implicit resolution. It is used here to compute the C
variable and update the Tmax
auxiliary variable.
@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; }
The implicit system is defined in the code block introduced by the @Integrator
code block.
This code block is called at each iteration of the implicit system, after the stress \({\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\) at the middle of the time step has been computed and stored in the sig
variable (this is handled by the StandardElasticity
brick).
@Integrator {
We first define some useful constants.
constexpr const stress eeps = 1.e-12;
const stress seps = eeps * young;
constexpr const auto id = Stensor::Id();
constexpr const auto esolver = StressStensor::FSESJACOBIEIGENSOLVER;
Then the negative part of the stress \({\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}\) and its derivative with respect to the stress \({\displaystyle \frac{\displaystyle \partial {\left<{\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}\right>_{-}}}{\displaystyle \partial {\left.\underline{\sigma}\right|_{t+\theta\,\Delta\,t}}}}\) are computed using the computeIsotropicFunctionAndDerivative
method of the stensor
class.
Notes
A scalar function \(f\) defines an isotropic function of a symmetric tensor \(\underline{a}\) as follows:
\[ f{\underline{a}} = \sum_{I=1}^{3}{\left<a_{I}\right>_{-}}\,\underline{n}_{I} \]
where \(a_{I}\) are the eigenvalues of \(\underline{a}\) and \(\underline{n}_{I}\) its eigentensors.
This function requires two scalar functions defining the isotropic function and its derivative. Those functions are called np
and dnp
here:
// negative part
const auto np = [](const stress x) -> stress { return x <= 0 ? x : 0; };
const auto dnp = [&seps](const stress x) -> real {
return std::abs(x) < seps ? 0.5 : ((x < 0) ? 1 : 0);
};
Then we call the computeIsotropicFunctionAndDerivative
method:
auto sn = StressStensor{};
auto dsn = Stensor4{};
std::tie(sn, dsn) =
template computeIsotropicFunctionAndDerivative<esolver>(np, dnp,
sig.0.1); seps *
The previous expression can be simplified in C++-17
:
const auto [sn,dsn] =
template computeIsotropicFunctionAndDerivative<esolver>(np, dnp,
sig.0.1); seps *
About the choice of the eigensolver
Here we used a the Jacobi algorithm as implemented in [2]. This algorithm is lower but more accurate than the default one (see here for an indepth comparison of available algorithms.
The next step is to compute \({\left.\eta\right|_{t+\theta\,\Delta\,t}}\) and its derivative \({\displaystyle \frac{\displaystyle \partial {\left.\eta\right|_{t+\theta\,\Delta\,t}}}{\displaystyle \partial {\left.sig\right|_{t+\theta\,\Delta\,t}}}}\) (See Equation (6)):
//----- Multiaxial correction coefficient
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 deta_dsig = gamma * dcm_dsig;
We are now able to compute the increment \(\Delta\,\underline{\varepsilon}^{\mathrm{LITS}}\) and its derivative:
const auto se = eval((1 + nulits) * sn - nulits * trace(sn) * id);
delits = C * eta * se;const auto dse_dsig = (1 + nulits) * dsn - nulits * ((id ^ id) * dsn);
const auto ddelits_dsig = C * (se ^ deta_dsig) + C * eta * dse_dsig;
Finally, we form the implicit equation associated with the elastic strain and its derivative:
feel += delits;
dfeel_ddeel += theta * ddelits_dsig * D;// end of @Integrator }
Note that we used the fact that feel
is implicitly initialized to \(\Delta\,\underline{\varepsilon}^{\mathrm{el}}-\Delta\,\underline{\varepsilon}^{\mathrm{to}}\).
The @UpdateAuxiliaryStateVariables
code block is called after the implicit resolution, after the computation of the stress at the end of the time step and after the update of the state variables (here the elastic strain).
We use the last estimate of the its increment to update the elits
variable:
@UpdateAuxiliaryStateVariables {
e_lits += delits; }