Evaluation of Prediction Accuracy for Anisotropic Yield Functions Using Cruciform Hole Expansion Test

: To evaluate the prediction accuracy of the anisotropic yield function, we propose an original cruciform hole expansion test. Displacements on two axes were applied to the cruciform specimens with a hole in the center. The thickness strain in the region near the hole was compared to the simulation results. Because this forming test is free of friction and bending, it is an appropriate method to assess the anisotropic yield function without the influences of friction or the Bauschinger effect, or the need to consider the stress-strain curve in high-strain region. Hill1948, YLD2000-2D, and spline yield function which is an improved version of the Vegter model were selected, and 6000 series aluminum alloy sheets (A6116-T4) were used in this study. The parameter identification method of the spline yield function also proposed in this paper using the pseudo plane strain tensile test and optimization software. As a result, the spline yield function has better predictive accuracy than the conventional anisotropic yield functions Hill1948 and YLD2000-2D.


Introduction
High-tensile-strength steel sheets and aluminum alloy sheets are used for weight reduction of car bodies. Because of the difficulty in forming these materials, preliminary evaluations by numerical simulation are necessary. Hence, the requirements for the prediction accuracy of simulations of sheet metal forming are high, and further improvement is essential. An anisotropic yield function is indispensable for highly accurate simulation. Within this paper, we evaluate the prediction accuracy of anisotropic yield functions Hill1948 [1], YLD2000-2D [2], and the spline yield function [3].
To evaluate the prediction accuracy of the anisotropic yield function, a hole-expansion test has been carried out in many studies. For examples in [4][5][6][7], the thickness strain distributions along the hole edge and its surroundings were compared with the simulation results. This approach is convenient and effective for assessing the validity of anisotropic modeling, because the anisotropy of the material affects the thickness strain distribution near the edge of the hole. However, the results rarely agree well with the experimental results. It is difficult to validate anisotropic modeling by only considering the change in thickness distribution. In general, this hole-expansion test is relatively simple, but includes contacts between the material and the dies, as well as bending and back-bending deformation. Thus, it is necessary to consider not only the anisotropy of the material, but also friction, the stress-strain relation in the high-strain region, and the Bauschinger effect. However, it is not possible to separate these effects within the calculations. The Coulomb friction model is a simple and popular friction model; however, friction during the forming process is affected by surface pressure, sliding speed, and distance. There are friction models that consider these additional variables [8][9][10]; however, none take anisotropy into account. In [11], the influence of anisotropic hardening was investigated using a hole-expansion test. In addition to the anisotropic hardening, Suzuki et al. [12] considered the Bauschinger effect using the Yoshida-Uemori model [13] and the stress-strain curve in the high-strain region, which was obtained from a plane bending test using a biaxial tensile machine. However, because no advanced friction model was considered only qualitative comparisons could be conducted.
In this study, we calibrate and validate the spline yield function, and compare the results with those of other widely used yield functions for the sheet material A6116 aluminum alloy. To validate the yield function, biaxial tensile tests were conducted using cruciform specimens with a hole in the center. Because this forming test does not include friction and bending, it is possible to a quantitatively compare the experimental and simulation results.

Anisotropic Yield Functions
Several important anisotropic yield functions have been proposed. The input variables of these functions differ; hence, the types of material tests required also differ among functions. To minimize simulation costs, simple functions with relatively high precision are most desirable. However, for functions with high complexity, it is often desirable to increase the predictive accuracy, despite additional costs and material tests. The features of the important anisotropic yield functions are described below.
1. Hill 1948 [1]: this is the most popular anisotropic yield function especially in the car manufacturing industry because of its simplicity. For calibration of this yield function, the results of uniaxial tensile tests in three directions are required. Usually, parameters are identified from the flow stress in the rolling direction (RD) and Lankford coefficients (r-values) in three directions. 2. YLD2000-2D [2]: this is a widely used yield function and is implemented in the several commercial finite element (FE) simulation software tools. In the case of face-centered cubic lattice metal, such as aluminum alloy, it is recommended that the order is set to 8 [14]. To calibrate this yield function, the flow stresses and r-values at 0°, 45°, and 90° with respect to the RD, the equi-biaxial stress value, and the strain increment direction at the equi-biaxial stress state are used. 3. Vegter [15]: this yield function is defined in the principal stress space. In addition to the uniaxial tensile tests in five or more directions and a material test of the equi-biaxial stress state, planestrain tensile tests are required in five or more directions. The Vegter model can consider the results of many material tests, but it requires substantial effort to determine all input variables, and it is difficult to conduct pure plane-strain tensile test.
Considering the feasibility of material testing and of industrial application, Hill 1948, YLD2000-2D, and the spline yield function, which is an improved version of the Vegter model and described in the next section, are compared in this paper.

Spline Yield Function
The spline yield function is based on the Vegter model, and is defined in the principal stress space and enhanced by making its structure flexible. For interpolation function, spline functions are used instead of Fourier functions of the Vegter model. In previous studies [16,17], higher ordered or multi segmented yield loci than the Vegter model were used. The structure of a spline yield function is mathematically identical to a spline surface used in a computer aided design (CAD) system. In contrast with the Vegter model, the spline yield function does not require a pure plane strain stress state. This is one of its important advantages over the Vegter model. In the next section, we describe the calibration method of the spline yield function's parameters using a pseudo plane strain tensile test and inverse analysis.
We consider the yield function with rotational symmetry with respect to the origin. We consider five points on the yield function, including the uniaxial stresses and the equi-biaxial stress in the first quadrant, as shown in Figure 1. This is a spline curve composed of four quadratic Bézier curve segments.
is determined as an intersection of two tangent lines of and and do not need to be the exact plane-strain stress state. ( ) = ( , ) which is stress point on the normalized yield locus in the principal stress space of 1 and 2 is described as follows: where, u is the parameter of the Bézier curve (0 ≦ ≦ 1). The direction of the principal stress, 1 and 2 are orthogonal to each other, and the smaller angle from RD is set to θ, and the other one is set to θ + 90°. To describe all the stress states, the angle θ between the direction of principal stress 1 and the RD must be continuously defined in the range of 0° to 90° in consideration of the symmetry condition due to the orthotropic anisotropy. In this study, five directions of 1 are defined discretely at 22.5° intervals, and ( , ) is defined continuously by interpolation using a spline basis function as follow: where ( ) means the "k"th ( ), if k = 1, the angle of σ1 respect to the rolling direction is 22.5°. ( ) is a one-dimensional cubic spline basis function and is given by the following boundary conditions ( = 0, 1, 2, 3,4).

4
• 90 ⋄ = (3) Figure 2. The convexity condition of the yield function and a relationship between the principal stress and the stress components in the global coordinate system is described in the Appendix A.
, , and and their normal vectors determined from the uniaxial tensile tests and bulge test, and the parameters pk and qk were introduced to identify , and their normal vectors in Figure 1. In Figure 3, is the intersection point of the tangent lines at and , is the intersection point of the lines from to and to . The parameters p0, p4 and q0, q4 are defined as follows: where, is the normal vector at . In Figure 3, is in the RD and is 90° to the RD. The parameters p1, p3, q1, and q3, as well as p0, p4, q0 and q4 are determined for at 22.5° to the RD and at 112.5° to the RD which is the same as at 67.5° to the RD because of symmetry. The parameters p2 and q2 are also determined for at 45° to the RD in the same way. The spline yield function was implemented using a user subroutine of LS-DYNA, which is a commercial code with dynamic explicit method. These parameters are identified by inverse analysis described in the next section.

Calibration of Spline Yield Function
To identify the material parameters for three yield functions, three different types of material experiments were carried out. The thickness of the A6116 specimen was 0.93 mm and T4 solution heat treatment was conducted. The specimens for the uniaxial tensile tests and the pseudo planestrain tensile tests were processed by wire cutting.

Uniaxial Tensile Tests
Uniaxial tensile tests are conducted for 5 specimens with tensile directions of 0°, 22.5°, 45°, 67.5° and 90° from the RD. Three valid repetitions were carried out for each experiment. The material properties obtained from these tests are summarized in Table 1. The normalized stress was determined when the plastic work in each direction was equal to the reference plastic work in the 0° tensile direction. The effective plastic strain was set to 0.05. The stress-strain curve in the 0° direction, which was extended from 0.2 by the Voce approximation, is shown in Figure 4. The anisotropy of the sheet metal is characterized by r-value, and the distribution of the anisotropic behavior of the material can be seen in Figure 5. Within this study, the r-value was calculated as width plastic strain divided by thickness plastic strain, the -value was defined as the ratio of the width plastic strain increment to the plastic thickness strain increment. An accurate calculation was undertaken for the tangential direction of the equi-plastic work surface in the uniaxial stress state by using the -values.

Bulge Tests
To measure the equi-biaxial stress and its strain increment ratio corresponding to the plastic strain of the uniaxial tensile test, hydraulic bulge tests were conducted. The diameter of the die was 123 mm. Three repetitions were carried out for each of experiment. The digital image correlation (DIC) system GOM ARAMIS was used to measure the curvature radius, , and the strain, or , at the apex of the bulge specimen. The equi-biaxial stress state at the apex of the specimen was assumed to be = = . For the determination of , the following equation was used.
where p is the oil pressure during the bulge test, and t is the sheet thickness at the apex of the specimen, which is calculated assuming constant volume during plastic deformation, and using the strain, and . In Table 2, the normalized stress and strain increment ratio are shown for a reference effective plastic strain of 0.05.

Pseudo Plane Strain Tensile Tests
Pseudo plane-strain tensile tests were carried out to determine the parameters for the spline yield function, which are described in Section 4.4. The specimen geometry, which is shown in Figure  6, with a width of 30 mm and length of 60 mm allows the same clamping system to be used for the uniaxial tensile tests. This geometry creates an approximately plan-strain stress state in the center of specimen; however, the strain is not homogenous over the width of the specimen. Similar to the uniaxial tensile test, three valid repetitions in the five tensile directions were carried out, and the load and strain were measured. ARAMIS was used to measure the strain at the center point, , of the specimen and the position , which is 5 mm in the tensile direction from the center (Figure 7). Figure  7 shows the maximum principal strain distribution of the specimen at 0° to RD. In order to compare with the simulation results later, the average strain was introduced as an axial representative strain. The average strain, , was defined as = { ( ) + 2 ( )}/3, where the tensile direction is x and the specimen width direction is y. Post processing parameters of DIC are used 0.2 mm as distance of points, 10 mm as facet size. Figure 8 shows the relationship between the nominal stress and the average strain, and Figure 9 shows the relationship between the strain in the width direction at and the average strain. These results were applied for inverse analysis using the FE analysis software LS-DYNA.   Inverse analysis using FE simulation was applied to determine pk and qk from the results of the pseudo plane strain tensile test. The optimization software LS-OPT was used for this inverse analysis, and a mean square error method was adopted, which had two variables and two objective functions. The input variables were pk and qk and a relationship between a tensile force and the average strain (defined in Section 3) was used for one of the objective functions. The other objective function was the relationship between the strain at P2 in the width direction and the average strain. A FE model for the pseudo plane strain tensile test was determined, as shown in Figure 10. The mesh size was investigated as far as necessary strain, which maximum average strain is 5.0%. The 1 mm mesh size and 0.5 mm mesh size were compared, and the difference of the load force was less than 3% when the average strain was 5.0%. The 1mm mesh was adopted in consideration of simulation cost. The solutions of this optimization, which are parameters of spline yield function, are shown in Table 3. Furthermore, comparisons of the experimental results in 0° direction and the results of the FE simulation after an optimization of the objective functions are shown in Figure 11.  Table 3. Parameters pk and qk for spline yield function.  The following conditions are fully satisfied in the spline yield function: The flow stresses and rvalues in 5 directions; the equi-biaxial stress and its strain increment direction; and the stresses near the plane strain in the five directions and their normal vectors. By contrast, the Vegter model could not be used, because this requires the pure plane-strain tensile test, so it is difficult to determine in the same material test.

Comparison of the Equi-Plastic Work Surface Shape
The normalized equi-plastic work surfaces of the three yield functions, identified by aforementioned material tests, are shown in Figures 12-14. Figure 12 also shows the results obtain using Hill 1948 and YLD2000-2D for comparison. A description of the parameter identification is given in the next section. The horizontal axis of Figure 12 corresponds to the RD and the vertical axis is 90° from the RD. The horizontal axis of Figure 13 is 22.5° with respect to the RD and the vertical axis is 112.5° (equal to 67.5° because of symmetry). The horizontal axis of Figure 14 is 45°, the vertical axis is 135° (equal to 45° because of symmetry). In Figures 13 and 14, only the spline yield function is shown. because other yield functions are difficult to show in the principal stress space except when the horizontal axis is the RD. The above-mentioned equi-plastic work surfaces were normalized with respect to the plastic work of the uniaxial tensile test for the 0° tensile direction. The reference plastic strain was 0.05. Depending on angle of the axis in the principal stress space, the shape of spline yield function changes substantially for this material.

Experiment
To evaluate the prediction accuracy of the spline yield function, we carried out in-plane cruciform tests, which are free of friction and bending using the biaxial tensile test machine [18], as shown in Figure 15. The cylinders of both axes can be controlled by load or displacement. The strain was measured by the LIMESS DIC system, using four cameras installed above the specimen. The cruciform specimens were manufactured by milling and the face of the cross section was fine-grained. To allow the best comparison of experimental and simulation data, the tests were carried out by displacement control because the displacement can be directly used as a boundary condition in the simulation model.
In total, six sets of experiments with three values of α and two displacement ratios were conducted (Table 4). Thereby, the angle between the horizontal arm of the specimen and the RD is ( Figure 16). Three valid repetitions were carried out. For a displacement ratio of 2:1, velocities of 4 and 2 mm min −1 were used. For the ratio 1:1, the velocity was set to 2 mm min −1 in both axes. The frequency of DIC system was 5 frames s −1 . The maximum and minimum principal true strain were determined using the software istra4D along the specified circle, Subsequently, the thickness strain was calculated assuming constant volume during plastic deformation. Post processing parameters of DIC are used 0.1 mm as distance of points, 10 mm as facet size.

Numerical Analysis
LS-DYNA was used for this FE analysis. In FE model, a fine mesh size of less than 1 mm was used along the center hole. The displacement constraints were given as boundary condition, and the assumed strain shell elements [19] were used. To determine the anisotropic parameters for Hill1948 yield function, three r-values at 0°, 45°, and 90° direction relative to the RD in Table 1 were used. The order m of YLD2000-2D was set to 8, the anisotropic parameters α1 to α8 [2] were numerically determined using the following experimental results: the flow stresses and the r-values at 0°, 45°, 90°, equi-biaxial stress and strain increment ratio rb in Table 2. Anisotropic parameters α1 to α8 and m are shown in Table 5.

Comparison of Numerical and Experimental Results
The total force of the step, where the maximum effective plastic strain of the simulation results of the spline yield function was 0.15, was taken as a reference load force to match the conditions of the simulation and experiment. The step corresponding to this load force was selected from the experiment and the simulation results of the other models. The effective plastic strain of 0.15 was  Figure 17. For numerical comparison, circles are defined for evaluation as shown in Figure 18. These positions were determined in the initial state before applying the displacement. The horizontal axis of the graph for numerical comparison ( Figure 18) is β, which is in the range from 0° to 360°. Comparisons of numerical and experimental results of thickness distributions are shown in Figures 19-27. The sum of the squared errors between experimental averages and simulation results for each case are also shown to quantify the predictive accuracy; The smaller sum of squared errors means better prediction.         Table 6, which support the above conclusion. The YLD2000-2D function generally has high accuracy, but this is Furthermore, the positions of the predicted maximum strain at approximately 90° and 270° in experiment D deviate from those in the experiment. There were no peaks at 90° and 270° of circle 1 in experiment F; however, the calculation result of YLD2000-2D had different upper and lower peaks from the experiments. The prediction accuracy of the Hill 1948 yield function is poor in the experiment A, B, and D. However, the results at 45° are in good agreement, which we attribute to the equal load on both axes, which is caused by symmetry of orthotropic anisotropy for the rolled material. Table 6. Comparison of sum of squared errors among the yield functions.

Spline Yield Function
Average of sum of squared errors between experimental averages and simulation results for 3 circles of 6 conditions in Table 4 0.03032 0.00114 0.00063

Concluding Remarks
The prediction accuracy of anisotropic yield functions was evaluated using an in plane biaxial forming tensile test. For the investigation of our proposed yield function proposed and conventional models, we used the aluminum alloy sheets of A6116-T4. The key messages of this work are summarized as follows: 1 A method for determining the parameters of the spline yield function has been developed, which includes a pseudo plane strain stress test by using inverse analysis. 2 For the evaluation of the prediction accuracy of the yield function, an in-plane cruciform tensile test was conducted. This kind of test is appropriate, because it is frictionless and free of bending. 3 The experimental and simulation results reveal high prediction accuracy of the spline yield function. Hence, the validity of the method for parameter determination is verified. 4 The spline yield function is the most accurate model for all considered conditions of the evaluation.
The good agreement of the spline yield function can be ascribed to the high number of input variables (i.e., the uniaxial stresses in five directions, the equi-biaxial stress, the stresses near the plane-strain stress state in five directions, and strain increment directions at all described stress points). There are some differences in the pseudo plane-strain tensile strength depending on the tensile directions. Only the spline yield function can consider these differences, which leads to the best accurate yield function for this material. The comparison of prediction accuracy does not simply determine the superiority or inferiority of the yield function. It should be recognized that each has different numbers of anisotropic parameters. The YLD2000-2D also shows high accuracy considering that the plane-strain tensile test is not needed. Thus, one fewer test yet high accuracy is a benefit of the YLD2000-2D. In further study, other anisotropic yield functions will be evaluated in this cruciform tensile test. thank Adam Brotchie, PhD, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.