Effect of Anisotropic Yield Functions on Prediction of Critical Process Window and Deformation Behavior for Hydrodynamic Deep Drawing of Aluminum Alloys

: Owing to the reduction of rupture instability and the avoidance of wrinkle defect, the hydrodynamic deep drawing (HDD) process is gradually becoming attractive for fabricating lightweight and complicated products. Meanwhile, since metallic materials present anisotropic deformation behavior, it is necessary to select an appropriate constitutive model for the prediction of plastic deformation behavior of applied material with high precision. In the present research, several anisotropic yield criteria, namely, Hill’48, Yld2000-2d, and BBC2005, were implemented to investigate the effects of yield functions on the prediction accuracy of the critical process window and deformation behavior for the HDD process of 2024 and 5754 aluminum alloys. Material constants in the yield criteria were determined by applying uniaxial and equi-biaxial tension tests and optimizing an error-function using the Levenberg–Marquardt algorithm. Furthermore, the process window diagram was computed utilizing the stress analytical model combined material properties with workpiece geometrical features. Numerical simulation results of predicted material anisotropic parameters, process window, and HDD deformation for aluminum alloys were compared with the experimental data. Through the comparison of diverse yield criteria based on materials’ anisotropic coefficients, critical process window prediction, earing profile, and thickness distribution, it was revealed that the Yld2000-2d and the BBC2005 yield criteria can offer more precise models of material behavior in planar anisotropy properties for the HDD process of 2024 and 5754 aluminum alloys.


Introduction
Aluminum alloy is continuously being applied in the aerospace and automotive industry sectors because of their lightweight characteristic and excellent combination of strength and corrosion resistance. Complex aluminum alloy sheet metal parts with complicated shapes, large curvature variations, and deep cavities are constantly applied in the design of advanced aerospace structural parts. The deep drawing process is one of the significant manufacturing processes capable of fabricating thin-walled components which can not only reduce the cost but also enhance the performance of the products [1]. However, the formability of aluminum alloy is inadequate for the low ductility at room temperature and the high degree of springback. As a result, many aluminum alloy sheet metal parts cannot be fabricated using a one-step deep drawing process due to the complex structure and poor formability. For these reasons, a multi-step deep drawing process or partial forming with a welding method are applied to manufacture the complicated sheet metal parts, which results in defects such as excessive process redundancy, low forming accuracy, and poor forming quality. The utilization of flexible media in the deep drawing process has received increasing attention due to the multitude of its advantages such as improving sheet forming performance, forming efficiency, and convenient in mold changing. Hence, to surmount the difficulty in forming thin-walled aluminum alloy components, a better choice would be the hydrodynamic deep drawing (HDD) process [2,3].
The hydrodynamic deep drawing is a process that uses liquid as force transmission medium instead of conventional rigid die and makes the blank fit the punch under the pressure of liquid. A schematic of the HDD process is shown in Figure 1. During the forming period of the HDD process, the workpiece is pressed into the chamber through the punch. Meanwhile, a cushion of pressurized viscous fluid from the chamber is generated to support the noncontact region of the workpiece. Due to the fact of such external support, the provided through-thickness compressive stress contributes to delay the onset of tensile instability as well as to reduce the occurrence of wrinkles. Moreover, the workpiece bulges backward under the effect of pre-bulging pressure before the punch contacts the blank, which is a valid approach to avoiding the origination of the wrinkle on the unsupported region [4] as shown in Figure 2. The HDD process takes possession of many virtues in comparison with the conventional deep drawing (CDD) process including high dimensional accuracy, desirable surface quality, improved cold formability, less springback, and shortened manufacturing cycle [5,6]. Above all, the HDD process has superior performance on the avoidance of wrinkle phenomenon due to the controllable fluid pressure. Obviously, the cavity pressure of the chamber liquid is one of the most crucial parameters in the HDD process. As a result, determining the critical pressure and process window plays an important role in forming thin-walled workpieces without ruptures and wrinkle defects. In recent years, several researchers have devoted their investigations to the prediction of the critical pressure in the HDD process. Meng et al. [7] performed a series of finite element analyses for the HDD process of aluminum alloy rectangular boxes with wide flange. Meanwhile, the effects of cavity pressure on forming quality and precision were explored according to theoretical analysis and experiments, and the process window of cavity pressure was established through stress analysis of the corner and straight regions. Bagherzadeh et al. [8] developed analytical models to investigate stress distribution and instability conditions in hydro-mechanical deep drawing (HMDD) of cylindrical bimetallic cups. It was demonstrated that the fluid pressure window for successful part forming could be rapidly predicted with reasonable accuracy by the analytic model compared to lengthy and costly finite element analysis or experimental trial and error. Wang et al. [9] computed the process window for the HDD process of a composite conical part with double concave features. The optimal pressure loading loci were promoted based on the feature of workpiece.
Numerical simulation significantly contributes to reducing the overall production costs and inaccuracy in sheet metal forming industries by shortening both the time required in the R&D stage and the time needed for a corresponding production implementation [10,11]. Furthermore, the increasing application of numerical simulation in the field of sheet metal forming is beneficial to solve problems in the manufacturing of qualified formed parts [12,13]. Hence, the finite element method (FEM) has become an indispensable analysis tool in the manufacturing processes design. Liu et al. [14] proposed an analytical approach to obtain a proper prediction of liquid pressure for the hydroforming of a curved surface shell by the combination of numerical simulation and theoretical models. Moreover, the effects of different parameters for the proper liquid pressure were discussed using multiple linear regression analysis. Lang et al. [15] optimized the forming parameters in FEM of sheet hydroforming based on the identification of parameters for constitutive models by inverse modeling. Based on the optimized parameters, the sheet hydroforming process can be analyzed more accurately to improve the robust design. Gorji et al. [16] investigated the forming of conicalcylindrical cups in the HDD process using numerical simulation and experiments. It was also illustrated that for the pressure path with a specified maximum amount, the workpiece was formed adequately with minimum sheet thickness reduction.  The consistency of the decisions built on the FEM basis is highly dependent on the attained degree of the physical parameters and numerical accuracy of the simulation [17][18][19][20][21][22]. Many researchers involved in this domain focused their efforts on improving both the quality of the theoretical models implemented in the simulation programs and the efficiency of their applied computational strategies [23,24]. In this regard, the constitutive model providing an accurate description of the plastic anisotropy of HDD process has been the subject of numerous studies. Hashemi et al. [25] proposed a process window diagram to predict the minimum and maximum critical pressure in HDD of conical cups using the Hill'48 yield criterion numerically. Azodi et al. [26] developed analytical models based on the Barlat'89 and Hill's non-quadratic yield criteria in the HMDD process to predict the maximum permissible fluid pressure with assuming plane strain tensile failure. Jalil et al. [3] developed theoretical models using the Barlat'89 yield criterion to analyze the critical bursting pressure in HDD of single-and double-layer sheets. According to the abovementioned research, it is generally believed that the Hill'48 and the Barlat'89 yield criteria can roughly describe the anisotropic plastic property of sheet metal, such as carbon steel and stainless steel, yet is unreasonable for the HDD process of aluminum alloys. Meanwhile, limited study has been reported for the applicable implementation of advanced anisotropic yield criteria for aluminum alloys in sheet hydroforming applications.
In this paper, the applicability of diverse anisotropic yield criteria, Hill'48, Yld2000-2d, and BBC2005, in theoretical analysis and numerical simulation of the HDD process of 2024 and 5754 aluminum alloys was carried out. A comparison between the numerical and experimental results is presented to evaluate the effectiveness of the selected yield criterion as well as the effect of the yield criterion on the prediction accuracy of process window and deformation behavior for HDD process of aluminum alloys in terms of forming force variation, thickness distribution, and earing phenomenon.

Anisotropic Yield Functions
The yield criterion expresses a relationship among the stress components in the transition from the elastic to the plastic regime which is a part of the most crucial criteria used in numerical simulation [27,28]. Due to the multi-pass rolling and heat treatment in the production process, aluminum alloy sheet presents remarkable anisotropic characteristics. The anisotropy of aluminum alloy has certain effects on the fracture position, earing, and forming limit of the part fabricating by HDD process. Therefore, the premise of accurately describing the deformation of aluminum alloy in the HDD process is to fully consider the influence of sheet anisotropy.
It can be mentioned that the most frequently applied yield criteria for the simulation of sheet metal forming process, Hill's quadratic functions, have been fitted to equally many test results. It is noticeable that a constitutive model making full use of the available eight material parameters obtained from uniaxial and biaxial tests would be deeply in demand. Hence, one four-parameter anisotropic yield criterion namely, Hill'48, two eight-parameter anisotropic yield criteria namely, Yld2000-2d and BBC2005, were employed in this research as well as the numerical approach to evaluate material anisotropy and predict deformation for HDD process of aluminum alloys.

Hill'48 Yield Function
In 1948, Hill [29] proposed a constitutive formula for the plastic yielding and deformation of anisotropic metals at a macroscopic level. The quadratic yield criterion is given by Equation (1): where f is the yield function. σxx, σyy, and σzz are the stresses in the rolling, transverse, and thickness directions, respectively. τxy, τyz, and τzx are the shear stresses in the xy, yz, and zx planes, respectively. F, G, H, L, M, and N are the constants that describe the anisotropy of the material. Under the assumption of plane stress condition ( σxx = τyz = τzx = 0 ), the yield criterion can be simplified as follow in Equation (2): In this study, the identification of the anisotropy parameters for Hill'48 model was carried out by applying two different methods. The classical approach uses the r values from three uniaxial tension tests (0°, 45°, and 90° to the rolling direction), which is labeled Hill'48-R. The second approach resorts to implementing the yield stresses from three uniaxial tension tests (0°, 45°, and 90° to the rolling direction), which is labeled Hill'48-S.

Yld2000-2d Yield Function
To alleviate the drawback of previous yield criteria in the description of yielding behavior of aluminum alloys, Barlat et al. [30] proposed the Yld2000-2d yield function which has eight anisotropy coefficients including tension yield stresses (i.e., σ0, σ45, σ90) and r values (i.e., r0, r45, r90) along the rolling direction, 45°, and the transverse directions as well as the yield stress σb and rb under the balanced biaxial tension condition.
The equation for the yield locus of Yld2000-2d is given by Equation (3): with where m is the Barlat exponent relevant to the crystal structure of the material. In this research, m was 8 because of the FCC structure of the aluminum alloy. φ was the sum of the two anisotropic equations, i.e., φ' and φ". X'i and X"j ( i, j = 1, 2 ) are the principal values of the matrixes X' and X" which is expressed by Equation (5): The elements of X' and X" are achieved from the following linear transformation of the Cauchy stress in Equation (6): Here, L' and L" are defined as follows in Equations (7)- (8): where 1~ 8 are eight anisotropy coefficients. In the model, the coefficients L ij ' and L"ij can be described by the set of coefficients 1~ 8 . All these coefficients are defined independently. The r values and the yield stresses from uniaxial tensile tests of the specimens taken along 0°, 90°, and 45° to the rolling direction were taken into account. Besides, the equi-biaxial yield stress and r value were required for determining the remaining coefficients of the yield model.

BBC2005 Yield Function
The BBC2005 yield function is a plane stress yield criterion developed to describe the deformation behavior of the orthotropic sheet metal [31]. The yield surface function of BBC2005 yield criterion is defined as Equation (9): where σ ij is the plane stress tensor. Y > 0 is an arbitrary reference yield stress. σ is the BBC2005 equivalent stress, given by Equation (10): where a and b are undetermined coefficients. The value of the integer exponent k can be adopted according to the crystallographic structure of the sheet metal: k = 3 for BCC alloys and k = 4 for FCC alloys. Γ, Λ, and Ψ are functions counting the planar components of the stress tensor which are given by Equation (11): where κ 1 , κ 2 , κ 3 , κ 4 , κ 5 , and κ 6 are yield criterion constants.

The Critical Process Window of HDD Process
In the HDD process, we proposed that the curve of chamber pressure with forming height is described as the cavity pressure loading locus. Only when a reasonable cavity pressure loading locus is designed within a certain range can the crack and wrinkle defects be avoided, and this area is defined as "critical process window" for the HDD process.
To analyze the critical process window of the sheet metal HDD process, several assumptions for calculating the critical cavity pressure are put forward as follows [32]: (1) The volume of the workpiece is constant throughout the whole process; (2) The radial and tangential directions are considered as principals; (3) A nonlinear hardening plastic behavior of the material is assumed based on the power law. Hollomon's equation expressed as Equation (12): σ̄= Kε̄n (12) where ε , K, and n are the equivalent strain, strain hardening coefficient, and strain hardening exponent, respectively. (4) The blank anisotropy can be depicted by the mean anisotropy coefficient ( r̅ ) as shown in Equation (13): where r 0 , r 45 , and r 90 are Lankford coefficients describing the anisotropy. The workpiece is divided into three regions according to the stress state including the corner flange area between the blank holder and the die (Region A), the curved region in contact with the pressurized fluid (Region B), and the area tightly compressed onto the surface of the punch due to the high pressure in chamber (Region C) as depicted in Figure 3. In Figure 3, for an axisymmetric radial element in the flange area, the transformation of the thickness of the corner flange was assumed to be neglected. Considering the radial and tangential directions as principal directions in such an element, the equilibrium formula for the flange area of the radial direction is shown in Equation (14): where σr and σθ are the radial stress and circumferential stress. r is the radius of material element. t is the material thickness. µ is the friction coefficient between the workpiece and the blank holder. p is the fluid pressure under the flange along the vertical direction. Rb is the current radius of the outer flange. rp is the profile radius of the punch shoulder. Rp is the radius of the punch. ρ is the radius of the part curvature around the die corner. hp is the current cup height. The stress state of Region B is similar to the flange area, while there is no friction force in this region due to the workpiece being completely separated from the die. The equilibrium formula in the normal direction of the workpiece in Figure 3 is depicted in Equation (15):

Analytical Modeling Applying the Hill'48 Yield Criterion
The associated flow rule was adopted which is expressed as follow in Equation (16): where λ is the plastic multiplier and f is the plastic potential, defined as a scalar function described by the yield criterion. By applying the Hill'48 yield criterion, the flow rule is written in Equation (17): where dε1, dε2, and dε2 are the plastic strain incremental components along the principal directions.
Considering the radial and tangential directions as principal directions, Equation (17) results in: Ignoring the thickness variations (dε r = dε θ ), the effective strain can be obtained as Equation (19): By substituting Equations (18) and (19) into Equation (14) and integrating, the radial stress of the workpiece in Region A is: where r (1) is the current radius in Region A, and r 0 is the initial radius of a supposed point in the flange region that moved to the current point with radius r (1) .
Adopting the volume constancy condition, the relationship between r (1) and r 0 (1) in Region A can be defined as the following in Equation (21) with the condition The current flange radius is determined as Equation (22): With regard to the boundary condition, the radial stress of the workpiece based on the Hill'48 yield criterion in Region B is calculated by Equation (23): where r (2) is the current radius in Region B.
The relation between the current and initial radius of the material element in Region B is defined as following Equation (24) where = arccos R p r t ρ

Analytical Modeling Applying the Yld2000-2d Yield Criterion
According to the plane strain assumption, the effective stress as a function with the ratio of tangential to radial stress ( γ = σ 2 σ 1 ⁄ ) can be written as Equation (25): where Substituting Equation (25) into the flow rule equation results in Equation (27): with Neglecting the thickness strain ( dε 1 = −dε 2 ), the strain work is written as Equation (30): The effective strain can be obtained as Equation (31): By substituting Equations (30) and (31) into Equation (14) and integrating, the radial stress of the workpiece in Region A is: Considering the plane stress assumption, γ can be obtained by solving m(γ) n(γ) = −1 ⁄ .
Therefore, the radial stress of the workpiece element based on the Yld2000-2d yield criterion in Region B is:

Analytical Modeling Applying the BBC2005 Yield Criterion
The effective stress of the BBC2005 yield function can be written as Equation (34): where Substituting Equation (9) into the flow rule equation results in Equation (36): The functions j(γ) and k(γ) are defined as Equation (38): The strain work solution of the BBC2005 is similar to the Yld2000-2d one; it is obtained as: By substituting Equation (39) into Equation (14) and integrating, the radial stress of the workpiece in Region A is: Therefore, the radial stress of the workpiece element based on the BBC2005 yield criterion in Region B is in Equation (41):

Critical Cavity Pressure
During the HDD process, the punch force (Fp) exerted on the workpiece and the equilibrium equation in the vertical direction of the workpiece are expressed as follows in Equations (42) and (43): The maximum tensile stress of the unsupported area (Region B) should not exceed the tensile strength of used material is the principle to calculate the upper critical pressure of critical process window. Therefore, the necking condition occurs around the punch corner and the forming force curve experiences the maximum point (dFp=0). When the cavity pressure is exceedingly higher, the fracture may appear around the die radius area. The specific procedure for calculating the upper critical pressure is referred to in a previous study [7].
The workpiece is separated from the die corner under the effect of cavity pressure to decrease the friction between the blank and the die. According to the critical geometric condition (ρ = r d ) and equilibrium equation, i.e., Equation (42), the lower critical cavity pressure can be presented as Equation (44): Calculation procedures for the critical process window is represented by the flow chart in Figure  4. All the calculation steps were executed utilizing Newton's method in MATLAB R2018b numerical computing language. The initial values of hp and ρ were set as 0. Furthermore, the values of Δh and Δρ were specified to be 0.05 and 0.1 mm, respectively.

Used Materials
In this study, two different types of aluminum alloy sheets, including 2024-O with a thickness of 1.2 mm and 5754-O with a thickness of 0.9 mm, were implemented to investigate the sheet metal HDD process. The chemical compositions of the 2024-O and 5754-O aluminum alloys are revealed in Table 1.
To determine the mechanical properties and the experimental data needed to obtain the anisotropy coefficients, uniaxial tension tests were performed at seven different angles including 0°, 15°, 30°, 45°, 60°, 75°, and 90° from the rolling direction. An MTS-E43.104 electronic universal testing machine was used to conduct the tensile tests with a tensile strain rate of 1.0 × 10 −3 s −1 . To calibrate the Yld2000-2d and BBC2005 yield functions, equi-biaxial tension tests were also performed on a planar biaxial tension testing machine using cruciform specimens [33]. The experiments were performed three times and the average value was taken for development of the yield criteria. The material characteristic coefficients obtained from experiments are presented in Table 2.

Experimental Procedure
The hydrodynamic deep drawing experiments were conducted on double-action specified HDD equipment with a forming capacity of 3500 kN and the maximum blank holder force of 2000 kN. The fluid pressure in die cavity and blank holder force can be regulated in real time according to various designed loading paths controlled by proportional pressure valves, and the maximum cavity pressure can reach 100 MPa. Figure 5. illustrates the die sets of the sheet metal HDD process and the composition of experimental equipment including the test machine, hydraulic control system and computer measurement and control system. The tool geometry is specified in Error! Reference source not found.. During the HDD process, pre-bulging pressure was applied on the blank, and a backward bulging was produced. Subsequently, the blank was drawn into the die chamber with the movement of the punch to a desired depth at a constant velocity of 15 mm/min. Meanwhile, the cavity pressure was adjusted by the proportional valve to press the workpiece onto the punch surface and avoid the defects of rupturing and wrinkling.

Calibration of Anisotropy Coefficient
The anisotropy coefficients of two type Hill'48 yield functions can be calibrated with experimental data. In addition, the eight measured material data were applied to calculate the coefficients of Yld2000-2d yield function. The anisotropy coefficients for BBC2005 function were determined implementing MATLAB R2018b software and optimizing an error function with the Levenberg-Marquardt algorithm. The materials' anisotropy coefficients for yield criteria are summarized in Error! Reference source not found..

Numerical Simulation
In this study, numerical simulations of the sheet metal HDD process were performed utilizing ABAQUS 6.13 software. In addition, for comparison purposes, the Hill'48 yield function, available in the library of ABAQUS software, the Yld2000-2d and BBC2005 yield functions, developed by VUMAT subroutine, were implemented in separate simulations. Due to the orthotropic material properties of aluminum alloy sheets, only a quarter section of the workpiece with the corresponding symmetry boundary conditions was considered. The tooling components were modeled as analytical rigid bodies with the four-node shell. A total of 2560 shell elements with reduced integration (ABAQUS S4R) and five section points through the sheet thickness were applied for the simulation. The contact condition in the tangential direction was governed by the Coulomb friction model. The friction coefficients at the blank/die, blank/blank holder, and blank/punch face interfaces were assumed to be 0.05, 0.08, and 0.1, respectively. Figure 6. demonstrates the dimensions of tooling in the sheet metal HDD process simulation.

Effect of Yield Functions on the Prediction of Material Anisotropy
A comparison between experimental consequences and predicted results obtained from diverse anisotropic yield criteria was performed to evaluate the flexibility of the proposed yield functions. The comparison was concentrated on the following performance aspects including prediction of the yield locus geometry and description of the planar distribution both for uniaxial yield stress and uniaxial coefficient of plastic anisotropy (r value).
In Figure 7, the yield locus projected on the zero-shear stress plane for Hill'48-R, Hill'48-S, Yld2000-2d, and BBC2005 anisotropic yield functions of AA2024 and AA5754 are illustrated together with experimental results. The tensile stresses were normalized with the uniaxial stress at rolling direction. For both the AA2024 and the AA5754 cases, predictions of the yield surface provided by BBC2005 and Yld2000-2d models were very similar, while the performances of the Hill'48 models signify differences to the others, being a little smoother in the region around equi-biaxial tension for both cases. This was because the Hill'48 yield function was less flexible than the others due to the fact that it does not include the equi-biaxial stress in its model.
The normalized yield stresses under uniaxial tension along various loading directions were predicted by diverse yield functions for the AA2024 and AA5754 sheet in Figure 8. and compared with the experimental data points. It was found that the results of the yield functions (i.e., Hill'48-S, Yld2000-2d, and BBC2005 models) matched well with that of the experimental data, while the Hill'48-R model underestimated the uniaxial yield stress for AA2024 and overestimated for the AA5754's in all directions. Moreover, it was apparent that the BBC2005 and Yld2000-2d criteria had better performances in the prediction of the yield stress for both cases. Since the Hill'48-S model was obtained based on the yield stresses, the predicted distribution fit exactly the experimental distribution. Whereas, the predictions based on Hill'48-R criterion were in poorer agreement with the experimental data for both materials.  The uniaxial anisotropy coefficient distributions in the plane of the AA2024 and AA5754 sheets, predicted utilizing different yield criteria, are presented in Figure 9. For the AA2024 case, as demonstrated in Figure 9a, the planar distribution of the r value predicted by Hill'48-S yield criterion was very inaccurate. In contrast, the variation of the r value described by Hill'48-R, Yld2000-2d and BBC2005 models closely followed the experimental results of AA2024 and AA5754. The BBC2005 function provided a slightly better prediction for the r value than the Yld2000-2d function. Concerning AA5754, it was revealed from Figure 9b that the highest r value was at the transverse direction, while the lowest value was at the rolling direction. It can be seen that the predictions for the r value from the BBC2005 yield function matched the experimental results for every orientation very well, with a slight deviation at 75°. The prediction of the r value for Hill'48-S yield function at every 15° from the rolling direction was not as good as the other yield functions because the coefficients for the Hill'48-S model were designed to fit yield stresses at 0°, 45°, and 90°.
The virtual uniaxial tension tests showed that the r values varied in the range of 0.545 to 0.939 and 0.707 to 0.956, while normalized yield stress varied in the range of 0.948 to 1.000 and 0.991 to 1.044 for both materials, i.e., both aluminum alloys revealed strong deformation anisotropy but weak strength anisotropy. Furthermore, the yield stress and r value predicted by the Yld2000-2d and BBC2005 models had better conformity with the experimental results in comparison with the predicted results using the Hill'48 model. It was concluded that the more parameters there were in the Yld2000-2d and BBC2005 yield criteria, the more precisely predicted yield stresses and r values could be obtained. To have a comprehensive evaluation tool of the anisotropic models, a global accuracy index was developed as follow in Equation (45): where ξ Y is the accuracy index associated to the prediction of the yield locus shape in the plane of the principal stresses; ξ S is the accuracy index associated to the prediction of the planar distribution of the uniaxial yield stress; ξ R is the accuracy index associated to the prediction of the planar distribution of the uniaxial coefficient of plastic anisotropy. ξ Y , ξ S , and ξ R are computed using the following formulas:     where di is the squared distance from an experimental point to the yield locus predicted by the yield criterion under testing; n is the total number of the available experimental points; σ θ i exp is the experimental uniaxial yield stress corresponding to the direction defined by the angle θi (measured from the rolling direction); σ θ i t is the predicted uniaxial yield stress associated to the same direction; r θ i exp is the experimental anisotropy coefficient corresponding to the direction defined by the angle θi (measured from the rolling direction); r θ i t is the predicted anisotropy coefficient corresponding to the same direction. The values of the individual accuracy indexes are shown in Figure 10. The best overall performance corresponded to the lowest value of the global yield criterion accuracy index ξ. It was revealed that the values of the anisotropic accuracy index of the Yld2000-2d and BBC2005 yield criteria for both materials were lower than the values for the two Hill'48 yield criteria. This demonstrates the higher reliability of the Yld2000-2d and BBC2005 yield functions rather than the other ones for both aluminum alloys.

The Adaptability of Yield Functions to Predict the Critical Process Window
The critical process window of cavity pressure can be acquired by calculating the lower and upper critical pressures utilizing the analytical procedure, illustrated in Figure 4, together with the material properties and the geometrical dimensions of the workpiece. The relationship between the reasonable cavity pressure and the punch stroke during the HDD process is a significant challenge. On the one hand, the material is bent to the radius of curvature of the punch much faster than the allowed ductility of the material when the pressure is higher than the upper critical pressure curve, which might lead to the rupture of the workpiece around the die shoulder. On the other hand, the material flow resistance is increased, and the sheet is unable to be entirely separated from the die radius if the cavity pressure is smaller than the lower critical pressure. In turn, the sharp thinning around the punch radius is increased which also results in the rupture of the part. Consequently, the safe area between the upper and lower critical profiles can provide reasonable cavity pressure versus punch stroke locus to ensure a flawless part. The critical upper and lower pressure curves predicted based on the Hill'48-R, Hill'48-S, Yld2000-2d, and BBC2005 yield criteria for both materials are illustrated in Figure 11. in which critical pressures in respect to the punch stroke are drawn. It was revealed that the critical pressure loading window predicted by the Hill'48-S, Yld2000-2d, and BBC2005 yield criteria was approximately the same, while the safety area in the critical process window predicted by the Hill48-R yield criteria was greater than the predictions of the other three criteria's results. To explore the adaptability of yield functions to predict the critical process window of cylindrical part, several loading loci were designed for hydroforming experiment for AA2024 and AA5754, respectively, as revealed in Figure 12a,b. Moreover, the maximum wall thickness reduction ratio for the fabricated workpieces are measured to evaluate the critical process window of cavity pressure, as depicted in Figure 13a,b. It can be observed that excessive or insufficient cavity pressure induced the severe thinning around the punch radius region in the results of Locus A and Locus B in Figure  12. and Locus A in Figure 13. For the loading Loci A and B, ruptures appeared because the cavity pressures were too low to entirely separate the blank from the die orifice at the beginning of hydroforming process. On the other hand, for the loading locus J, apparent thickness reduction and fracture are occurred before the blank thoroughly coated with the punch due to the excessive pressure, which is not conducive to the subsequent forming. It is demonstrated that the Locus A in Figure 10 is lower than the lower critical loci of all the yield criteria in Figure 9a. Although the Locus B is in the safe region predicted by Hill'48-R model, it crosses the lower limit of the loading path predicted by the other three yield criteria. Moreover, the Locus J is higher than the upper limit of all yield criteria in Figure 9a. Although the path I is in the safe area predicted by Hill'48-R, it has an intersection with the lower critical pressure locus predicted by the other three yield criteria. It can be seen from Figure  11 that the same situation also occurred in the investigation of AA5754.   Figure 14. compares the experimental critical pressure loading loci with those predicted by the four yield criteria for both aluminum alloys. It can be clearly observed that the critical pressure curves predicted by the Hill'48-R, Hill'48-S, Yld2000-2d and BBC2005 models all had proper conformity with the experimental results, but the Hill'48-R model predicted a larger safe region than the Yld2000-2d and BBC2005 models. To be more specific, the predictions of the Hill'48-S, Yld2000-2d, and BBC2005 criteria were much closer to the experimental safe region for both the AA2024 and AA5754 alloys. Therefore, it can be proved that the Hill48-S, Yld2000-2d, and BBC2005 models were more accurate than the Hill48-R model in predicting the critical process window of the aluminum alloy sheet metal HDD process.

Effect of Yield Functions on the Prediction Accuracy of Deformation Behavior
The comparison between experimental and numerical punch force evolution during the hydroforming process of AA2024 and AA5754 are displayed in Figure 15. The punch force was a global variable and, therefore, was roughly insensitive to the yield criteria with the range between the lower and upper punch force predictions reaching approximately 20% for both materials.
For the 2024 aluminum alloy, the experimental evolution was accurately predicted by the Hill'48-S, Yld2000-2d and BBC2005 models. However, the Hill'48-R yield criterion underestimated the necessary punch force required for further drawing after 10 mm, although the slope was correctly predicted. Concerning the AA5754, the punch force evolution obtained with the BBC2005 yield model, which seem to be closer to the experimental one, was located among the results achieved with the Hill'48 yield criteria, utilizing the two approaches for material parameters identification. Moreover, the Hill'48-S and Yld2000-2d models led to similar results in terms of punch force evolution. The punch force predicted by the Hill'48-R model was higher than applying other three models.
The earing profile is defined by the height of the drawn cup measured along the circumferential direction of the cup using the angle with the rolling direction. A comparison of the experimental earing profile with the computationally predicted earing profiles is featured in Figure 16. It can be stated that the earing profile predicted by the BBC2005 and Yld2000-2d models was in rather good qualitative agreement with experiments for AA2024 and AA5754. Although the two calibrated forms of Hill'48 yield criterion could predict the locations of the peaks and valleys, the earing profiles were overestimated by Hill'48-R, underestimated by Hill'48-S for AA2024, and overestimated by Hill'48-S for AA5754. The earing profile was poorly predicted by the Hill'48 yield model. On the other hand, the Yld2000-2d function predicted the earing profile but not as precisely as the earing predicted by the BBC2005 yield model which is evident from Figure 16.  The capability and accuracy of the selected yield functions in the prediction of thickness variations in the HDD process were investigated. The wall thickness variations obtained by simulation at the rolling and the transverse directions were measured and compared with the corresponding experimental results in three different regions of AA2024 and AA5754 as plotted in Figure 17 and Figure 18, respectively. Based on the simulation results for AA2024, shown in Figure  17, predictions of thickness variation at the rolling and the transverse directions were nearly the same for all yield criteria and approximately coincided with the experimental results. Moreover, it was observed that the simulation results based on the Yld2000-2d and BBC2005 yield functions were in better agreement with the experimental results for both aluminum alloys than those based on the other three yield functions. Concerning the consequences for the AA5754 alloy in Figure 18, the wall thickness in Region C was overestimated by the Hill'48-R function in both directions.

Conclusions
In this research, three anisotropic yield criteria including Hill'48, Yld2000-2d, and BBC2005 were implemented to describe the anisotropic behavior in the hydrodynamic deep drawing process of aluminum alloys. Theoretical models based on different yield criteria were developed to determine the critical pressure at instability points and the critical process window in the HDD process of cylindrical cups. Furthermore, the accuracy of the predicted materials' anisotropic coefficients was evaluated by the experimental results. The effect of the anisotropic yield criterion on the prediction of the critical process window and material deformation behavior in the HDD process was investigated. The conclusions have been drawn as follows: (1) The Yld2000-2d and BBC2005 yield formulas performed better in predicting the anisotropic coefficients, yield locus, and directional flow stresses for the AA2024 and AA5754 alloys. However, the Hill'48 yield function provided relatively poor capability in the prediction of the r value and uniaxial yield stress due to the fact that it underestimated the biaxial yield stress for aluminum alloys.
(2) Although all of the four anisotropic yield models had proper conformity with the experimental results of the critical process window, the predictions using the Yld2000-2d and BBC2005 yield criteria were much closer to the experimental safe region for both the AA2024 and AA5754 alloys.
(3) The FEM results based on the Yld2000-2d and BBC2005 yield functions were in closer agreement with experimental data in terms of forming load, cup earing height, and wall thickness distribution for both AA2024 and AA5754 alloys. The flexible yield functions, like BBC2005 and Yld2000-2d, with more anisotropic parameters provided better predictability in modeling the critical process window and anisotropic deformation in the HDD process of aluminum alloys.