B. SÖHNGEN AND K. WILLNER
Institute of Applied Mechanics 
FriedrichAlexanderUniversität ErlangenNürnberg 
Egerlandstr. 5, 91058 Erlangen, Germany 
email: benjamin.soehngen@fau.de, kai.willner@fau.de 
web page: www.ltm.fau.de 
In this work an anisotropic material model at finite strains with nonlinear mixed (isotropic and kinematic) hardening is used for the identification of the hardening parameters of sheet steel. The algorithmic system is thereby reduced to a single equation return mapping. For the identification, a cruciform specimen is loaded biaxially in an alternating shear test to provoke the kinematic hardening behavior and prevent the sheet from buckling. The material parameters are found through an optimization strategy by comparing the deformation field from the experiment to that from a finite element (FE) simulation. The resulting cost function is minimized by means of a gradientbased 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 sheetbulk 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 everincreasing complexity of the constitutive models demand for appropriate identification methods by simultaneously decreasing experimental effort.
In this contribution a material law capable of describing elastoplastic material behavior utilizing a Hilltype yield function in combination with a mixed hardening law is used to identify parameters for the dualphase steel DP600.
In this work, scalar quantities are represented by lower case letters . Boldfaced characters denote secondorder tensors; in particular denotes the secondorder identity tensor . Fourth and sixthorder tensors are defined by calligraphic symbols . The dyadic product between two secondorder tensors yield a tensor of fourthorder and can be expressed in index notation as . The fourthorder tensor is a linear map from to which is given as in this work (accordingly for sixthorder tensors). The fourthorder 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 secondorder tensor . The Euclidean norm of a secondorder tensor is defined as where the double dot product reads .
The formulation of large strain plasticity is based on the algorithm presented by Miehe2001 [Miehe2001]. It consists of the calculation of stresses and elasticity moduli with regard to a SethHill strain tensor. For the application to elastoplasticity, the algorithm can be interpreted as a pre and postprecessing of the stresses and elastoplastic material tangent originating from a “small strain” algorithm based on the Hencky (logarithmic) strain tensor

(1) 
where denotes the right CauchyGreen tensor, as shown by Miehe2002 [Miehe2002]. Referring to as the workconjugate stress measure to the logarithmic strain and as the consistent material tangent, the transformation to the Lagrangian stress and module is performed by

(2) 
via the fourth and sixthorder projection tensors and that are defined as

(3) 
The first and second derivative of the logarithmic strain tensor can be obtained by usage of a spectral decomposition of the right CauchyGreen 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 Hilltype anisotropic plasticity. The yield function therefore reads

(4) 
where the tensor norm is given by

(5) 
with the relative stress tensor being defined by the stress tensor and the backstress tensor

(6) 
The anisotropic plasticity is modeled by the fourthorder 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) secondorder tensors for stress and strainlike 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 secondorder tensors or the tensor product between a fourth and a secondorder tensor do not yield the same result in Voigt notation, and , but using the representation according to Mandel, and . With the choice of Mandel notation the Hill tensor can be displayed as

(7) 
with the anisotropy parameters which are based on the Hill parameters

(8) 
Remark The tensor has the property of a deviatoric projection tensor [Miehe2002] and specifically for reduces to the deviatoric identity tensor which yields the classical plasticity for the tensor norm, .
Considering isotropic hardening, the yield stress in eqn:yield_function is represented by the exponential law according to Hockett1975 [Hockett1975],

(9) 
where and are stresslike quantities and and describe the curvature of the yield stress. As the algorithm is not limited to this formulation, other well known hardening laws like the ones from Ludwik, Voce or Swift could be used as well.
The shift of the yield surface in stress space is considered by Armstrong1966 [Armstrong1966] nonlinear kinematic hardening through the evolution of the backstress by

(10) 
with the material parameter representing the kinematic hardening modulus and the rate of saturation. For usage in a finite element code the evolution law is discretized in time by an implicit Euler scheme, resulting in

(11) 
Declaring the flow rule to be of an associative type, the flow vector reads and the evolution of the plastic strain tensor is defined through

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

(13) 
where indicates a trial value, the updated stress then reads

(14) 
Combining eqn:stress_update with eqn:kinematic_hardening_analy into the relative stress gives

(15) 
where the definitions for the stress tensor and the fourthorder tensor ,

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

(17) 
hence being solely a function of the unknown plastic multiplier . By substitution into the definition of the yield function (eqn:yield_function) a scalar equation with only one unknown is derived. To solve the (nonlinear) equation for a NewtonRaphson scheme is applied, giving

(18) 
with the iteration counter. The derivative of the yield function w.r.t. the plastic multiplier can be displayed after straightforward usage of the chain rule as

(19) 
with the auxiliary variables

(20) 
In the algorithmic implementation, the related tangent modulus plays the major role for achieving (quadratic) convergence in the NewtonRaphson 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 elastoplastic tangent modulus. To develop the algorithmic tangent the system of equations (eqn:plastic_strain_tensor,eqn:kinematic_hardening_analy,eqn:relative_stress,eqn:stress_update,eqn:yield_function) is linearized, giving

(21) 
Solving the linearized system for leads to the desired tangent operator

(22) 
with the fourthorder auxiliary tensors

It should be noted that the consistent tangent modulus is calculated at the end of the stress update, hence all quantities are evaluated at the updated state and the index has been suppressed to improve readability.
The continuum tangent modulus is defined by the rate relation

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

(26) 
which  inserted back into eqn:continuum_tangent_rate_relation along with the consistency condition and the rate of the plastic strain tensor (eqn:plastic_strain_tensor)  yields the continuum tangent operator as

(27) 
Calculating the limit case for the consistent tangent operator leads, after some rearrangement and by using , to the continuum tangent operator as expected, . It is emphasized that although the chosen plasticity model is associative, the kinematic hardening law is of a nonassociative type. This leads to a nonsymmetric elastoplastic 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 nonsymmetric systems. The loss of symmetry thereby only refers to the major symmetry of , leading to a nonsymmetric representation in Mandel notation . The algorithmic equations are listed in tab:algorithm for convenience.
~ ~ Geometric preprocessing: , 
~ ~ Plasticity algorithm: 
~ ~ ~ Trial stress: 
~ ~ ~ Evaluate yield function: , 
~ : if then

~ ~ Geometric postprocessing: , 
Remark As already mentioned above, all equations can be represented in some compressed notation. Every fourthorder 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 electromechanical actuators of which the two of each axis are operated in master/slave control to keep the specimen centered and inhibit rigid body motion. To prevent buckling, the loading along the two axes is applied in an alternating fashion rather than by pure tension/compression. Therefore, the applied force is equal in amount but opposite in direction, resulting in an alternating shear loading.
(a) Dimensions of the cruciform specimen in  
Figure 1: Finite element mesh of specimen (dots mark optimization nodes) 
The captured images are subsequently analyzed with the digital image correlation (DIC) software Aramis to generate the displacements. For the numerical computation a FEmesh consisting of 3080 eightnoded 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 equibiaxial 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 aforementioned alternating shear loading of the specimen is depicted in fig:specimen_load. Starting off with the initial unloaded configuration the speckle patterns for the three peak loading values (, and ) are shown. The deformation state can be judged by the elliptic deformation of the originally circular bore in the middle. For the loading a linear force rate of is chosen.The identification process is based on a comparison of the displacements obtained from experiment and FE simulation. The minimization function can be formulated as

(28) 
with denoting the nodal displacements from the FEcalculation and the displacements coming from the experimental fullfield 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 FEmesh. The discretization introduced by DIC is thereby about ten times finer than that of the FEmesh, 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 gradientbased and gradientfree 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 gradientbased approaches should be chosen over gradientfree 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 elastoplastic 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 wellsuited and in the area of parameter optimization widely used (see e.g. [Chaparro_2008]) algorithm is the LevenbergMarquardt 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 [roundmode=places,roundprecision=1,outputexponentmarker=]1.e6. The parameters for the anisotropy of the yield surface for the considered material have been identified in [Landkammerinpress] and are listed in tab:DP600_Hillparameter.
For that identification the same cruciform specimen geometry was loaded with equibiaxial tension. In the elastic region the material is assumed to behave isotropically with an elastic modulus of and a Poisson's ratio of . The optimization of the mixed hardening parameters is based on three different initial values sets that are given in tab:DP600_kinematic_hard along with the mean and standard deviation (SD) of the three optimization runs.
in MPa  in MPa  
Set 1  
Set 2  
Set 3  
Mean of optimization  
SD 
Evolution of isotropic hardening w.r.t. the equivalent plastic strain  Evolution of backstress in loading direction w.r.t. the equivalent plastic strain 
(a) Evolution of isotropic hardening w.r.t. the equivalent plastic strain  
Figure 4: Evolution of backstress in loading direction w.r.t. the equivalent plastic strain 
In this contribution a single equation return mapping algorithm for large strain anisotropic plasticity with isotropic and kinematic hardening is presented. The identification of the mixed hardening parameters for the dualphase steel DP600 with a LevenbergMarquardt 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 FEpackages 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  SheetBulk Metal Forming” (www.tr73.de).
Published on 24/05/17
Submitted on 17/05/17
Licence: Other