Identification of Anisotropic Yield Functions Using FEMU and an Information-Rich Tensile Specimen

To fully exploit the predictive accuracy of advanced anisotropic yield functions, a large number of classical mechanical tests is required for calibration purposes. The Finite Element Model Updating (FEMU) technique enables to simultaneously extract multiple anisotropic parameters when fed with heterogeneous strain fields obtained from a single information-rich experiment. This inverse approach has the potential to mitigate the experimental calibration effort by resorting to a single, yet more complex experiment augmented with Digital Image Correlation. In this paper, we inversely identify the sought anisotropic parameters of two selected yield functions for a low carbon steel sheet based on the previously designed information-rich tensile specimen. The experimentally acquired strain field data is used to inversely identify the Hill48 yield criterion and the Yld2000-2d yield function, respectively. The results are compared with conventional calibration methods for both anisotropic yield functions. The inverse identification is then thoroughly studied using virtual experiments enabling to disentangle the effect of the material model error and the strain reconstruction error (DIC), respectively. It is shown that the material model error dominates the inverse identification of the Hill48 yield criterion. The reduced material model error for the Yld2000-2d yield function enables obtain inversely identified anisotropic parameters that are closer to the reference parameters. The paper clearly shows the importance of the predictive accuracy of the selected anisotropic yield function when applying inverse identification. Keywords: Anisotropic yield criteria; Material parameters identification; Heterogeneous mechanical tests; Inverse identification; DIC.


Introduction
Sheet metals typically exhibit an orthotropic plastic material behaviour related to the rolling process, causing different plastic behaviour along different directions [1]. The application of Finite Element Analysis (FEA) modelling to accurately simulate sheet metals forming behaviour heavily relies on the predictive accuracy of the adopted anisotropic yield function. Numerous anisotropic plastic constitutive models exist to precisely reproduce the anisotropic plastic material behaviour [2]. The more advanced models involve multiple anisotropic parameters requiring an extensive experimental campaign to accurately calibrate the model using conventional material tests. Strong proof of concept [1,[3][4][5][6][7][8][9] is available showing that a heterogeneous mechanical test along with an inverse identification strategy can reduce the experimental calibration effort. Assuming that the selected anisotropic yield function enables to describe the actual plastic material response, a well-designed heterogeneous mechanical test "activates" all anisotropic parameters leading to a robust and accurate inverse identification process. The strategy for such a strain heterogeneity design methodology has been proposed by Souto et al. [10]. It needs to be noted that the approach consists of accepting the shortcomings of the selected anisotropic yield function. In that case, a heterogeneous mechanical test should be designed in accordance with the forming process to be simulated afterwards [3]. This paper aims at a better understanding of the role of the selected anisotropic yield function in the inverse identification using FEMU [3]. The starting point of this paper is a previously designed tensile specimen that maximizes the strain heterogeneity in the Area Of Interest (AOI) [10]. The material under investigation is a low carbon steel sheet that was fully characterized by Coppieters et al. [11] using conventional material tests yielding the reference values for the anisotropic parameters. The anisotropic parameters of both selected anisotropic yield criteria are inversely identified based on the heterogeneous experiment and compared with the reference. The assessment is based on the residual cost function values at reference and identified anisotropic parameters. The material model error is disentangled from the strain reconstruction error by resorting to virtual experimentation (FEDEFmodule [12] in MatchID [13]) that enables to accurately mimic the complete DIC measurement chain. The paper is organized as follows. In Section 2, the selected anisotropic yield function and the reference parameters acquired from conventional material testing are introduced. In addition, the heterogeneous mechanical test is briefly discussed along with the FE model used in the FEMU process. Section 3 embarks on the experimental set up, the principles of FEMU and the metrics to quantify the identification quality. Section 4 presents the identification results and thoroughly investigates the role of the selected anisotropic yield function using virtual experimentation.

Material model
The test material used in this study is 1.2 mm thick industrial cold-rolled steel sheet. The plastic anisotropy of this material was characterized using an extensive experimental campaign reported by Coppieters et al. [11]. The isotropic elastic behavior was modelled using a Young's modulus E of 219 GPa and Poisson coefficient v of 0.3. The material is assumed to be plastically orthotropic, and two selected plane stress anisotropic yield functions are studied, namely the Hill48 yield criterion and the Yld2000-2d yield function. For the Hill48 criterion, the equivalent stress can be expressed as: with F, G, H and N being the anisotropic parameters. For Yld2000-2d [14], the equivalent stress can be calculated as: where the exponent "M" is a material parameter and typically based on the micro-structure and ' i X , '' i X are the principal values of two stress tensors ' X , '' X . The model involves eight independent anisotropic parameters, i α (with 1,...,8 i = ). The selected anisotropic yield functions were integrated within the FEMU framework using the Unified Material Model Driver for Plasticity (UMMDp) [15]. Standard mechanical tests (uniaxial, biaxial tensile tests and multi-axial tube expansion tests) were used to calibrate the selected anisotropic yield functions and more details can be found in Coppieters et al. [11]. Table 1 shows the reference parameters obtained through the conventional test campaign. It must be noted that the evolution of the equivalent stress eq σ as a function of the equivalent plastic strain pl eq ε is considered to be known and determined from a uniaxial tensile test in the Rolling Direction (RD). The latter data is used to calibrate the strain hardening behavior using Swift's hardening law: The identified model parameters, i.e. deformation resistance K , the initial strain 0 ε and the hardening exponent n can be found in Table 1.

FEA model
The heterogeneous specimen shape was designed using shape optimization driven by the maximization of the strain heterogeneity. More details can be found in [16]. It must be noted that slight geometrical deviations were observed between numerical specimen design and the manufactured specimen. Therefore, the final numerical geometry is a digitization of the manufactured specimen. Fig. 1(a) illustrates the geometry and Fig. 1(b) show local seed size, element formulation (S4R), boundary conditions and material orientation, respectively. Two material orientations are studied in this paper. In the first specimen, denoted as Test C TD (see Fig.1 (b)), the tensile direction was aligned with the Transverse Direction (TD). In the second specimen, denoted as Test C R45, the tensile direction is aligned with 45D (i.e. 45 o with respect to the RD direction). As shown in Fig. 1, an AOI of 22 mm by 30 mm was used.

Experimental Set-up
The experiments (Test C TD and Test C R45) were conducted on a regular tensile machine with a load capacity of 10 kN using a constant crosshead speed of 5 mm/min ensuring quasi-static conditions. Stereo DIC was used to measure the displacement fields within the AOI using MatchID [13]. The DIC settings for the two tests are obtained by a DIC performance analysis as described in [17]. The detailed DIC settings for the two specimens are listed in Table 2. The stereo DIC system was synchronized with the load cell so that, for each image, the system records the corresponding tensile force. The experimental set up is illustrated in Fig. 2.   Fig. 3 shows the experimentally acquired tensile force as a function of the testing time. It can be inferred that both tests yield virtually the same global force response. Note that although force-time curves for the two-orientation tests are similar, they correspond to different deformation fields.

Fig. 3
Force curves with respect to time step for Test C TD and Test C R45.

Key Engineering Materials Vol. 926
To investigate whether the identified parameters depend on the applied load level (e.g. second load step) and how much data (e.g. the number of load steps within the maximum of four load steps) is used in the identification, the following four different load level are considered (see Table 3) within the grey shaded box of Fig. 3. Each of the selected load steps corresponds to a similar maximum equivalent plastic strain within the AOI. As shown in Fig.3, the maximum equivalent plastic strain-pl eq ε for the first load step and for the fourth load step is around 0.09 and 0.20, respectively. The most straightforward FEMU approach is to apply the measured tensile forces at the edges of the FE model, i.e. a force-driven simulation approach, as shown in Fig. 1(b). The experimental data associated with the four load steps are employed to inversely identify the plane stress Hill48 yield criterion as well as Yld2000-2d using FEMU.

FEMU and identification quality assessment
The basic principle of FEMU is to minimize the discrepancy between the experimental and numerical strain response by tuning the unknown anisotropic parameters of a chosen yield function used in the FE model. The Levenberg-Marquardt [18] algorithm is used to minimize ( ) C θ . Given the forcedriven simulation approach, only the strain discrepancy is considered by a least squares cost function that reads as: where θ is the vector of unknown anisotropic parameters and n the number of data points in the DIC measurement at load step j . The m RMS (root mean square) value can be obtained by Eq. (5), which was used to normalize each strain component thereby eliminating the magnitude difference in Eq. (4). The subscripts "exp" and "num" indicate experimental and numerical response, respectively. The parameters identification quality can be expressed by the Average Relative Error (A.R.E) for all identified parameters as: FEMU was fed with experimental data to identify the parameters of Hill48 model. The identified parameters of the Hill48 yield criterion obtained for Test C TD and Test C R45 are reported in Table  4 and Table 5, respectively. The table shows the results for the identification using different individual load step experimental data (e.g., Second load step means only the second load step experimental data was fed in FEMU) as well as all load steps experimental data (four load steps means four load steps experimental data were fed in FEMU). Coppieters et al. [11] showed that the investigated material exhibits differential work hardening in the range 0 < pl eq ε < 0.03. Isotropic hardening (i.e. the shape of the yield function is constant) is valid beyond this strain range. Given that the considered load steps generate strain levels that can be associated with isotropic hardening, it can be also inferred that the effect of considering single or multiple (4) load steps is marginal. From Table 4 and Table 5, it can be inferred that the identification of the parameter H is relatively stable and very accurately identified (<5%). Nevertheless, the overall identified parameter quality (A.R.E) is unacceptably high. The reason for this lies in the choice of the reference anisotropic parameters. Indeed, it was shown by Coppieters et al. [11] that the plastic behaviour of the test material cannot be accurately captured by the r-based (i.e., based on the Lankford coefficients) Hill48 yield criterion. This can be exemplified by computing the cost function when using the reference anisotropic parameters. For example, at second load step, this yields a cost function value of 30.65 which is considerably larger than the cost function after inverse identification, i.e. 9.89. This is basically a consequence of adopting an anisotropic yield function that cannot capture the actual material behaviour. The FEMU approach merely fits then the adopted model to the available experimental data by minimizing the discrepancy between the numerical and experimental strain field(s).

Table 4 Comparison of reference parameters values and identified parameters values for Test C TD
at different load step. 1 Cost function computed at reference parameters of Hill48 with real experimental data at second load step. 2 Cost function computed at reference parameters of Hill48 with virtual experimental data at second load step.
The remaining discrepancy is the sum of several error sources [3,8]: • FEA discretization error: the FEA data needs to be interpolated to the DIC data points. This error mainly depends on the mesh size and here we conducted a mesh convergence study leading to a local element size of 0.3 mm.
• Geometry and the boundary conditions error, which is due to the approximation of the geometry and the boundary conditions in the FEA model. The geometry was derived from the manufactured sample, hence deemed sufficiently accurate. Given that the grips are far from the AOI, the force-driven approach is deemed to be accurate.

•
Mapping errors between the FEA model and the experiment. Before the FEMU procedure, The coordinate systems was carefully aligned. • Measurement errors (e.g. adopted DIC settings). The optimal DIC settings were obtained from performance analysis via MatchID module. Moreover, the influence of DIC measurement becomes negligible when large plastic deformation is of interest [8].
Given the above, it can be concluded that the largest uncertainty is probably caused by the numerical material model. This was also a finding in the work of Cooreman et al. [3]. Nevertheless, it is until now unclear too which extent the DIC error plays a role of importance. To disentangle the role of a material model error and the DIC error, we resort to virtual (or synthetic) DIC experiments. Such synthetic data is generated based on an FE model that adopts a material constitutive model that is considered the ground truth. The displacements generated by this FE model are then used to numerically deform a speckle pattern. The data is then processed using a DIC code to retrieve the displacement fields and strain fields.  Table 4. According to previous research [3-4, 8, 14], material orientation is of fundamental importance when it comes to identifying the plastic anisotropy of steel sheet. To investigate the material orientation influence for identified parameters, the same procedures were repeated for Test C R45.
The results are summarized in Table 5. It can be seen that the same order of magnitude is found for A.R.E, yet it appears to be more stable and independent of the considered load step. This could be interpreted as if the parameters obtained with the Test C R45 are more repeatable than Test C TD. Additionally, the identified parameters of N is consistent (around 1.6), which can probably ascribed to the higher sensitivity of N in Test C R45. The virtual experiment is also conducted at second load step for Test C R45, and the DIC error is 42.19 as shown in Table 5. It is clear that the DIC error in this experiment is larger than Test C TD. This is probably related to the quality of the speckle pattern and DIC settings. For Test C R45, the cost function value at reference parameters for actual experimental data at second load step is 116.13, indicating that the r -based Hill48 criterion fails to reproduce the experiment. Indeed, the estimated model error is here 73.94 which is significantly larger than for Test C TD (20.55). Table 5 Comparison of reference parameters values and identified parameters values for Test C R45 at different load step. 1 Cost function computed at reference parameters of Hill48 with real experimental data at second load step. 2 Cost function computed at reference parameters of Hill48 with virtual experimental data at second load step.

Achievements and Trends in Material Forming
Although it was previously shown a numerical study of Zhang et al. [16] that Test C R45 exhibits a better identifiability than Test C TD, this could not be confirmed experimentally in this study due to the large material model error associated with the Hill48 yield criterion. Indeed, for a simple model such as the Hill48 model, it is clearly shown that the identified parameters largely depend on the adopted heterogeneous experiment. Similar conclusions were repeatedly reported [3,8,16], e.g.
Coppieters et al. [20] argued that the validity of the inversely identified Hill48 yield function is confined to the experiment used in the FEMU procedure. The test configurations of Test C TD and Test C R45 were designed to maximize the strain heterogeneity. This leads to an information-rich experiment which increases the probability of a good identifiability [16]. An important caveat, however, is that the plastic deformation in the test becomes more complex putting more demands on the accuracy of the selected material model (here the anisotropic yield function). Inverse identification of an "inferior" material model using a heterogeneous experiment, designed by maximizing the strain heterogeneity, is discouraged as it leads to a model calibration that is only valid for the considered heterogeneous experiment. In general, if the material models error is significant, Cooreman et al. [3] recommended to inversely identify parameters via designed mechanical test in accordance with the material forming which will be simulated afterwards. However, to the author's best knowledge, there has been no progress in designing a heterogeneous material test based on the dominating material states during a forming operation. The novelty in this section is that a unique method to disentangle the material model error from the DIC error in a heterogeneous experiment is proposed. Finally, our findings in this section suggest that more advanced constitutive models are required in conjunction with the adopted heterogeneous experiment. Hence, a more complex numerical model, namely, Yld2000-2d, was applied in the next section.

Yld2000-2d identification
Exactly the same FEMU process with identical experimental data was employed to identify the parameters of the Yld2000-2d model for Test C TD and Test C R45. For a fair comparison with the reference parameters of Yld2000-2d model (calibrated at pl eq ε =0.24, as shown in Table 1), the experimental data from the fourth load step was used. Consequently, the anisotropic parameters of Yld2000-2d were identified at fourth load step for the Test C TD and Test C R45, respectively. Similar to Hill48, assuming eq σ as the yield stress in rolling direction [6], the following relation between the model parameters should be met:    4 shows the comparison of the stress state distribution as well as its corresponding density contour plots using reference parameters (e.g., Fig.4 (a) and Fig.4 (c)). The reference yield loci and the inversely identified loci are shown Fig.4 (b) and Fig.4 (d) at different values of 0 xy σ σ . It can be inferred that the inversely identified yield loci fit very well the stress states probed in Test C TD and Test C R45, as shown in Fig. 4. Due to the constraint expressed by Eq. (8), the two inversely identified yield loci coincide at uniaxial tension in the RD (aligned with x direction). Nevertheless, the inversely identified yield loci are inaccurate in regions (see dashed black rectangle in Fig. 4(a) and Fig. 4(c)) that do not contain experimental data. In that sense, Yld2000-2d seems to be merely fitted to the available experimental data. However, there are a number of observations that enable to arrive at a more advanced identification approach. First of all, it can be seen that with the increase of 0 xy σ σ , the yield locus for Test C TD shows a larger discrepancy with the reference yield locus than the yield locus for Test C R45, which is illustrated by the dashed blue rectangles (see Fig. 4(b) and Fig. 4(d)). This indicates that a proper selection of the material orientation within the heterogeneous experiment can enhance the identifiability. As shown by Güner et al. [6], the experimentally measured equibiaxial stress value can be included in the cost function to improve the identification quality. Yet, this solution increases the experimental work load again. The work of Zhang et al. [22] suggests that, based on an identifiability analysis, the identification error associated with 7 α and 8 α is caused by parameter interactions. An identification strategy that enables to mitigate the effect of parameter Key Engineering Materials Vol. 926 interaction is currently investigated and will be published in a forthcoming paper. Finally, it must be noted that the problem can be solved if a complex tensile-load specimen is also used to add biaxial tension, as proposed by Macek et al. [23].

Conclusions
The FEMU technique combined with DIC is employed to inversely identify the Hill48 and Yld2000-2d function using a heterogeneous tensile specimen designed for maximum strain heterogeneity. From the analysis, two main conclusions can be drawn: • A method based on virtual experimentation is proposed to disentangle the DIC error from the material model error. It is shown that for the material under investigation and the adopted heterogeneous experiment, the material model error associated with the Hill48 yield criterion is significantly larger than the error associated with the Yld2000-2d yield function.
• When using an optimized heterogeneous experiment that maximizes the strain heterogeneity, an accurate anisotropic yield function needs to be selected to enable proper identification. In case the material model error is significant, it is recommended to use a heterogeneous experiment that generates similar conditions as in forming process that will be simulated afterwards.