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