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$^{*}$ 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

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

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

${\boldsymbol {E}}^{0}={\frac {1}{2}}\ln {\boldsymbol {C}}\,,$

(1)

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

${\boldsymbol {S}}={\boldsymbol {T}}:{\mathcal {P}}\quad {\hbox{and}}\quad {\mathcal {C}}_{\mathrm {L} }^{\mathrm {ep} }={\mathcal {P}}^{\mathsf {T}}:{\mathcal {C}}^{\mathrm {ep} }:{\mathcal {P}}+{\boldsymbol {T}}:{\mathcal {L}}$

(2)

via the fourth and sixth-order projection tensors ${\textstyle {\mathcal {P}}}$ and ${\textstyle {\mathcal {L}}}$ that are defined as

${\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}}}}\,.$

(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. ${\textstyle {\boldsymbol {C}}}$, as shown in [Miehe2001].

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

$\varphi ={\mathop {}}\!\left\vert \left\vert {\boldsymbol {\xi }}\right\vert \right\vert _{\mathcal {H}}-{\mathop {}}\!{\sqrt {\frac {2}{3}}}{\mathop {}}\!\sigma _{\hbox{Y}}$

(4)

where the tensor norm ${\textstyle {\mathop {}}\!\left\vert \left\vert {\boldsymbol {\xi }}\right\vert \right\vert _{\mathcal {H}}}$ is given by

${\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$

(5)

with the relative stress tensor ${\textstyle {\boldsymbol {\xi }}}$ being defined by the stress tensor ${\textstyle {\boldsymbol {T}}}$ and the backstress tensor ${\textstyle {\boldsymbol {B}}}$

${\boldsymbol {\xi }}={\boldsymbol {T}}-{\boldsymbol {B}}\,.$

(6)

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

${\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}}$

(7)

with the anisotropy parameters ${\textstyle a_{i}}$ which are based on the Hill parameters ${\textstyle h_{i}}$

${\begin{array}{lll}a_{1}={\dfrac {1}{h_{2}^{2}}}+{\dfrac {1}{h_{3}^{2}}}-{\dfrac {1}{h_{1}^{2}}}\quad &a_{2}={\dfrac {1}{h_{3}^{2}}}+{\dfrac {1}{h_{1}^{2}}}-{\dfrac {1}{h_{2}^{2}}}\quad &a_{3}={\dfrac {1}{h_{1}^{2}}}+{\dfrac {1}{h_{2}^{2}}}-{\dfrac {1}{h_{3}^{2}}}\\a_{4}={\dfrac {1}{h_{6}^{2}}}\quad &a_{5}={\dfrac {1}{h_{5}^{2}}}\quad &a_{6}={\dfrac {1}{h_{4}^{2}}}\,.\end{array}}$

(8)

Remark The tensor ${\textstyle {\mathcal {P}}}$ has the property ${\textstyle {\mathcal {P}}:{\boldsymbol {I}}={\boldsymbol {0}}}$ of a deviatoric projection tensor [Miehe2002] and specifically for ${\textstyle h_{i}=1}$ reduces to the deviatoric identity tensor ${\textstyle {\mathcal {I}}^{\mathrm {dev} }={\mathcal {I}}^{\mathrm {sym} }-{\frac {1}{3}}{\boldsymbol {I}}\otimes {\boldsymbol {I}}}$ which yields the classical ${\textstyle \mathrm {J_{2}} }$-plasticity for the tensor norm, ${\textstyle {\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 }$ .

Considering isotropic hardening, the yield stress in eqn:yield_function is represented by the exponential law according to Hockett1975 [Hockett1975],

${\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)\,,$

(9)

where ${\textstyle \sigma _{0}}$ and ${\textstyle \sigma _{\infty }}$ are stress-like quantities and ${\textstyle a}$ and ${\textstyle b}$ 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

${\dot {\boldsymbol {B}}}={\dot {\gamma }}\left[k_{1}{\boldsymbol {N}}-k_{2}{\boldsymbol {B}}\right]\,,$

(10)

with the material parameter ${\textstyle k_{1}}$ representing the kinematic hardening modulus and ${\textstyle k_{2}}$ 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

${\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]\,.$

(11)

Declaring the flow rule to be of an associative type, the flow vector reads ${\textstyle {\boldsymbol {N}}=\partial {\text{/}}\varphi {\partial {\boldsymbol {T}}}}$ and the evolution of the plastic strain tensor is defined through

${\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}$

(12)

2.3 Single equation return-mapping

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

${\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]\,,$

(13)

where ${\textstyle {\widetilde {\left[\cdot \right]}}}$ indicates a trial value, the updated stress then reads

${\boldsymbol {T}}_{n+1}={\widetilde {\boldsymbol {T}}}_{n+1}-{\mathop {}}\!\Delta \gamma _{n+1}{\mathcal {C}}^{\mathrm {el} }:{\boldsymbol {N}}_{n+1}\,.$

(14)

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

${\boldsymbol {\xi }}_{n+1}={\boldsymbol {T}}_{n+1}-{\boldsymbol {B}}_{n+1}$ $={\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}\,,$

(15)

where the definitions for the stress tensor ${\textstyle {\widehat {\boldsymbol {\xi }}}_{n+1}}$ and the fourth-order tensor ${\textstyle {\mathcal {B}}_{n+1}}$,

${\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} }$

(16)

are introduced. Since in the converged step ${\textstyle \varphi =0}$ and thus ${\textstyle {\mathop {}}\!\left\vert \left\vert {\boldsymbol {\xi }}\right\vert \right\vert _{\mathcal {H}}={\sqrt {\text{2/3}}}\,{\mathop {}}\!\sigma _{\hbox{Y}}}$ holds, eqn:relative_stress can be rewritten as

${\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}\,,$

(17)

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

${\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}\,,$

(18)

with ${\textstyle k}$ 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

${\frac {\partial \varphi }{\partial {\mathop {}}\!\Delta \gamma }}={\boldsymbol {N}}_{n+1}:\left[{\mathcal {A}}_{n+1}\right]^{-1}:\left\{k_{2}\kappa _{n+1}^{2}{\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.$ $\left.\left.-\eta _{n+1}k_{1}k_{2}\kappa _{n+1}^{2}{\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 }}\,,$

(19)

with the auxiliary variables

$\kappa ={\frac {1}{1+k_{2}{\mathop {}}\!\Delta \gamma }}\quad {\hbox{and}}\quad \eta ={\frac {{\mathop {}}\!\Delta \gamma }{{\sqrt {2/3}}\,\sigma _{\mathrm {Y} }}}\,.$

(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

${\mathop {}}\!\mathrm {d} {\boldsymbol {\xi }}={\mathop {}}\!\mathrm {d} {\boldsymbol {T}}-{\mathop {}}\!\mathrm {d} {\boldsymbol {B}}$ ${\mathop {}}\!\mathrm {d} {\boldsymbol {T}}={\mathcal {C}}^{\mathrm {el} }:\left[{\mathop {}}\!\mathrm {d} {\boldsymbol {E}}-{\mathop {}}\!\mathrm {d} {\boldsymbol {E}}^{\mathrm {pl} }\right]$ ${\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]$ ${\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$ ${\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\,.$

(21)

Solving the linearized system for ${\textstyle {\mathop {}}\!\mathrm {d} {\boldsymbol {T}}/{\mathop {}}\!\mathrm {d} {\boldsymbol {E}}}$ leads to the desired tangent operator

${\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}$

(22)

with the fourth-order auxiliary tensors

${\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} }$ (23) ${\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}}\,.$ (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 ${\textstyle \left[\cdot \right]_{n+1}}$ has been suppressed to improve readability.

The continuum tangent modulus is defined by the rate relation

${\dot {\boldsymbol {T}}}={\mathcal {C}}_{\mathrm {c} }^{\mathrm {ep} }:{\dot {\boldsymbol {E}}}={\mathcal {C}}^{\mathrm {el} }:\left[{\dot {\boldsymbol {E}}}-{\dot {\boldsymbol {E}}}^{\mathrm {pl} }\right]\,,$

(25)

where ${\textstyle \left[\cdot \right]_{\mathrm {c} }}$ denotes the continuum operator. Under usage of the consistency condition ${\textstyle {\dot {\varphi }}=0}$ an expression for the rate of the backstress can be derived as

${\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}}}\,,$

(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

${\mathcal {C}}_{\mathrm {c} }^{\mathrm {ep} }={\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}\,.$

(27)

Calculating the limit case ${\textstyle {\mathop {}}\!\Delta \gamma \to 0}$ for the consistent tangent operator leads, after some rearrangement and by using ${\textstyle {\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 }}}$, to the continuum tangent operator as expected, ${\textstyle \lim \limits _{{\mathop {}}\!\Delta \gamma \to 0}{\mathcal {C}}^{\mathrm {ep} }={\mathcal {C}}_{\mathrm {c} }^{\mathrm {ep} }}$. 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 ${\textstyle {\mathcal {C}}^{\mathrm {ep} }}$, leading to a non-symmetric representation in Mandel notation ${\textstyle {\mathcal {{}^{\hbox{M}}C}}^{\mathrm {ep} }\neq [{\mathcal {{}^{\hbox{M}}C}}^{\mathrm {ep} }]^{\mathsf {T}}}$. The algorithmic equations are listed in tab:algorithm for convenience.

Table. 1 Single-equation return-mapping algorithm

~ $\blacksquare$~ Geometric preprocessing: ${\textstyle {\boldsymbol {C}}={\boldsymbol {F}}^{\mathsf {T}}{\boldsymbol {F}}}$ , ${\textstyle {\boldsymbol {E}}_{n+1}^{0}={\frac {1}{2}}\ln({\boldsymbol {C}})}$

~ $\blacksquare$~ Plasticity algorithm:

~ ~  ~ Trial stress: ${\textstyle {\widetilde {\boldsymbol {T}}}_{n+1}={\mathcal {C}}^{\mathrm {el} }:\left[{\boldsymbol {E}}_{n+1}-{\boldsymbol {E}}_{n}^{\mathrm {pl} }\right]}$

~ ~  ~ Evaluate yield function: ${\textstyle \varphi =\lVert {\widetilde {\boldsymbol {\xi }}}_{n+1}\rVert _{\mathcal {H}}-{\mathop {}}\!{\sqrt {\frac {2}{3}}}{\mathop {}}\!\sigma _{\hbox{Y}}({\mathop {}}\!{\bar {\varepsilon }}_{n}^{\mathrm {pl} })}$ , ${\textstyle {\widetilde {\boldsymbol {\xi }}}_{n+1}={\widetilde {\boldsymbol {T}}}_{n+1}-{\boldsymbol {B}}_{n}}$

~   : if $\varphi \leq 0.0{\text{d}}0$ then Elastic step Set: ${\boldsymbol {T}}_{n+1}={\widetilde {\boldsymbol {T}}}_{n+1}$ , ${\boldsymbol {B}}_{n+1}={\boldsymbol {B}}_{n}$ , ${\mathcal {C}}_{n+1}^{\mathrm {ep} }={\mathcal {C}}^{\mathrm {el} }$ else end Plastic step Set: $k=0$ , ${\mathop {}}\!\Delta \gamma _{n+1}^{0}=0.0{\text{d}}0$ , ${\boldsymbol {\xi }}_{n+1}^{0}={\widetilde {\boldsymbol {\xi }}}_{n+1}$ , $\varepsilon _{\mathrm {tol} }=1.0{\text{d}}-8$ while $|\varphi ^{k}|>\varepsilon _{\mathrm {tol} }$ do $\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}}$ ${\mathop {}}\!\Delta \gamma _{n+1}^{k+1}={\mathop {}}\!\Delta \gamma _{n+1}^{k}-\left[\left.{\frac {\partial \varphi }{\partial {\mathop {}}\!\Delta \gamma }}\right\vert ^{k}\right]^{-1}\varphi ^{k}$ ${\widehat {\boldsymbol {\xi }}}^{\,k+1}={\widetilde {\boldsymbol {T}}}_{n+1}-\kappa ^{k+1}{\boldsymbol {B}}_{n}$ , ${\boldsymbol {\xi }}_{n+1}^{k+1}=\left[{\mathcal {A}}\right]^{-1}:{\widehat {\boldsymbol {\xi }}}^{\,k+1}$ ${{\mathop {}}\!{\bar {\varepsilon }}^{\mathrm {pl} }}_{n+1}^{k+1}={\mathop {}}\!{\bar {\varepsilon }}_{n}^{\mathrm {pl} }+{\mathop {}}\!{\sqrt {\frac {2}{3}}}{\mathop {}}\!\Delta \gamma _{n+1}^{k+1}$ , ${\mathop {}}\!\sigma _{\hbox{Y}}^{k+1}={\mathop {}}\!\sigma _{\hbox{Y}}({{\mathop {}}\!{\bar {\varepsilon }}^{\mathrm {pl} }}_{n+1}^{k+1})$ $\varphi ^{k+1}=\left\vert \left\vert {\boldsymbol {\xi }}_{n+1}^{k+1}\right\vert \right\vert _{\mathcal {H}}-{\mathop {}}\!{\sqrt {\frac {2}{3}}}{\mathop {}}\!\sigma _{\hbox{Y}}^{k+1}$ $k=k+1$ end ${\boldsymbol {B}}_{n+1}=\kappa _{n+1}\left[k_{1}{\boldsymbol {N}}_{n+1}+{\mathop {}}\!\Delta \gamma _{n+1}{\boldsymbol {B}}_{n}\right]$ , ${\boldsymbol {T}}_{n+1}={\boldsymbol {\xi }}_{n+1}+{\boldsymbol {B}}_{n+1}$ ${\mathcal {C}}_{n+1}^{\mathrm {ep} }$ from eqn:Materialtangente_kin_1,eqn:Materialtangente_kin_2,eqn:Materialtangente_kin_3

~ $\blacksquare$~ Geometric postprocessing: ${\textstyle {\boldsymbol {S}}={\boldsymbol {T}}:{\mathcal {P}}}$ , ${\textstyle {\mathcal {C}}_{\mathrm {L} }^{\mathrm {ep} }={\mathcal {P}}^{\mathsf {T}}:{\mathcal {C}}^{\mathrm {ep} }:{\mathcal {P}}+{\boldsymbol {T}}:{\mathcal {L}}}$

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.

(a) Dimensions of the cruciform specimen in $\mathrm {mm}$
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 (${\textstyle \varepsilon _{1}=-0.5\varepsilon _{2}}$), pure shear (${\textstyle \varepsilon _{1}=-\varepsilon _{2}}$), uniaxial tension (${\textstyle \varepsilon _{1}=-2\varepsilon _{2}}$), plane strain (${\textstyle \varepsilon _{2}=0}$) and equi-biaxial tension (${\textstyle \varepsilon _{1}=\varepsilon _{2}}$).

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 (${\textstyle |F_{\mathrm {x/y} }|=8{\mbox{kN}}}$, ${\textstyle |F_{\mathrm {x/y} }|=9{\mbox{kN}}}$ and ${\textstyle |F_{\mathrm {x/y} }|=1x10^{1}{\mbox{kN}}}$) 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 ${\textstyle |{\dot {F}}_{\mathrm {x/y} }|=1x10^{2}{\mbox{N/s}}}$ is chosen.

File:Draft Söhngen 975815040-picture-f834aa.png	File:Draft Söhngen 975815040-picture-f6bba8.png	File:Draft Söhngen 975815040-picture-dff45b.png	File:Draft Söhngen 975815040-picture-11401c.png
File:Draft Söhngen 975815040-picture-e80114.png
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

$\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}\,,$

(28)

with ${\textstyle {\boldsymbol {x}}_{\mathrm {Sim} }}$ denoting the nodal displacements from the FE-calculation and ${\textstyle {\boldsymbol {x}}_{\mathrm {Exp} }}$ 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 ${\textstyle {\boldsymbol {P}}}$ 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 ${\textstyle {\boldsymbol {\alpha }}}$ contains the parameters to be optimized, in the present case it consists of the four parameters describing the mixed hardening behavior, ${\textstyle {\boldsymbol {\alpha }}=\left[\sigma _{0},\sigma _{\infty },k_{1},k_{2}\right]^{\mathsf {T}}}$. 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 ${\textstyle f({\boldsymbol {\alpha }})}$, 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 ${\textstyle {\boldsymbol {\alpha }}}$ 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=$\mathrm {d}$]1.e-6. The parameters for the anisotropy of the yield surface ${\textstyle h_{i}}$ for the considered material have been identified in [Landkammerinpress] and are listed in tab:DP600_Hillparameter.

Table. 2 Hill parameters for DP600 [Landkammerinpress]

${\textstyle h_{1}}$	$h_{2}$	$h_{3}$	$h_{4}$	$h_{5}$	$h_{6}$
${\textstyle 0.9759}$	$1.0039$	$0.9714$	$1.0571$	$1.0000$	$1.0000$

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 ${\textstyle E=210000}$ and a Poisson's ratio of ${\textstyle \nu =0.3}$. 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.

Table. 3 Initial and identified parameters of the mixed hardening formulation for DP600

	$\sigma _{0}$in MPa	$\sigma _{\infty }$in MPa	$k_{1}$	$k_{2}$
Set 1	$3.5x10^{2}$	$1.5x10^{3}$	$1x10^{3}$	$1x10^{1}$
Set 2	$3.5x10^{2}$	$1.5x10^{3}$	$5x10^{3}$	$1x10^{2}$
Set 3	$4x10^{2}$	$3x10^{3}$	$1x10^{3}$	$1x10^{2}$
Mean of optimization	$2.6322x10^{2}$	$8.217x10^{2}$	$1.74349x10^{4}$	$1.74102x10^{2}$
SD	$3.30969x10^{-2}$	$2.85087x10^{-2}$	$9.03331$	$5.51514x10^{-2}$

Observing a relative SD of ${\textstyle 5x10^{-4}}$ 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 (${\textstyle {\mathop {}}\!{\bar {\varepsilon }}^{\mathrm {pl} }=\int _{0}^{t}{\dot {\gamma }}{\mathop {}}\!\mathrm {d} t}$). The evolution of the backstress in loading direction w.r.t. ${\textstyle {\mathop {}}\!{\bar {\varepsilon }}^{\mathrm {pl} }}$ is depicted in fig:evolution_backstress. As can be seen the kinematic hardening modulus ${\textstyle k_{1}}$ leads to an initially steep incline while the parameter ${\textstyle k_{2}}$ causes the increase in backstress to saturate with advancing plastic strain for the identified parameters.

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

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