Research on Forming Limit Stress Diagram of Advanced High Strength Dual-Phase Steel Sheets

The Forming Limit Stress Diagram (FLSD) can accurately describe the forming process of high-strength steel. However, obtaining FLSD is relatively difficult. In order to predict fracture in advanced high-strength dual-phase (DP) steels, limit maximum and limit minimum principal strains of sheet were obtained through multiple sets of test and simulation. Two material parameters, strength coefficient K and hardening exponent n are introduced into the FLSD function which is established by the strain-stress transformation function. The function shows that the k-value determines the value of the maximum principal stress, while the n-value affects the curvature of the curve. Verification of correctness by testing and simulation to within 10% accuracy. This paper explores a new approach to FLSD research based on material properties, which can expand the application scope of FLSD.


Introduction
As excellent automotive lightweight material, Advanced high-strength steels (AHSS), especially advanced high-strength dual-phase (DP) steels are widely used in body structures with high-strength requirements. Usually, as with other metallic materials, the formability of AHSS sheet is assessed by forming limit stress diagram (FLD). FLD is represented by a curve in the principal strain coordinate system, which divides the state space of the two principal strains (all stress states from uniaxial tension to biaxial tension) in the sheet into two regions: a safety zone and a fracture zone.
However, although widely used, the current FLD is unsatisfactory to describe AHSS forming. Due to the local plastic deformation of the edge, the property of the AHSS sheet processed by the shear process has changed significantly. Their uneven hardening phenomenon of the shear edge area is more prominent, which greatly improves the probability of edge cracking of the material in the subsequent forming process. Wang K et al. [1] found that when the pre-strain is larger than the the limit strain which can be obtained under linear path, the FLD will become very unstable due to the path influence. Another typical problem that occurs during stamping of AHSS is the shear fracture [2,3]. Fu Q T et al. [4] investigated the strength coefficient K, hardening index n, and stress triaxiality associated with this shear fracture.
In order to avoid the springback in the forming process and improve the size precision of forming, the smaller bending radius is employed in designing of the mold, meanwhile increase the blank holder force during the forming as much as possible without creating fractures. With the decrease of the bending angle of the mold, part of the material close to the bending angle of the mold is easy to produce a shear fracture. There is no significant necking and almost no thinning at the shear fracture. Both of these phenomena are hard to predict in FLD. Numerous studies have found that FLD has disadvantages of strainpath dependent, which makes this method become ineffective in the analysis of complex The M-K theory considers the thickness of the sheet metal to be inhomogeneous and that this inhomogeneity is amplified during stretching until fracture. The theory uses a model of a sheet with a groove to carry out the study of a sheet subjected to biaxial stretching. Where zone A is the uniform deformation zone, zone B is the groove zone, and 1 σ , 2 σ are the maximum and minimum principal stresses in Figure 1.
This theory was based on some assumptions: (1) Plastic deformation of the sheet due to biaxial tensile stress, but the volume remains unchanged.
where 1 dε , 2 dε , 3 dε are the plastic strains of the sheet in three directions.
(2) The second principal strain increment in zone A and zone B are equal.
where 2A dε , 2B dε are second principal strain increment in zone A and zone B.
(3) The forces in zone A and zone B in the 1st principal direction are in balance where 1A σ , 1B σ are the first principal stresses of zone A and zone B; A t , B t are the thickness of zone A and zone B.
(4) When simple loading is applied to Zone A, each principal stress and principal strain increase in proportion, and the ratio is constant: where d x σ , d x ε are the incremental principal stress and incremental principal stain, x σ , x ε are the principal stress and principal stain.
Combining M-K theory with material constitutive model and corresponding yield criterion, the maximum and minimum limit principal strains 1 ε and 2 ε can be obtained by iterative calculation, and the sheet FLD can be obtained. The M-K theory considers the thickness of the sheet metal to be inhomogeneous and that this inhomogeneity is amplified during stretching until fracture. The theory uses a model of a sheet with a groove to carry out the study of a sheet subjected to biaxial stretching. Where zone A is the uniform deformation zone, zone B is the groove zone, and σ 1 , σ 2 are the maximum and minimum principal stresses in Figure 1.
This theory was based on some assumptions: (1) Plastic deformation of the sheet due to biaxial tensile stress, but the volume remains unchanged.
where dε 1 , dε 2 , dε 3 are the plastic strains of the sheet in three directions.
(2) The second principal strain increment in zone A and zone B are equal.
where dε 2A , dε 2B are second principal strain increment in zone A and zone B.
(3) The forces in zone A and zone B in the 1st principal direction are in balance where σ 1A , σ 1B are the first principal stresses of zone A and zone B; t A , t B are the thickness of zone A and zone B.
(4) When simple loading is applied to Zone A, each principal stress and principal strain increase in proportion, and the ratio is constant: where dσ x , dε x are the incremental principal stress and incremental principal stain, σ x , ε x are the principal stress and principal stain. Combining M-K theory with material constitutive model and corresponding yield criterion, the maximum and minimum limit principal strains ε 1 and ε 2 can be obtained by iterative calculation, and the sheet FLD can be obtained.

Obtaining Stress State
The limit stress value of metal sheet is difficult to measure in the production. Based on the relationship between stress and strain, the transition from FLD to FLSD can be realized. Panich [15] compared the FLSD of 980 steel calculated from different yield criteria and found that the FLSD data obtained differed significantly according to the different yield criteria. Although the FLSD derived from different theories is not unique, the established FLSD failure criterion will be close to the original FLD results, as long as the same plasticity theory and yield criteria as the FLD calculation model are used in the numerical simulation.
The stress-strain equations in this paper are in the plane-stress (σ 3 = 0) conditions, σ 1 , σ 2 are the maximum and minimum principal stresses.
According to the literature [16,17], the transformation relationship from strain state to stress state is shown in the Equations (5) and (6): where ε 1i and ε 2i are the pre strains in two directions, ε 1j and ε 2j are the final strains in two directions, and σ(ε) is the hardening function, ε(ε) is effective strain function, α ε 2 is a function of material parameters. In this paper, the Swift hardening law Equation (7) is used to fit the stress-strain curve of DP steel with uniform elongation.
where σ is stress, K is the strength coefficient, n is the hardening exponent, ε 0 is prestrain, ε p is plastic strain. Based on the non preload test conducted in this study, there is no pre strain. Equations (5) and (6) are transformed into Equations (8) and (9).
By different yield criteria, the FLSD obtained is also different. In this paper, based on the Hill'48 yield criteria selected and the Material properties of the four kinds of sheets, the FLSD data are converted. The equivalent stress σ is solved by the function of the stress σ 1 and σ 2 .
The equivalent strain ε is The proportionality coefficient of equivalent strain to maximum principal stress is The relationship between α and β can be defined by The plastic strain ratio r, is the ratio of the true width strain to the true thickness strain at a particular value of length strain. where ε t is the true thickness strain, ε b is the true width strain. The r-value is provided by the cooperating company, as shown in Table 1. The stress-strain transformation relationship can be solved by combining the above functions. The schematic diagram of the transformation process is shown in Figure 2.

1+ + r rβ
The plastic strain ratio r, is the ratio of the true width strain to the true thickness strain at a particular value of length strain. b t ε r ε  (15) where t ε is the true thickness strain, b ε is the true width strain.
The r-value is provided by the cooperating company, as shown in Table 1. The stress-strain transformation relationship can be solved by combining the above functions. The schematic diagram of the transformation process is shown in Figure 2.

Experimental Design
The material used in this paper was DP590, DP780, DP980, and DP1180 sheets with a nominal thickness of 1 mm. The microstructure of DP steel mainly consists of two phases: ferrite and martensite. The composition of these four steels was provided by the producer, Shanghai Baosteel Group, and is shown in Table 2. The mechanical properties of the four steels were obtained after metal tensile test, as shown in Table 3, where E is the modulus of elasticity, ν is Poisson's ratio, K is strength coefficient and n is the hardening exponent. Plotting FLD requires obtaining the maximum and minimum principal strain data of the sheet, which can be obtained by Nakazima dome tests on specimens of different sizes [18]. In this paper, experiments are designed with reference to GB/T 15825.8-2008, Sheet metal formability and test methods-Part 8: Guidelines for the determination of forming-limit diagrams. The Nakazima test in this paper was carried out using a YZ32-160S hydraulic press by Zibo Aoheng Hydraulic Machinery Co. (Zibo, China). The Stamping test machine is mainly composed of blank, punch, die and blank holder. The structural principle is shown in Figure 3. During the Nakazima dome test, the blank holder moves downwards to provide sufficient pressure to secure the sheet, and then the punch moves upwards to cause deformation of the sheet until the sheet fractures.

DP980
216,000 0.29 1180 0.08 DP1180 218,000 0.29 1300 0.06 Plotting FLD requires obtaining the maximum and minimum principal strain data of the sheet, which can be obtained by Nakazima dome tests on specimens of different sizes [18]. In this paper, experiments are designed with reference to GB/T 15825.8-2008, Sheet metal formability and test methods-Part 8: Guidelines for the determination of forminglimit diagrams. The Nakazima test in this paper was carried out using a YZ32-160S hydraulic press by Zibo Aoheng Hydraulic Machinery Co. (Zibo, China). The Stamping test machine is mainly composed of blank, punch, die and blank holder. The structural principle is shown in Figure 3. During the Nakazima dome test, the blank holder moves downwards to provide sufficient pressure to secure the sheet, and then the punch moves upwards to cause deformation of the sheet until the sheet fractures. In order to obtain FLD data of DP steel of different grades, three groups of specimens with full length of 200 mm and width varies were designed ( Figure 4). The specific dimensions of specimens are shown in Table 4. In order to obtain FLD data of DP steel of different grades, three groups of specimens with full length of 200 mm and width varies were designed ( Figure 4). The specific dimensions of specimens are shown in Table 4.    L  180  180  180  180  180  180  180  180  180  L1  40  40  40  40  40  ----L2  20  20  20  20  20  ----W1  40  60  80  100  120  120  120  160  180  W2  30  50  70  90  110  120  140  160  180  W3  20  40  60  80  100  120  140  160  180 Before the Nakazima test, in order to measure the strain of specimens, an electronic marking machine was used to etch an array composed of 1 mm circular grids on each specimen through electrochemical methods ( Figure 5). The maximum and minimum principal strains were calculated from the deformation variables of the grid before and after the test. Assume that the circular grid is elliptical after distortion, and the maximum axis is recorded as d1 and the minimum axis as d2, and the d1 and d2 are approximately regarded as two principal strain directions on a point in the surface of the specimen sheet. The diagram of grid circle distortion is shown in Figure 6. After the test, the process stress and strain data and the real stress and strain data are calculated according to Equations (16) and (17).   L  180  180  180  180  180  180  180  180  180  L1  40  40  40  40  40  ----L2  20  20  20  20  20  ----W1  40  60  80  100  120  120  120  160  180  W2  30  50  70  90  110  120  140  160  180  W3  20  40  60  80  100  120  140  160  180 Before the Nakazima test, in order to measure the strain of specimens, an electronic marking machine was used to etch an array composed of 1 mm circular grids on each specimen through electrochemical methods ( Figure 5). The maximum and minimum principal strains were calculated from the deformation variables of the grid before and after the test. Assume that the circular grid is elliptical after distortion, and the maximum axis is recorded as d 1 and the minimum axis as d 2 , and the d 1 and d 2 are approximately regarded as two principal strain directions on a point in the surface of the specimen sheet. The diagram of grid circle distortion is shown in Figure 6. After the test, the process stress and strain data and the real stress and strain data are calculated according to Equations (16) and (17).
L2 20 20 20 20 20  ----W1  40  60  80  100  120  120  120  160  180  W2  30  50  70  90  110  120  140  160  180  W3  20  40  60  80  100  120  140  160  180 Before the Nakazima test, in order to measure the strain of specimens, an electronic marking machine was used to etch an array composed of 1 mm circular grids on each specimen through electrochemical methods ( Figure 5). The maximum and minimum principal strains were calculated from the deformation variables of the grid before and after the test. Assume that the circular grid is elliptical after distortion, and the maximum axis is recorded as d1 and the minimum axis as d2, and the d1 and d2 are approximately regarded as two principal strain directions on a point in the surface of the specimen sheet. The diagram of grid circle distortion is shown in Figure 6. After the test, the process stress and strain data and the real stress and strain data are calculated according to Equations (16) and (17).   L2  20  20  20  20  20  ----W1  40  60  80  100  120  120  120  160  180  W2  30  50  70  90  110  120  140  160  180  W3  20  40  60  80  100  120  140  160  180 Before the Nakazima test, in order to measure the strain of specimens, an electronic marking machine was used to etch an array composed of 1 mm circular grids on each specimen through electrochemical methods ( Figure 5). The maximum and minimum principal strains were calculated from the deformation variables of the grid before and after the test. Assume that the circular grid is elliptical after distortion, and the maximum axis is recorded as d1 and the minimum axis as d2, and the d1 and d2 are approximately regarded as two principal strain directions on a point in the surface of the specimen sheet. The diagram of grid circle distortion is shown in Figure 6. After the test, the process stress and strain data and the real stress and strain data are calculated according to Equations (16) and (17).

Experimental Results and Discussion
Specimens after the Nakazima test are shown in Figure 7. By observing the fracture specimens in Figure 7, it is found that the places where the fracture surfaces of specimens of different sizes occur are roughly the same. As the tensile and yield strengths of the steel increase, greater blank holder force is required to compress the sheets. The blank holder force increases from about 100 kN for DP590 to about 300 kN for DP1180, while reducing the forming depth and causing earlier failure. Therefore, steel sheet with high tensile strength and yield strength is difficult to be stamped and forming. The fracture surfaces of specimens of most sizes produce arc fracture tracks near the center of the specimen.
The grid distortion of failure critical area is measured and the data are processed to obtain the maximum and minimum principal strains of each grid. Limit strain data after processing is shown in Figure 8. It is found from the limit strain data of nakazima test that the value of the limit principal strain of sheet forming increases with the increase of the sheet strength, and the value range of the minimum principal strain also expands with the increase of the sheet strength. The trend of the FLC on the left and right sides of the four DP steel is the same, showing a decreasing trend in the left branch as a whole, and an increasing trend in the right branch. specimens in Figure 7, it is found that the places where the fracture surfaces of specimens of different sizes occur are roughly the same. As the tensile and yield strengths of the steel increase, greater blank holder force is required to compress the sheets. The blank holder force increases from about 100 kN for DP590 to about 300 kN for DP1180, while reducing the forming depth and causing earlier failure. Therefore, steel sheet with high tensile strength and yield strength is difficult to be stamped and forming. The fracture surfaces of specimens of most sizes produce arc fracture tracks near the center of the specimen. The grid distortion of failure critical area is measured and the data are processed to obtain the maximum and minimum principal strains of each grid. Limit strain data after processing is shown in Figure 8. It is found from the limit strain data of nakazima test that the value of the limit principal strain of sheet forming increases with the increase of the sheet strength, and the value range of the minimum principal strain also expands with the

FLD of DP Steel Sheets
Strain data of the FLD from the Nakazima test are not enough to draw an accurate FLC. In order to obtain more strain data, the finite element model of Nakazima test was established in Abaqus. The fracture failure maximum principal strain of some specimens is obtained through simulation to enrich the test data. The accuracy of the simulation model is verified by comparing with fracture results of the Nakazima test.
Similar to actual experiments, the simulation model is also composed of blank, punch, die and blank holder. The blank is a deformable shell, and other parts are set as rigids. The thickness of blank is set to 1 mm. After mesh size analysis, 2 mm size mesh can save the calculation time and ensure the accuracy of the calculation, the analysis result is shown in the Figure 9. The mesh size of the sheet is 2.0 × 2.0 mm and mass scaling factor was set to 20 so as to reduce the total simulation time. The element type of the blank is S4R.

FLD of DP Steel Sheets
Strain data of the FLD from the Nakazima test are not enough to draw an accurate FLC. In order to obtain more strain data, the finite element model of Nakazima test was established in Abaqus. The fracture failure maximum principal strain of some specimens is obtained through simulation to enrich the test data. The accuracy of the simulation model is verified by comparing with fracture results of the Nakazima test.
Similar to actual experiments, the simulation model is also composed of blank, punch, die and blank holder. The blank is a deformable shell, and other parts are set as rigids. The thickness of blank is set to 1 mm. After mesh size analysis, 2 mm size mesh can save the calculation time and ensure the accuracy of the calculation, the analysis result is shown in the Figure 9. The mesh size of the sheet is 2.0 × 2.0 mm and mass scaling factor was set to 20 so as to reduce the total simulation time. The element type of the blank is S4R.
punch, die and blank holder. The blank is a deformable shell, and other parts are set as rigids. The thickness of blank is set to 1 mm. After mesh size analysis, 2 mm size mesh can save the calculation time and ensure the accuracy of the calculation, the analysis result is shown in the Figure 9. The mesh size of the sheet is 2.0 × 2.0 mm and mass scaling factor was set to 20 so as to reduce the total simulation time. The element type of the blank is S4R. The simulation uses a user subroutine of Abaqus, VUMAT. VUMAT is used to define the constitutive model of materials which is written in fortran. When the subroutine is called, it is provided with state at the start of the increment (stress, solution-dependent The simulation uses a user subroutine of Abaqus, VUMAT. VUMAT is used to define the constitutive model of materials which is written in fortran. When the subroutine is called, it is provided with state at the start of the increment (stress, solution-dependent state variables) of blocks of material points. The simulation program applies the Hill (1948) orthotropic yield criterion and Swift hardening law Equation (7) for the calculation.
The Hill (1948) orthotropic yield condition reads where σ is the Hill (1948) equivalent stress, σ x represents the component of the Cauchy stress tensor, the six constants, F, G, H, L, M, N are anisotropic parameters. Anisotropy parameters of DP steel sheet in Hill'48 anisotropy yield criterion selected in this paper are calculated by the plastic strain ratio r in the previous paper, and their relationship is shown in Equation (19). L, M, and N are difficult to distinguish, and it is generally considered that L = M = N. Anisotropic parameters of these four steels is shown in Table 5. Establishing contact and friction on the surfaces of the rigids and the deformable shell, and the friction coefficient of the whole model is set to 0.15. Limit all degrees of freedom of the die during the simulation. The first step compresses the blank downwards with the blank holder, and the second step deforms the blank with the punch upwards. The finite element model of the Nakazima test is illustrated in Figure 10. Establishing contact and friction on the surfaces of the rigids and the deformable shell, and the friction coefficient of the whole model is set to 0.15. Limit all degrees of freedom of the die during the simulation. The first step compresses the blank downwards with the blank holder, and the second step deforms the blank with the punch upwards. The finite element model of the Nakazima test is illustrated in Figure 10. The bulge fracture simulation was carried out for DP590, DP780, DP980 and DP1180 respectively, and the strain nephograms at the fracture time were obtained. The simulated strain nephograms of some sheets are shown in Figure 11. At the same time, the maximum and minimum limit principal strain data of the bulge specimens of different sizes are recorded and sorted out. The bulge fracture simulation was carried out for DP590, DP780, DP980 and DP1180 respectively, and the strain nephograms at the fracture time were obtained. The simulated strain nephograms of some sheets are shown in Figure 11. At the same time, the maximum and minimum limit principal strain data of the bulge specimens of different sizes are recorded and sorted out. After obtaining the limit strain data and FLD0 points of each sheet, binary linear fitting is carried out. The expression after fitting is shown in Equation (20), and the goodness of fit is 0.968. Combine the test data points with the simulation data points to form FLD, as shown in Figure 12 where n is the hardening exponent, t representing the thickness of the steel sheet (only steel sheets with thickness less than 3 mm are considered in this paper). After obtaining the limit strain data and FLD 0 points of each sheet, binary linear fitting is carried out. The expression after fitting is shown in Equation (20), and the goodness of fit is 0.968. Combine the test data points with the simulation data points to form FLD, as shown in Figure 12  where n is the hardening exponent, t representing the thickness of the steel sheet (only steel sheets with thickness less than 3 mm are considered in this paper).
where n is the hardening exponent, d 1 and d 2 are parameters to be fitted. Their fitted values are shown in Table 6.
(20). Curves on both sides of the lowest point are fitted based on Equation (21). 0 =1.656 +0.032 0.025 where n is the hardening exponent, t representing the thickness of the steel sheet (only steel sheets with thickness less than 3 mm are considered in this paper).
where n is the hardening exponent, d1 and d2 are parameters to be fitted. Their fitted values are shown in Table 6.

Forming Limit Stress Diagram
The limit stress data of DP590, DP780, DP980 and DP1180 are determined by combining the conversion function with the performance parameters of the four kinds of sheets. To fit the FLSD formula, the FLSD fitting curves of DP590, DP780, DP980 and DP1180 are established, and the limit stress data points are processed. Observe the curves in the Figure 13, the trend is relatively gentle compared with FLD. The limit principal stress ranking between each type of DP steel is consistent with the strength ranking of the steel. The curves are obtained by multiple linear regression fitting in combination with the performance parameters of the sheets, the coefficient of the quadratic term of curves is related to the hardening exponent n. as shown in Equation (22). The expression established after fitting is as follows: where σ 1 , σ 2 are the maximum and minimum limit principal stresses, K is the strength coefficient, and n is the hardening exponent. Due to fluctuations in the stress data obtained through experiments, there are inevitable errors in the fitting formula. The RMSE (root mean squared error) is a standard for intuitively testing the accuracy of regression equations.
where n is number of items, Y i is original or observed value, f (x i ) is value from fitting formula.
where 1  , 2  are the maximum and minimum limit principal stresses, K is the strength coefficient, and n is the hardening exponent.
where n is number of items, i Y is original or observed value,   i f x is value from fitting formula. The RMSE calculated from the FLSD of DP590, DP780, DP980 and DP1180 are 6.64, 10.50, 11.70, 12.91. This means that the fitting formula in this study needs to consider the existence of small errors when applied.
According to the expression observation curve and the strain hardening function, when solving the equivalent stress σ, its value is exponentially related to the strength coefficient K and hardening exponent n. Therefore, the limit stress value in the FLSD curve increases with the increase of K and n. The K and the strength of the sheet have the same trend of change in the performance, while the n decreases with the increase of the strength of the sheet. The function shows that the k-value determines the value of the maximum principal stress, while the n-value affects the curvature of the curve. It can be seen from the Figure 13 that the limit stress value of DP1180 sheet is greater than that of DP590, so the limit stress value in the sheet generally increases with the increase of sheet strength, and the influence of strength coefficient K on FLSD is greater than that of hardening exponent n.

Nakazima Test Simulation
The FLSD fitting formulas obtained in the previous chapter can guide sheet metal forming by computer simulation. The settings of the simulation process are basically the same as above. FLSD can be conveniently applied to form a simulation as a fracture criterion by using VUMAT subprogram in Abaqus. The simulation model uses shell elements, which define two principal stresses in the plane direction and the principal stresses in the normal direction is zero. Call the two principal stresses of the shell element through the vumat subroutine, and divide the two principal stresses into the maximum principal stress and the minimum principal stress based on their magnitude relationship. Once the numerical relationship between the maximum principal stress and the minimum principal stress of a element exceeds the limit value determined by the FLSD function, the element will be deleted to simulate failure.
Taking the bulging test of specimens with a width of 60 as an example, Figure 14 shows the test results and Figure 15 shows the simulation results. The first diagram of each simulation test shows the maximum principal stress distribution at the moment before fracture, while the second shows the principal stress distribution at the moment of fracture. A comparison of test and simulation fracture depth is shown in the Table 7. pal stress of a element exceeds the limit value determined by the FLSD function, the element will be deleted to simulate failure.
Taking the bulging test of specimens with a width of 60 as an example, Figure 14 shows the test results and Figure 15 shows the simulation results. The first diagram of each simulation test shows the maximum principal stress distribution at the moment before fracture, while the second shows the principal stress distribution at the moment of fracture. A comparison of test and simulation fracture depth is shown in the Table 7.  Taking the bulging test of specimens with a width of 60 as an example, Figure 14 shows the test results and Figure 15 shows the simulation results. The first diagram of each simulation test shows the maximum principal stress distribution at the moment before fracture, while the second shows the principal stress distribution at the moment of fracture. A comparison of test and simulation fracture depth is shown in the Table 7.    Using FLSD as the failure criterion for the simulation model of the Nakazima dome test, the calculation results are basically consistent with the FLD and also consistent with the experimental results, which proves the reliability of FLSD as a fracture criterion.

Box Forming Test
In order to verify the fitted forming limit curve, the box forming test was carried out by YZ32-160S hydraulic press. The test equipment was roughly the same as the bulging tests. Punch and blank holder are replaced with box die. In order to study the effect of different fillet on the forming effect, a box die with four sizes of fillet was experimentally designed. The size of box die with unequal fillet is shown in Figure 16. The rectangular sheet of DP590, DP780, DP980, and DP1180 were tested by the press with the size of 190 mm in length and 180 mm in width. A simulation model was designed based on the mold size the simulation model is shown in Figure 17. Simulating the forming process of the box forming experiment, and the limit stress value is obtained. FLSD damage data is imported into Abaqus, and the accuracy of FLSD data of each sheet is verified by outputting the maximum principal stress through damage evolution settings.  Figure 17.Simulating the forming process of the box forming experiment, and the limit stress value is obtained. FLSD damage data is imported into Abaqus, and the accuracy of FLSD data of each sheet is verified by outputting the maximum principal stress through damage evolution settings.  The test and simulation results of the box forming and the simulated principal stress values are shown in Figure 18. By observing the FLSD simulation results of the four types of sheet DP590, DP780, DP980 and DP1180, and comparing maximum and minimum principal stress values obtained from the stamping forming test and simulation of each type of sheet box, it is found that the maximum principal stress gap of the four types of sheet is large and the minimum principal stress gap is small. Maximum principal stress of DP590 is kept around 1000 MPa; the maximum principal stress simulation value of DP780  Figure 17.Simulating the forming process of the box forming experiment, and the limit stress value is obtained. FLSD damage data is imported into Abaqus, and the accuracy of FLSD data of each sheet is verified by outputting the maximum principal stress through damage evolution settings.  The test and simulation results of the box forming and the simulated principal stress values are shown in Figure 18. By observing the FLSD simulation results of the four types of sheet DP590, DP780, DP980 and DP1180, and comparing maximum and minimum principal stress values obtained from the stamping forming test and simulation of each type of sheet box, it is found that the maximum principal stress gap of the four types of sheet is large and the minimum principal stress gap is small. Maximum principal stress of DP590 is kept around 1000 MPa; the maximum principal stress simulation value of DP780 The test and simulation results of the box forming and the simulated principal stress values are shown in Figure 18. By observing the FLSD simulation results of the four types of sheet DP590, DP780, DP980 and DP1180, and comparing maximum and minimum principal stress values obtained from the stamping forming test and simulation of each type of sheet box, it is found that the maximum principal stress gap of the four types of sheet is large and the minimum principal stress gap is small. Maximum principal stress of DP590 is kept around 1000 MPa; the maximum principal stress simulation value of DP780 is about 1100 MPa; the maximum principal stress of DP980 fluctuates around at 1200 MPa; the maximum principal stress of DP1180 reaches 1300 MPa. The maximum principal stress and the minimum principal stress values of the four kinds of DP steel sheet show an increasing trend, but the growth range of the minimum principal stress is small with the increase of the sheet strength. The comparison of box forming test and simulated stamping depth is shown in the Table 8. The test and simulation results of the box forming and the simulated principal stress values are shown in Figure 18. By observing the FLSD simulation results of the four types of sheet DP590, DP780, DP980 and DP1180, and comparing maximum and minimum principal stress values obtained from the stamping forming test and simulation of each type of sheet box, it is found that the maximum principal stress gap of the four types of sheet is large and the minimum principal stress gap is small. Maximum principal stress of DP590 is kept around 1000 MPa; the maximum principal stress simulation value of DP780 is about 1100 MPa; the maximum principal stress of DP980 fluctuates around at 1200 MPa; the maximum principal stress of DP1180 reaches 1300 MPa. The maximum principal stress and the minimum principal stress values of the four kinds of DP steel sheet show an increasing trend, but the growth range of the minimum principal stress is small with the increase of the sheet strength. The comparison of box forming test and simulated stamping depth is shown in the Table 8.    Through comparative test and simulation, the stamping depth value and test error of the simulation model are less than 10% when the blank holder force value is the same, which effectively shows that the FLSD values of sheets are consistent with the established FLSD data, and proves the accuracy of the established FLSD.

Conclusions
In this paper, the Nakazima test of DP590, DP780, DP980 and DP1180 is designed, and the required forming limit data points of four materials are obtained. The numerical simulation of Nakazima test is carried out to improve the accuracy of FLD fitting. Finally, the FLD and its unified expression of four kinds of dual phase steel sheets are established. Based on the stress-strain relationship of M-K theory, the conversion between FLD and FLSD is realized, and the FLSD of four kinds of dual phase steel sheets is established. By studying the correlation with the sheet property parameters, the unified FLSD expression of four kinds of dual phase steel sheets is established. The main conclusions of this research are: (1) The limit stress data of the four kinds of DP steel (DP590, DP780, DP980, DP1180) were obtained and FLSD was plotted through the Nakazima dome test and simulation. The differences in limit stresses of the four steels and the effect of such differences on forming were compared. It was verified by numerical simulation that the FLSD obtained in this paper is within 10% error from experiment and can be used as a fracture criterion to accurately simulate the fracture of DP steel sheets. (2) The influence strength coefficient K and hardening exponent n of DP steel of FLSD function is analyzed, and these parameters are introduced into the FLSD function. The K and the strength of the sheet have the same trend of change in the performance, while the n decreases with the increase of the strength of the sheet. The function shows that the k-value determines the value of the maximum principal stress, while the n-value affects the curvature of the curve. A high fitting FLSD calculation method is obtained. The FLSD function proposed in this paper can be extended to the fracture calculation of DP steel with other strength grade, reducing the difficulty of FLSD application.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

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