Identification of nonlinear kinematic hardening parameters for sheet metal from biaxial loading tests

==Identification of Nonlinear Kinematic Hardening Parameters for Sheet  Metal from Biaxial Loading Tests==

'''B. SÖHNGEN<math>^{*}</math> AND K. WILLNER<math>^{*}</math>'''

{|
|-
|<math>^{*}</math>Institute of Applied Mechanics
|-
| Friedrich-Alexander-Universität Erlangen-Nürnberg
|-
| Egerlandstr. 5, 91058 Erlangen, Germany
|-
| e-mail: [mailto:benjamin.soehngen@fau.de benjamin.soehngen@fau.de], [mailto:kai.willner@fau.de kai.willner@fau.de]
|-
| web page: [http://www.ltm.fau.de www.ltm.fau.de]
|}

==Abstract==

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

==1 INTRODUCTION==

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.

==2 MATERIAL MODEL==

===2.1 Notation===

In this work, scalar quantities are represented by lower case letters <math display="inline">a,b, \dots </math> . Boldfaced characters <math display="inline">\boldsymbol{A}, \boldsymbol{B}, \dots </math> denote second-order tensors;  in particular <math display="inline">\boldsymbol{I}</math> denotes the second-order identity tensor <math display="inline">\boldsymbol{A} \cdot \boldsymbol{I} = \boldsymbol{A}</math>. Fourth- and sixth-order tensors are defined by calligraphic symbols <math display="inline">\mathcal{A}, \mathcal{B}, \dots </math>. The dyadic product between two second-order tensors yield a tensor of fourth-order <math display="inline">\mathcal{C}=\boldsymbol{A} \otimes \boldsymbol{B}</math> and can be  expressed in index notation as <math display="inline">C_{ijkl} = A_{ij} B_{kl}</math>. The fourth-order tensor <math display="inline">\mathcal{C}</math> is a linear map from <math display="inline">\boldsymbol{A}</math> to <math display="inline">\boldsymbol{B}</math> which is given as <math display="inline">\boldsymbol{B} = \mathcal{C} : \boldsymbol{A}</math> in this work (accordingly for sixth-order tensors). The fourth-order symmetric identity tensor is defined through the Kronecker delta as <math display="inline">I^{\mathrm{sym}} = \frac{1}{2} ( \delta _{ik} \delta _{jl} + \delta _{jk} \delta _{il} )</math> and has the property that it projects a symmetric tensor on itself <math display="inline">\boldsymbol{A}^{\mathrm{sym}} = \mathcal{I}^{\mathrm{sym}} : \boldsymbol{A}^{\mathrm{sym}}</math> for any symmetric second-order tensor <math display="inline">\boldsymbol{A}^{\mathrm{sym}} = \frac{1}{2}( \boldsymbol{A} + \boldsymbol{A}^\mathsf{T})</math>.  The Euclidean norm of a second-order tensor is defined as <math display="inline">\left\vert \left\vert \boldsymbol{A} \right\vert \right\vert = \sqrt{\boldsymbol{A}:\boldsymbol{A}}</math> where the double dot product reads <math display="inline">a = \boldsymbol{A} : \boldsymbol{B} = A_{ij}B_{ij}</math>.

===2.2 Model formulation===

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}^0 = \frac{1}{2} \ln \boldsymbol{C} \, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\boldsymbol{C} = \boldsymbol{F}^\mathsf{T}\boldsymbol{F}</math> denotes the right Cauchy-Green tensor, as shown by Miehe2002 [Miehe2002].  Referring to <math display="inline">\boldsymbol{T}</math> as the work-conjugate stress measure to the logarithmic strain <math display="inline">\boldsymbol{E}^0</math>  and <math display="inline">\mathcal{C}^{\mathrm{ep}}</math> as the consistent material tangent, the transformation to the Lagrangian stress and module is performed by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{S} = \boldsymbol{T} : \mathcal{P}  \quad   \hbox{and}  \quad   \mathcal{C}^{\mathrm{ep}}_{\mathrm{L}} = \mathcal{P}^\mathsf{T}: \mathcal{C}^{\mathrm{ep}} : \mathcal{P} + \boldsymbol{T} : \mathcal{L} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

via the fourth and sixth-order projection tensors <math display="inline">\mathcal{P}</math> and <math display="inline">\mathcal{L}</math> that are defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{P} = 2 \frac{ \partial  \boldsymbol{E}^0 }{\partial \boldsymbol{C}}  \quad   \hbox{and}  \quad   \mathcal{L} = 4 \frac{ \partial ^2  \boldsymbol{E}^0 }{\partial \boldsymbol{C} \partial \boldsymbol{C}}  \, . </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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. <math display="inline">\boldsymbol{C}</math>, as shown in [Miehe2001].

The model regards nonlinear isotropic and kinematic hardening as well as Hill-type anisotropic plasticity. The yield function therefore reads

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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)
|}

where the tensor norm <math display="inline">\mathop{}\!\left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}}</math> is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathop{}\!\left\vert \left\vert \boldsymbol{\xi } \right\vert \right\vert_{\mathcal{H}}= \mathop{}\!\left\vert \left\vert \boldsymbol{\xi }^\mathsf{T}: \mathcal{H} : \boldsymbol{\xi } \right\vert \right\vert </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

with the relative stress tensor <math display="inline">\boldsymbol{\xi }</math> being defined by the stress tensor <math display="inline">\boldsymbol{T}</math> and the backstress tensor <math display="inline">\boldsymbol{B}</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\xi } = \boldsymbol{T} - \boldsymbol{B} \, . </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

The anisotropic plasticity is modeled by the fourth-order Hill tensor <math display="inline">\mathcal{H}</math> 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, <math display="inline">\boldsymbol{A}:\boldsymbol{B} \neq \boldsymbol{{}^{\hbox{V}}A}\cdot \boldsymbol{{}^{\hbox{V}}B}</math> and <math display="inline">\mathcal{A}:\boldsymbol{B} \neq \mathcal{{}^{\hbox{V}}A} \cdot \boldsymbol{{}^{\hbox{V}}B}</math>, but using the representation according to Mandel, <math display="inline">\boldsymbol{A}:\boldsymbol{B} = \boldsymbol{{}^{\hbox{M}}A} \cdot \boldsymbol{{}^{\hbox{M}}B}</math> and <math display="inline">\mathcal{A}:\boldsymbol{B} = \mathcal{{}^{\hbox{M}}A} \cdot \boldsymbol{{}^{\hbox{M}}B}</math>. With the choice of Mandel notation the Hill tensor can be displayed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| 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 \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)
|}

with the anisotropy parameters <math display="inline">a_i</math> which are based on the Hill parameters <math display="inline">h_i</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{lll}a_1 = \dfrac{1}{h^2_2} + \dfrac{1}{h^2_3} - \dfrac{1}{h^2_1} \quad &a_2 = \dfrac{1}{h^2_3} + \dfrac{1}{h^2_1} - \dfrac{1}{h^2_2} \quad &a_3 = \dfrac{1}{h^2_1} + \dfrac{1}{h^2_2} - \dfrac{1}{h^2_3}  \\ a_4 = \dfrac{1}{h^2_6} \quad &a_5 = \dfrac{1}{h^2_5} \quad   &a_6 = \dfrac{1}{h^2_4} \, . \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

'''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&nbsp;[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&nbsp;[Hockett1975],

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathop{}\!\sigma _{\hbox{Y}}(\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}) = \sigma _\infty + [ \sigma _0 - \sigma _\infty ] \exp \left(a \, {\mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}^b \right)  \, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

where <math display="inline">\sigma _0</math> and <math display="inline">\sigma _\infty </math> are stress-like quantities and <math display="inline">a</math> and <math display="inline">b</math> 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&nbsp;[Armstrong1966] nonlinear kinematic hardening through the evolution of the backstress by

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\boldsymbol{B}} = \dot \gamma \left[k_1 \boldsymbol{N} - k_2 \boldsymbol{B} \right]\, ,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

with the material parameter <math display="inline">k_1</math> representing the kinematic hardening modulus and <math display="inline">k_2</math> 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

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{B}_{n+1} = \frac{1}{\mathop{}\!\Delta \gamma _{n+1} + k_2} \left[k_1 \boldsymbol{N}_{n+1} +  \mathop{}\!\Delta \gamma _{n+1} \boldsymbol{B}_n \right]  \, .   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

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>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\boldsymbol{E}}^{\mathrm{pl}} = \dot \gamma \boldsymbol{N}     \, ,     \quad      \hbox{resp.}     \quad      \boldsymbol{E}_{n+1}^{\mathrm{pl}} = \boldsymbol{E}_{n}^{\mathrm{pl}}  + \mathop{}\!\Delta \gamma _{n+1} \boldsymbol{N}_{n+1}      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

===2.3 Single equation return-mapping===

Using a radial return-mapping algorithm (see e.g.&nbsp;[Neto:1252726]) and starting off with the trial stress as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\widetilde{\boldsymbol{T}}_{n+1} = \mathcal{C}^{\mathrm{el}} : \widetilde{\boldsymbol{E}}_{n+1}^{\mathrm{el}} = \mathcal{C}^{\mathrm{el}} : \left[\boldsymbol{E}_{n+1} - \boldsymbol{E}_{n}^{\mathrm{pl}} \right]\, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">\widetilde{\left[\cdot \right]}</math> indicates a trial value, the updated stress then reads

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{T}_{n+1} = \widetilde{\boldsymbol{T}}_{n+1} - \mathop{}\!\Delta \gamma _{n+1}  \mathcal{C}^{\mathrm{el}} : \boldsymbol{N}_{n+1} \, .   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Combining eqn:stress_update with eqn:kinematic_hardening_analy into the relative stress gives

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\xi }_{n+1} = \boldsymbol{T}_{n+1} - \boldsymbol{B}_{n+1}   </math>
|-
| style="text-align: center;" | <math>   =   \widehat{\boldsymbol{\xi }}_{n+1} - \frac{\mathop{}\!\Delta \gamma _{n+1}}{\left\vert \left\vert \boldsymbol{\xi }_{n+1} \right\vert \right\vert_{\mathcal{H}}}  \mathcal{B}_{n+1} : \mathcal{H}  : \boldsymbol{\xi }_{n+1}   \, ,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where the definitions for the stress tensor <math display="inline">\widehat{\boldsymbol{\xi }}_{n+1}</math> and the fourth-order tensor <math display="inline">\mathcal{B}_{n+1}</math>,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\widehat{\boldsymbol{\xi }}_{n+1} = \widetilde{\boldsymbol{T}}_{n+1} - \frac{1}{1+k_2 \mathop{}\!\Delta \gamma _{n+1}} \boldsymbol{B}_{n}  \quad   \hbox{ and }  \quad   \mathcal{B}_{n+1} = \mathcal{C}^{\mathrm{el}} + \frac{k_1}{1+k_2 \mathop{}\!\Delta \gamma _{n+1}} \mathop{}\!\mathcal{I}^{\mathrm{\,sym}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

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;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\xi }_{n+1}   = \left[\mathop{}\!\mathcal{I}^{\mathrm{\,sym}} + \frac{\mathop{}\!\Delta \gamma _{n+1}}{\mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\sigma _{\hbox{Y}}}  \mathcal{B}_{n+1} : \mathcal{H} \right]^{-1} : \widehat{\boldsymbol{\xi }}_{n+1}   =  \left[\mathcal{A}_{n+1} \right]^{-1} : \widehat{\boldsymbol{\xi }}_{n+1}   \, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

hence being solely a function of the unknown plastic multiplier <math display="inline">\mathop{}\!\Delta \gamma </math>. 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 <math display="inline">\mathop{}\!\Delta \gamma </math>  a Newton-Raphson scheme is applied, giving

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| 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 \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)
|}

with <math display="inline">k</math> 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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

with the auxiliary variables

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\kappa = \frac{1}{1 + k_2 \mathop{}\!\Delta \gamma }  \quad   \hbox{and}  \quad   \eta = \frac{\mathop{}\!\Delta \gamma }{\sqrt{2/3} \, \sigma _{\mathrm{Y}}}  \, . </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

===2.4 Consistent elasto-plastic tangent modulus===

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathop{}\!\mathrm{d}\boldsymbol{\xi } = \mathop{}\!\mathrm{d}\boldsymbol{T} - \mathop{}\!\mathrm{d}\boldsymbol{B}   </math>
|-
| style="text-align: center;" | <math>   \mathop{}\!\mathrm{d}\boldsymbol{T} = \mathcal{C}^{\mathrm{el}} : \left[\mathop{}\!\mathrm{d}\boldsymbol{E} - \mathop{}\!\mathrm{d}\boldsymbol{E}^{\mathrm{pl}} \right]   </math>
|-
| style="text-align: center;" | <math>   \mathop{}\!\mathrm{d}\boldsymbol{B} = \kappa \left[\left[k_1 \boldsymbol{N} - k_2 \boldsymbol{B} \right]\mathop{}\!\mathrm{d}\mathop{}\!\Delta \gamma + k_1 \mathop{}\!\Delta \gamma \frac{\partial \boldsymbol{N}}{\partial \boldsymbol{\xi }} : \mathop{}\!\mathrm{d}\boldsymbol{\xi } \right]   </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}\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)
|}

Solving the linearized system for <math display="inline">\mathop{}\!\mathrm{d}\boldsymbol{T} / \mathop{}\!\mathrm{d}\boldsymbol{E}</math> leads to the desired tangent operator

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C}^{\mathrm{ep}}  = \frac{\mathop{}\!\mathrm{d}\boldsymbol{T}}{\mathop{}\!\mathrm{d}{\boldsymbol{E}}}   =   \mathcal{C}_1 :         \left[   \left[\mathcal{C}^{\mathrm{el}} \right]^{-1} :    \mathcal{C}_1    +    \mathcal{C}_2      \right]^{-1}          </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

with the fourth-order auxiliary tensors

<span id="eq-23"></span>
<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C}_1 =   - \mathop{}\!\sqrt{\frac{2}{3}}k_1 \kappa \mathop{}\!\Delta \gamma \mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \Delta \gamma }\frac{\partial \boldsymbol{N}}{\partial \boldsymbol{\xi }}    + k_2 \kappa \boldsymbol{B} \otimes \boldsymbol{N}   - k_1 \kappa \boldsymbol{N} \otimes \boldsymbol{N}   - \mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \Delta \gamma }\mathop{}\!\mathcal{I}^{\mathrm{\,sym}}            </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|-
| style="text-align: center;" | <math>   \mathcal{C}_2 =        -\mathop{}\!\sqrt{\frac{2}{3}}\mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \Delta \gamma }\mathop{}\!\Delta \gamma \frac{\partial \boldsymbol{N}}{\partial \boldsymbol{\xi }}   -  \boldsymbol{N} \otimes \boldsymbol{N}      \, .   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}
|}

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 <math display="inline">\left[\cdot \right]_{n+1}</math> has been suppressed to improve readability.

The continuum tangent modulus is defined by the rate relation

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\boldsymbol{T}} = \mathcal{C}^{\mathrm{ep}}_{\mathrm{c}} : \dot{\boldsymbol{E}}  =   \mathcal{C}^{\mathrm{el}} : \left[\dot{\boldsymbol{E}} - \dot{\boldsymbol{E}}^{\mathrm{pl}} \right]  \, ,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

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;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\boldsymbol{B}} = \frac{ \left[k_1 \boldsymbol{N} - k_2 \boldsymbol{B} \right]: \boldsymbol{N} }{ \frac{2}{3} \mathop{}\!\frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}}+ \left[k_1 \boldsymbol{N} - k_2 \boldsymbol{B} \right]: \boldsymbol{N} } \dot{\boldsymbol{T}}   \, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C}^{\mathrm{ep}}_{\mathrm{c}}  = \frac{\dot{\boldsymbol{T}}}{\dot{\boldsymbol{E}}}   = \left[\left[\mathcal{C}^{\mathrm{el}} \right]^{-1} + \frac{1}{\frac{2}{3} \frac{\partial \sigma _{\mathrm{Y}}}{\partial \mathop{}\!\bar{\varepsilon }^{\mathrm{pl}}} +  \left[k_1 \boldsymbol{N} - k_2 \boldsymbol{B} \right]: \boldsymbol{N} }   \boldsymbol{N} \otimes \boldsymbol{N} \right]^{-1}   \, . </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

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%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Single-equation return-mapping algorithm
|-
|   ~&nbsp;<math>\blacksquare </math>~&nbsp;Geometric preprocessing: <math display="inline">\boldsymbol{C} = \boldsymbol{F}^\mathsf{T}\boldsymbol{F}</math> , <math display="inline">\boldsymbol{E}^{0}_{n+1} = \frac{1}{2} \ln (\boldsymbol{C})</math>  
|-
| ~&nbsp;<math>\blacksquare </math>~&nbsp;Plasticity algorithm:
|-
| ~&nbsp;~&nbsp;&nbsp;~&nbsp;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> 
|-
| ~&nbsp;~&nbsp;&nbsp;~&nbsp;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>
|-
| ~&nbsp;&nbsp; :  '''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>
:: 
:'''else'''
::
:'''end'''
:   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>|\varphi ^k| > \varepsilon _{\mathrm{tol}}</math> '''do'''
:: <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 \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>{\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>\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>
::  
:'''end'''
: <math>\boldsymbol{B}_{n+1} = \kappa _{n+1} \left[k_1 \boldsymbol{N}_{n+1} +  \mathop{}\!\Delta \gamma _{n+1} \boldsymbol{B}_n \right]</math> , <math>\boldsymbol{T}_{n+1} = \boldsymbol{\xi }_{n+1} + \boldsymbol{B}_{n+1}</math>
: <math>\mathcal{C}^{\mathrm{ep}}_{n+1}</math> from eqn:Materialtangente_kin_1,eqn:Materialtangente_kin_2,eqn:Materialtangente_kin_3
:     


|-
| ~&nbsp;<math>\blacksquare </math>~&nbsp;Geometric postprocessing: <math display="inline">\boldsymbol{S} = \boldsymbol{T} : \mathcal{P}</math> , <math display="inline">\mathcal{C}^{\mathrm{ep}}_{\mathrm{L}} = \mathcal{P}^\mathsf{T}: \mathcal{C}^{\mathrm{ep}} : \mathcal{P} + \boldsymbol{T} : \mathcal{L}</math>

|}

'''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.

==3 PARAMETER IDENTIFICATION==

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.

===3.1 Experimental and numerical setup===

The considered material is a DP600 sheet steel with 2mm thickness. In the experimental investigation a cruciform specimen with the dimensions depicted in  fig:cruciform_dimensions is used. To obtain the deformation field a speckle pattern which is captured with a stereo camera setup is applied to the surface of the specimen.   <div id='img-1a'></div>
<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Söhngen_975815040-biax-mit6mmBohrung_edit2.png|480px|Dimensions of the cruciform specimen in \SImm]]
|[[Image:draft_Söhngen_975815040-FEmesh.png|400px|Finite element mesh of specimen (dots mark optimization nodes)]]
|- style="text-align: center; font-size: 75%;"
| (a) Dimensions of the cruciform specimen in <math>\mathrm{mm}</math>
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''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&nbsp;fig:cruciform_strain_states major vs. minor strain is plotted for the optimization nodes along with the strain states for  uniaxial compression (<math display="inline">\varepsilon _1 = -0.5 \varepsilon _2</math>), pure shear (<math display="inline">\varepsilon _1 = - \varepsilon _2</math>), uniaxial tension (<math display="inline">\varepsilon _1 = -2 \varepsilon _2</math>),  plane strain (<math display="inline">\varepsilon _2 = 0</math>) and equi-biaxial tension (<math display="inline">\varepsilon _1 = \varepsilon _2</math>).  <div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Söhngen_975815040-picture-663b19.png|600px|Distribution of major vs. minor strain in the specimen during loading]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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&nbsp;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%;"
|-
|[[Image:draft_Söhngen_975815040-picture-f834aa.png|144px|]]
|[[Image:draft_Söhngen_975815040-picture-f6bba8.png|144px|]]
|[[Image:draft_Söhngen_975815040-picture-dff45b.png|144px|]]
|[[Image:draft_Söhngen_975815040-picture-11401c.png|144px|]]
|-
| colspan="4"|[[Image:draft_Söhngen_975815040-picture-e80114.png|600px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="4" | '''Figure 3:''' 
|}

===3.2 Inverse Identification===

The identification process is based on a comparison of the displacements obtained from experiment and FE simulation. The minimization function can be formulated as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol{\alpha }} f(\boldsymbol{\alpha }) = \min _{\boldsymbol{\alpha }} \frac{1}{2} \left\vert \left\vert \boldsymbol{x}_{\mathrm{Sim}} (\boldsymbol{\alpha }) - \boldsymbol{P} \boldsymbol{x}_{\mathrm{Exp}} \right\vert \right\vert ^{2} _{2}   \, , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

with <math display="inline">\boldsymbol{x}_{\mathrm{Sim}}</math> denoting the nodal displacements from the FE-calculation and <math display="inline">\boldsymbol{x}_{\mathrm{Exp}}</math> 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 <math display="inline">\boldsymbol{P}</math> 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 <math display="inline">\boldsymbol{\alpha }</math> contains the parameters to be optimized, in the present case it consists of the   four parameters describing the mixed hardening behavior, <math display="inline">\boldsymbol{\alpha } = \left[\sigma _0,\sigma _\infty ,k_1,k_2 \right]^\mathsf{T}</math>. 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 <math display="inline">f(\boldsymbol{\alpha })</math>, 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.&nbsp;[Schmaltz2013].  The application of a finite difference procedure to evaluate the  sensitivity of the cost function w.r.t. the design parameters <math display="inline">\boldsymbol{\alpha }</math> seems a necessary evil. An attempt to an analytical derivation for elasto-plastic material behaviour  was made by Cooreman_2007&nbsp;[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=<math>\mathrm{d}</math>]1.e-6. The parameters for the anisotropy of the yield surface <math display="inline">h_i</math> for the considered material have been identified in [Landkammerinpress] and are  listed in tab:DP600_Hillparameter.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Hill parameters for DP600 [Landkammerinpress]
|-
| <math display="inline">h_1</math> 
| <math>h_2</math>
| <math>h_3</math>
| <math>h_4</math>
| <math>h_5</math>
| <math>h_6</math>
|-
| <math display="inline">0.9759</math> 
| <math>1.0039</math>
| <math>0.9714</math>
| <math>1.0571</math>
| <math>1.0000</math>
| <math>1.0000</math>

|}

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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ 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 _\infty </math>in MPa 
| <math>k_1</math>
| <math>k_2</math>
|-
|  Set 1 
| <math>3.5x10^{2}</math>
| <math>1.5x10^{3}</math>
| <math>1x10^{3}</math>
| <math>1x10^{1}</math>
|-
| Set 2 
| <math>3.5x10^{2}</math>
| <math>1.5x10^{3}</math>
| <math>5x10^{3}</math>
| <math>1x10^{2}</math>
|-
| Set 3 
| <math>4x10^{2}</math>
| <math>3x10^{3}</math>
| <math>1x10^{3}</math>
| <math>1x10^{2}</math>
|-
|    Mean of optimization 
| <math>2.6322x10^{2}</math>
| <math>8.217x10^{2}</math>
| <math>1.74349x10^{4}</math>
| <math>1.74102x10^{2}</math>
|-
| SD 
| <math>3.30969x10^{-2}</math>
| <math>2.85087x10^{-2}</math>
| <math>9.03331</math>
| <math>5.51514x10^{-2}</math>

|}

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>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Söhngen_975815040-picture-a4f73d.png|600px|Evolution of isotropic hardening w.r.t. the equivalent plastic strain]]
|[[Image:draft_Söhngen_975815040-picture-a47046.png|400px|Evolution of backstress in loading direction w.r.t. the equivalent plastic strain]]
|- style="text-align: center; font-size: 75%;"
| (a) Evolution of isotropic hardening w.r.t. the equivalent plastic strain
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Evolution of backstress in loading direction w.r.t. the equivalent plastic strain
|}

==4 CONCLUSION AND OUTLOOK==

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).