Semi-analytic parameter identification for complex yield functions

This paper presents a new parameter identification scheme for complex yield criteria of sheet metals using experimentally determined load and strain distributions and combining analytical stress analysis on sectional strain data with common inverse analysis. The approach is suitable for specimen with a non-homogenous strain distribution in the specimen and reduces computing time compared to common methods like iteratively updated FEM considerably. It serves to identify a set of material parameters in a fast and precise manner, describing the material


Introduction
For sheet metal forming industry, the numerical simulation of manufacturing processes is a competitive factor as reliable prediction allows the industry to shorten cycles of process development and to reduce time to production launch. Forming simulation supports engineers and toolmakers in the prediction of forming behaviour, formability and final part properties during the process design. In order to comply with the progressive demands in reducing the discrepancy between numerical prediction and experimental results, an improvement of the phenomenological modelling of forming behaviour combined with efficient experimental testing is of high relevance [1].

State of the art
In metal forming, the plastic material behaviour is usually modelled by a flow curve, describing the true stress response as a function of true strain for uniaxial loading. In combination with a yield function the material behaviour under multiaxial loading is described.
To overcome the drawbacks of isotropic yield functions as Tresca [2] or von Mises [3] complex yield criteria were developed in order to describe anisotropic behaviour more accurately. For many steel grades Hill's quadratic yield criterion [4] is sufficient. An advantage of this criterion is the possibility of direct parameter calculation based on experimental results from only three standard tensile tests. However, Hill's criterion is not appropriate for materials like high strength steels or aluminium that have a distinctive anisotropy. Therefore, more complex yield functions like Yld2000-2D [5], BBC2000 [6] and Vegter [7] were developed in order to describe the yield surface shape more accurately and improve the results of numerical simulations. As these yield functions need additional points in stress space, extensive efforts for experimental testing and parameter identification are required [1]. Additional testing data of biaxial, plain strain and shear specimen is commonly used to model an accurate material response with these complex yield criteria.
Compared to standard experiments like tensile tests, the challenge for most of the additional experiments is the appearance of non-homogenous strain and stress distributions, which bias the results of experimental analysis. For this reason, some experiments are specifically designed to improve the results. In shear test configurations [8], the geometric dimension of specimen is adapted, so that rotational effects are reduced and the quasi-homogenous centre of the shear zone is enlarged. Within these tests, the shear stress is often approximated analytically by using the testing force and assuming the stress distribution to be uniform over the shear zone. But even with special specimen design, the shear stress distribution is quite inhomogeneous and therefore leads to imprecise results for the stress states. Hence, a direct analytical calculation of material parameters based on global measurement of complex specimen is not appropriate for an accurate description of the yield surface.
To cope with the problems associated with nonhomogenous strain distribution appearing in experimental testing, several approaches for parameter identification were proposed in the past. For instance, it is possible to transform the crystal anisotropy into a global anisotropy using crystal plasticity methods to obtain a phenomenological description of the yield surface [9]. Using this technique, inhomogeneity on grain scale is used to predict the macroscopic mechanical behaviour. A drawback of this method is, that only defects in the polycrystalline microstructure, but no further macroscopic defects of the sheet metal can be considered.
In contrast, identification methods based on kinematic full-field measurements and inverse analysis offer the possibility to model the material behaviour based on the observed macroscopic phenomena [10]. Within these approaches, a testing machine coupled to an optical strain measurement system records the experimental results as input for the identification scheme. Subsequently, possible strategies specially developed for parameter identification of plastic material behaviour are briefly described.
The most popular approach for parameter identification is finite element method updated (FEMU), using iterative finite element analysis within an inverse analysis [11][12][13][14][15]. The method uses the discrepancy between numerical results and experimentally determined loads and strain fields to define an objective function. Using an optimisation algorithm, numerical results with varying parameters are iteratively computed by finite element analysis until the objective function is minimised. The approach was utilised in diverse studies to determine material parameters for sheet metal forming processes. For instance, inverse analysis with iterative FEA was used by Güner et al. [11] to identify the yield function parameters of Yld2000-2d [5] using the FE-code Abaqus-Explicit and the Levenberg-Marquardt Algorithm. In a subsequent study of Aydın et al. [12], a parameter identification for the evolution of the yield surface Yld2000-2d [5] was performed using the optimisation tool Ls-Opt and the FE-code Ls-Dyna with experimental data of standard tensile tests and a plane strain tension test. Further studies calibrate the yield surface on experiments designed for the validation of specific forming conditions in order to get a processoriented material response. Bäck et al. [13] presented a parameter identification strategy based on evolutionary algorithms and its application in calibration of diverse yield criteria with specially designed validation tools for stretch forming conditions. Using specifically designed experiments in combination with FEMU, it is also possible to identify material parameters for yield surface and flow curve simultaneously with one single test [14]. Furthermore, the FEMU approach can be used for the identification of material parameters for kinematic hardening models. Yin et al. [15] proposed an in-plane torsion test and used Abaqus-Explicit in combination with a trust region method to identify kinematic hardening parameters. The presented studies show, that FEMU can be used as a computationally expensive, but very powerful tool for parameter identification.
In order to increase the computational efficiency, the application of the virtual fields method (VFM) for parameter identification in elasto-plastic problems was proposed by Avril et al. [16]. The algorithm uses fullfield surface measurements of the deformation to calculate a stress field with the predictor-corrector method of Sutton et al. [17]. Stress components are iteratively computed with varied and updated material parameters until the principle of virtual work used as optimisation criterion is satisfied. In the work of Avril et al. [16], four constitutive parameters for elastic and plastic material behaviour were identified on one single specimen assuming a von Mises yield criterion.
According to the study, the approach was able to reduce the computing time by a factor of 100 compared to the FEMU approach. This promising approach was adjusted in a subsequent study for anisotropic yield criteria [18]. In the study, two specially designed tests were used to retrieve six constitutive parameters of the Hill yield criterion [4] and the Swift hardening law [19].
FEMU and VFM show good results for adjusting a complex yield surface to experimental results with arbitrary geometry. But as finite element simulations are costly, the widely-used iterative FEM shows only poor efficiency regarding the computation time. In contrast, the virtual field method presents an approach with considerably reduced computational efforts. Nevertheless, the potential in stress analysis and identification is not thoroughly used yet. The direct analytical calculation of stress states is still used in practice due to its simplicity and computational efficiency compared to the presented approaches. Hence, there is a need for improving the efficiency and accuracy of (semi-)analytical yield surface characterisation regarding computing time and reduced input of experimental results.

Objective
To overcome the presented drawbacks of current parameter identification schemes, the development of an alternative identification strategy which significantly reduces computational costs is necessary. Within this strategy, field data, commonly recorded in standard setups of testing machines coupled with optical strain measurement, shall be used in an efficient way to identify material parameters. Unexploited potential in the combination of inverse analysis and direct analytical evaluation over cross-section area with simplified stress analysis shall be used in order to create a highly performant identification strategy. Subsequently, an alternative semi-analytical approach enabling fast parameter identification for complex yield functions is presented.

Approach
The proposed approach can be used for the determination of yield function parameters for experimental analysis of arbitrary sheet specimen containing non-homogenous distribution in the deformation zone. As no further finite element analysis is performed and a simplified computational procedure is used, a fast identification is possible. Prerequisites for the identification are the strain history distribution of one or more intersections through the plastic region, the initial sheet thickness, a flow curve approximation in uniaxial direction and a first estimation for the parameter set of the yield function. The approach can be subdivided into three parts, namely a simplified stress analysis, an internal load calculation and an inverse analysis. As mentioned in section 2, a direct analytical stress calculation using global forces and the cross-section area is only possible if a uniform stress distribution can be assumed. In contrast, the calculation of resulting global forces by stress integration over the section area is straightforward for every kind of stress distribution.
Using this relation, we can analyse arbitrary sections with heterogeneous stress distribution on a specimen (see Fig. 1). For correct assumption of the stress distribution, the static equilibrium of forces has to be satisfied. Using stress analysis on surface measurements [17,20] analogue to VFM, this relationship can be used in inverse analysis. Subsequently, the algorithm of the new approach is briefly described. Before the identification, one or more cross-sections have to be defined as curves over the plastic region. Thereby, the time-and space-resolved strain field is reduced to time-resolved strain distributions over the predefined intersections. The strain distribution history subsequently is used as input data for the inverse analysis to identify yield function parameters.
Based on the strain history, the corresponding stress components are calculated. In contrast to VFM, that uses the elasto-plastic stress computation of Sutton et al. [17], a simplified direct stress analysis on incremental plastic strain is applied. Therefore, the idea of direct stress analysis on plastic strain measurements [20] is modified and formulated for general anisotropic plane stress yield functions. Within the framework, a stress response is computed for each point on the intersection curve using the incremental strain stress relationship and the given strain history. Furthermore, it is important to mention, that the stress analysis is performed in the coordinate system of the material orientation. The computational procedure presented in Fig. 2 is given below.
Assuming the total strain to be plastic, first the incremental plastic strain d dߝ p is calculated based on the optically measured strain. Subsequently, the location on the yield surface is identified by finding the point on the yield surface with a normal vector n collinear to the incremental plastic strain d dߝ p according to the associated flow rule. In a next step, the subsequent yield locus is computed for isotropic work hardening. Therefore, the yield surface is incrementally evolved until plastic work of the provided uniaxial flow curve equals the plastic work of the current stress point. Plastic work of this point is calculated using the measured strain and the stress corresponding to the current yield surface. The resulting accumulative strain ߝeq p is saved for consideration in the following time step as new initial yield surface. Furthermore, for the following procedure, the stress components are transformed to the global coordinate system, as the stress results are computed for the local material orientation system. Subsequently, internal loads are calculated in the global coordinate system by numerical integration of the identified stress components over the intersectional area (see Fig. 1). The intersectional area is defined piecewise by an interpolated curve and updated thickness.
The presented computational procedure of stress analysis and internal load calculation is performed for each time step and used within inverse analysis for the parameter identification. An optimisation algorithm [21] is performed with iteratively varying yield surface parameters until an objective function Φ is minimised: The objective function Φ satisfies the equilibrium of forces and is calculated as the weighted total mean square error of computed internal loads F int compared to experimental determined external loads F ext . Depending on the material model and the number of parameters to optimise, further constraining experimental data like the Lankford-coefficients has to be considered in the objective function. Result of the procedure is an unambiguous parameter set fitted to the material response of the examined intersections on the specimen.

Virtual testing
For validation of the proposed parameter identification, the scheme is implemented in a software prototype and tested with reference data. Instead of using real experimental data, virtual experiments generated by FEA are used to develop the computational procedure. The quasi-ideal results of the fictitious experiment can be used to identify and redesign the effects modelled in the FE reference. Exemplary for the computational procedure, a Swift power law [19] for the flow curve ߪܻ (2) and a Hill quadratic yield criterion ‫ܨ‬ (3) for orthotropic materials [4] are used to model the hardening behaviour: Analogue to real experiments, strain data of an intersectional line in the plastic region computed for the virtual experiment is used as input for the identification scheme. Instead of using facets like in optical measurements the nodal FEA results are used for the following calculations. Compared to real experiments, the crucial advantage of using the FE reference is the possibility to observe internal variables as stress states and the certainty, that only numerical errors but no measurement mistakes will occur. This leads to easier testing and development of the procedure, as deviations between reference and new calculations can be measured.

Shear test
Subsequently, a shear test [22] is used for the validation of the new approach. Therefore a virtual specimen is generated by the commercial finite element code Marc 2013.1.0. using an Updated Lagrange procedure and implicit time integration. Fig. 3 shows the testing setup including geometry and mesh used for the validation. The specimen has a sheet thickness of 1.0 mm and a shearing zone with a length of 4.72 mm defined by the slots in the specimen. A bilinear thin shell formulation (type 139) is used. The mesh is refined in the shearing zone with an initial element size of 0.25 mm x 0.25 mm. The material orientation (R.D.) is collinear to the X-Axis. The material behaviour of the virtual specimen is modelled with the same Swift power law and an equivalent Hill yield criterion. The assumed material parameters for the virtual specimen are given in table 1  and table 2 for a mild steel (EN 1.0338 / DC04).
After FE simulation, the results are used to extract strain data over an intersectional curve in the shearing zone of the virtual specimen. Fig. 4 shows the deformed virtual specimen and the intersection through the shearing zone. Subsequently, the history of the corresponding nonhomogenous strain distributions over the shearing zone (Fig. 5) is used as input for a validation of the computational procedure.

Evaluation of the new approach
Using this intersectional data of the virtual specimen as reference, the suitability of the new approach is evaluated in the following section. The evaluation focuses on stress analysis and internal load and briefly describes the inverse analysis for parameter identification.

Stress analysis
The simplified stress analysis proposed within the new approach is the key method which enables the computational procedure's performance as no further FEA is needed. As the computed stress is directly related to the internal force used as optimisation criterion, an accurate and reliable stress calculation is of major interest. Based on the strain distribution history of the virtual experiment, the stress components ‫,ݔݔߪ‬ ‫,ݕݕߪ‬ ‫ݕݔ߬‬ are determined with the new approach and compared with the FE reference (Fig. 6). In accordance to the FE reference, stress estimations for actual A 4 parameters present good results for the inner part of the intersection. In the outer part, close to the edge, the deviations are up to 70 MPa higher compared to the FE reference (Fig. 6b). These deviations can be traced to simplifications of the nodal stress interpolation in the FE reference. Nodal results are usually interpolated from the integration points in direct neighbourhood. Therefore, the interpolated strain doesn't have to match with the interpolated stress regarding the nodal results. As the results for strain and rotation differ stronger among neighbouring integration points in the edge region, the calculation on nodal results in edge zones becomes less suitable than in the quasihomogenous centre region.
To investigate the impact of this simplification, several elements were isolated (blue elements in Fig. 7). Then, the stress components on the integration points were estimated and compared with the FE reference (see Fig. 8). The accuracy of the stress analysis is significantly increased for a calculation on integration points. Mentionable deviations only occur on the first and last intersection point and are below 10 MPa compared to the FE reference. It is worth to mention, that the plot cannot be directly compared to the distribution on nodal results of Fig. 6, as it represents the results for an intersection translated to the integration points.
As the main deviation only appears in the relatively small edge regions and is homogenised by integration over the intersectional area later, the calculation on nodal results is sufficient for our further investigations, if a small total deviation for the internal load calculation is accepted. Furthermore, a comparatively poor stress calculation is observed for the first time steps, which can be traced to the simplification of only considering the plastic response. As elastic strains are evanescent small, this simplification is sufficient as long as the results are excluded from parameter identification for the first time steps with predominant elastic response.

Sensitivity of stress analysis
Subsequently, the variation of the Hill function parameters is examined in order to investigate their sensitivity on stress computation. Therefore, reference parameters of table 2 are varied by simply changing the A 4 parameter to a value of 1.15385. In the following sections, this new parameter set is called the initial estimation. By changing the A 4 parameter, the normal on the yield surface deviates drastically from its original direction and work hardening is estimated differently. In Fig. 9 the dashed curves shows the stress components calculated on a wrongly estimated initial parameter A 4 . The wrong estimation leads to biased results for shear and normal stress. The strongly diverging value of the shear stress ‫ݕݔ߬‬ is induced by changing the yield surface with an enlarged A 4 parameter. In contrast a smaller A 4 value leads to shear stress values higher than the actual reference.

Internal load calculation
After stress analysis, internal loads can be calculated in the global coordinate system using the nodal coordinates. As stress analysis is performed with respect to the material orientation, every stress state has to be transformed from the system of material orientation to the global system before integration. In addition to the initial orientation, the rotation of the material during the deformation has to be considered. Reasons for this rotation are mainly related to rigid body rotation due to tilting of the specimen and to shear rotation induced by the rotation of local orthogonal element system of the used 3D shell during shear deformation. Fig. 7 shows the material orientation (red arrows) in comparison to the initially collinear global X-axis. Considering these rotations, an internal load in the global system was successfully computed for each time step of the virtual experiment. Fig. 10 shows the results for forces in X and Y direction based on the stress results calculated for the actual parameters in section 5.1. The comparison of internal and external loads shows only small deviations.
For the estimated load in Y-direction a deviation of 2% can be observed (Fig. 10b). These deviations are mainly caused by inaccuracies of computed stress components. The second relevant aspect is the interpolation of the nodal results for rotation, analogue to section 5.1. But, as deviations mainly occur in small edge areas, their impact is not high for the internal load calculation, as they are homogenized after integration over the intersection area.

Inverse analysis
In order to investigate the suitability of the new approach, the implementation of the inverse analysis is tested using the virtual shear test and an initial estimation for the parameter A 4 . The remaining parameters A 1 , A 2 and A 3 are  constrained to the values given in table 2.
Starting with a falsely estimated A 4 parameter as initial estimation, we can observe a wrong result for the internal load induced by incorrect stress estimation. For our assumed initial estimation of section 5.3, this leads to a computed force in X direction which is about 150 N smaller than the external load. In contrast, smaller A 4 values lead to force values higher than the external load.
During inverse analysis, the parameters are iteratively adapted until the internal load converges the external load. Fig. 11 shows the development of convergence for A 4 during the parameter identification. After parameter optimisation, stress components ( Fig.  9 dotted line) and internal loads (Fig.10) approximately result in a same response as the FE-reference. The dashed line in Fig. 11 shows the development of A 4 for an initial value of 1.15385. The final value of 0.86064 is reached at iteration step 30. An equivalent result for inverse analysis is computed for an initial value of 0.69231 (dotted line). The currently not runtime optimized implementation needs less than 5 minutes for this identification on a standard desktop computer. Compared to the FEMU approach, the computation time is reduced by a factor of more than 100.
As presented in Fig. 11, a small deviation of the final compared to the actual A 4 value is observed. As the optimisation considers forces in X and Y direction, biased force calculations for the Y direction like presented in Fig. 10b lead to this divergent results. By introducing weight factors and only considering the load in X direction, an A 4 value of 0.88668 with a smaller force deviation in X direction is calculated by the algorithm. This result shows, that accurate stress estimation, transformation and integration is essential for the parameter identification, as the computational procedure is very sensitive to small deviations of loads in the Y direction.

Summary and outlook
In this study, a new approach for efficient and performant parameter identification for complex yield functions is proposed. Compared to other approaches as finite element method updated, the numerical and experimental efforts for the evaluation of complex specimen with nonhomogenous strain distribution can be reduced considerably.
Within the approach, a simplified stress analysis using the plastic strain stress relationship was used. To reduce computation time, the analysis is performed on the strain distribution of a user defined intersection in the plastic region. By simplified stress integration, a corresponding internal load is computed and used with measured testing forces in an inverse analysis to identify the yield function parameters for balanced forces. For testing, the computation procedure was developed and validated for a Hill yield criterion with fictitious virtual experiments generated by FEA. Using a parameter set describing the actual yield surface, computed stress components and corresponding internal loads only show little deviation compared to the reference data. An inverse analysis with one single varied parameter showed that the procedure is suitable for parameter identification.
With our approach, we introduce a performant tool for parameter identification for complex yield functions. Especially for yield functions needing a higher amount of sampling data to describe distinctive anisotropy, great potential is given. Further developments will focus on the implementation of complex yield functions like Yld2000-2D and on a transfer to experimental analysis of real specimen. Future studies should also consider phenomena like combined isotropic-kinematic hardening. Additionally, the scheme for refined stress computation based inside the approach could be used for other purposes using in-process analysis.