m (Bsoehngen moved page Draft Söhngen 975815040 to Söhngen 2017a) |
|||
(3 intermediate revisions by 2 users not shown) | |||
Line 84: | Line 84: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\ | + | | style="text-align: center;" | <math>\varphi = \mathop{}\!\left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}}- \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\sigma _{\hbox{Y}} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (4) | | style="width: 5px;text-align: right;white-space: nowrap;" | (4) | ||
Line 120: | Line 120: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\mathcal{{}^{\hbox{M}}H} = \frac{1}{3} \begin{bmatrix}a_2+a_3& -a_3& -a_2& 0& 0& 0& \\ -a_3& a_1+a_3& -a_1& 0& 0& 0& \\ -a_2& -a_1& a_2+a_1& 0& 0& 0& \\ 0& 0& 0& 2 \cdot \ | + | | style="text-align: center;" | <math>\mathcal{{}^{\hbox{M}}H} = \frac{1}{3} \begin{bmatrix}a_2+a_3& -a_3& -a_2& 0& 0& 0& \\ -a_3& a_1+a_3& -a_1& 0& 0& 0& \\ -a_2& -a_1& a_2+a_1& 0& 0& 0& \\ 0& 0& 0& 2 \cdot \text{3/2}\, a_6& 0& 0 \\ 0& 0& 0& 0& 2 \cdot \text{3/2}\,a_5& 0 \\ 0& 0& 0& 0& 0& 2 \cdot \text{3/2}\, a_4 \end{bmatrix} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (7) | | style="width: 5px;text-align: right;white-space: nowrap;" | (7) | ||
Line 137: | Line 137: | ||
|} | |} | ||
− | '''Remark''' The tensor <math display="inline">\mathcal{P}</math> has the property <math display="inline">\mathcal{P}:\boldsymbol{I} = \boldsymbol{0}</math> of a deviatoric projection tensor [Miehe2002] and specifically for <math display="inline">h_i = 1</math> reduces to the deviatoric identity tensor <math display="inline">\mathcal{I}^{\mathrm{dev}} = \mathcal{I}^{\mathrm{sym}} - \frac{1}{3} \boldsymbol{I} \otimes \boldsymbol{I}</math> which yields the classical <math display="inline">\mathrm{J_2}</math>-plasticity for the tensor norm, <math display="inline">\sqrt{\ | + | '''Remark''' The tensor <math display="inline">\mathcal{P}</math> has the property <math display="inline">\mathcal{P}:\boldsymbol{I} = \boldsymbol{0}</math> of a deviatoric projection tensor [Miehe2002] and specifically for <math display="inline">h_i = 1</math> reduces to the deviatoric identity tensor <math display="inline">\mathcal{I}^{\mathrm{dev}} = \mathcal{I}^{\mathrm{sym}} - \frac{1}{3} \boldsymbol{I} \otimes \boldsymbol{I}</math> which yields the classical <math display="inline">\mathrm{J_2}</math>-plasticity for the tensor norm, <math display="inline">\sqrt{\text{3/2}} \left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}\vert _{h_i = 1}} = \sqrt{3\, \mathrm{J_2}(\boldsymbol{\xi }) } = \sqrt{\text{3/2}} \left\vert \left\vert \mathrm{dev}(\boldsymbol{\xi }) \right\vert \right\vert</math> . |
Considering isotropic hardening, the yield stress in eqn:yield_function is represented by the exponential law according to Hockett1975 [Hockett1975], | Considering isotropic hardening, the yield stress in eqn:yield_function is represented by the exponential law according to Hockett1975 [Hockett1975], | ||
Line 179: | Line 179: | ||
|} | |} | ||
− | Declaring the flow rule to be of an associative type, the flow vector reads <math display="inline">\boldsymbol{N} = | + | Declaring the flow rule to be of an associative type, the flow vector reads <math display="inline">\boldsymbol{N} = \partial \text{/} \varphi {\partial \boldsymbol{T}}</math> and the evolution of the plastic strain tensor is defined through |
<span id="eq-12"></span> | <span id="eq-12"></span> | ||
Line 246: | Line 246: | ||
|} | |} | ||
− | are introduced. Since in the converged step <math display="inline">\ | + | are introduced. Since in the converged step <math display="inline">\varphi = 0</math> and thus <math display="inline">\mathop{}\!\left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}}= \sqrt{\text{2/3}} \, \mathop{}\!\sigma _{\hbox{Y}}</math> holds, eqn:relative_stress can be rewritten as |
{| class="formulaSCP" style="width: 100%; text-align: left;" | {| class="formulaSCP" style="width: 100%; text-align: left;" | ||
Line 265: | Line 265: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\mathop{}\!\Delta \gamma ^{k+1} = \mathop{}\!\Delta \gamma ^k - \left[\left.\frac{\partial \ | + | | style="text-align: center;" | <math>\mathop{}\!\Delta \gamma ^{k+1} = \mathop{}\!\Delta \gamma ^k - \left[\left.\frac{\partial \varphi }{\partial \mathop{}\!\Delta \gamma }\right\vert_{\mathop{}\!\Delta \gamma ^k} \right]^{-1} \varphi ^k \, , </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (18) | | style="width: 5px;text-align: right;white-space: nowrap;" | (18) | ||
Line 277: | Line 277: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\frac{\partial \ | + | | style="text-align: center;" | <math>\frac{\partial \varphi }{\partial \mathop{}\!\Delta \gamma } = \boldsymbol{N}_{n+1} : \left[\mathcal{A}_{n+1} \right]^{-1} : \left\{k_2 \kappa ^2_{n+1} \boldsymbol{B}_n - \left\{\left[\frac{1}{\sqrt{{2}/{3}} \sigma _{\mathrm{Y}}} - \frac{\eta _{n+1}}{\sigma _{\mathrm{Y}}} \frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\Delta \gamma } \right]\mathcal{B}_{n+1} \right.\right. </math> |
|- | |- | ||
| style="text-align: center;" | <math> \left.\left. - \eta _{n+1} k_1 k_2 \kappa ^2_{n+1} \mathop{}\!\mathcal{I}^{\mathrm{\,sym}}\right\}: \mathcal{P} : \boldsymbol{\xi }_{n+1} \right\}- \sqrt{\frac{2}{3}} \frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\Delta \gamma } \, , </math> | | style="text-align: center;" | <math> \left.\left. - \eta _{n+1} k_1 k_2 \kappa ^2_{n+1} \mathop{}\!\mathcal{I}^{\mathrm{\,sym}}\right\}: \mathcal{P} : \boldsymbol{\xi }_{n+1} \right\}- \sqrt{\frac{2}{3}} \frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\Delta \gamma } \, , </math> | ||
Line 313: | Line 313: | ||
| style="text-align: center;" | <math> \mathop{}\!\mathrm{d}\boldsymbol{E}^{\mathrm{pl}} = \mathop{}\!\Delta \gamma \frac{\partial \boldsymbol{N}}{\partial \boldsymbol{\xi }} : \mathop{}\!\mathrm{d}\boldsymbol{\xi } + \boldsymbol{N} \mathop{}\!\mathrm{d}\mathop{}\!\Delta \gamma </math> | | style="text-align: center;" | <math> \mathop{}\!\mathrm{d}\boldsymbol{E}^{\mathrm{pl}} = \mathop{}\!\Delta \gamma \frac{\partial \boldsymbol{N}}{\partial \boldsymbol{\xi }} : \mathop{}\!\mathrm{d}\boldsymbol{\xi } + \boldsymbol{N} \mathop{}\!\mathrm{d}\mathop{}\!\Delta \gamma </math> | ||
|- | |- | ||
− | | style="text-align: center;" | <math> \mathop{}\!\mathrm{d}\ | + | | style="text-align: center;" | <math> \mathop{}\!\mathrm{d}\varphi = \mathop{}\!\frac{1}{\left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}}}\mathcal{H} : \boldsymbol{\xi } : \mathop{}\!\mathrm{d}\boldsymbol{\xi } - \sqrt{\frac{2}{3}} \frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\Delta \gamma } \mathop{}\!\mathrm{d}\mathop{}\!\Delta \gamma = 0 \, . </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (21) | | style="width: 5px;text-align: right;white-space: nowrap;" | (21) | ||
Line 363: | Line 363: | ||
|} | |} | ||
− | where <math display="inline">\left[\cdot \right]_{\mathrm{c}}</math> denotes the ''continuum'' operator. Under usage of the consistency condition <math display="inline">\dot{\ | + | where <math display="inline">\left[\cdot \right]_{\mathrm{c}}</math> denotes the ''continuum'' operator. Under usage of the consistency condition <math display="inline">\dot{\varphi }=0</math> an expression for the rate of the backstress can be derived as |
{| class="formulaSCP" style="width: 100%; text-align: left;" | {| class="formulaSCP" style="width: 100%; text-align: left;" | ||
Line 389: | Line 389: | ||
Calculating the limit case <math display="inline">\mathop{}\!\Delta \gamma \to 0</math> for the consistent tangent operator leads, after some rearrangement and by using <math display="inline">\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}= \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \Delta \gamma }</math>, to the continuum tangent operator as expected, <math display="inline">\lim \limits _{\mathop{}\!\Delta \gamma \to 0} \mathcal{C}^{\mathrm{ep}} = \mathcal{C}^{\mathrm{ep}}_{\mathrm{c}}</math>. It is emphasized that although the chosen plasticity model is associative, the kinematic hardening law is of a non-associative type. This leads to a non-symmetric elasto-plastic tangent modulus that needs to be taken into consideration when solving the global equilibrium, as the overall stiffness matrix loses its symmetry as well and therefore needs a solver that can handle non-symmetric systems. The loss of symmetry thereby only refers to the major symmetry of <math display="inline">\mathcal{C}^{\mathrm{ep}}</math>, leading to a non-symmetric representation in Mandel notation <math display="inline">\mathcal{{}^{\hbox{M}}C}^{\mathrm{ep}} \neq [\mathcal{{}^{\hbox{M}}C}^{\mathrm{ep}}]^\mathsf{T}</math>. The algorithmic equations are listed in tab:algorithm for convenience. | Calculating the limit case <math display="inline">\mathop{}\!\Delta \gamma \to 0</math> for the consistent tangent operator leads, after some rearrangement and by using <math display="inline">\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}= \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \Delta \gamma }</math>, to the continuum tangent operator as expected, <math display="inline">\lim \limits _{\mathop{}\!\Delta \gamma \to 0} \mathcal{C}^{\mathrm{ep}} = \mathcal{C}^{\mathrm{ep}}_{\mathrm{c}}</math>. It is emphasized that although the chosen plasticity model is associative, the kinematic hardening law is of a non-associative type. This leads to a non-symmetric elasto-plastic tangent modulus that needs to be taken into consideration when solving the global equilibrium, as the overall stiffness matrix loses its symmetry as well and therefore needs a solver that can handle non-symmetric systems. The loss of symmetry thereby only refers to the major symmetry of <math display="inline">\mathcal{C}^{\mathrm{ep}}</math>, leading to a non-symmetric representation in Mandel notation <math display="inline">\mathcal{{}^{\hbox{M}}C}^{\mathrm{ep}} \neq [\mathcal{{}^{\hbox{M}}C}^{\mathrm{ep}}]^\mathsf{T}</math>. The algorithmic equations are listed in tab:algorithm for convenience. | ||
− | |||
{| class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;" | {| class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;" | ||
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Single-equation return-mapping algorithm | |+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Single-equation return-mapping algorithm | ||
Line 399: | Line 398: | ||
| ~ ~ ~ Trial stress: <math display="inline">\widetilde{\boldsymbol{T}}_{n+1} = \mathcal{C}^{\mathrm{el}} : \left[\boldsymbol{E}_{n+1} - \boldsymbol{E}_{n}^{\mathrm{pl}} \right]</math> | | ~ ~ ~ Trial stress: <math display="inline">\widetilde{\boldsymbol{T}}_{n+1} = \mathcal{C}^{\mathrm{el}} : \left[\boldsymbol{E}_{n+1} - \boldsymbol{E}_{n}^{\mathrm{pl}} \right]</math> | ||
|- | |- | ||
− | | ~ ~ ~ Evaluate yield function: <math display="inline">\ | + | | ~ ~ ~ Evaluate yield function: <math display="inline">\varphi = \lVert \widetilde{\boldsymbol{\xi }}_{n+1} \rVert _{\mathcal{H}} - \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\sigma _{\hbox{Y}}(\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}_n)</math> , <math display="inline">\widetilde{\boldsymbol{\xi }}_{n+1} = \widetilde{\boldsymbol{T}}_{n+1} - \boldsymbol{B}_n</math> |
|- | |- | ||
− | | ~ : '''if''' <math>\ | + | | ~ : '''if''' <math>\varphi \leq 0.0\text{d}0</math> '''then''' |
:: Elastic step Set: <math>\boldsymbol{T}_{n+1} = \widetilde{\boldsymbol{T}}_{n+1}</math> , <math>\boldsymbol{B}_{n+1} = \boldsymbol{B}_{n}</math> , <math>\mathcal{C}^{\mathrm{ep}}_{n+1} = \mathcal{C}^{\mathrm{el}}</math> | :: Elastic step Set: <math>\boldsymbol{T}_{n+1} = \widetilde{\boldsymbol{T}}_{n+1}</math> , <math>\boldsymbol{B}_{n+1} = \boldsymbol{B}_{n}</math> , <math>\mathcal{C}^{\mathrm{ep}}_{n+1} = \mathcal{C}^{\mathrm{el}}</math> | ||
:: | :: | ||
Line 407: | Line 406: | ||
:: | :: | ||
:'''end''' | :'''end''' | ||
− | : Plastic step Set: <math>k=0</math> , <math>\mathop{}\!\Delta \gamma ^0_{n+1} = | + | : Plastic step Set: <math>k=0</math> , <math>\mathop{}\!\Delta \gamma ^0_{n+1} = 0.0\text{d}0</math> , <math>\boldsymbol{\xi }^0_{n+1} = \widetilde{\boldsymbol{\xi }}_{n+1}</math> , <math>\varepsilon _{\mathrm{tol}}=1.0\text{d}-8</math> |
− | : '''while''' <math>|\ | + | : '''while''' <math>|\varphi ^k| > \varepsilon _{\mathrm{tol}}</math> '''do''' |
− | :: <math>\left.\frac{\partial \ | + | :: <math>\left.\frac{\partial \varphi }{\partial \mathop{}\!\Delta \gamma } \right\vert^{k} = \left.\frac{\partial \varphi }{\partial \mathop{}\!\Delta \gamma } \right\vert_{\mathop{}\!\Delta \gamma ^k}</math> |
− | :: <math>\mathop{}\!\Delta \gamma ^{k+1}_{n+1} = \mathop{}\!\Delta \gamma ^k_{n+1} - \left[\left.\frac{\partial \ | + | :: <math>\mathop{}\!\Delta \gamma ^{k+1}_{n+1} = \mathop{}\!\Delta \gamma ^k_{n+1} - \left[\left.\frac{\partial \varphi }{\partial \mathop{}\!\Delta \gamma } \right\vert^{k} \right]^{-1} \varphi ^k</math> |
:: <math>\widehat{\boldsymbol{\xi }}^{\,k+1} = \widetilde{\boldsymbol{T}}_{n+1} - \kappa ^{k+1} \boldsymbol{B}_{n}</math> , <math>\boldsymbol{\xi }^{k+1}_{n+1} = \left[\mathcal{A} \right]^{-1} : \widehat{\boldsymbol{\xi }}^{\,k+1}</math> | :: <math>\widehat{\boldsymbol{\xi }}^{\,k+1} = \widetilde{\boldsymbol{T}}_{n+1} - \kappa ^{k+1} \boldsymbol{B}_{n}</math> , <math>\boldsymbol{\xi }^{k+1}_{n+1} = \left[\mathcal{A} \right]^{-1} : \widehat{\boldsymbol{\xi }}^{\,k+1}</math> | ||
:: <math>{\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}^{k+1}_{n+1} = \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}_{n} + \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\Delta \gamma ^{k+1}_{n+1}</math> , <math>\mathop{}\!\sigma _{\hbox{Y}}^{k+1} = \mathop{}\!\sigma _{\hbox{Y}}({\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}^{k+1}_{n+1})</math> | :: <math>{\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}^{k+1}_{n+1} = \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}_{n} + \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\Delta \gamma ^{k+1}_{n+1}</math> , <math>\mathop{}\!\sigma _{\hbox{Y}}^{k+1} = \mathop{}\!\sigma _{\hbox{Y}}({\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}^{k+1}_{n+1})</math> | ||
− | :: <math>\ | + | :: <math>\varphi ^{k+1} = \left\vert \left\vert \boldsymbol{\xi }^{k+1}_{n+1} \right\vert \right\vert_{\mathcal{H}} - \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\sigma _{\hbox{Y}}^{k+1}</math> |
:: <math>k = k + 1</math> | :: <math>k = k + 1</math> | ||
:: | :: | ||
Line 442: | Line 441: | ||
|[[Image:draft_Söhngen_975815040-FEmesh.png|400px|Finite element mesh of specimen (dots mark optimization nodes)]] | |[[Image:draft_Söhngen_975815040-FEmesh.png|400px|Finite element mesh of specimen (dots mark optimization nodes)]] | ||
|- style="text-align: center; font-size: 75%;" | |- style="text-align: center; font-size: 75%;" | ||
− | | (a) Dimensions of the cruciform specimen in <math>\ | + | | (a) Dimensions of the cruciform specimen in <math>\mathrm{mm}</math> |
|- style="text-align: center; font-size: 75%;" | |- style="text-align: center; font-size: 75%;" | ||
| colspan="2" | '''Figure 1:''' Finite element mesh of specimen (dots mark optimization nodes) | | colspan="2" | '''Figure 1:''' Finite element mesh of specimen (dots mark optimization nodes) | ||
Line 457: | Line 456: | ||
Though the specimen exhibits mainly pure shear as expected it can be seen that also other strain states are induced. | Though the specimen exhibits mainly pure shear as expected it can be seen that also other strain states are induced. | ||
− | The afore-mentioned alternating shear loading of the specimen is depicted in fig:specimen_load. Starting off with the initial unloaded configuration the speckle patterns for the three peak loading values (<math display="inline">|F_{\mathrm{x/y}}|=\ | + | The afore-mentioned alternating shear loading of the specimen is depicted in fig:specimen_load. Starting off with the initial unloaded configuration the speckle patterns for the three peak loading values (<math display="inline">|F_{\mathrm{x/y}}|=8 \mbox{kN}</math>, <math display="inline">|F_{\mathrm{x/y}}|=9 \mbox{kN}</math> and <math display="inline">|F_{\mathrm{x/y}}|=1x10^{1} \mbox{kN}</math>) are shown. The deformation state can be judged by the elliptic deformation of the originally circular bore in the middle. For the loading a linear force rate of <math display="inline">|\dot{F}_{\mathrm{x/y}}| = 1x10^{2} \mbox{N/s}</math> is chosen. <div id='img-3'></div> |
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
Line 510: | Line 509: | ||
|} | |} | ||
− | For that identification the same cruciform specimen geometry was loaded with equi-biaxial tension. In the elastic region the material is assumed to behave isotropically with an elastic modulus of <math display="inline">E= | + | For that identification the same cruciform specimen geometry was loaded with equi-biaxial tension. In the elastic region the material is assumed to behave isotropically with an elastic modulus of <math display="inline">E=210000</math> and a Poisson's ratio of <math display="inline">\nu = 0.3</math>. The optimization of the mixed hardening parameters is based on three different initial values sets that are given in tab:DP600_kinematic_hard along with the mean and standard deviation (SD) of the three optimization runs. |
Line 516: | Line 515: | ||
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Initial and identified parameters of the mixed hardening formulation for DP600 | |+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Initial and identified parameters of the mixed hardening formulation for DP600 | ||
|- | |- | ||
− | | | + | | |
| <math>\sigma _0</math>in MPa | | <math>\sigma _0</math>in MPa | ||
| <math>\sigma _\infty </math>in MPa | | <math>\sigma _\infty </math>in MPa | ||
Line 522: | Line 521: | ||
| <math>k_2</math> | | <math>k_2</math> | ||
|- | |- | ||
− | | | + | | Set 1 |
− | | <math> | + | | <math>3.5x10^{2}</math> |
− | | <math> | + | | <math>1.5x10^{3}</math> |
− | | <math> | + | | <math>1x10^{3}</math> |
− | | <math> | + | | <math>1x10^{1}</math> |
|- | |- | ||
| Set 2 | | Set 2 | ||
− | | <math> | + | | <math>3.5x10^{2}</math> |
− | | <math> | + | | <math>1.5x10^{3}</math> |
− | | <math> | + | | <math>5x10^{3}</math> |
− | | <math> | + | | <math>1x10^{2}</math> |
|- | |- | ||
| Set 3 | | Set 3 | ||
− | | <math> | + | | <math>4x10^{2}</math> |
− | | <math> | + | | <math>3x10^{3}</math> |
− | | <math> | + | | <math>1x10^{3}</math> |
− | | <math> | + | | <math>1x10^{2}</math> |
|- | |- | ||
− | | | + | | Mean of optimization |
− | | <math> | + | | <math>2.6322x10^{2}</math> |
− | | <math> | + | | <math>8.217x10^{2}</math> |
− | | <math> | + | | <math>1.74349x10^{4}</math> |
− | | <math> | + | | <math>1.74102x10^{2}</math> |
|- | |- | ||
| SD | | SD | ||
− | | <math> | + | | <math>3.30969x10^{-2}</math> |
− | | <math> | + | | <math>2.85087x10^{-2}</math> |
− | | <math> | + | | <math>9.03331</math> |
− | | <math> | + | | <math>5.51514x10^{-2}</math> |
|} | |} | ||
− | Observing a relative SD of <math display="inline"> | + | Observing a relative SD of <math display="inline">5x10^{-4}</math> and less confirms the significance of the mean values and hence the global nature of the minimum. To fulfill the above- mentioned convergence criteria the Levenberg-Marquardt optimization of the three sets needs between 14 (Set 2) and 25 (Set 3) iterations. Thereby, the number of iterations refer to the optimization algorithm itself and does not represent the total number of FE simulations which is greater due to the evaluation of the finite differences. In fig:isotropic_hardening the isotropic hardening function for the initial values and the mean of the identified parameters is plotted vs. the equivalent plastic strain (<math display="inline">\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}=\int _0^t \dot{\gamma } \mathop{}\!\mathrm{d}t</math>). The evolution of the backstress in loading direction w.r.t. <math display="inline">\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}</math> is depicted in fig:evolution_backstress. As can be seen the kinematic hardening modulus <math display="inline">k_1</math> leads to an initially steep incline while the parameter <math display="inline">k_2</math> causes the increase in backstress to saturate with advancing plastic strain for the identified parameters. <div id='img-4a'></div> |
<div id='img-4'></div> | <div id='img-4'></div> | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" |
B. SÖHNGEN AND K. WILLNER
Institute of Applied Mechanics |
Friedrich-Alexander-Universität Erlangen-Nürnberg |
Egerlandstr. 5, 91058 Erlangen, Germany |
e-mail: benjamin.soehngen@fau.de, kai.willner@fau.de |
web page: www.ltm.fau.de |
In this work an anisotropic material model at finite strains with nonlinear mixed (isotropic and kinematic) hardening is used for the identification of the hardening parameters of sheet steel. The algorithmic system is thereby reduced to a single equation return mapping. For the identification, a cruciform specimen is loaded biaxially in an alternating shear test to provoke the kinematic hardening behavior and prevent the sheet from buckling. The material parameters are found through an optimization strategy by comparing the deformation field from the experiment to that from a finite element (FE) simulation. The resulting cost function is minimized by means of a gradient-based method.
keywords Kinematic Hardening, Anisotropic Plasticity, Parameter Identification, Biaxial Loading, Sheet Metal
Minimizing production costs and overall weight of products is a major goal for manufacturers in the 21st century. In the context of sheet metals an ongoing development is the incorporation of functional elements through bulk forming operations by means of the technology of sheet-bulk metal forming [Merklein2012]. Here, both sheet and bulk forming is applied to thin sheets to generate complex shape elements and reduce waste.
The accurate rendering of forming processes by the finite element method (FEM) requires suitable material models and proper identification of therein specified parameters. This can be a challenging task as the ever-increasing complexity of the constitutive models demand for appropriate identification methods by simultaneously decreasing experimental effort.
In this contribution a material law capable of describing elasto-plastic material behavior utilizing a Hill-type yield function in combination with a mixed hardening law is used to identify parameters for the dual-phase steel DP600.
In this work, scalar quantities are represented by lower case letters . Boldfaced characters denote second-order tensors; in particular denotes the second-order identity tensor . Fourth- and sixth-order tensors are defined by calligraphic symbols . The dyadic product between two second-order tensors yield a tensor of fourth-order and can be expressed in index notation as . The fourth-order tensor is a linear map from to which is given as in this work (accordingly for sixth-order tensors). The fourth-order symmetric identity tensor is defined through the Kronecker delta as and has the property that it projects a symmetric tensor on itself for any symmetric second-order tensor . The Euclidean norm of a second-order tensor is defined as where the double dot product reads .
The formulation of large strain plasticity is based on the algorithm presented by Miehe2001 [Miehe2001]. It consists of the calculation of stresses and elasticity moduli with regard to a Seth-Hill strain tensor. For the application to elasto-plasticity, the algorithm can be interpreted as a pre- and postprecessing of the stresses and elasto-plastic material tangent originating from a “small strain” algorithm based on the Hencky (logarithmic) strain tensor
|
(1) |
where denotes the right Cauchy-Green tensor, as shown by Miehe2002 [Miehe2002]. Referring to as the work-conjugate stress measure to the logarithmic strain and as the consistent material tangent, the transformation to the Lagrangian stress and module is performed by
|
(2) |
via the fourth and sixth-order projection tensors and that are defined as
|
(3) |
The first and second derivative of the logarithmic strain tensor can be obtained by usage of a spectral decomposition of the right Cauchy-Green tensor and subsequent differentiation of the eigenvalues and -vectors w.r.t. , as shown in [Miehe2001].
The model regards nonlinear isotropic and kinematic hardening as well as Hill-type anisotropic plasticity. The yield function therefore reads
|
(4) |
where the tensor norm is given by
|
(5) |
with the relative stress tensor being defined by the stress tensor and the backstress tensor
|
(6) |
The anisotropic plasticity is modeled by the fourth-order Hill tensor which shows minor (as well as major) symmetry properties and can therefore be represented in a contracted notation. Besides the popular Voigt notation there are others that show certain advantages over the former. The main disadvantages of Voigt's notation are the different representation of (symmetric) second-order tensors for stress- and strain-like quantities (coefficient of 2 for shear strain) and that certain tensor operations will not be transferred to the reduced notation. E.g. the scalar product between two symmetric second-order tensors or the tensor product between a fourth- and a second-order tensor do not yield the same result in Voigt notation, and , but using the representation according to Mandel, and . With the choice of Mandel notation the Hill tensor can be displayed as
|
(7) |
with the anisotropy parameters which are based on the Hill parameters
|
(8) |
Remark The tensor has the property of a deviatoric projection tensor [Miehe2002] and specifically for reduces to the deviatoric identity tensor which yields the classical -plasticity for the tensor norm, .
Considering isotropic hardening, the yield stress in eqn:yield_function is represented by the exponential law according to Hockett1975 [Hockett1975],
|
(9) |
where and are stress-like quantities and and describe the curvature of the yield stress. As the algorithm is not limited to this formulation, other well known hardening laws like the ones from Ludwik, Voce or Swift could be used as well.
The shift of the yield surface in stress space is considered by Armstrong1966 [Armstrong1966] nonlinear kinematic hardening through the evolution of the backstress by
|
(10) |
with the material parameter representing the kinematic hardening modulus and the rate of saturation. For usage in a finite element code the evolution law is discretized in time by an implicit Euler scheme, resulting in
|
(11) |
Declaring the flow rule to be of an associative type, the flow vector reads and the evolution of the plastic strain tensor is defined through
|
(12) |
Using a radial return-mapping algorithm (see e.g. [Neto:1252726]) and starting off with the trial stress as
|
(13) |
where indicates a trial value, the updated stress then reads
|
(14) |
Combining eqn:stress_update with eqn:kinematic_hardening_analy into the relative stress gives
|
(15) |
where the definitions for the stress tensor and the fourth-order tensor ,
|
(16) |
are introduced. Since in the converged step and thus holds, eqn:relative_stress can be rewritten as
|
(17) |
hence being solely a function of the unknown plastic multiplier . By substitution into the definition of the yield function (eqn:yield_function) a scalar equation with only one unknown is derived. To solve the (nonlinear) equation for a Newton-Raphson scheme is applied, giving
|
(18) |
with the iteration counter. The derivative of the yield function w.r.t. the plastic multiplier can be displayed after straightforward usage of the chain rule as
|
(19) |
with the auxiliary variables
|
(20) |
In the algorithmic implementation, the related tangent modulus plays the major role for achieving (quadratic) convergence in the Newton-Raphson iteration for the global equilibrium. It enters the calculation of the overall stiffness and needs to be consistent with the algorithm for stress calculation, hence the name of a consistent elasto-plastic tangent modulus. To develop the algorithmic tangent the system of equations (eqn:plastic_strain_tensor,eqn:kinematic_hardening_analy,eqn:relative_stress,eqn:stress_update,eqn:yield_function) is linearized, giving
|
(21) |
Solving the linearized system for leads to the desired tangent operator
|
(22) |
with the fourth-order auxiliary tensors
|
It should be noted that the consistent tangent modulus is calculated at the end of the stress update, hence all quantities are evaluated at the updated state and the index has been suppressed to improve readability.
The continuum tangent modulus is defined by the rate relation
|
(25) |
where denotes the continuum operator. Under usage of the consistency condition an expression for the rate of the backstress can be derived as
|
(26) |
which - inserted back into eqn:continuum_tangent_rate_relation along with the consistency condition and the rate of the plastic strain tensor (eqn:plastic_strain_tensor) - yields the continuum tangent operator as
|
(27) |
Calculating the limit case for the consistent tangent operator leads, after some rearrangement and by using , to the continuum tangent operator as expected, . It is emphasized that although the chosen plasticity model is associative, the kinematic hardening law is of a non-associative type. This leads to a non-symmetric elasto-plastic tangent modulus that needs to be taken into consideration when solving the global equilibrium, as the overall stiffness matrix loses its symmetry as well and therefore needs a solver that can handle non-symmetric systems. The loss of symmetry thereby only refers to the major symmetry of , leading to a non-symmetric representation in Mandel notation . The algorithmic equations are listed in tab:algorithm for convenience.
~ ~ Geometric preprocessing: , |
~ ~ Plasticity algorithm: |
~ ~ ~ Trial stress: |
~ ~ ~ Evaluate yield function: , |
~ : if then
|
~ ~ Geometric postprocessing: , |
Remark As already mentioned above, all equations can be represented in some compressed notation. Every fourth-order tensor shows at least minor symmetry and by making use of e.g. Mandel representation, the calculation time on local Gauss point level can be reduced substantially.
The experimental investigation is carried out on a biaxial testing machine with four electro-mechanical actuators of which the two of each axis are operated in master/slave control to keep the specimen centered and inhibit rigid body motion. To prevent buckling, the loading along the two axes is applied in an alternating fashion rather than by pure tension/compression. Therefore, the applied force is equal in amount but opposite in direction, resulting in an alternating shear loading.
![]() |
![]() |
(a) Dimensions of the cruciform specimen in | |
Figure 1: Finite element mesh of specimen (dots mark optimization nodes) |
The captured images are subsequently analyzed with the digital image correlation (DIC) software Aramis to generate the displacements. For the numerical computation a FE-mesh consisting of 3080 eight-noded 3D continuum elements (two in thickness direction) with linear shape functions is used. Assuming a uniform load distribution generated by the clamping, the respective edge nodes are coupled through nodal ties and the measured forces are taken as boundary conditions. To assess the residual of experimental and numerical displacements only a subset of the FE nodes are considered. These optimization nodes are shown in fig:cruciform_mesh along with the discretization.
To encourage the heterogeneity of the deformation state a bore of 6mm in diameter is introduced at the center of the specimen, resulting in various strain states. In fig:cruciform_strain_states major vs. minor strain is plotted for the optimization nodes along with the strain states for uniaxial compression (), pure shear (), uniaxial tension (), plane strain () and equi-biaxial tension ().Distribution of major vs. minor strain in the specimen during loading |
Figure 2: Distribution of major vs. minor strain in the specimen during loading |
Though the specimen exhibits mainly pure shear as expected it can be seen that also other strain states are induced.
The afore-mentioned alternating shear loading of the specimen is depicted in fig:specimen_load. Starting off with the initial unloaded configuration the speckle patterns for the three peak loading values (, and ) are shown. The deformation state can be judged by the elliptic deformation of the originally circular bore in the middle. For the loading a linear force rate of is chosen.The identification process is based on a comparison of the displacements obtained from experiment and FE simulation. The minimization function can be formulated as
|
(28) |
with denoting the nodal displacements from the FE-calculation and the displacements coming from the experimental full-field measurement. As the discretization of the numerical analysis and that of the image correlation do in general not coincide, the operator is used to map the experimentally obtained deformation values on the nodal positions of the FE-mesh. The discretization introduced by DIC is thereby about ten times finer than that of the FE-mesh, hence the possible error derived from interpolation should be marginal. The vector contains the parameters to be optimized, in the present case it consists of the four parameters describing the mixed hardening behavior, . As a remark, the identification of the mixed hardening parameters demand a strategy where both phenomena are considered simultaneously. The parameters for isotropic hardening may not be taken from e.g. a uniaxial tensile test.
To minimize the cost function , an appropriate optimization algorithm has to be chosen. In general, one distinguishes between gradient-based and gradient-free methods, the latter having the advantage of not needing derivatives of the minimization functions. Nevertheless, in the context of material parameter identification, experience shows that gradient-based approaches should be chosen over gradient-free methods. Not only having an increased rate of convergence by nature, the amount of function evaluations can be significantly lower, despite using a finite difference scheme to approximate the gradient, see as well e.g. [Schmaltz2013]. The application of a finite difference procedure to evaluate the sensitivity of the cost function w.r.t. the design parameters seems a necessary evil. An attempt to an analytical derivation for elasto-plastic material behaviour was made by Cooreman_2007 [Cooreman_2007]. The derivation however is simplified substantially and only applicable to isotropic plasticity with a strictly homogeneous strain field (simple tensile test). Moreover, the sensitivity of every node has to be calculated in an iterative manner, leading to the question about performance gain over finite differences, which keeps unanswered.
A well-suited and in the area of parameter optimization widely used (see e.g. [Chaparro_2008]) algorithm is the Levenberg-Marquardt method. It shows quadratic convergence (see e.g. [Nocedal2006]) and is used with a forward finite differences scheme in this work. The minimum of the cost function is assumed to be found once either the relative step size of the design parameters or of the objective function itself reaches a threshold of [round-mode=places,round-precision=1,output-exponent-marker=]1.e-6. The parameters for the anisotropy of the yield surface for the considered material have been identified in [Landkammerinpress] and are listed in tab:DP600_Hillparameter.
For that identification the same cruciform specimen geometry was loaded with equi-biaxial tension. In the elastic region the material is assumed to behave isotropically with an elastic modulus of and a Poisson's ratio of . The optimization of the mixed hardening parameters is based on three different initial values sets that are given in tab:DP600_kinematic_hard along with the mean and standard deviation (SD) of the three optimization runs.
in MPa | in MPa | |||
Set 1 | ||||
Set 2 | ||||
Set 3 | ||||
Mean of optimization | ||||
SD |
Evolution of isotropic hardening w.r.t. the equivalent plastic strain | Evolution of backstress in loading direction w.r.t. the equivalent plastic strain |
(a) Evolution of isotropic hardening w.r.t. the equivalent plastic strain | |
Figure 4: Evolution of backstress in loading direction w.r.t. the equivalent plastic strain |
In this contribution a single equation return mapping algorithm for large strain anisotropic plasticity with isotropic and kinematic hardening is presented. The identification of the mixed hardening parameters for the dual-phase steel DP600 with a Levenberg-Marquardt approach is performed on a cruciform specimen, which is loaded in an alternating shear test on a biaxial testing machine. To assure the obtained minimum of the cost function to be of a global nature, the optimization is run with three distinct initial values sets, whereby all lead to similar material parameters with only small deviations. The material model is written in Fortran and can be used (with rather minor adjustments) in commercial FE-packages like Marc or Abaqus through the provided user subroutines hypela2 and UMAT, respectively.
To further enhance the identification of material parameters with the proposed method, uncertainties in the experimental investigation will be considered in a next step. This includes in particular the determination of the deformation field through digital image correlation where the uncertainty in dependence of the deformation state can be included in the minimization function by weighting factors.
Acknowledgments. This work is supported by the German Research Foundation (DFG) within the Transregional Collaborative Research Centre TCRC 73: “Manufacturing of complex functional components with variants by using a new sheet metal forming process - Sheet-Bulk Metal Forming” (www.tr-73.de).
Published on 24/05/17
Submitted on 17/05/17
Licence: Other