1% Implementing an isotropic elliptic yield criterion of the Green type 2% Thomas Helfer 3% 30/11/2017 4 5\newcommand{\bts}[1]{{\left.#1\right|_{t}}} 6\newcommand{\mts}[1]{{\left.#1\right|_{t+\theta\,\Delta\,t}}} 7\newcommand{\ets}[1]{{\left.#1\right|_{t+\Delta\,t}}} 8\newcommand{\trace}[1]{{\mathrm{tr}\paren{#1}}} 9\newcommand{\tenseur}[1]{\underline{#1}} 10\newcommand{\tenseurq}[1]{\underline{\underline{\mathbf{#1}}}} 11\newcommand{\tns}[1]{{\underset{\tilde{}}{\mathbf{#1}}}} 12\newcommand{\transpose}[1]{{#1^{\mathop{T}}}} 13 14\newcommand{\tepsilonto}{\tenseur{\varepsilon}^{\mathrm{to}}} 15\newcommand{\tepsilonel}{\tenseur{\varepsilon}^{\mathrm{el}}} 16\newcommand{\tepsilonp}{\tenseur{\varepsilon}^{\mathrm{p}}} 17\newcommand{\tdepsilonp}{\tenseur{\dot{\varepsilon}}^{\mathrm{p}}} 18\newcommand{\tsigma}{\underline{\sigma}} 19\newcommand{\sigmaeq}{\sigma_{\mathrm{eq}}} 20\newcommand{\trace}[1]{{\mathrm{tr}\paren{#1}}} 21\newcommand{\Frac}[2]{{{\displaystyle \frac{\displaystyle #1}{\displaystyle #2}}}} 22\newcommand{\deriv}[2]{{\displaystyle \frac{\displaystyle \partial #1}{\displaystyle \partial #2}}} 23\newcommand{\sderiv}[2]{{\displaystyle \frac{\displaystyle \partial^{2} #1}{\displaystyle \partial #2^{2}}}} 24\newcommand{\dtot}{{{\mathrm{d}}}} 25\newcommand{\paren}[1]{{\left(#1\right)}} 26 27This article shows how to implement an isotropic elliptic yield 28criterion of the Green type in `MFront`. Such an example illustrates 29the usage of `StandardElasticity` brick (see 30[this page](BehaviourBricks.html)). This page is inspired by the paper 31of Fritzen et al. (See @fritzen_computational_2013). 32 33The whole implementation is available 34[here](./gallery/plasticity/GreenPerfectPlasticity.mfront). 35 36# Description of the behaviour 37 38The behaviour is described by a standard decomposition of the strain 39\(\tepsilonto\) in an elastic and a plastic component, respectively 40denoted \(\tepsilonel\) and \(\tepsilonp\): 41 42\[ 43\tepsilonto=\tepsilonel+\tepsilonp 44\] 45 46## Elastic behaviour 47 48The stress \(\tsigma\) is related to the the elastic strain 49\(\tepsilonel\) by a the standard Hooke law expressed using the Lamé 50coefficients: 51\[ 52\tsigma = \lambda\,\trace{\tepsilonel}\tenseur{I}+2\,\mu\,\tepsilonel 53\] 54 55Cette relation est directement prise en charge par la brique 56`StandardElasticity`. 57 58## Plastic flow 59 60The behaviour is based on the definition of an equivalent stress 61\(\sigmaeq\) defined as follows: 62\[ 63\sigmaeq=\sqrt{\Frac{3}{2}\,C\,\tenseur{s}\,\colon\,\tenseur{s}+F\,\trace{\tsigma}^{2}} 64\] 65where \(\tenseur{s}\) is the deviatoric stress tensor: 66\[ 67\tenseur{s}=\tsigma-\Frac{1}{3}\,\trace{\tsigma}\,\tenseur{I} 68\] 69This equivalent stress is an homogeneous function of the stress of 70degree 1. 71 72The plastic part of the behaviour is described by the following yield 73surface: 74\[ 75f\paren{\tsigma} = \sigmaeq-\sigma_{0} 76\] 77where \(\sigma_{0}\) is the yield stress. 78 79Following the principle of maximum plastic dissipation, the plastic 80flow is assumed to be associated, so the flow direction 81\(\tenseur{n}\) is given by \(\deriv{f}{\tsigma}\): 82 83\[ 84\tenseur{n} = \deriv{f}{\tsigma} = \Frac{1}{\sigmaeq}\,\paren{\Frac{3}{2}\,C\,\tenseur{s}+F\,\trace{\tsigma}\,\tenseur{I}} 85\] 86 87Since the equivalent stress is an homogeneous function of degree 1, 88the equivalent plastic strain \(p\) is well defined and the plastic 89strain rate can be written: 90\[ 91\tdepsilonp=\dot{p}\,\tenseur{n} 92\] 93 94# Integration algorithm 95 96The previous constitutive equations will be integrated using a 97standard implicit scheme. 98 99The first thing that the integration will determine is whether a 100plastic flow occurs during the time step. To do this, a plastic 101prediction of the stress is made. If this prediction is inside the 102yield surface, the step is purely elastic and the integration is 103trivial. On the other hand, the material must lie on the yield surface 104at the end of the time step. 105 106## State variables 107 108The state variable considered will be: 109 110- the elastic strain \(\tepsilonel\), which is automatically declared 111 by the `StandardElasticity` brick as `eel`. 112- the equivalent plastic strain \(p\). 113 114## Plastic loading case 115 116### Implicit system 117 118Assuming a plastic loading, the system of equations to be solved is: 119\[ 120\left\{ 121\begin{aligned} 122 \Delta\,\tepsilonel-\Delta\,\tepsilonto+\Delta\,p\,\mts{\tenseur{n}} &= 0 \\ 123 f\paren{\mts{\sigmaeq}} &= 0 \\ 124\end{aligned} 125\right. 126\] 127 128where \(\mts{X}\) is the value of \(X\) at \(t+\theta\,\Delta\,t\), 129\(\theta\) being a numerical parameter. In the following, the first 130(tensorial) equation is noted \(f_{\tepsilonel}\) and the second 131(scalar) equation is noted \(f_{p}\). 132 133In practice, it is physically sound to make satisfy exactly the yield 134condition at the end of the time step (otherwise, stress extrapolation 135can lead to stress state outside the yield surface and spurious 136oscillations can also be observed). This leads to the choice 137\(\theta=1\). 138 139### Computation of the jacobian 140 141The jacobian \(J\) of the implicit system can be decomposed by blocks: 142\[ 143J= 144\begin{pmatrix} 145\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel} & \deriv{f_{\tepsilonel}}{\Delta\,p} & \\\\ 146\deriv{f_{p}}{\Delta\,\tepsilonel} & \deriv{f_{p}}{\Delta\,p} \\ 147\end{pmatrix} 148\] 149 150The expression of the previous terms is given by: 151 152\[ 153\left\{ 154\begin{aligned} 155\deriv{f_{\tepsilonel}}{\Delta\,\tepsilonel} &= \tenseur{I} + dp\,\deriv{\mts{\tenseur{n}}}{\Delta\,\tepsilonel}\\ 156\deriv{f_{\tepsilonel}}{\Delta\,p} &= \mts{\tenseur{n}}\\ 157\deriv{f_{p}}{\Delta\,\tepsilonel} &= -\theta\,\mts{\tenseur{n}}\,\colon\,\mts{\tenseurq{D}}\\ 158\deriv{f_{p}}{\Delta\,p} &= 0 159\end{aligned} 160\right. 161\] 162 163\[ 164\begin{aligned} 165\deriv{\mts{\tenseur{n}}}{\Delta\,\tepsilonel} 166&=\deriv{\mts{\tenseur{n}}}{\mts{\sigma}}\,\deriv{\mts{\sigma}}{\mts{\tepsilonel}}\,\deriv{\mts{\tepsilonel}}{\Delta\,\tepsilonel}\\ 167&=\Frac{\theta}{\sigmaeq}\,\paren{\Frac{3}{2}\,C\,\paren{\tenseurq{I}-\Frac{1}{3}\,\tenseur{I}\,\otimes\,{\tenseur{I}}}+F\,\tenseur{I}\,\otimes\,{\tenseur{I}}-\mts{\tenseur{n}}\,\otimes\,\mts{\tenseur{n}}}\,\colon\,\mts{\tenseurq{D}}\\ 168&=\Frac{\theta}{\sigmaeq}\,\paren{\Frac{3}{2}\,C\,\,\tenseurq{I}+\paren{F-\paren{\Frac{C}{2}}}\,\tenseur{I}\,\otimes\,\tenseur{I}-\mts{\tenseur{n}}\,\otimes\,\mts{\tenseur{n}}}\,\colon\,\mts{\tenseurq{D}}\\ 169\end{aligned} 170\] 171 172This expression could be further simplified. 173 174\paren{\mts{\lambda}\,\tenseur{I}\,\otimes\,\tenseur{I}+2\,\mts{\mu}\,\tenseurq{I}} 175 176## Elastic loading case 177 178Assuming an elastic loading, the system of equations to be solved is 179trivially: 180\[ 181\left\{ 182\begin{aligned} 183 \Delta\,\tepsilonel-\Delta\,\tepsilonto &= 0 \\ 184 \Delta\,p &= 0 \\ 185\end{aligned} 186\right. 187\] 188 189The jacobian associated with this system is the identity matrix. 190 191# Implementation 192 193In the plastic loading case, the system of equations to be solved is 194*non-linear*. We choose the Newton-Raphson algorithm with an 195analytical jacobian. Compared to other algorithm available in `MFront` 196(Runge-Kutta, Broyden, Newton-Raphson with numerical jacobian, etc..), 197this algorithm is generally the *most efficient* in pratice. 198 199## Preamble 200 201The implementation begins with the choice of the `Implicit` domain 202specific language (dsl): 203 204~~~~{.cpp} 205@DSL Implicit; 206~~~~ 207 208Note that this dsl automatically declares the elastic strain `eel` as 209a state variable. 210 211As discussed before, we explicit state that a fully implicit 212integration will be used by default: 213 214~~~~{.cpp} 215@Theta 1; 216~~~~ 217 218This value can be changed at runtime by modifying the `theta` 219parameter. 220 221The stopping criterion is chosen low, to ensure the quality of the 222consistent tangent operator (the default value, \(10^{-8}\) is enough 223to ensure a good estimation of the state variable evolution, but is 224not enough to have a proper estimation of the consistent tangent 225operator): 226 227~~~~{.cpp} 228@Epsilon 1e-14; 229~~~~ 230 231We then declare that we want to support all the modelling hypotheses: 232 233~~~~{.cpp} 234@ModellingHypotheses {".+"}; 235~~~~ 236 237The support of plane stress modelling hypotheses are handled by the 238`StandardElasticity` brick and will not be discussed here. 239 240### Usage of the `StandardElasticity` brick 241 242To implement this behaviour, we will use the `StandardElasticity` 243brick which provides: 244 245- Automatic computation of the stress tensor at various stages of the 246 behaviour integration. 247- Automatic computation of the consistent tangent operator. 248- Automatic support for plane stress and generalized plane stress 249 modelling hypotheses (The axial strain is defined as an additional 250 state variable and the associated equation in the implicit system is 251 added to enforce the plane stess condition). 252- Automatic addition of the standard terms associated with the elastic 253 strain state variable. 254 255This behaviour brick is fully described [here](BehaviourBricks.html). 256 257The usage of the `StandardElasticity` is declared as follows: 258 259~~~~{.cpp} 260@Brick StandardElasticity; 261~~~~ 262 263### Elastic stiffness tensor 264 265The elastic stiffness tensor \(D\) is defined using 266`@ComputeStiffnessTensor` keyword by giving the elastic material 267properties as constants: 268 269~~~~{.cpp} 270@ComputeStiffnessTensor<UnAltered> {150e9,0.3}; 271~~~~ 272 273This computed stiffness is stored in a variable `D`. A second variable 274`D_tdt` is also introduced. As the material properties are constants, 275`D_tdt` is an alias to `D`. 276 277The elastic material properties can be changed at runtime time by 278modifying the following parameters: `YoungModulus` and `PoissonRatio`. 279 280Rather than constants, one can also use correlations implemented in 281seperate `MFront` files. This allows to take into account the 282dependency of the material properties with the temperature for 283example. In this case, the variable `D` contains the stiffness tensor 284at \(t+\theta\,\Delta\,t\) and the variable `D_tdt` contains the 285stiffness tensor at \(t+\Delta\,t\). 286 287Another possibility is to use the `@RequireStiffnessTensor` 288keyword. In this case, the elastic material properties must be 289computed by the calling solver at the end of the time step (and 290furnished to the mechanical behaviours through hidden material 291properties). 292 293### Variable declarations 294 295#### State variables 296 297As stated earlier, the state variable `eel` is automatically declared 298by the `Implicit` dsl. 299 300The equivalent plastic strain state variable `p` is declared as: 301 302~~~~{.cpp} 303@StateVariable real p; 304~~~~ 305 306We then associate the appropriate glossary name to this variable: 307 308~~~~{.cpp} 309p.setGlossaryName("EquivalentPlasticStrain"); 310~~~~ 311 312#### Parameters 313 314The definition of yield surface introduces three material coefficients 315\(C\), \(F\) and \(\sigma_{0}\) that we declare as a parameter: 316 317~~~~{.cpp} 318@Parameter C = 0.8; 319C.setEntryName("GreenYieldCriterion_C"); 320@Parameter F = 0.2; 321F.setEntryName("GreenYieldCriterion_F"); 322@Parameter s0 = 150e6; 323s0.setGlossaryName("YieldStress"); 324~~~~ 325 326The `YieldStress` is an entry of the glossary (see 327[here](glossary.html)). 328 329#### Local variable 330 331To select the implicit system associated either with the elastic or 332plastic loading case, we introduce a boolean local variable `b`. 333 334~~~~{.cpp} 335@LocalVariable bool b; // if true, plastic loading 336~~~~ 337 338## Initialization of the local variables, determination of the loading case 339 340Before solving the implicit system, the code block introduced by the 341`@InitLocalVariables` keyword is executed. For this behaviour, this 342block will select either the elastic or plastic loading case. 343 344We first make an *elastic prediction* of the stress at the end of the 345time step. We use the `computeElasticPrediction` method introduced by 346the `StandardElasticity` brick. This method takes into account the 347modelling hypothesis, which is mandatory for plane stress modelling 348hypotheses. We then make an elastic prediction of the Hill equivalent 349stress and check whether or not this elastic prediction is inside the 350elastic domain. The latter information is stored in the boolean value 351`b` which will be `false` (no plastic loading) if the loading is 352elastic. 353 354~~~~{.cpp} 355@InitLocalVariables{ 356 const auto sig_el = computeElasticPrediction(); 357 const auto s_el = deviator(sig_el); 358 const auto tr_el = trace(sig_el); 359 const auto seq = sqrt(3*C*(s_el|s_el)/2+F*tr_el*tr_el); 360 b = seq-s0 > 0; 361} 362~~~~ 363 364## Implicit system and jacobian 365 366Finally, we describe how the implicit system and the computation of 367the jacobian is written in a code block introduced by the 368`@Integrator` keyword. 369 370We use the following facts: 371 372- The equations of implicit system are initialized to the state 373 variables increments (i.e. `feel` is initialized to `deel`). 374- The jacobian \(J\) is initialized to the identity 375 (i.e. `dfeel_ddeel` is initialized to the identity tensor). 376- The increment of the total strain is automatically deduced from 377 `feel` by the `StandardElasticity` brick. 378 379Apart from those facts, the code is an almost direct translation of 380the mathematical expression described in previous sections and boils 381down to the following lines of code: 382 383~~~~{.cpp} 384@Integrator{ 385constexpr const auto id = Stensor::Id(); 386 constexpr const auto id4 = Stensor4::Id(); 387 if(b){ 388 const auto hC = C/2; 389 const auto s = deviator(sig); 390 const auto tr = trace(sig); 391 const auto seq = sqrt(3*hC*(s|s)+F*tr*tr); 392 const auto iseq = 1/(max(seq,real(1.e-10*young))); 393 const auto n = eval(iseq*(3*hC*s+F*tr*id)); 394 // elasticity 395 feel += dp*n; 396 dfeel_ddeel += theta*dp*iseq*(3*hC*id4+(F-hC)*(id^id)-(n^n))*D; 397 dfeel_ddp = n; 398 // plasticity 399 fp = (seq-s0)/young; 400 dfp_ddp = strain(0); 401 dfp_ddeel = theta*(n|D)/young; 402 } 403} 404~~~~ 405 406# References 407 408<!-- Local IspellDict: english --> 409