1% Coupling the Fichant-La Borderie damage behaviour and Torelli' load induced strain behaviour
2% Thomas Helfer
3% 16/12/2020
4
5\newcommand{\absvalue}[1]{{\left|#1\right|}}
6\newcommand{\ppos}[1]{{\left<#1\right>_{+}}}
7\newcommand{\pneg}[1]{{\left<#1\right>_{-}}}
8\newcommand{\paren}[1]{{\left(#1\right)}}
9\newcommand{\tenseur}[1]{\underline{#1}}
10\newcommand{\tenseurq}[1]{\underline{\underline{\mathbf{#1}}}}
11\newcommand{\D}{\ets{\tenseurq{D}}}
12\newcommand{\tepsilonto}{\tenseur{\varepsilon}^{\mathrm{to}}}
13\newcommand{\tepsilonel}{\tenseur{\varepsilon}^{\mathrm{el}}}
14\newcommand{\tepsilonelp}{\ppos{\tenseur{\varepsilon}^{\mathrm{el}}}}
15\newcommand{\eeq}{\varepsilon^{\mathrm{el}}_{\mathrm{eq}}}
16\newcommand{\lits}{\mathrm{LITS}}
17\newcommand{\telits}{\tenseur{\varepsilon}^{\lits}}
18\newcommand{\tdelits}{\tenseur{\dot{\varepsilon}}^{\lits}}
19
20\newcommand{\tsigmae}{\underline{\sigma}^{\mathrm{eff}}}
21\newcommand{\tsigma}{\underline{\sigma}}
22\newcommand{\trace}[1]{{\mathrm{tr}\paren{#1}}}
23\newcommand{\sigmaH}{\sigma_{H}}
24\newcommand{\Frac}[2]{{{\displaystyle \frac{\displaystyle #1}{\displaystyle #2}}}}
25\newcommand{\deriv}[2]{{\displaystyle \frac{\displaystyle \partial #1}{\displaystyle \partial #2}}}
26\newcommand{\derivtot}[2]{{\displaystyle \frac{\displaystyle \mathrm{d} #1}{\displaystyle \mathrm{d} #2}}}
27\newcommand{\sigmaeq}{\sigma_{\mathrm{eq}}}
28\newcommand{\bts}[1]{{\left.#1\right|_{t}}}
29\newcommand{\mts}[1]{{\left.#1\right|_{t+\theta\,\Delta\,t}}}
30\newcommand{\ets}[1]{{\left.#1\right|_{t+\Delta\,t}}}
31
32This document describes how to build a behaviour of concrete describing
33the damage evolution following Fichant-La Borderie and the load induced
34thermal strain (`LITS`) following Torelli et al (See
35[@fichant_endommagement_1996;@gangnant_etude_2016;@torelli_confinement-dependent_2018;@draup_development_2019]
36for details).
37
38We strongly advice readers to refer to the following documents before
39continuing:
40
41- The implementation of the Fichant-La Borderie damage behaviour is
42  detailled [here](FichantLaBorderieDamageBehaviour.html).
43- The implementation of the the load induced thermal strain is detailled
44  [here](LoadInducedThermalStrainBehaviourTorelli2018.html).
45
46# Description of the coupled behaviour
47
48The behaviour is described by a standard split of the strain
49\(\tepsilonto\) in an elastic and a viscoplastic parts, respectively
50denoted \(\tepsilonel\) and \(\telits\):
51
52\[
53\tepsilonto=\tepsilonel+\telits
54\]{#eq:eel}
55
56## Damage behaviour
57
58The stress \(\tsigma\) is related to the the elastic strain
59\(\tepsilonel\) by the Fichant-La Borderie damage behaviour:
60
61\[
62\tsigma = f_{FLB}\paren{d,\tepsilonel}
63\]{#eq:FLB}
64
65where \(d\) is a scalar damage variable. See [this
66page](FichantLaBorderieDamageBehaviour.html) for details.
67
68### Damage evolution
69
70The damage evolution is a function of the equivalent elastic strain
71\(\eeq\), which will be defined hereafter:
72
73- Under a given threshold \(\epsilon_{0}\), no damage occurs.
74- Once this threshold is reached, the damage \(d\) is defined as
75  follows:
76\[
77d=1-\Frac{\epsilon_{0}}{\eeq}\,\exp\paren{B_{t}\,\paren{\epsilon_{0}-\eeq}}
78\]{#eq:d}
79
80However, the latter expression must be modified to take into account the
81irreversibility of the damage evolution. Let \(\bts{d}\) be the value of
82the damage at the beginning of the time step and \(\ets{d}\) its value
83at the end of the time step, \(\ets{d}\) is determined as follows:
84
85\[
86\ets{d}=\max\paren{\bts{d},1-\Frac{\epsilon_{0}}{\ets{\eeq}}\,\exp\paren{B_{t}\,\paren{\epsilon_{0}-\ets{\eeq}}}}
87\]
88
89where \(\ets{\eeq}\) is the value of the equivalent elastic strain at
90the end of the time step.
91
92### Equivalent strain
93
94The equivalent elastic strain \(\eeq\) is defined as:
95
96\[
97\eeq=\sqrt{\tepsilonelp\,\colon\,\tepsilonelp}=
98\sqrt{\sum_{I=1}^{3}\ppos{\varepsilon^{\mathrm{el}}_{I}}^2}
99\]
100
101where \(\tepsilonelp\) is the positive part of the elastic strain and
102\(\left(\varepsilon^{\mathrm{el}}_{I}\right)_{i \in [1:3]}\) are its
103eigenvalues.
104
105### Computation of the stress
106
107The stress tensor \(\tsigma\) is then computed as follows:
108\[
109\ets{\tsigma}=\paren{1-\ets{d}}\,\ppos{\ets{\tsigmae}}+\paren{1-\ets{d}^{a}}\,\pneg{\ets{\tsigmae}}
110\]
111
112where the effective stress tensor \(\tsigmae\) is computed using the
113standard Hooke law:
114
115\[
116\ets{\tsigmae}=\D\,\colon\,\ets{\tepsilonel}
117\]
118
119The behaviour is assumed isotropic, so that the stiffness tensor \(\D\)
120can be related to the first and second Lamé coefficients, denoted
121respectively \(\lambda\) and \(\mu\), as follows:
122
123\[
124\D=\ets{\lambda}\,\tenseur{I}\,\otimes\,\tenseur{I}+2\,\ets{\mu}\,\tenseurq{I}
125\]
126
127
128## Load induced thermal strain
129
130The evolution of the load induced thermal strain \(\telits\) is given by:
131
132\[
133\tdelits=
134\Frac{\eta\paren{\tsigma}\,\beta\paren{T}}{\sigma_{u0}}
135\left[
136  -\nu_{\lits}\trace{\pneg{\tsigma}}\,\tenseur{I}
137  +\paren{1-\nu_{\lits}}\,\pneg{\tsigma}
138\right]\,\dot{T}_{e}
139\]{#eq:delits}
140
141where:
142
143- \(\pneg{\sigma}\) is the negative part of the stress tensor
144  \(\tsigma\). This negative part is computed as follows:
145  \[
146    \pneg{\sigma} = \sum_{I=1}^{3}\pneg{\sigma_{I}}\,\tenseur{n}_{I}
147  \]
148  where \(\sigma_{I}\) is the Ith eigenvalue of \(\tsigma\) and
149  \(\tenseur{n}_{I}\) the associated eigen tensor.
150- \(\eta\paren{\tsigma} = 1+C_{m}\paren{\tsigma}\,\gamma\). \(\eta\)
151  describes the effect of the stress triaxiality.
152- \(C_{m}\paren{\tsigma}=\Frac{-\trace{\pneg{\tsigma}}}{\sqrt{\pneg{\tsigma}\,\colon\,\pneg{\tsigma}}}\).
153  \(C_{m}\) is the triaxiality index linked to the principal stresses.
154- \(\sigma_{u0}\), \(\nu_{\lits}\) and \(\gamma\) are material
155  parameters. More precisely \(\sigma_{u0}\) is the compressive strength
156  of the material and \(\gamma\) is a material parameter calibrated to
157  experimental data for the appropriate temperature range.
158- \(\beta\paren{T}\) is a four order polynomial of the temperature
159  describing uniaxial temperature-LITS behaviour observed in the
160  literature:
161  \[
162    \beta\paren{T}=\beta_{0}+\beta_{1}\,T+\beta_{2}\,T^{2}+\beta_{3}\,T^{3}+\beta_{4}\,T^{4}
163  \]
164
165\(\dot{T}_{e}\) is the effective temperature rate defined as follows:
166
167\[
168\dot{T}_{e}=
169\left\{
170\begin{aligned}
171\dot{T} &\quad\text{if}\quad T=T_{\mathrm{max}}\\
1720 &\quad\text{if}\quad T<T_{\mathrm{max}}\\
173\end{aligned}
174\right.
175\]
176
177where \(T_{\mathrm{max}}\) is the maximal temperature seen by the
178material over time.
179
180# Implicit algorithm
181
182Following the implementation of Torelli' LITS behaviour (see
183[here](LoadInducedThermalStrainBehaviourTorelli2018.html) for details),
184the implicit scheme can be reduced to one non linear equation whose only
185unknown is the elastic strain increment:
186
187\[
188f_{\tepsilonel}=\Delta\,\tepsilonel+\Delta\,\telits\paren{\Delta\,\tepsilonel}-\Delta\,\tepsilonto=\tenseur{0}
189\]{#eq:feel}
190
191Computing the jacobian \(\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}\)
192is straightforward:
193
194\[
195\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel}=\tenseurq{I}+\deriv{\Delta\,\telits}{\mts{\tsigma}}\,\colon\,\deriv{\mts{\tsigma}}{\Delta\,\tepsilonel}
196\]
197
198The derivative \(\deriv{\Delta\,\telits}{\mts{\tsigma}}\) has been
199computed on [this
200page](LoadInducedThermalStrainBehaviourTorelli2018.html) while the
201derivative \(\deriv{\mts{\tsigma}}{\Delta\,\tepsilonel}\) is computed on
202[this page](FichantLaBorderieDamageBehaviour.html).
203
204# Implementation
205
206## Choice of the domain specific language
207
208The domain specific language suitable for implementing an implicit
209scheme is called `Implicit`. The `@DSL` keyword allows specifying the
210domain specific language to be used:
211
212~~~~{.cxx}
213@DSL Implicit;
214~~~~
215
216## Choice of the behaviour name
217
218The `@Behaviour` allows giving a name the behaviour:
219
220~~~~{.cxx}
221@Behaviour FichantLaBorderieDamageAndTorelliLoadInducedThermalStrainBehaviour;
222~~~~
223
224Defining the name of behaviour is required.
225
226
227## Metadata
228
229The `@Author`, `@Date` and `@Description` keywords allows adding
230metadata to the behaviour (respectively the name of the author of the
231implementation, the date at which the implementation was written, a
232short description), as follows:
233
234~~~~{.cxx}
235@Author Thomas Helfer;
236@Date 15 / 12 / 2019;
237@Description {
238  "Implementation of a behaviour combining Fichant-La Borderie' "
239  "description of damage and Torelli'LITS"
240}
241~~~~
242
243## Numerical parameters
244
245We specify two numerical parameters impacting the implicit scheme:
246
247- the \(\theta\) parameter (default value is \(1/2\)). In pratice, we
248  set to \(1\) to indicate that the resolution is purely implicit, but
249  we won't use this value in the implementation.
250- the stopping criteria (default value is \(10^{-8}\)). This may be
251  sufficient in most cases, but not for the computation of the
252  consistent tangent operator.
253
254~~~~{.cxx}
255@Epsilon 1.e-14;
256@Theta 1;
257~~~~
258
259## Thermoelastic properties, computation of the thermal strain
260
261The `@ElasticMaterialProperties` keyword allows defining the material
262properties. When those are defined, the `Implicit` `DSL` automatically
263defines the first and second Lamé coefficients in local variables
264called:
265
266- `lambda` and `mu` at time \(t+\theta\,\Delta\,t\)
267- `lambda_tdt` and `mu_tdt` at time \(t+\Delta\,t\)
268
269~~~~{.cxx}
270@ElasticMaterialProperties{47000, 0.25};
271~~~~
272
273The `@ComputeThermalExpansion` keyword allows computing the thermal
274strain. It is followed by the mean thermal expansion coefficient.
275
276~~~~{.cxx}
277@ComputeThermalExpansion 10e-6;
278~~~~
279
280~~~~{.cxx}
281// parameters of the Fichant-La Borderie part
282@Parameter Bt = 3690.070983;
283@Parameter e0 = 1.03e-04;
284@Parameter a = 3;
285~~~~
286
287~~~~{.cxx}
288// parameters of the LITS part
289@Parameter gamma = 1.5;
290@Parameter sigmultimate = 50.;
291@Parameter nulits = 0.50;
292@Parameter tcrit = 0.;
293@Parameter b[5] = {2.7031065533E-05, -1.0209170592E-06, 6.1200423753E-9,  //
294                   -1.2632648735E-11, 6.9158539621E-15};
295~~~~
296
297
298~~~~{.cxx}
299@AuxiliaryStateVariable real d;
300d.setGlossaryName("Damage");
301@AuxiliaryStateVariable StrainStensor elits;
302elits.setEntryName("LoadInducedThermaStrain");
303@AuxiliaryStateVariable temperature Tmax;
304Tmax.setEntryName("MaximalValueOfTheTemperature");
305~~~~
306
307~~~~{.cxx}
308//! LITS increment
309@LocalVariable StrainStensor delits;
310//! Creep coefficient
311@LocalVariable real C;
312~~~~
313
314~~~~{.cxx}
315@LocalVariable StiffnessTensor dsig_ddeel;
316@LocalVariable real d_p;
317@LocalVariable StiffnessTensor De;
318~~~~
319
320~~~~{.cxx}
321@InitLocalVariables {
322  Tmax = max(max(tcrit, Tmax), T);
323  const auto T_ = T + theta * dT;
324  const auto beta = b[0] + T_ * (b[1] + T_ * (b[2] + T_ * (b[3] + T_ * b[4])));
325  const auto dTe = max(T + dT - max(tcrit, Tmax), temperature(0));
326  C = (beta / (-sigmultimate)) * dTe;
327  De = lambda_tdt * Stensor4::IxI() + 2 * mu_tdt * Stensor4::Id();
328} // end of @InitLocalVariables
329~~~~
330
331~~~~{.cxx}
332@Integrator {
333  constexpr const auto id = Stensor::Id();
334  constexpr const auto esolver = StressStensor::FSESJACOBIEIGENSOLVER;
335  constexpr const stress eeps = 1.e-12;
336  const stress seps = 1.e-12 * young;
337~~~~
338
339~~~~{.cxx}
340  // positive part
341  const auto pp = [](const real x) { return x > 0 ? x : 0; };
342  // derivative of the positive part
343  const auto dpp = [&eeps](const real x) { return std::abs(x) < eeps ? 0.5 : ((x < 0) ? 0 : 1); };
344  // square of the posititve part
345  auto square_ppos = [](const strain& v) { return v > 0 ? v * v : 0; };
346  // elastic strain at the midle of the time step
347  const auto e = eval(eel + deel);
348  // eigen values and eigen tensors of the elastic strain
349  auto e_vp = tvector<3u, strain>{};
350  auto m = tmatrix<3u, 3u, strain>{};
351  e.template computeEigenVectors<esolver>(e_vp, m);
352  // update the damage
353  const auto e_eq = sqrt(square_ppos(e_vp[0]) + square_ppos(e_vp[1]) + square_ppos(e_vp[2]));
354  // effective stress at t+dt
355  const auto Cd = (e0 / e_eq) * exp(Bt * (e0 - e_eq));
356  d_p = (e_eq > e0) ? 1 - Cd : 0;
357  const auto bp = d_p > d;
358  const auto de = bp ? d_p : d;
359  // derivative with respect to the damage
360  const auto dde_ddeel = [&]() -> Stensor {
361    if (!bp) {
362      return Stensor(real(0));
363    }
364    // positive part of the total strain
365    const auto ep = StrainStensor::computeIsotropicFunction(pp, e_vp, m);
366    // derivative of the damage
367    const auto dde_deq = Cd * (Bt + 1 / e_eq);
368    const auto dep_ddeel =
369        StrainStensor::computeIsotropicFunctionDerivative(pp, dpp, e_vp, m, eeps * 0.1);
370    const auto deq_dep = ep / e_eq;
371    return dde_deq * deq_dep * dep_ddeel;
372  }();
373  // function of the damage to simplify expressions
374  const auto de_a = pow(de, a);
375  const auto fpd = (1 - de);
376  const auto fpn = (1 - de_a);
377  // effective stress at the end of the time step
378  const auto l_tr_e = lambda_tdt * trace(e);
379  const auto s = eval(l_tr_e * id + 2 * mu_tdt * e);
380  const auto s_vp = tvector<3u>{l_tr_e + 2 * mu_tdt * e_vp[0],  //
381                                l_tr_e + 2 * mu_tdt * e_vp[1],  //
382                                l_tr_e + 2 * mu_tdt * e_vp[2]};
383  const auto sp = StressStensor::computeIsotropicFunction(pp, s_vp, m);
384  const auto dsp = StressStensor::computeIsotropicFunctionDerivative(pp, dpp, s_vp, m, seps * 0.1);
385  const auto sn = eval(s - sp);
386  const auto dsn = eval(Stensor4::Id() - dsp);
387  // final stress
388  sig = fpd * sp + fpn * sn;
389  // derivative of the stress
390  dsig_ddeel = ((fpd - fpn) * dsp + fpn * Stensor4::Id()) * De;
391  if (bp) {
392    const auto ide = 1 / max(eeps, de);
393    const auto dfpd_dd = -1;
394    const auto dfpn_dd = -a * de_a * ide;
395    dsig_ddeel += dfpd_dd * (sp ^ dde_ddeel) + dfpn_dd * (sn ^ dde_ddeel);
396  }
397~~~~
398
399~~~~{.cxx}
400  // LITS part
401  const auto sn_eq = sqrt(sn | sn);
402  const auto isn_eq = 1 / max(seps, sn_eq);
403  const auto cm = -trace(sn) * isn_eq;
404  const auto dcm_dsig = eval((-isn_eq * id + trace(sn) * power<3>(isn_eq) * sn) | dsn);
405  const auto eta = 1 + (cm - 1) * gamma;
406  const auto se = eval((1 + nulits) * sn - nulits * trace(sn) * id);
407  delits = C * eta * se;
408  const auto deta_dsig = gamma * dcm_dsig;
409  const auto dse_dsig = (1 + nulits) * dsn - nulits * ((id ^ id) * dsn);
410  const auto ddelits_dsig = C * (se ^ deta_dsig) + C * eta * dse_dsig;
411~~~~
412
413
414~~~~{.cxx}
415  // elasticity
416  feel += delits;
417  dfeel_ddeel += ddelits_dsig * dsig_ddeel;
418} // end of @Integrator
419~~~~
420
421~~~~{.cxx}
422@UpdateAuxiliaryStateVariables {
423  d = max(d, d_p);
424  elits += delits;
425}
426~~~~
427
428~~~~{.cxx}
429@TangentOperator {
430  if (smt == ELASTIC) {
431    Dt = De;
432  } else if (smt==CONSISTENTTANGENTOPERATOR){
433    Stensor4 ddeel_ddeto;
434    getPartialJacobianInvert(ddeel_ddeto);
435    Dt = dsig_ddeel * ddeel_ddeto;
436  }
437}
438~~~~
439
440# References
441
442