Experimental determination of the failure surface for DP980 high-strength metal sheets considering stress triaxiality and Lode angle

To meet requirements for reduced fuel consumption of cars, the use of components made of sheets of high-strength steel instead of conventional steel has been on the rise. However, low ultimate elongation of high-strength steel often causes problems during plastic deformation and more research is needed to improve failure predictability. Ductile failure models available in commercial finite element analysis (FEA) packages require the material’s tensile strength and failure strain for failure prediction. For stress states that are more complex than the uniaxial case, accurate prediction of the how, when, and where failure occurs has been become problematic and it has been investigated by numerous researchers. In this study, we investigate the prediction of failure in DP980 sheets under triaxle stress states. We first determine the shapes of specimens using certain triaxial stress states, such as pure shear, uniaxial tension, biaxial deformation, which are induced by corresponding tensile tests. When failure occurs, equivalent strain at the failure locus is obtained by means of digital image correlation (DIC) and then plotted against triaxiality and Lode angle, based on which the triaxiality failure diagram (TFD) is established to implement in the FEA program of LS-DYNA. Validation is made by comparing the numerical results with burring test data. Good agreement was found for failure locus and strain distribution at the time of failure.


Introduction
In the context of fuel consumption, weight reduction plays an ever-increasing role in the car industry, which means that conventional steels are being replaced by high-strength steel (HSS) sheets showing elevated tensile strength. However, the higher tensile strength of such HSS sheets is usually at the expense of ultimate elongation, leading to frequent failure under plastic deformation. This issue calls for in-depth research, both numerical and experimental, on how, when, and where failure occurs and how this can be predicted accurately. Research on high-strength metal sheets has been actively pursued [1,2], and Fig. 1 gives an impression on which high-strength sheets have been developed and contributed to weight reduction [3]. For connection with other structures, burrs must be removed from sheets in which holes were drilled. It is well known that formability changes with the quality of the surface of the hole [4]. In particular, as shown in Fig. 2, high-strength metal sheets are prone to tearing during burring, owing to their low ultimate elongation. Commonly, the failure limit of metal sheets is defined as the elongation at failure in tensile testing, and failure is determined based on the forming limit diagram (FLD), which is a plot of major against minor strains [6]. FLD can be obtained from dome height tests by applying the circle grid analysis (CGA) method, which geometrically considers necking and failure at the time of fracture [7,8]. During sheet metal forming, FLD becomes inapplicable (or at least limited in its applicability) once the deformation path departs from its original one or if the stress state is not (purely) tensile. As formability is here determined based on the assumption that the strain ratio is constant for stress states in the region between uniaxial tension and biaxial tension. Further, the failure strain at the failure locus cannot be deduced accurately from the deformed grid geometry after the dome height test. These limitations associated with the strain-based FLD led to the development of a stress-based FLD [9] and triaxiality failure diagram (TFD) to evaluate the formability for a large range of stress states [10,11].
Commercial FEA packages are equipped with diverse ductile damage models, which can be applied to the problem of failure in sheet metal forming. As an initial approach, Cockcroft and Latham [12] assumed that fracture occurs when the plastic work associated with the maximum tensile stress reaches a limit. Oh et al. [13] modified the equation proposed by Cockcroft and Latham by normalizing it. Brozzo et al. [14] implemented the influences of tensile stress and mean stress and corrected the equation by fitting it to experimental data and numerical results. Oyane [15] introduced a yield condition for materials with voids and plastic potential and assumed that fracture occurs once the volumetric strain of the void reaches a limit. The damage models mentioned above require tensile strength and fracture strain obtained from the tensile tests. As became clear in Choi et al. [5], accurate prediction of shear failure phenomena observed in forming of high-strength steels is problematic with the above-mentioned models (Fig. 2).
Consequently, the triaxiality parameter η was introduced as a parameter representing the degree of stress triaxiality and fracture was then predicted based on η-ε f plots. However, according to the study by Bai and Wierzbicki [16], the fracture strain is not only a function of η but also the function of Lode angle θ, which is the angle between the principal stresses on the yield surface.
In this study, we derive a failure surface that considers η and θ for the HSS (DP980). Tensile tests were performed with various specimen shapes to consider all necessary stress states to determine η. Failure simulation is performed in the FEA program of LS-DYNA [17] and among the many damage models available, the GISSMO (Generalized Incremental Stress State damage MOdel) [18,19] was chosen, as the GISSMO matches the numerical and experimental results regarding the failure behavior under all stress states. Moreover, the GISSMO has the advantage that it can be made mesh sizeindependent by means of a scale factor. For any stress state, the failure locus can be predicted numerically through reverse engineering analysis by fitting the load-displacement curve obtained by FEA to tensile test data. Based on this, the triaxial stress surface is determined and burring tests are carried out for validation. Finally, the load-punch stroke-curve obtained though burring tests are compared with numerical results obtained using the FLD and TFD.

Theory
Among the many approaches, the strain-based FLD is applied often for sheet formability prediction. This is because the measurement of stresses is practically impossible and strain measurement is relatively straightforward during stamping. The forming limit curve (FLC) established from forming limit strains helps to determine the region of deformation and the probability of failure. The experimental and theoretical establishment of FLCs has been studied in abundance.
Theoretical approaches focused on the material's loss of stability and the establishment of failure criteria. Theoretical approaches for formability evaluation are relatively advanced but mainly experimental methods are used. The most common methods are certainly the Marciniak [21] and Nakazima [22] tests. In the Nakazima test, a semi-spherical punch is employed to stretch a sheet. By varying the specimen width, diverse deformation paths between uniaxial and biaxial can be established.   Figure 3 shows the FLC and the FFL (fracture forming line) for AA7075-O specimens with thickness of 1.6 mm obtained by the Nakazima test. However, the Nakazima method assumes a constant deformation path (ρ = ε 2 / ε 1 ), which limits its applicability. As shown in Fig. 3, beyond the FLD, the deformation path departs from its previous path and failure advances without significant change of the minor strain [20]. Further, precise determination of the failure point becomes problematic because the failure strain is calculated based on the specimen shape after fracture. To overcome these deficiencies, stress-based failure limit is determined here and to evaluate deformability for a large range of stress states, including not only tension but also shear, failure limit diagrams accounting for stress triaxiality have to be applied.

Experimental setup
As mentioned in the "Theory" section, by using specimens of different widths in the Nakazima test, we can determine formability for diverse stress states. However, this does not include shear and, what is more, friction between punch and specimen may affect test results.
By contrast with the Nakazima test, friction does not have an impact in tensile testing and by varying the specimen shape, we can finally induce any tensile stress state, including shear, and obtain corresponding forming limits. In this study, to determine failure by tensile testing only, we chose seven representative deformation modes. Suitable specimen shapes had previously been determined (by FEA) so that during the tensile test, failure occurs in the center of the specimen and that η remains constant throughout the test. Specimen geometries that give representative η values are shown in Fig. 4. As tensile tester, the universal testing machine (UTM) Shimadzu AG-X was employed; Table 1 lists parameters of the UTM and the test settings.
Each test was performed at least three times per specimen to ensure reproducibility. DP980 sheets of 1.0-mm thickness with rolling direction (RD) were used as a test material. DP steel contains generally 0.06-0.15 wt% C, 1.5-3 wt% Mn, Si and Nb, Cr, Mo. Its overall behavior, however, is controlled by volume fraction, morphology (size, aspect ratio, interconnectivity, etc.), grain size and carbon content [23]. The chemical composition of the DP980 used in this study is given in with strain gauges-as common in standard tensile test methods [24,25]-digital image correlation (DIC) was used here to determine local failure strains. A speckle pattern is applied on the specimen surface to generate contrast and from the difference between the patterns before and after the test, local strains are measured. Figure 5 shows the whole test rig with the DIC device and the speckle pattern on a uniaxial tensile specimen. During the tensile test, two stereo ultra-high-speed cameras captured 900 frames per second. Detailed camera parameters are provided in Table 3. Based on these recordings, the strain is finally computed in the commercial DIC software GOM ARAMIS [26]. Based on DIC, Fig. 6 shows the equivalent strain distributions in shear, uniaxial tension, and large-R and small-R specimens right before fracture. With an extensometer, we could get strain values in the failure region only up to a point of uniform elongation, where necking occurs (around 10%). However, local elongations obtained with DIC reached up to 46% in the fracture region. For establishing a stress-strain relation in the plastic range from tensile test data, i.e., an equivalent stress-equivalent strain relation, we need some mathematical formulation that can be implemented in an FE program. Generally, it is difficult to obtain true-stress-true strain data because from the point of necking up to failure, the cross-sectional area of the tensile specimen does not reduce at a constant aspect ratio.
Thus, we need to be able to obtain the maximum uniform elongation. For the region beyond uniform elongation, a stress-strain relation can be established by fitting equivalent stress-equivalent plastic strain data to the Hollomon [σ eq = Kε eq n ] or Swift equations [σ eq = K(ε 0 + ε eq ) n ], where K, ε 0 , and n are strength coefficient, initial strain and strain hardening coefficient, respectively. Here, we applied the Swift equation for DP980 (Fig. 7); coefficients are listed in Table 4.

Triaxiality failure diagram (TFD)
The TFD is a plot of triaxiality η vs. equivalent strain ε eq , where η is a scalar and defined as the ratio of equivalent stress σ eq to mean stress σ m . In conjunction with the von Mises yield condition and isotropic hardening, η can be expressed as   Here, σ 1 , σ 2 , and σ 3 are principal stresses in directions 1 (longitudinal), 2 (width), and 3 (thickness), respectively. η is a stress-based expression representative of the load path and independent of the deformation path. Generally, during sheet deformation, the stress state is in-plane (σ 3 = 0) so that triaxiality can be calculated readily from Eq. (1) by inserting σ 1 and σ 2 . For instance, for shear σ 1 = σ 2 = 0 and thus η = 0; for uniaxial tension σ 1 = σ eq , σ 2 = 0 and thus η = 1/3; and for plane strain, η = 0.577. By plotting the failure strains for all considered stress states against the corresponding η values, we can establish the failure prediction curve. Obviously, failure prediction becomes more accurate, when the more combinations of σ 1 and σ 2 were considered in tensile tests. Figure 8 represents the TFD and the instability curve for diverse stress states (η) obtained in DP980. Independent of specimen and stress state, necking occurs at the maximum load; therefore, this point can be defined as an initiation of instability. Failure (circles) and instability initiation points (squares) are the equivalent strains derived by means of DIC from images at the time of failure and maximum load, respectively. Considering only the tensile region (η > 1/3), the failure strain decreases monotonically with increasing η. This is a typical triaxial stress state-dependent failure behavior and it is in accordance with the failure models of Johnson-Cook [27], Gurson [28], and others. However, for shear (η = 0), the failure strain is lower than those from uniaxial tension but larger than those from other stress states (large-R, notch). This is not in agreement with the exponential failure models by Johnson-Cook [27] and Gurson [28], where the failure strain decreases monotonically with increasing η. This means that while accurate presentation of the region between tension and plane strain is possible, an error occurs when extending to shear, biaxial tension, and compression. Thus, formulating failure for triaxial stress states requires a bi-polar expression or a combination of exponential expressions. To do so, we adopted the failure model by B&W [16], as mentioned in the "Introduction" section, and express the fracture strain ε critical p , as a fitting function of η and the Lode angle θ

Experimental TFD for DP980
where C 1~C3 , D 1~D6 are material parameters. For plane stress (σ 3 = 0), the Lode angle θ can be expressed as θ = 0 means shear or plane strain state, whereas θ = 1, − 1 mean uniaxial tension, and uniaxial compression states for rodshaped specimens. On the one hand, instability initiates earliest in uniaxial tension. Among all specimens considered here, the lowest equivalent strain that leads to necking is ε eq = 0.12. Since the equivalent strain is ε eq = 0.6 at failure, the phase of diffuse necking and local necking is extremely long. It may result in a huge error in numerical analysis when we simply assume that failure occurs rapidly after necking and the necking strain in uniaxial tension is the failure strain. On the other hand, for shear and plane strain, instability initiation nearly coincides with failure so that the assumption that necking and failure occur at the same time is valid (Fig. 8). This is backed by the observation that in sheet drawing, the thickness cannot be maintained uniformly and once necking initiates, the sheet tears immediately. Thus, application of the FLD for formability evaluation is even more problematic when the deformation path varies during forming, even for tensile-dominated sheet drawing analysis.
Fitting the experimental data in Fig. 8 to the B&W model, we get the relation between triaxiality, Lode angle, and failure strain as depicted in Fig. 9. The six B&W model parameters given in Table 5 were determined from in total seven test data sets by using functions of at least second order. The highest failure strain is obtained for compression (η = − 1/3). When we look in positive η-direction, the two curve sections are V-shaped, with the first-the region between uniaxial compression and uniaxial tension (− 1/3 ≤ η ≤ 1/3)-hitting its low at η = 0 (shear) whereupon it increases toward uniaxial tension (η = 1/3) and the second-the region between uniaxial tension and biaxial tension (1/3 ≤ η ≤ 2/3)-having its minimum at plane strain. If we projected the curves on the ε f -η-plane, we would get an asymmetric W-shaped curve.

Burring test
To validate the accuracy of the TFD established in the "Experimental TFD for DP980" section, we employed the burring test device as shown in Fig. 10. As specimens,    DP980 sheets (of thickness 1.0 mm) equal to those for obtaining the TFD were used. The rectangular specimens of 95 mm in both length and width are chosen with a hole diameter of 5 mm in the center. The specimen is positioned between die and blank holder and a cylindrical punch concentric with the hole in the specimen is moved upwards. With increasing upward movement, the hole in the specimen widens and burrs form it. To capture the specimen deformation inside the mold during the test, we employed, similar to the method used for tensile testing, two stereo cameras making 1000 fps. From the many full field pictures captured during the test, we can extract local strains at the time of fracture and locus accurately. Figure 11 shows images of a deformed specimen at the time of fracture. In experiments failure is observed in four spots along the circumference. Figure 11a, b provide top and side views, respectively, and Fig. 11c shows the equivalent strain distribution from the DIC techniques. The height of the deformed specimen is 5.2 mm and the local equivalent strain at the point of failure is 0.42, which is larger than the value of 0.35 that the previously obtained TFD provided for plane strain (η = 0.577) and lower than the maximum value of 0.6 for uniaxial tension. Thus, failure occurs depending on triaxiality. The deformation path, i.e., the minor strain to major strain ratio ρ = ε 2 /ε 1 , can be determined by tracing the strain at failure locus and it may change with the hole diameter. Here, ρ is − 1.7 (= 0.34/− 0.2) at failure locus; this means that the failure locus is subjected to a stress state between shear and uniaxial tension. Therefore, it can be predicted that tearing sets in when the total strain reaches values between the failure strains in the TFD for shear and uniaxial tension.

Simulation of burring test
Transferring the specimen shape and test conditions from the experiment, we established an FE model consisting of specimen and mold in LS-DYNA [17] (Fig. 12). The specimen was modeled as a disk-shaped blank with a hole of diameter = 5 mm. The mesh consists of 3-and 4-node shell elements. The minimum element size is 0.25 mm and no re-meshing or mesh refinement algorithm was used to minimize the error associated with mesh size. Due to symmetry reasons, only a quarter model was established. Material properties of DP980 as given in Table 4 (which were obtained from tensile tests) are assigned to the specimen. The piecewise linear plasticity model is implemented in LS-DYNA (*MAT_24, *MAT_ PIECEWISE_LINEAR_PLASTICITY) [17] as flow stress model. Damage modeling and element removal were performed according to the GISSMO [18,19]. (The corresponding keyword in LS-DYNY is *MAT_ ADD_EROSION.) According to the GISSMO damage model, damage initiates when the stress starts to decrease beyond a certain limit on the true stress-true strain curve. Damage initiation occurs after plasticization and then damage accumulates gradually until failure occurs when the damage parameter D becomes a certain value. D = 0 denotes the region prior to necking and D = 1 denotes failure, we write D Here, ε f , ε eq , and n D are equivalent fracture strain, equivalent strain, and damage coefficient, respectively.
Damage accumulates according to the following rule where, Δε eq is the equivalent plastic strain increment. Further, an instability factor F is introduced, and it is assumed that F grows with ongoing deformation. For an increment in F, we write Here, where D c is defined as the onset of instability, i.e., the D value where F reaches unity, and ε eq,loc is the triaxiality dependent equivalent plastic instability strain. Softening caused by void nucleation and growth is accounted for by reducing the stress in an element according to Here, σ r is the reduced flow stress when plasticization is considered, σ is the stress that would be induced if there was no plasticization and m is the fading exponent, governing material softening. Figure 13 schematically shows the influences of parameters on damage progress according to the GISSMO, and the stress-strain curve. The GISSMO parameters used in this study are listed in Table 6. Here, to obtain failure strains independent of mesh density, we applied a scale factor that had been determined previously. Accordingly, the results become mesh-independent (Fig. 14). To validate the failure surface established in this study, we compare punch stroke-load data obtained by the burring test explained above ( Fig. 10) with those obtained with FE model in Fig. 12  (Fig. 15a). Black circles represent experimental data while red and blue curves represent FE results disregarding and considering the GISSMO, respectively. It becomes clear that, with GISSMO, numerical and experimental punch stroke-load curves nearly coincide, while without GISSMO, the numerical curve departs from the experimental curve with the initiation of necking. Even the equivalent strain distributions around the failure locus are quite similar in addition to the validation of the failure surface established in this study (Fig. 15b).

Conclusion
Summarizing this study, we can derive the following conclusions.
(1) For accurate failure prediction in sheet metal forming for stress states other than tensile, it is necessary to establish a failure curve that considers stress triaxiality; FLD itself is insufficient under general triaxial stress states. (2) We determined diverse specimen shapes using which we can generate shear and tensile stress states (excluding compression). (Local) failure strains were experimentally obtained by means of DIC. The failure surface, plot of failure strain ε f against triaxiality η, and Lode angle θ were established. Different ε f were measured for different θ. (3) Using the GISSMO parameters, we could implement 3D failure surface using B&W damage model in FEA. Loadpunch stroke curves, failure loci, and failure strains obtained by real burring test and corresponding FE analysis were compared, and good agreement was found only with the GISSMO, and thereby validating the failure surface obtained for DP980.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.