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


Institute of Applied Mechanics
Friedrich-Alexander-Universität Erlangen-Nürnberg
Egerlandstr. 5, 91058 Erlangen, Germany
web page:


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.


2.1 Notation

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 .

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


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


via the fourth and sixth-order projection tensors and that are defined as


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


where the tensor norm is given by


with the relative stress tensor being defined by the stress tensor and the backstress tensor


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


with the anisotropy parameters which are based on the Hill parameters


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],


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


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


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


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


where indicates a trial value, the updated stress then reads


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


where the definitions for the stress tensor and the fourth-order tensor ,


are introduced. Since in the converged step and thus holds, eqn:relative_stress can be rewritten as


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


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


with the auxiliary variables


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


Solving the linearized system for leads to the desired tangent operator


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


where denotes the continuum operator. Under usage of the consistency condition an expression for the rate of the backstress can be derived as


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


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.

Table. 1 Single-equation return-mapping algorithm
~ Geometric preprocessing: ,
~ Plasticity algorithm:
~ ~  ~ Trial stress:
~ ~  ~ Evaluate yield function: ,
~   : if then
Elastic step Set: , ,
Plastic step Set: , , ,
while do
from eqn:Materialtangente_kin_1,eqn:Materialtangente_kin_2,eqn:Materialtangente_kin_3

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

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.
Dimensions of the cruciform specimen in \SImm Finite element mesh of specimen (dots mark optimization nodes)
(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.
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


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.

Table. 2 Hill parameters for DP600 [Landkammerinpress]

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.

Table. 3 Initial and identified parameters of the mixed hardening formulation for DP600
in MPa in MPa
Set 1
Set 2
Set 3
Mean of optimization
Observing a relative SD of 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 (). The evolution of the backstress in loading direction w.r.t. is depicted in fig:evolution_backstress. As can be seen the kinematic hardening modulus leads to an initially steep incline while the parameter 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


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” (

Back to Top

Document information

Published on 24/05/17
Submitted on 17/05/17

Licence: Other

Document Score


Views 101
Recommendations 0

Share this document