Modeling plastic anisotropy evolution of AISI 304 steel sheets by a polynomial yield function

In this study, a numerical model for the evolution of plastic anisotropy is investigated for the purpose of stamping method design by Finite Element (FE) analysis and proved experimentally via process simulations of a cold-rolled austenitic stainless steel (AISI 304) sheet. The plastic anisotropy of the sheets is described with a fourth-order homogenous polynomial yield function and this modelling approach is enhanced by plastic strain dependent material coefficients. Tensile tests of coupon specimens taken along the different directions from rolling direction, and flow strength and deformation anisotropies are described with the planar variations of yield stress and plastic strain ratio computed at four plastic strain levels (0.002, 0.02, 0.05 and 0.18). A new numerical approach is, then, applied to identify polynomial coefficients ensuring an orthotropic positive-definite, convex yield surface with a well-defined stress gradient at every loading point on plane stress subspace. The developed computational model is implemented into general purpose explicit FE analysis software Ls-Dyna by a user-defined material model subroutine (UMAT) and applied in the stamping simulation of AISI 304 steel rectangular cups for the house-hold applications. The computed thickness distributions and the flange geometries were compared with measurements and it was observed that the best predictions were done with material parameters at %5 plastic strain level.


Introduction
A proper description of material anisotropy has a great influence on the prediction accuracy of sheet metal stamping simulations commonly based on the finite element (FE) analysis as a computational methodology, and it is well known that the plane stress orthotropic yield functions have an apparent impact on these issues among other numerical parameters. The plane stress orthotropic yield functions contain several modeling parameters utilizing specific material deformation paths, usually along proportional deformation paths, for the initial yield stresses and Lankford parameters (r values or plastic strain ratio) along different directions. However, experimental studies proved that these mechanical properties are not constant and change during plastic deformation. Hu [1] observed that the plastic strain ratios of steels which have strong textured structure decrease with the increase of plastic strain. Hu [2] then investigated this relationship for randomtextured, deep drawing and poor drawability sheets and determined that the degree of texture in the tested material affects the alteration of the anisotropy with the plastic deformation. Truszkowski [3] determined the r values of hot rolled nickel and aluminum sheets at different degrees of straining and concluded that fiber texture appearing during the deformation influences the material anisotropy.
Savoie et al. [4] carried out experiments on 5XXX series aluminum alloy and determined that the variation of the r value is associated with Luders bands. Safaei and De Waele [5] investigated the evolution of anisotropy for an interstitial-free steel sheet (IF300) and observed the variations in r-values of material at different plastic work values. Taghvaipour et al. [6] performed experiments on Ti alloys and determined that r-value of the material is constant in the elastic region, however it changes rapidly at the beginning of yielding. These results lead to a modeling need to account this variation of material anisotropy to be taken into account in FE analyses. Several studies have been carried out on the evolution of anisotropy with plastic strain in the literature. Cvitanic et al. [7] modeled the anisotropic behavior of DC06 sheet using Hill48 [8] and Karafillis-Boyce (1993) [9] yield criteria. They defined the coefficients of these criteria based on the plastic strain and predicted iso-error maps successfully. Zamiri and Pourboghrat [10] included the variation of r values with plastic strain to Hill48 criterion and applied it to the bulge test of niobium sheet. They correctly predicted the strain localization of the formed part with this evolutionary yield function. Unlike the previous study [10], Volk et al. [11] considered the alteration of both yield stress ratios and r-values with plastic strain and described Hill48 coefficients depend on plastic strain. They carried out the bending simulations of DP800 sheet and improved the springback prediction accuracy with this material model. Safaei et al. [12] proposed the non-associated Yld2000-2d [13] model to represent the anisotropy of DC06 material and defined model parameters depend on the plastic strain. They showed that the model could successfully predict the evolution of r-values and yield stresses during plastic deformation. Lian et al. [14] proposed non-associated Hill48 plasticity model depend on the anisotropy evolution and used this model to predict forming limit diagram (FLD) of ferritic stainless steel. They successfully validated the model with experimental results. Aretz [15] examined the effect of anisotropy evolution on the prediction of localized necking and showed that it has a significant impact on the prediction accuracy of FLD. Bandyopadhyay et al. [16] performed the FE simulations with and without considering the anisotropy evolutions for limiting dome height and cup drawing tests. They showed that the prediction accuracy of FE simulations could be improved by considering the effect of anisotropy evolution. Kuwabara et al. [17] identified the coefficients of Yld2000-2d material model depend on plastic strain and applied this yield function to the hole expansion forming simulation of aluminum alloys. They declared that FE simulations performed with this material model better correlated with experiments compared to the isotropic hardening models. Yoon et al. [18] investigated the effect of the evolution of anisotropy on the earing prediction in cylindrical cup drawing. They used two anisotropic yield functions (Yld2000-2d and Cazacu-Plunkett-Barlat [19]) to define the anisotropic behavior of AA5042-H2 sheet material and also represented the coefficients of these yield criteria as functions of equivalent plastic strain. They indicated that Cazacu-Plunkett-Barlat criterion based on plastic strain could successfully predicted the cup profile and height. Similar studies were carried out by Wang et al. [20] and Choi et al. [21]. Cai et al. [22] carried out an extensive study related to the variation of anisotropy with plastic strain. They defined both the exponent and the anisotropy coefficients of the Yld2000-2d yield criterion as a function of plastic strain. Researchers performed the numerical simulation of hydro-bulging and cup drawing tests and showed that the developed model improved the prediction accuracy of FE simulations.
Numerous studies have been carried out with these evolutionary models, but it is seen from the literature that generally quadratic and non-quadratic yield functions have been used, while the usage of the polynomial yield functions is limited in these studies. Firstly, Gotoh [23] investigated a polynomial yield criterion in the literature and could successfully define the anisotropic behavior of aluminum alloys. However, he didn't take notice of convexity in his criterion. Therefore, this yield criterion couldn't gain popularity in metal forming. Cazacu and Barlat [24] developed the sixth degree polynomial criterion and applied to aluminum alloys. They obtained successful results, however this criterion has complex coefficient identification procedure which is based on nonlinear formulas. Hu suggested fourth-order polynomial models for plane stress [25] and three dimensional stress states [26]. He obtained successful results for highly anisotropic aluminum alloys, but some oscillations could be observed in both stress and strain ratios with these models [27]. Yoshida et al. developed sixth-order polynomial criterion and represented the plastic behaviors of IF steel and high strength steel sheets [28]. Then they investigated the evolution of anisotropy of AA6022-T43 with this model and validated the prediction results at each plastic strain level [29]. Although applications of the polynomial type models have been limited in the literature, as a matter of fact, polynomial type functions have more direct formulation and they are simpler when compared the other models. In addition to that their gradients can be obtained simply due to the polynomial structure of these models and this is an important advantage for implementation into a FE code. Therefore, in this study the fourth-order polynomial material model was used and implemented in the dynamic-explicit FE code Ls-Dyna by writing a user defined material subroutine (UMAT). The aim of this study was to investigate the effect of plastic strain on the evolution of anisotropy and AISI 304 stainless steel was selected as test material in this study. For this purpose, the coefficients of the polynomial yield function were estimated at four different plastic strain values (%0.2, %2, %5 and %18) by using the results obtained from the uniaxial tensile tests carried out along different directions and the in-plane variations of yield stress ratios and r values were predicted at same plastic strain values. Then, deep drawing of an industrial part was considered as an application and FE simulations of the process were performed with these coefficients. Finally, the numerical results were compared with experiments and the effect of plastic strain on material anisotropy was investigated. As stated above, AISI 304 stainless steel sheet was preferred in this study. This material has wide usage area in different industries such as nuclear, defense, food, marine industries etc. [30]. 304 stainless steel is a metastable material that austenite transforms to martensite by plastic deformation and this provides advantage in mechanical properties of this material. This transformation delays crack initiation, necking or fracture and it improves formability of the material [31].
This study divided into six sections. In Sect. 2, the experimental studies which is performed to define material anisotropy is presented. In Sect. 3, the fourth-order polynomial yield criterion and its identification procedure are introduced. In Sect. 4, material anisotropy is assessed with polynomial yield criterion at different plastic strain levels and its application on deep drawing simulations of an industrial part is described. In Sect. 5, the numerical and experimental results are compared and the prediction capability of the yield criterion is evaluated. In Sect. 6, the main findings are summarized and the study is concluded.

Experimental studies
Experimental works in this study are performed in two phases. Firstly, uniaxial tensile tests are conducted with coupon specimens taken from three orientations with respect to rolling direction, and flow stress and anisotropy coefficients are determined. Next, rectangular cup drawing tests are performed using a 160-ton-capacity doubleaction hydraulic press with the actual stamping tooling for production.

Uniaxial tensile test
Uniaxial tensile tests were carried out using specimens machined according to ASTM E8M standard on a 100kN capacity Shimadzu universal tensile test machine with a quasi-static strain rate of 0.008/s. The specimens were tested along the three directions (rolling-RD, diagonal-DD and transverse-TD) to evaluate the anisotropy of the material and experiments were repeated three times for each direction. Flow curves in three directions and biaxial flow curve are shown in Fig. 1. Information about biaxial flow curve of the material was taken from Yadav [32].
Lankford coefficient is an indicator of plastic anisotropy and this parameter depends on the plastic strain. Therefore, the variations of Lankford coefficients with plastic strain were investigated to determine the evolution of anisotropy with plastic deformation. Figure 2 shows instantaneous variation of r values with increasing plastic strain p for three directions [33]. It can be seen from Figs. 1 and 2 that both yield stresses and r values of the material are sensitive to the plastic strain and direction.

Rectangular cup drawing test
After performing the uniaxial tensile tests, rectangular cup drawing tests were carried out to evaluate the plastic deformation of the material in multiaxial stress state. Deep drawing of a rectangular part was taken into account as application in the article, since non-uniform material flow is observed and different types of forming defects (tearing and wrinkling) could be occurred at different locations of the part. Experiments were performed on a 160  ton capacity hydraulic press. The process consists of two stages: clamping and drawing. In the clamping stage, the die moves downward, and the blank is clamped between the die and the binder by applying the binder force (BF). In the drawing stage, the die and the binder move down and the blank is forced into the die cavity by the stationary punch. Mineral oil was applied on blank-die and blankbinder surfaces. In the experiments, the part was formed with 20 mm/s die velocity and 340kN BF was applied during 80 mm die travel height. Tools and the drawn part are shown in Fig. 3.

The fourth-order polynomial yield function
In this study, the fourth-order polynomial type material model (Poly4) was used. The model was proposed by Soare [34]. The yield criterion (F) can be given by.
in which, eq and 0 denote the equivalent stress and yield stress, respectively. eq for plane stress state can be expressed as follows: where a 1 , a 2 , a 3 …a 9 are coefficients describing the material anisotropy, 11 , 22 and 12 represent the normal and shear stresses, respectively. Positivity and convexity conditions have taken into consideration in the identification procedure and upper and lower bounds on coefficients are derived for satisfying of these conditions. procedure for the coefficients of Poly4 yield criterion. denotes the yield stress ratio along direction with respect to the rolling direction ( 0 ) and it is determined by dividing 0 = ∕ 0 . In this study, a separate method was applied from Soare in the determining the coefficients of a 6 and a 8 . Summary of the identification procedure are given below: Step 1: The coefficients a 1 , a 2 , a 4 and a 5 are calculated analytically and the equations are given below: Step 2: Positivity and convexity conditions for b are checked and then the coefficient a 3 is determined. The inequalities for b and the equation for a 3 are given below: where M 1 and M 2 are the roots of second order equation depend on the coefficient a 3 . Detailed information for positivity and convexity conditions can be found in [34].  Fig. 3 Tools and the drawn part [42] Step 3: Positivity and convexity conditions for r 45 are checked and then the coefficient a 9 is determined by using the expressions which are given below: Step 4: The coefficients a 6 and a 8 are determined by optimizing the difference between theoretical and experimental values (error or distance function).
w h e re ω (i) 1 a n d ω (i) 2 a re we i g h t s fo r a n d r , r e s p e c t i v e l y.
Step 5: Interval for the coefficients a 6 and a 8 is checked by using the inequalities (11) to satisfy positivity and convexity conditions. In this step, a different method was applied from Soare's identification procedure. In our method, bounding domain of the coefficients a 6 and a 8 was divided into 100 pieces and the error function was calculated for each pair a 6 , a 8 by using Eq. (11). A minimum error value was selected (10 -30 ) for determination of a 6 and a 8 coefficients. When an error value which is lower than a minimum value is found, corresponding the coefficients a 6 and a 8 are accepted as optimum coefficients.
Step 6: The coefficient a 7 is determined by using Eq. (12): 15 0 and 75 0 were taken as input in the determination of Poly4 coefficients and these angles were determined according to Eq. (13) and Eq. (14) in this study.
The flow chart which summarizes the identification procedure of Poly4 yield criterion is given in Fig. 4.
The representation of the angular variations of the stress and strain ratios is an indicator in the evaluation of the prediction capability of a material model. Therefore, directional variations of stress ratio and Lankford

Application studies
In this section, the evolution of anisotropy with plastic strain was investigated and the study was conducted on two applications. Firstly, the parameters of the yield function were determined at different plastic strain levels and   Fig. 5 Material orthotropic frame (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22) and loading coordinate system (xx-yy) they were used to predict the directional variations of the yield stress ratios and r values. Then, rectangular cup drawing simulations were carried out with these coefficients to evaluate the variation of the prediction capability with plastic strain.

Prediction of the directional variations of yield stress ratios and r values
It has been observed from Fig. 1 that the proportionality between the flow curves (RD, DD and TD directions) distorted during plastic deformation. It is expected from this result that the coefficients of the yield function and correspondingly the shape of the yield surface change with plastic strain. Therefore, yield stress ratios and r values at different plastic strains p were determined and parameters of the yield function were estimated at specific plastic strain values. The values of p were taken to be 0.002, 0.02, 0.05 and 0.18 in the study. Experimental results and calculated Poly4 coefficients at the selected plastic strain values were given in Tables 1 and 2, respectively. Planar variations of yield stress ratio and r value at specific plastic strain values were predicted according to Poly4 yield criterion and the predicted results were compared with the experimental results. Comparisons of theoretical and experimental results for yield stress ratio and r values were shown in Fig. 6 and Fig. 7, respectively.
As it is seen from the Figs. 6 and 7 that Poly4 yield criterion accurately predicted both the angular variation of the yield stress ratios and of the Lankford coefficients at selected plastic strain levels. After determination of Poly4 coefficients for each plastic strain level, these coefficients were substituted into Eq. (2) and then yield function was obtained by using Eq. (1). In these equations, normal and shear stresses were normalized by dividing 0 . Then normalized yield surface for material was predicted and contours of the yield surface were plotted at selected plastic strain levels to observe the alteration in the shape of the yield surface with plastic strain (Fig. 8).
It is seen from Fig. 8 that the shape of yield surface varies with increase in plastic strain and this variation is more pronounced; in particular, in the biaxial tensile and biaxial compression regions.

Finite element analysis
Deep drawing of a rectangular cup is analyzed using general-purpose Ls-Dyna FE solver and, since the Poly4 yield criterion is not an option for the plasticity capabilities in this software, an orthotropic rate-independent plasticity model based on Poly4 yield surface is numerically implemented using UMAT subroutine of Ls-Dyna.
The implementation is composed of time-integration of plasticity equations at each integration point of FE mesh [35]. UMAT subroutine provides the strain increments and requires the calculation of the corresponding stress increments as input, on the return, back to main routine. In this context, a backward Euler discretization of the differential equations is employed in order to obtain algebraic incremental plasticity relations [36]. Due to nonlinearity of resulting algebraic equations, an iterative solution scheme is applied by a successive substitution and updating the backstress and yield function during iterations. Following the convergence of plastic strain increment, stress and strain tensors at the end of increment are also updated [30]. After implementation of the yield criterion, FE model of the process was created. Quarter model was created due to the symmetry. Shape and dimensions of the tools (die, punch and binder) are shown in Fig. 9. The blank was meshed with 2 mm quadrilateral shell elements and full-integrated element formulation with five integration points through the thickness direction were used in the model.
The blank was modeled as an elastoplastic object, while tools were modeled as rigid bodies. Forming-one-way-surface-to-surface contact algorithm was used in the FE model. Friction coefficient was assumed to be 0.125 for blank-punch interface due to dry condition, while this parameter was assumed to be 0.05 for blankdie and blank-binder surfaces due to lubricated condition [37]. Swift's hardening law was used in the definition of the hardening curve [38].
where K is the strength coefficient, 0 is the initial plastic strain and n is the strain hardening exponent. The parameters of the Swift hardening rule were identified from (25) = K 0 + p n experimental flow curves by curve fitting technique. Nonlinear least squares and trust-region algorithm were used. The identified parameters for each direction and comparison of predicted flow curve by Swift model with experimental data were given in Table 3 and Fig. 10, respectively. Die displacement was defined in the FE model. 85 kN BF was used in analyses due to the quarter model. It was started from zero and reached to 85 kN within 0.1 s, while die is started to move in 0.115 s in order to suppress oscillation and thus reduce the kinetic energy in the blank. Curves related to proses parameters were shown in Fig. 11.

Results and discussion
The computed results from FE analyses were compared with experiments. Thickness distributions and the flange profile results were considered in the comparison. Thicknesses of the formed parts were measured along three directions by using micrometer. Comparisons of predicted and measured thickness distributions for each direction were shown in Figs. 12, 13 and 14. It is seen from the figures that significant variations were observed in the punch radius region. Because, in this region material is exposed to stretching and it can be observed from the Fig. 8 that the shape of the yield loci is distorted especially in biaxial tensile and compression regions. Similar results were observed by Kuwabara et al. [39], Hill et al. [40] and Hill and Hutchinson [41]. In addition to that it could be understood from the figures that the computed thickness distributions along three directions from FE simulation performed using the material data at %5 plastic strain value matched with the experiment better compared to other three simulations. More thinning was predicted from FE simulations carried out with experimental data at % 0.2 and % 2 plastic strain values. Therefore, the predicted results from these two plastic  Flange profile of a drawn cup is not fully symmetrical due to non-uniform material flow in deep drawing process. Plastic anisotropy is one of the primary reasons caused to this phenomenon. Therefore, flange profile was considered for evaluation the effect of anisotropy in this study. After the deep drawing experiments, the formed parts were scanned and flange geometries of the parts were determined experimentally. Figure 15 shows comparison of experimental flange profile with FE predicted flange profiles. It was observed that both of the FE computed flange profiles are compatible with the experiment. However, only slight differences between predicted flange profiles and experiment were observed along the RD and TD. More draw-in was noticed in FE prediction with %5 plastic strain, while less draw-in was noticed in FE prediction with %18 plastic strain. The amounts of percentage error between FE and experimental results along the RD and TD directions are given in Table 4.
In the present work, a novel coefficient identification program is applied. As previously stated, our method is based on the error analysis in a bounded region. In the literature, numerical optimization methods are applied to determine the coefficients and these methods require gradient information of objective and constraint functions. However, functions are nonlinear and determination of the gradients are computationally expensive and impractical. Our method doesn't require any gradient information and therefore solution could be obtained easily.

Conclusions
In this study, the evolution of anisotropy with plastic strain was investigated for AISI 304 stainless steel. Anisotropy of the material was defined with the fourth-order polynomial material model and the model was incorporated into dynamic explicit FE code Ls-Dyna by UMAT subroutines. In order to evaluate the evolution of the anisotropy, the parameters of the yield criterion were estimated at different plastic strain levels and they were used to predict the angular variations of the yield stress ratios and r values. Then rectangular cup drawing simulations were performed with these coefficients together with developed user material subroutine UMAT. The predicted thickness   Experimental and numerical results obtained from this study could be summarized as follows: • From the uniaxial tensile test results, it was observed that the initial proportionality of the hardening curves and Lankford coefficients of the material change with plastic strain. This indicates that material anisotropy evolves during plastic deformation. • The optimized parameters of Poly4 yield criterion at four different equivalent plastic strains were used in the prediction of the directional variations of the yield stress ratios and Lankford coefficients and excellent agreement was obtained between the experiment and computational results. This result shows that Poly4 yield criterion has high predictive capability at various levels of plastic deformation. • Yield surfaces for the material were predicted by using Poly4 yield criterion at selected plastic strain values and it was realized that the shape of the yield loci evolves with plastic deformation. In this case, shape distortion of the yield loci is clearly observed in the biaxial tensile and biaxial compression regions, while little variation occurs in the uniaxial tensile and uniaxial compression regions. This result reveals distortional hardening phenomenon. • Rectangular cup drawing simulations were performed with the Poly4 coefficients estimated from different plastic strain values and the computed results were compared with experiments. The comparisons showed that the most successful predictions were done with FE simulation using material data obtained from %5 plastic strain value. This means that the material properties at the initial yield point can't completely represent the material behavior during the forming.
In the present work, the parameters of the yield function were determined at four different plastic strain levels and prediction capability of the function is evaluated. In the future, instantaneous variation of the model parameters with plastic strain should be considered and a polynomial constitutive model based on plastic strain should be developed to improve the prediction accuracy of FE simulations.