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