On the Use of Homogeneous Polynomial Yield Functions in Sheet Metal Forming Analysis

Sheet metal forming techniques are a major class of stamping and manufacturing processes of numerous parts such as doors, hoods, and fenders in the automotive and related supplier industries. Due to series of rolling processes employed in the sheet production phase, automotive sheet metals, typically, exhibit a significant variation in the mechanical properties especially in strength and an accurate description of their so-called plastic anisotropy and deformation behaviors are essential in the stamping process and methods engineering studies. One key gradi-ent of any engineering plasticity modeling is to use an anisotropic yield criterion to be employed in an industrial content. In literature, several orthotropic yield functions have been proposed for these objectives and usually contain complex and nonlinear formulations leading to several difficulties in obtaining positive and convex functions. In recent years, homogenous polynomial type yield functions have taken a special attention due to their simple, flexible, and generalizable structure. Furthermore, the calculation of their first and second derivatives are quite straight-forward, and this provides an important advantage in the implementation of these models into a finite element (FE) software. Therefore, this study focuses on the plasticity descriptions of homogeneous second, fourth and sixth order polynomials and the FE implementation of these yield functions. Finally, their performance in FE simulation of sheet metal cup drawing processes are presented in detail.


Introduction
Sheet materials represent significant anisotropic behavior due to their thermomechanical process history. Anisotropy states the variation of the mechanical properties with direction. This material property is determined from tensile test and it is calculated by dividing width plastic strain increments to thickness plastic strain increments. From this definition, it is seen that anisotropy indicates the resistance to the thinning. Therefore, it can be said that increasing anisotropy values improves the deep drawability of the material. Two approaches are applied in the description of the anisotropy. The first approach is the phenomenological approach in which global material behavior is determined according to the average behavior of all grains. The second approach is crystal plasticity which investigates the behavior of one grain to determine the material behavior.
In the phenomenological plasticity approach, the transition from the elastic deformation to plastic deformation is defined with yield functions [1]. A yield function establishes the relationship between principal stresses and yield stress of the material. Plastic flow occurs when the yield function reaches a critical value which is the yield stress of the material. Therefore, yield condition actually indicates a state of equilibrium and it can be defined by the following equation: where σ and σ y denote the equivalent and yield stresses, respectively. Eq. (1) defines a surface in three dimensional stress space and it is called as yield surface. According to Drucker's postulate [2], this surface must be closed, convex, and smooth in order to establish a relationship between plastic strain increments and stresses. In the literature, Tresca and von Mises are well known and have been most commonly used yield criteria. However, these yield criteria are isotropic and they could not give satisfactory results for sheet metal forming processes. Therefore, the usage of anisotropic yield functions is required for representation of sheet metal behavior and several anisotropic yield functions have been proposed by researchers. The first phenomenological anisotropic yield function was proposed by Hill in 1948 [3]. Hill added some coefficients to von Mises criterion to transform isotropic von Mises criterion into an anisotropic form. Hill's quadratic criterion could be used for both plane stress (2D) and general stress (3D) states. The criterion has four coefficients for 2D stress state, and it has six coefficients for 3D stress state. These coefficients could be obtained analytically according to stress or plastic strain ratios. Hill48 quadratic criterion has a simple form and useful coefficient identification procedure. However, this criterion could not simultaneously predict the variations of the stress and strain ratios within the sheet plane. Therefore, it could not successfully define the plastic behavior of highly anisotropic materials such as Al-Mg alloys, Ti alloys, etc. Different type yield criteria have been applied to accurately describe the anisotropic behavior of these materials. The most popular approach used to derive an anisotropic yield criterion is the linear transformation method. In this method, Cauchy stress tensor or the deviatoric stress tensor is transformed linearly, and an anisotropic yield function is obtained by substitution of this transformed tensor in an isotropic yield function [4]. Yld89 is one of the functions developed by this approach. Barlat and Lian [5] applied linear transformation method to Hosford 1972 [6] isotropic yield criterion and developed this anisotropic material model. The criterion has four coefficients and it could be used for only 2D stress state. Then, Barlat et al. [7] extended this yield criterion for 3D stress state and developed a criterion has six coefficients in 1991. However, these yield criteria could not accurately describe the anisotropic behavior of especially Al-Mg alloys. Another yield criterion based on linear transformation approach was developed by Karafillis and Boyce [8] in 1993. Karafillis and Boyce generalized Hosford's yield function and proposed an isotropic yield function. Then researchers applied to linear transformation approach and developed an anisotropic yield criterion. They applied their developed yield criterion for modeling of AA2008-T4 alloy and could successfully define the angular variations of both stress and plastic strain ratios of the material. Barlat et al. have inspired by this method and developed Yld2000 and Yld2004 yield criteria, respectively [9,10]. From these models, Yld2000 could only be used for plane stress condition, whereas the other could be used for both plane stress and general stress states. Yld2004 criterion has 18 coefficients and it could successfully describe in-plane variations of plastic properties of highly anisotropic aluminum alloys. These models are effective in the representation of the anisotropic behavior. However, their parameter identification procedures consist of complex nonlinear formulas and computation of the derivatives is difficult.
Another method which is applied to derive anisotropic yield function is the polynomial approach. Due to inability of quadratic Hill48 criterion, Hill suggested that the usage of general homogeneous polynomials as yield functions in 1950 [11]. In the literature, firstly Gotoh [12,13] applied this method and modeled the anisotropic behavior of commercial Al-killed steel and Cu-(1/4)H sheets with fourthorder polynomial yield function. Gotoh determined explicitly the coefficients of the polynomial function for these materials and successfully predicted the angular variations of the plastic properties. However, Gotoh did not take into account the convexity of the yield surface in the parameter identification. This deficiency was noticed by Soare et al. [14] and they proposed changes to Gotoh's identification procedure. This modification has contributed to the applicability of the polynomial criteria and important results have been obtained.
In the present work, polynomial yield criteria, their modeling capability and applications on the sheet metal forming simulations have been investigated. Article consists of four sections. In Section 2, the theoretical background of the developed polynomial yield functions are briefly explained. Then, applications of polynomial criteria and results are presented. In Section 4, the main conclusions and findings are summarized.

Homogeneous polynomial yield functions
It is seen from the literature that the second, the fourth, and the sixth-order homogeneous polynomials have been used as yield functions. Therefore, the general formulation of these functions are explained in this section.

Second-order polynomial yield function
Conventional quadratic Hill48 yield criterion can be defined as second-order polynomial yield function (P 2 ). The form of the criterion for plane stress state could be written as follows: (2) a 1 , a 2 , a 3 , and a 4 are function parameters and they can be determined based on stress or plastic strain ratios. The equations related to stress and plastic strain ratios are given in below. The coefficients determined with stress or strain based definition are distinguished by subscripts σ and R, respectively.

Fourth-order polynomial yield function
For plane stress state, the fourth-order polynomial yield function (P 4 ) is expressed as following: where a 1 , a 2 , a 3 … .a 9 are the material coefficients. In order to determine these nine coefficients, nine experimental data are required. Direct approach for coefficient determination can lead to oscillations in the predictions of the plastic strain or yield stress ratios. Therefore, Soare et al. [14] proposed a different coefficient identification procedure and derived upper and lower bounds on coefficients to obtain a convex and smooth yield surface. In this section, the coefficient identification procedure developed by Soare is explained: (i) Firstly, the first five coefficients are determined with explicit formulas are given below: where r 0 and r 90 indicate plastic strain ratios (r-values) along rolling and transverse directions, whereas σ 90 denotes yield stress ratio along transverse direction.
(ii) The coefficient a 3 is determined according to the Eq. (7).
where σ b indicates the biaxial yield stress ratio.
(iii) The coefficient a 9 is determined according to Eq. (8) where σ 45 and r 45 indicate the yield stress and plastic strain ratios along the diagonal direction.
(iv) The coefficients a 6 and a 8 are determined with the minimization of the error (distance) function given in Eq. (9).
where w 1 and w 2 are the weight coefficients for stress and plastic strain ratios at the interval angles. In this minimization problem, interval angles could be 15 0 -75 0 , 30 0 -60 0 or 22.5 0 -67.5 0 . After determination of the coefficients a 6 and a 8 , these coefficients are checked for positivity and convexity of the yield surface. In order to obtain convex and smooth yield surface, a 6 and a 8 must satisfy the following inequalities: v) The coefficient a 7 is determined with Eq. (11) Inequalities related to convexity and positivity conditions are given detailed in [14].

Applications of polynomial yield functions
Three validation studies are generally performed in the literature in order to evaluate the prediction capability of orthotropic yield criteria: These are the description of the planar variations of plastic properties, the prediction of the earing profile and number of ears in cup drawing test, and prediction of the thickness strain distributions along the different directions in a drawn part, respectively. Obtained results with polynomial yield functions are presented in below.

Description of the directional properties
Soare et al. [14] investigated the prediction capability of the polynomial yield functions. They described the anisotropic behavior of AA2090-T3 with P 4 and P 6 yield criteria. Figures 1 and 2 show the P 4 and P 6 predictions of the angular variation of plastic properties for AA2090-T3 alloy, respectively.
It is seen from Figures 1 and 2 that both criteria could simultaneously predict the angular variations of stress and plastic strain ratio. In addition to that the predictions of P 6 criterion were more successful than P4 criterion especially at interval angles. Sener et al. [15] investigated the evolution of anisotropic behavior of Al5754 with P 2 and P 4 yield criteria. They determined the coefficients of the yield functions at four different plastic strain levels and predicted the angular variations of yield stress and plastic strain ratios. Then, researchers compared the predicted results from yield criteria with experimental data for each plastic strain level. Figures 3 and 4 show the comparison results for P 2 and P 4 criteria, respectively.
It is seen from Figures 3 and 4 that P 2 criterion could only accurately predict the variation of r-values in the sheet plane, while P 4 criterion could predict both the angular variations of stress and strain ratios. This result is related to the identification procedures of the yield criteria. As it is declared in Section 2 that, P 2 criterion takes as input either stress or strain ratios. However, the coefficients of P 4 criterion   are calibrated with both stress and strain ratios. In addition to description of the planar anisotropy, researchers investigated the variation of the yield locus shape with plastic strain. Figure 5a and b show the variation of yield locus contours with plastic strain for P 2 and P 4 yield criteria, respectively.
It is seen from Figure 5 that the contours of the yield locus are changed with plastic strain and this evolution is more pronounced in P 4 criterion.

Prediction of the earing profile
Cup drawing is a test which is used for validation of an anisotropic yield criterion. If material has a strong anisotropy, the height of the formed cup is not uniform and a series of crests and valleys are observed around the cup perimeter. This waviness in the top edge of a cup is called as earing and four, six or eight ears could be occurred in a drawn cup depend on the degree of the anisotropy [16,17]. Soare et al. [14] investigated the prediction capability of polynomial yield functions on the cup drawing test. They implemented P 4 and P 6 yield criteria into FE code ABAQUS and performed FE analyses of the test. Researchers also studied the effect of element type on the predictions and they carried out simulations with shell and solid elements. After FE analyses, they predicted the number of ears, cup height, and compared the numerical results with the Yld96 criterion and experiment. Yld96 criterion was selected as reference by the researchers due to involving the same number of material coefficients of both criteria. Figures 6 and 7 show the geometry of the drawn cup and the comparison of the predicted cup profiles from P 4 and Yld96 yield criteria with experiment for AA2090-T3 alloy.
It is seen from Figure 7 that P 4 and Yld96 criteria could successfully predict cup heights, however the predictions of P 4 were closer to the experiment in the rolling direction. Both criteria predicted two extra ears along the transverse direction (90°a nd 270°). It was also observed that there are no significant differences between the predictions of P 4 -2D, and P 4 -3D models. Researchers also investigated the capability of P 6 criterion on earing prediction and compared the predictions with Yld2004 and experiment. These comparisons are shown in Figure 8.
From the comparisons, it is observed that P 6 criterion could accurately predict both the number of ears and cup height. Another observation in this study is related to Yld2004 and P 6 predictions. Both criteria gave similar results and this shows that P 6 has higher capability in the modeling of the anisotropy.

Prediction of thickness strains in rectangular cup drawing
Another study related to polynomial yield functions was carried out by Sener et al. [18]. They investigated the anisotropic behavior of AISI 304 stainless steel with P 4 yield criterion. Investigation was conducted on the uniaxial tensile test and a rectangular cup drawing process. Criterion could successfully describe stress anisotropy and r-value variations. Researchers implemented the criterion into explicit FE code Ls-Dyna by using user defined material subroutines and performed FE simulation of rectangular cup drawing process. They investigated the thickness distributions and flange geometry. Figures 9 and 10 show the comparisons of the numerical and experimental results in terms of the thickness distributions and flange geometry of the cup.
It is seen from the Figures 9 and 10 that the predicted thickness distributions and flange geometry matches well with the experimental results. Then, Sener et al. [19] expanded the study [18] and studied the variation of anisotropy during plastic Experimental and predicted cup profiles from the fourth-order polynomial and Yld96 criteria for AA2090-T3 [14].

Figure 6.
Drawn cup [14]. deformation experimentally and numerically. They carried out FE simulations of same industrial part at different plastic strain levels (0.2%, 2%, 5%, and 18%) and compared P 4 predictions with experimental data. Figure 11 shows the comparison of the predicted thickness distributions along the three directions with experiment.
It is seen from Figure 11 that different thickness predictions were obtained at different plastic strain levels. After the comparison of the predicted thickness results with experiment, researchers eliminated two strain levels and then they investigated the flange geometry results (Figure 12).     From the comparison of the predicted and experimental flange geometry results, it is seen that numerical results were matched well with the experiment.

Conclusions
In the present study, homogeneous anisotropic polynomial yield functions, their types, and application areas in the metal forming process were investigated. In the literature, generally anisotropic yield functions derived from linear transformation approach are used. These functions have high modeling capability and they could be used for different materials. However, yield functions based on linear transformation approach have some disadvantages. They have complex coefficient identification procedure and nonlinear formulas. Therefore, calculations of the first and second order gradients of these models are difficult and it causes to difficulties in the implementation of the models into FE codes. On the other hand, polynomial yield functions have a generalized, simple structure and derivatives of these functions could easily calculated.
It is seen from the studies carried out in the literature that researchers generally use the fourth and the sixth order polynomial functions to model of the anisotropic behavior of the materials. Based on the results obtained from the studies performed in the literature, the following conclusions could be drawn: a. Homogeneous polynomial yield functions have high modeling capability in the description of anisotropic behavior.
b. Homogeneous polynomial yield functions could be used for both plane stress and generalized stress state. This provides the flexibility to the polynomial yield criteria.
c. Sixth-order homogeneous polynomial yield function could predict six or more ears in a deep drawn cup.
d. Homogeneous polynomial yield functions could model body centered and face centered cubic materials without the need of any exponent related to crystallographic structure.
e. Apart from the linear transformation approach, polynomial models may not satisfy convexity requirements for each stress state. Therefore, the user should consider convexity conditions and has to investigate the model parameters in terms of convexity and positivity conditions.
f. The modeling capability of the fourth-order polynomial yield function is similar with Yld96 yield function, whereas predictions of the sixth-order polynomial yield function close to Yld2004-18p model.