Study of Anisotropic Behavior in Sheet Metal Forming

Since sheet metal exhibits significant anisotropy in processing and forming, which has a significant impact on its performance during processing, forming, and use, we explore the anisotropic behavior of materials in the forming process of sheet metal. The ability of the Yld2000-2d criterion to describe anisotropic behavior is analyzed, and its accuracy for characterization of the anisotropic behavior of metal plates is improved, based on which anisotropic behavior is predicted in three-dimensional space. Theoretical and experimental results on the anisotropy of sheet metal are compared, and two materials, 5754O aluminum alloy and DP980 steel plate, are tested and analyzed, and the anisotropic behaviors, such as three-point bending and cylindrical deep-drawing, are well predicted.


Introduction
Sheet metals that have undergone remanufacturing and multiple rolling will present fibrous tissue or preferential orientation caused by crystallization, which will lead to anisotropy, meaning the material exhibits different mechanical properties in different directions.This anisotropic behavior has obvious influences on deformation during the forming process, such as earring and fracture on straight walls in deep drawing [1][2][3].In general, the anisotropy of machined sheet metal and the resulting stresses and loose areas have a significant effect on its performance during forming and use, as in the following examples.(1) Mechanical property differences: anisotropy will lead to differences in the mechanical properties of the metal plate in different directions, such as strength, hardness, ductility, and so on.(2) Increased molding difficulty: anisotropy causes the deformation capacity in different directions to be different, which may lead to deformation in some directions in the molding process being difficult to achieve or requiring greater molding force.(3) Stress concentration: anisotropy makes it possible for stress concentration to occur when the metal plate is subjected to force, increasing the risk of local instability and fatigue damage.If the anisotropy is ignored, it may result in poorly formed or cracked aluminum sheets.If the anisotropic behavior during the sheet metal forming process can be predicted accurately, the forming quality can be greatly improved by utilizing the anisotropic characteristic effectively [4,5].
In order to accurately describe the anisotropic behavior of sheet metals, researchers all over the world have proposed many functions, such as Hill series yield functions [6][7][8], Hosford yield function [9], Barlat series yield function [10][11][12], and so on.Every yield function has its advantages and disadvantages.For example, the Hill48 yield function that is widely used in engineering applications is a simple function, the parameters of which are easy to solve, but its prediction of the anisotropic behavior of certain metals is not accurate enough [13].In 1972, Hosford [14] proposed a yield function that can accurately describe the behaviors of cubic center and face center materials, but it does not include shear stress, which is inconvenient for engineering applications.In 2003, Barlat et al. [15] proposed a plane stress anisotropic yield function, which has high accuracy for anisotropic materials, especially for aluminum alloys [16,17].The Yld2000-2d yield function has been widely used in the sheet metal forming field.At present, some special software already includes this yield function [18], while general software like ABAQUS6.4 needs to develop user subroutines to realize the application of the Yld2000-2d anisotropic yield function [19].There are also some scholars who have developed anisotropic models for the anisotropy of a certain type of material; for example, Li et al. [20] developed a model that captures the tensile asymmetry of magnesium alloys, and Ji et al. [21] developed an anisotropic model that captures the anisotropy of aluminum alloys with respect to the strain rate.A generalized evolutionary plasticity model that takes into account thermal effects on the evolution of flow behavior was proposed by Lian et al. [22] to describe the temperature dependence of the anisotropic plastic flow behavior of the studied materials.
To accurately predict the plastic deformation behavior of materials in finite element simulations, it is necessary to use more accurate and convenient models [23][24][25].In many cases, sheet metal experiences a three-dimensional stress state during the forming process, and the stress in the thickness direction will obviously affects the metal forming properties.Existing studies show that the stress in the thickness direction has obvious influences on sheet metal forming processes [26,27].Existing studies also show that when compressive stress exists in the thickness direction, the forming limit curves will rise, and when tensile stress exists in the thickness direction, the forming limit curves will become lower [28,29].In the sheet metal forming cases where the thickness stress is obvious, if the plane stress yield function is still used, the simulation accuracy might be decreased.Since some sheet metal forming simulations needs to consider the stress in the thickness direction, solid elements are more commonly adopted, and a 3D stress yield function should then be used.
Many three-dimensional stress anisotropic yield functions exist, among which the Hill48 and Yld2004-18p yield functions are widely accepted for engineering applications.Mu et al. [30] characterized the deformation behavior of DC06 steel plates more accurately by improving Hill48, and Rong et al. [31] used yld2004 to describe the thermal anisotropic behavior of AA7075 with high similarity.The Hill48 yield function, which has less flexibility, has low accuracy sometimes.It will lead to abnormal phenomena when describing aluminum alloy behaviors.The Yld2004-18p yield function has very high accuracy due to its high flexibility, which can describe the anisotropy of sheet metals very well.For example, the Yld2004 yield function can predict six and eight ears in deep drawing, while many other yield functions cannot [32].The Yld2004 3D stress yield function has high accuracy, but the formula is complicated, and many experimental results are needed to determine the parameters' values.
In many cases, the anisotropy of the sheet metal is relatively simple (for example, considering the anisotropic behaviors along three or four directions is enough), but 3D stress should be considered (i.e., the normal stress of sheet metal should not be neglected).Therefore, in order to better characterize the anisotropic behavior of sheet metal in practical machining applications, it is necessary to establish a yield function suitable for threedimensional stress states so that a more convenient and accurate model can be constructed.
In this study, the anisotropic mechanical properties of two commonly used metal sheets (5754O and DP980) were analyzed.The Yld2000-2d yield function proposed by Barlat et al. [15] is extended into three-dimensional stress space for materials that are hydrostatic-stress-independent; this is called the Yld2000-3d yield function, and with it, the anisotropic properties of aluminum alloy plates under three-dimensional stress states are described more precisely.A finite element simulation of a three-point bending model of 5754O aluminum alloy sheet was performed by inserting the model into ABAQUS through the UMAT (VUMAT) subroutine.The accuracy of the new model in capturing the anisotropy of the aluminum alloy sheet is verified by comparing the simulation results with experimental results.In addition, deep drawing tests and finite element simulations of DP980 plates were conducted to further validate the generalizability of the method for characterizing the anisotropy of metal plates.

Yld2000-2d Yield Criterion Anisotropy Characterization Capability
The Yld2000-2d yield function [15] is expressed as: where In Equation ( 1), ϕ is the sum of the two isotropic functions.In Equation ( 2), X ′ i and X ′′ j are the principal values of the matrices X ' and X " : 2 ) 2 ) 2 ) 2 ) (3) The elements of X ' and X " are converted from Cauchy stress: where where m is the material parameter, and α 1 − α 8 are the eight anisotropic parameters of the anisotropic yield function, which can be calculated with the material properties.When the eight parameters are equal to 1, the function will return to the isotropic function.

Prediction of Anisotropic Behavior of Materials in Three-Dimensional Stress State
As stated above, the stress components included in the Yld2000-2d yield function are σ xx , σ yy , and σ xy in the X-Y plane.In many cases, however, the normal stress and shear stress in the thickness direction cannot be neglected and have a significant effect on the sheet metal forming process.In this study, the Yld2000-2d yield function is modified, and the normal stress along the thickness direction is included firstly.That is, the Yld2000-2d yield function for the 2D case is extended into three-dimensional stress space.
As shown in Figure 1a, the material point is subjected to different stresses in different directions and is in equilibrium.Convert the 3D stress state of Figure 1a by subtracting a hydrostatic stress equal to the algebraic value of thickness direction stress σ z (adding a hydrostatic stress, the algebraic value of which is −σ z ).The converted stress state is shown in Figure 1b.According to the continuity hypothesis, for the hydrostatic-stress-independent materials, the hydrostatic stress will only cause elastic volume changes, and will not affect the plastic deformation rules.Therefore, the plastic deformation caused by the two stress states shown in Figure 1 will be completely the same.
hydrostatic stress equal to the algebraic value of thickness direction stress z  (adding a hydrostatic stress, the algebraic value of which is z - ).The converted stress state is shown in Figure 1b.According to the continuity hypothesis, for the hydrostatic-stressindependent materials, the hydrostatic stress will only cause elastic volume changes, and will not affect the plastic deformation rules.Therefore, the plastic deformation caused by the two stress states shown in Figure 1 will be completely the same.As for the normal stress, the above conversion expression from the three-dimensional to two-dimensional stress space is: Besides the three normal stresses Thus, the Yld2000-2d yield function is extended into three-dimensional stress space and is named "Yld2000-3d" in this study, the expression of which is where As for the normal stress, the above conversion expression from the three-dimensional to two-dimensional stress space is: Besides the three normal stresses σ x , σ y , and σ z , the three-dimensional stress components also include the three shear stress components σ xy , σ yz , and σ zx .The shear stresses σ yz and σ zx are handled with the method that was used for σ xy in the Yld2000-2d yield function.
Thus, the Yld2000-2d yield function is extended into three-dimensional stress space and is named "Yld2000-3d" in this study, the expression of which is where where The components of X ' and X " are converted from the Cauchy stress as follows: From Equations ( 13) and ( 14), it can be seen that in the yield function expression, the coefficients of the three shear stress components σ xy , σ yz , and σ zx are L ′ 44 , L ′ 55 , L ′ 66 and L ′′ 44 , L ′′ 55 , L ′′ 66 , which represent the shear anisotropic properties and can be expressed as where α 1 -α 12 are anisotropic material parameters.
Because the material property data in the XZ and YZ planes are difficult to obtain in basic sheet metal tests, the shear isotropy is assumed in this study, namely the anisotropic coefficients for σ xy , σ yz , σ zx are assumed to be the same.Liu and Raemy used similar methods in their yield functions [33,34].Therefore, in Equations ( 13) and ( 14), we have In other words, we have α 7 = α 9 = α 10 and α 8 = α 11 = α 12 .Surely, the shear anisotropy can be represented if corresponding material properties in the XZ and YZ planes can be obtained.
It is shown that the developed yield function still used the linear transformation method.Yoshida et al. [35] perform the same conversions on a six-order polynomial, and the convexity of a certain cross-section of the three-dimensional yield function is considered.In this study, the three-dimensional stress yield function developed can ensure strict convexity.It is easy to prove that its convexity is in consistency with the original Yld2000-2d yield function.Besides the uniaxial and the biaxial tensile mechanical properties, the thickness stress components σ zz , σ yz , σ zx are also included in the Yld2000-3d yield function.On the other hand, from Equations ( 13) and ( 14), it is easy to find that the Yld2000-3d yield function only considers the value of the thickness stress and does not consider the anisotropy of the thickness stress.Cazacu [36] found that in the sheet-forming processes, the existence of the thickness stress should be considered, but the anisotropy of the thickness direction stress does not need to be considered.
The determination of the parameters is an important step for the application of the developed Yld2000-3d yield function.Comparing it with Yld2000-2d yield function, no more parameters are included in the Yld2000-3d function, and the determining method in this study is exactly the same as that of the original Yld2000-2d yield function [16] (only some anisotropic properties under planes stress state are needed), so it will not be presented in detail.As a result, the anisotropy in the metal sheet-forming process can be captured more conveniently and accurately for practical engineering use.

Experimental and Material Data
Biaxial tensile testing of cruciform specimens was performed on an established biaxial tensile testing system.The test system is modified and improved from the biaxial loading test machine shown in Figure 2, which can realize the synergistic motion of four axes in two directions and can carry out biaxial tensile tests under different proportional loading paths.The characteristics of the test system are as follows: tensile speed 0.6-6 mm/min, force transducer accuracy ≤2%, chuck maximum stroke 400 mm, working environment temperature −10-50 • C. It can be seen that the bi-directional tensile test system meets the needs of the test in terms of requirements for precision, reliability, and accuracy.
Biaxial tensile testing of cruciform specimens was performed on an established biaxial tensile testing system.The test system is modified and improved from the biaxial loading test machine shown in Figure 2, which can realize the synergistic motion of four axes in two directions and can carry out biaxial tensile tests under different proportional loading paths.The characteristics of the test system are as follows: tensile speed 0.6-6 mm/min, force transducer accuracy ≤2%, chuck maximum stroke 400 mm, working environment temperature −10-50 °C.It can be seen that the bi-directional tensile test system meets the needs of the test in terms of requirements for precision, reliability, and accuracy.Biaxial tensile testing of cruciform specimens is usually performed using both loadcontrolled and strain-controlled methods.The study of yielding behavior is generally carried out in stress space, and the yield criterion is usually expressed in terms of stresses; hence, the load control approach is used in this paper.For the cruciform specimen, a 5754O (t = 1 mm) aluminum alloy plate with r < 1 was selected.The test specimens were Biaxial tensile testing of cruciform specimens is usually performed using both loadcontrolled and strain-controlled methods.The study of yielding behavior is generally carried out in stress space, and the yield criterion is usually expressed in terms of stresses; hence, the load control approach is used in this paper.For the cruciform specimen, a 5754O (t = 1 mm) aluminum alloy plate with r < 1 was selected.The test specimens were prepared by laser cutting, and the specific dimensions were obtained from Wang et al.'s paper [19].
Before the test, the extensometer for strain measurement and the force transducer for load measurement were calibrated, and the calibration errors were less than 1%.During the test, the main control direction tensile speed is 2 mm/min, the initial measurement area of the strain in the center area is 50 mm × 50 mm, and the strain value during deformation is measured by an extensometer installed on both surfaces of the specimen at a distance of 50 mm, respectively.The measurement range of the extensometer is 10 mm, and the pins of the extensometer are fixed to the specimen by clamping during the measurement.
The ratios of the rolling direction to the load perpendicular to the rolling direction used in the test were 4:0, 4:1, 4:2, 4:3, 4:4, 3:4, 2:4, 1:4, and 0:4, respectively.The test is stopped when the load curve decreases or when rupture of the specimen is observed, and the data collected in the test are the load on the cross arm of the specimen and the deformation in the center area.Each set of data was repeated three times, and its average value was finally solved to obtain the mechanical properties of 5754O, as shown in Tables 1 and 2.

Three-Dimensional Stress State Material Anisotropy
The Yld2000-2d yield function can accurately predict the mechanical behaviors of 5754O aluminum alloy [37,38] from existing research results.In this study, considering the normal stress and shear stress in the sheet thickness direction, 5754O aluminum alloy is adopted as an example to verify the accuracy of the developed Yld2000-3d yield function.The uniaxial and biaxial tensile properties of 5754O aluminum alloy are shown in Tables 1 and 2, where β represents the plastic strain direction under plane stress condition, and the r-value is the coefficient of anisotropy in the thickness direction.They have the following definitions: r = dε y /dε z (18) By the principal of volume invariance, dε x + dε y + dε z = 0, and the r-value and β can be transformed into each other computationally.
According to the parameter determination method of Yld2000-3d, which is the same as that of the Yld2000-2d yield function, the material constants of 5754O aluminum alloy sheets are obtained as: The yield surface of 5754O aluminum sheet represented by the Yld2000-3d yield function is shown in Figure 3. Figure 4a,b show the uniaxial tensile anisotropy of a 5754O aluminum sheet predicted by the Yld2000-3d yield function.From the theoretical and experimental results, it is shown that the Yld2000-3d (Yld2000-2d) yield function can accurately characterize the uniaxial and biaxial tensile properties of a 5754O aluminum alloy sheet.The tricomponent yield surfaces in the σ x − σ y , σ y − σ z , and σ z − σ x stress space based on the Yld2000-3d yield function are shown in Figures 5a, 6a and 7a, respectively.The results showed that the tricomponent yield surfaces of 5754O aluminum sheet in different stress spaces are very different from each other due to anisotropy. Figure 5b is the plastic strain direction in σ x − σ y stress space, where φ = arctan(σ y /σ x ) is the loading angle.It is shown that the plastic strain direction for different loading angle agrees very well with the experimental results, which further verifies the anisotropic description ability of the Yld2000-3d (Yld2000-2d) yield function.Figures 6b and 7b are the plastic strain direction in σ y − σ z and σ z − σ x stress space predicted by the Yld2000-3d yield function, respectively.the loading angle.It is shown that the plastic strain direction for different loading angle agrees very well with the experimental results, which further verifies the anisotropic description ability of the Yld2000-3d (Yld2000-2d) yield function.Figures 6b and 7b

5754O Aluminum Alloy Plate Anisotropy Prediction Verification
In this study, the complete implicit back-Euler integration mapping algorithm is adopted to implement the developed Yld2000-3d yield function in ABAQUS commercial software since it has been used successfully for the implementation of Yld2000-2d in ABAQUS [39].The specific stress transfer process in the subroutine is derived in detail in In the FEM simulations, the two-dimensional S4R shell elements and three-dimensional C3D8R solid elements are used based on the Yld2000-2d yield function and the Yld2000-3d yield function, respectively.The dies are established with discrete rigid elements, limiting the movement of the three-point-bending thin plate specimen in the xdirection (Ux) and the rotation in both directions of the y-axis and z-axis (URy, URz).The FEM model is shown in Figures 9 and 10.In the FEM simulations, the two-dimensional S4R shell elements and three-dimensional C3D8R solid elements are used based on the Yld2000-2d yield function and the Yld2000-3d yield function, respectively.The dies are established with discrete rigid elements, limiting the movement of the three-point-bending thin plate specimen in the x-direction (Ux) and the rotation in both directions of the y-axis and z-axis (URy, URz).The FEM model is shown in Figures 9 and 10.As shown in Figure 11, the spring-back result of the solid elements model (tw four layers) with the Yld2000-3d yield function agrees better with the experimental res than that of the shell elements (five or seven integral points) with the Yld2000-2d y function.This is because when using the solid elements in the simulations, the Yld2 3d three-dimensional stress yield function comprehensively considers the normal st and shear stress of the sheet thickness direction.The accuracy of the new model for As shown in Figure 11, the spring-back result of the solid elements model (tw four layers) with the Yld2000-3d yield function agrees better with the experimental re than that of the shell elements (five or seven integral points) with the Yld2000-2d function.This is because when using the solid elements in the simulations, the Yld2 3d three-dimensional stress yield function comprehensively considers the normal s and shear stress of the sheet thickness direction.The accuracy of the new model for As shown in Figure 11, the spring-back result of the solid elements model (two or four layers) with the Yld2000-3d yield function agrees better with the experimental results than that of the shell elements (five or seven integral points) with the Yld2000-2d yield function.This is because when using the solid elements in the simulations, the Yld2000-3d three-dimensional stress yield function comprehensively considers the normal stress and shear stress of the sheet thickness direction.The accuracy of the new model for capturing the anisotropy of the 5754O aluminum alloy in the 3D stress space is further verified.There is no obvious difference between the two simulated results based on four layers and two layers of solid element, or the two simulated results based on five and seven integral points shell elements.As shown in Figure 11, the spring-back result of the solid elements model (two or four layers) with the Yld2000-3d yield function agrees better with the experimental results than that of the shell elements (five or seven integral points) with the Yld2000-2d yield function.This is because when using the solid elements in the simulations, the Yld2000-3d three-dimensional stress yield function comprehensively considers the normal stress and shear stress of the sheet thickness direction.The accuracy of the new model for capturing the anisotropy of the 5754O aluminum alloy in the 3D stress space is further verified.There is no obvious difference between the two simulated results based on four layers and two layers of solid element, or the two simulated results based on five and seven integral points shell elements.

DP980 Steel Plate Anisotropy Prediction Verification
In this section, the ability of the extended 3D model (Yld2000-3d) to capture the anisotropy of the rest of the plates and processing is further tested by numerical simulation of cylindrical deep drawing tests on DP980 plates, in the hope that the study of 5754O aluminum alloy plates will provide a new idea for the description of the anisotropic behavior of metal sheet materials.The material properties and the cylindrical deep drawing test of DP980 steel sheets given by Cai et al. [40] were adopted.The dimensions of the tools for the deep drawing tests are specified as follows: the punch diameter is 50.00 mm, the punch profile radius is 5.00 mm, the opening diameter of the die is 53.64 mm, and the die profile radius is 13.00 mm.The thickness of the DP980 steel sheet used in the test is 1.2 mm.The blank holding forces (BHF) used in the deep drawing tests are 15 kN equally, and the Coulomb coefficient of friction is set to be 0.2 in the simulations.The stamping processing speed is 20 mm/min.Mechanical properties of DP980 steel alloy sheet are shown in Table 3.As before, in the FEM simulations, the two-dimensional S4R shell elements and threedimensional C3D8R solid elements are used.The dies are established with discrete rigid elements.The simulated profiles of the deep drawing test are shown in Figure 12.The comparison between the experimental and numerical predicted earring profiles is presented in Figure 13.The comparison between the experimental and numerical predicted earring profiles is presented in Figure 13.As shown in Figure 13, the simulated results based on theYld2000-3d yield function and three-dimensional C3D8R solid elements are in better agreement with the experimental results than other yield functions.Although the Hill's quadratic and Yld2000-2d yield function could predict the locations of the peaks and valleys, the amplitudes of the earring profiles are underestimated.Since no anisotropy is considered, the Mises yield function could not predict the earring phenomenon.The Yld2000-3d yield function could provide precise predictions of the peaks' heights, and relatively reasonable estimates of the valleys' heights, indicating that the model generalizes for the remaining plates in characterizing the anisotropy.

Conclusions
The uniaxial tensile test, biaxial tensile test, and three-point bending models of the aluminum alloy were sequentially conducted and compared with finite element simulations to verify the accuracy of the model extended to a three-dimensional stress state (Yld2000-3d), which captures the anisotropic mechanical properties of the 5754O aluminum alloy with greater precision.The model was also used to test the rest of the plates and different loading methods, further extending the generalization of the model to different materials and deformation processes.We will continue to test the model for more complex deformation modes or new materials in the following work.
(1) Considering that some actual sheet metal forming processes are under a three-dimensional stress condition, the Yld2000-2d anisotropic yield function is extended into three-dimensional space for hydrostatic-stress-independent materials; this extended function is named the Yld2000-3d yield function.Taking the 5754O aluminum sheet as an example, the yield surface and the plastic strain ratio diagram are drawn based on the Yld2000-3d three-dimensional stress yield function developed in this study.As shown in Figure 13, the simulated results based on theYld2000-3d yield function and three-dimensional C3D8R solid elements are in better agreement with the experimental results than other yield functions.Although the Hill's quadratic and Yld2000-2d yield function could predict the locations of the peaks and valleys, the amplitudes of the earring profiles are underestimated.Since no anisotropy is considered, the Mises yield function could not predict the earring phenomenon.The Yld2000-3d yield function could provide precise predictions of the peaks' heights, and relatively reasonable estimates of the valleys' heights, indicating that the model generalizes for the remaining plates in characterizing the anisotropy.

Conclusions
The uniaxial tensile test, biaxial tensile test, and three-point bending models of the aluminum alloy were sequentially conducted and compared with finite element simulations to verify the accuracy of the model extended to a three-dimensional stress state (Yld2000-3d), which captures the anisotropic mechanical properties of the 5754O aluminum alloy with greater precision.The model was also used to test the rest of the plates and different loading methods, further extending the generalization of the model to different materials and deformation processes.We will continue to test the model for more complex deformation modes or new materials in the following work.
(1) Considering that some actual sheet metal forming processes are under a threedimensional stress condition, the Yld2000-2d anisotropic yield function is extended into three-dimensional space for hydrostatic-stress-independent materials; this extended function is named the Yld2000-3d yield function.Taking the 5754O aluminum sheet as an example, the yield surface and the plastic strain ratio diagram are drawn based on the Yld2000-3d three-dimensional stress yield function developed in this study.The results show that the developed Yld2000-3d anisotropic yield function can accurately describe the uniaxial and biaxial mechanical properties of 5754O aluminum alloy sheets, and it can predict the mechanical behaviors in YZ and XZ planes.(2) The Yld2000-3d function developed in this study is implemented in ABAQUS software by using the UMAT (VUMAT) user subroutine.The three-points bending test of a 5754O aluminum alloy sheet is carried out and the corresponding FEM simulations are performed based on the Yld2000-2d and Yld2000-3d yield functions.It is shown that the simulated results of three-points bending based on the Yld2000-3d yield function with a solid element agree better with the experimental results than those based on the Yld2000-2d yield function with a shell element.(3) Further testing of other sheet metal forming processes was carried out by simulating deep drawing tests of DP980 steel plates using the Yld2000-3d, Yld2000-2d, Hill48, and Mises yield functions.The comparison between the simulated and experimental results of the deep drawing test shows that the Yld2000-3d yield function has the highest accuracy among them in predicting the earring profile.(4) Besides the advantages of the Yld2000-2d yield function, the stress components related to the thickness direction are also considered, and the Yld2000-3d yield function has high accuracy in describing the anisotropic behavior under a 3D stress state.Compared to the Yld2000-2d yield function, there is no additional parameter in the Yld2000-3d yield function, which makes it convenient in engineering applications.modulus; ρ is the hardening curve, and the subscript represents the analysis step; and σ and ε are the three-dimensional stress and strain, respectively, where During the stress updating process, with regard to a certain three-dimensional strain increment ∆ε, assume the updated stress as the elastic stress firstly, and we have: where the superscript t represents an attempting condition, and the plastic variable is the same as the previous value.
If the given elastic tolerance is satisfied with regard to the following condition, then the step "n + 1" is the elastic condition: If the above-mentioned condition is not satisfied, then this step is the elastic-plastic stage.Set the above attempting condition as the starting value of the plastic adjusting program.The new attempting stress should be on the new yield surface, that is: The consistency requirements need to be satisfied on the new yield surface, that is: Equation (A12) is the nonlinear equation to solve the equivalent plastic strain increment ∆ε p n+1 .The prediction-adjusting system based on Newton-Raphson method is adopted to solve ∆ε p n+1 .Each three-dimensional tensor in this system is updated with each update of ∆ε In this study, the solving method of the tangent modulus that is used in the literature is adopted, and the tangent modulus Cep that is suitable for the three-dimensional stress yield function is deduced:

Appendix B
The three-dimensional "Yld2000-3d" yield function is expressed as follows: The elements in X ' and X " are converted by Cauchy stress: The parameters are determined by the following method: The first derivative of the three-dimensional "Yld2000-3d" yield function is as follows:

Appendix C
Because the ∂ϕ ∂σ expression in B complicated, the ∂σ∂σ expression is lengthier.Therefore, the second derivate expression of each stress component is not given in detail.The overall expression of ∂ 2 ϕ ∂σ∂σ is given here.In order to express it conveniently, the improved yield function is simply converted as follows: 2 ) ) 2 ) The second derivate of the "Yld2000-3d" three-dimensional stress yield function is as follows:

Figure 1 .
Figure 1.Stress state analysis of the one material point: (a) three-dimensional stress; (b) equivalent two-dimensional stress.


, the three-dimensional stress components also include the three shear stress components , the method that was used for xy  in the Yld2000-2d yield function.

Figure 1 .
Figure 1.Stress state analysis of the one material point: (a) three-dimensional stress; (b) equivalent two-dimensional stress.

Figure 3 .
Figure 3.The normalized plastic work contours in

Figure 3 .
Figure 3.The normalized plastic work contours in σ x − σ y plane of 5754O aluminum alloy sheet based on different yield functions.Materials 2023, 16, x FOR PEER REVIEW 9 of 19

Figure 4 .Figure 5 .
Figure 4.The uniaxial tensile properties of 5754O aluminum alloy sheet predicted based on different yield functions: (a) stress prediction; (b) r value prediction.

Figure 4 .
Figure 4.The uniaxial tensile properties of 5754O aluminum alloy sheet predicted based on different yield functions: (a) stress prediction; (b) r value prediction.

Figure 4 .Figure 5 .
Figure 4.The uniaxial tensile properties of 5754O aluminum alloy sheet predicted based on different yield functions: (a) stress prediction; (b) r value prediction.

Figure 9 .Figure 10 .
Figure 9. FEM model established using S4R shell elements with five and seven integration poi

Figure 9 .Figure 9 .Figure 10 .
Figure 9. FEM model established using S4R shell elements with five and seven integration points.

Figure 10 .
Figure 10.FEM model established using C3D8R solid elements with two and four layers of elements: (a) two layers of elements; (b) four layers of elements.

Figure 10 .
Figure 10.FEM model established using C3D8R solid elements with two and four layers of elements: (a) two layers of elements; (b) four layers of elements.

Figure 11 .
Figure 11.The experimental and simulated three-point bending spring-back results.Figure 11.The experimental and simulated three-point bending spring-back results.

Figure 11 .
Figure 11.The experimental and simulated three-point bending spring-back results.Figure 11.The experimental and simulated three-point bending spring-back results.
before, in the FEM simulations, the two-dimensional S4R shell elements and three-dimensional C3D8R solid elements are used.The dies are established with discrete rigid elements.The simulated profiles of the deep drawing test are shown in Figure12.

Materials 2023 , 19 Figure 13 .
Figure 13.Experimental and simulation results of DP980 cylindrical deep drawing earring profile with different yield functions.

Figure 13 .
Figure 13.Experimental and simulation results of DP980 cylindrical deep drawing earring profile with different yield functions.

Table 1 .
Mechanical properties of 5754O aluminum alloy sheets.

Table 2 .
Biaxial mechanical properties of 5754O aluminum alloy sheets.

Table 3 .
Mechanical properties of DP980 steel alloy sheets.