Comparison of MMC and BW fracture models for constructing 3D ductile fracture envelope of AA6063 during cold forming

: The 3D ductile fracture envelopes of AA6063-T6 was developed to predict and prevent its fracture. Smooth round bar (SRB) tension tests were carried out to characterize the flow stress, and a series of experiments were conducted to characterize the ductile fracture firstly, such as notched round bar (NR) tension tests, compression tests and torsion tests. These tests cover a wide range of stress triaxiality (ST) and Lode parameter (LP) to calibrate the ductile fracture criterion. Plasticity modeling was performed, and the predicted results were compared with corresponding experimental data to verify the plasticity model after these experiments. Then the relationship between ductile fracture strain and ST with LP was constructed using the modified Mohr–Coulomb (MMC) model and Bai and Wierzbicki (BW) model to develop the 3D ductile fracture envelope. Finally, a new ductile damage model was proposed based on the 3D fracture envelope of AA6063. The final results show that the predicted results from the proposed ductile fracture model showed good agreement with experimental results.

AA6063 aluminum alloy is an Al-Mg-Si alloy with excellent plasticity and machinability. It is widely used in the automobile, construction and energy industry.
The plastic forming processes are commonly classified into hot working and cold working according to whether the processing temperature exceeds the recrystallization temperature of alloys [1] . Cold working is usually carried out at room temperature and leads to improvement in mechanical properties and decrease in plasticity due to increase in dislocation density. It is necessary to improve the mechanical properties of AA6063 by cold working. However, cold working is easy to cause plastic damage, so it is very important to investigate the mechanism of plastic damage in cold working process [2,3] .
The stress state of a point is usually represented by stress triaxiality (ST) and the Lode parameter (LP). From the microscopic point of view, the plastic damage is caused by the nucleation, growth, and aggregation of voids under a high ST. Under a lower ST, the shear band is the main leading cause of ductile fracture [4,5] . Based on this damage formation mechanism, three types of damage models are constructed: (1) failure criterion, (2) porous plasticity theory, and (3) phenomenological damage model [6] . For the first type of model, the damage is obtained by integrating the internal stress variables [7] . Although this is a simple model, it cannot predict the complex deformation path and large plastic deformation. For the second approach, the effect of voids on the von Mises yield surface is considered in the model, as proposed by Gurson [8] . Then, the GTN damage model was formed by the supplement of Tvergaard [9] and Needleman [10] . The advantage of this model is that it is more accurate for the damage prediction under a relatively low ST. On the other hand, its many parameters make calculations complicated and sensitive to mesh size.
For the third approach, the phenomenological damage model was first used to solve engineering problems. Phenomenological damage models can be divided into two types: coupled and uncoupled damage models. Lemaitre proposed a continuous damage mechanical model with damage variable D as internal variable, and established a coupled phenomenological damage model. [11] . On this basis, He where is ST, and ̅ is LP. Therefore, it is essential to construct a 3D fracture envelope in the space ( ̅ , , ̅ ) for uncoupled damage. Bai and Wierzbicki first constructed the 3D fracture envelope in the space ( ̅ , , ̅ ), which defines the fracture strain as a function of ST and LP [13] . The Mohr-Coulomb model is formed based on the maximum shear stress to predict shear damage. Bai and Wierzbicki use the relationship between principal stress and ST and LP to transform the Mohr-Coulomb model into the MMC model [14] .
The phenomenological damage model uses parameters that are calibrated from experimental fracture tests. For a precise calibration of a ductile fracture criterion, fracture experiments covering a wide range of stress states must be performed [15] . The fracture of materials with positive ST can be obtained by using SRB tension test, NR tension test, and plate tension test [16][17][18] . The fracture of materials with nearly zero ST can be obtained by the torsion and notched tube tension tests [19,20] . The fracture of materials with negative ST can be obtained by a cylindrical compression test. Because cracks are difficult to form by the cylindrical compression test, the notched cylindrical compression test is used to control the formation of cracks [21] .
The phenological damage model of aluminum alloys has been reported. Lou developed a weight function of an uncoupled shear ductile fracture criterion and applied it to model the ductile fracture of AA6082-T2 [22] .  [23] . Hossein et al. used the MMC model to AA6061-T6 to examine the damage evolution in the roll forming of U-channel sections [24] . However, a 3D fracture envelope has not yet been constructed for AA6063 alloy.
The paper aimed to develop an uncoupled phenomenological plastic damage model of AA6063 and compare it with other models in terms of accuracy to predict the damage of AA6063 during cold deformation. To obtain a wide range of stress triaxiality variation, tension, torsion, and compression tests were conducted to determine the relationships between the fracture strain and ST and LP. The relationship between the stress state and fracture strain is determined by finite element (FE) simulation. The 3D fracture envelope was calibrated by the GA optimization technique using the relationship between ST, LP, and fracture strain. The developed damage model was implemented into Abaqus 6.14 for the FE simulation. The accuracy of the model was verified by the specific force-displacement curve and the distribution of the damage field. Finally, the two models were compared in terms of fitting quality and damage prediction ability.

Material description
AA6063-T6, which is an Al-Mg-Si alloy with medium strength and heat treatment, was used in these experiments. It is widely used in the construction and transportation industry, and its chemical composition is shown in Table 1.  Fig. 2(a). To obtain the effect of different ST on ductile fracture, different notched specimens were designed according to Eq. (2) [13] : where η is the ST, R is the radius of the neck in the notched bar specimen, and a is the radius of the smallest cross section.  (2) were 20 mm, 10 mm, and 5 mm, respectively. The dimension and shape of these specimens are shown in Fig. 2(a). The gauge lengths of SRB and NR were both 25 mm. The SRB and NR tension tests were carried out with a speed of 0.25 mm/s until the specimen was broken.
Because the hydrostatic pressure of the material is close to 0 in the pure shear state, shear tests were carried out to obtain the ductile fracture at low ST [20] . Through the torsion test, the specimen reached a pure shear state. The shape and dimension of the torsion specimen are shown in Fig. 2   (a) SRB and NR tension test, (b) torsion test, and (c) compression test.

Experimental results
Macroscopic fracture surfaces were investigated to study the fracture mechanism of specimens. Fig. 3 shows the macroscopic fracture phenomena of different tension specimens. When the initial ST is low, the cracks grow along 45° direction, such as for SRB and NR20. This is because when the ST is low, the stress state of these specimens varies between uniaxial tension and balanced biaxial tension, and the maximum shear stress is approximately 45°. The maximum shear stress leads to the coalescence of voids in 45° direction [25] . The macroscopic fracture surface of the specimens with high initial ST (NR5, NR10) mainly presents the shape of "cup cone".
The cup-cone fracture has a macroscopic appearance with numerous dimples in the central interior area of the fracture surface and an oblique appearance at the outer surface. The fracture is caused by the dominant effect of the tensile stress [26] . There are many dimples in the inclined surfaces in Fig. 3, indicating that ductile fracture takes place due to the maximum tensile stress in NR5 and NR10.  increases when necking occurs [19,27] . With the increase in ST, the plasticity of the material decreases, and the rate of void occurrence increases. Fig. 4 shows that NR5 is the first to enter the softening stage and fracture the largest ST. The slowest to enter the softening stage and fracture is the SRB with the minimum ST.  that the bearing capacity of the cylinder with the height-to-diameter ratio of 1 is stronger than that of the specimen with the height-to-diameter ratio of 1.5. Due to the friction between the specimen and test equipment, the specimen has a severe nonuniform deformation, resulting in a drum shape of the specimen. Due to the increase in the cross-sectional area and the inhibition of void nucleation under negative ST, the load-bearing capacity of the specimen increases; thus, the forcedisplacement curve continuously increases.

Stress-strain relationship
We assumed that the material is isotropic and conforms to the von Mises yield criterion. Firstly, the engineering stress-strain curve in Fig. 7 was calculated according to the SRB force-displacement curve in Fig. 4. According to the stressstrain curve in Fig. 7, Young's modulus of elasticity was 14611 MPa, from which the true stress-strain curve in Fig. 7 was calculated until the ultimate tension strength, thereby excluding the effects of damage when fitting the hardening model in Section 4. Therefore, the true stress-strain curve in Fig. 7 keeps increasing. The commercial software Abaqus 6.14 was used to simulate the SRB tension test.
The geometric model is shown in Fig. 8(a). To simplify the calculation, only the part of the specimen within the gauge distance was simulated. The mapped mesh was created using an 8-node linear hexahedron element with a reduced integration size of 0.5 mm. The corrected curve of the equivalent stress and equivalent plastic strain (stress-strain relationship) was determined using the true curve in a multi-linear form by the trial and error method. It obtained the correct equivalent stress and equivalent plastic strain relationship ( Fig. 8(b)) by comparing the force-displacement curves from the experimental and simulation results until they matched, as shown in Fig.   8(a).

Plasticity model
The von Mises yield criterion is used to describe the flow criterion of a metal.
The strain hardening behavior of metals is described by the power hardening law: where A, B, and n are material constants obtained by fitting the experimental data in Fig. 8(b). In addition, due to the needs of the ductile fracture model in Section 5, the Swift hardening law was used to fit the experimental data as follows [28] : where K is a material constant, 0 is a prestrain-like material constant, and is the strain hardening exponent.
Through the results of the SRB tension test, the plasticity parameters required for the plastic model in Table 2 were obtained by fitting the data in Fig. 8(b). The fracture experiments were simulated by using these parameters. These FE simulations, described in the following, were done using Abaqus/Explicit with the power hardening law. All simulations were replicated in the same way using C3D8R 8-node bilinear axisymmetric quadrilateral and hourglass control. The critical elements were usually selected where the equivalent plastic strain is larger, and the damage is easy to occur. Based on previous findings [29,30] , the selection of critical elements is shown in Fig. 9. Since the real mesh size is very small, the grid division in Fig. 9 is used as a schematic diagram to clearly indicate the location of critical elements. After the simulation, the stress states of the critical elements selected on the model were extracted from the FE model to obtain the stress states of the specimens in other tests.  The force-displacement curves of SRB from simulations compared with the corresponding experiments are shown in Fig. 8(a). The force-displacement curve obtained from the experiment is consistent with the simulation results. In addition, necking can be observed in the simulation result. The variation of ST and LP with equivalent plastic strain extracted from the critical element is shown in Fig. 10, where LP remains -1 while ST increases. This indicates that the critical element is continuously elongated along the axial direction [18] .  The simulation of the compression test was performed as an axisymmetric case.
Punches and sets were regarded as a discrete rigid body. The force-displacement curves from simulations compared with experiments are shown in Fig. 12(a). The force-displacement curves obtained from the experiments agree well with the simulation results for both height-to-diameter ratios. The location of critical elements for obtaining stress state variables was in the middle of the axes, as shown in Fig. 9 (e-f). The variation of ST and LP with equivalent plastic strain is shown in Fig. 12(b), which indicates that ST remains at a specific amplitude and suppresses crack initiation. The torque-rad curve obtained from torsion test is shown in Fig. 13(a). The torque-rad curve obtained from the experiments is consistent with the simulation results. The location of critical elements was on the outer surface in the model, as shown in Fig. 9 (g). The variation of ST and LP with plastic strain is shown in Fig.   13(b). It was found that ST and LP approach zero, which leads to the failure of the specimen by shear stress.

Ductile fracture criterion
ST is dimensionless variable that describes the stress state of a point in a continuous medium by the ratio of stress invariants and is defined as follows: (5) where and ̅ are the mean stress and equivalent stress. Another important parameter is LP, which is related to the third stress invariant as follows: where 2 and 3 represent the second and third stress invariants.
In this paper, two phenomenological models were used to predict the ductile fracture of AA6063 in cold deformation. The damage prediction ability of two different ductile fracture criterion was verified for AA6063 under different loading conditions. The first model was the MMC criterion, which was derived from the Mohr-Coulomb criterion. The expression of the fracture strain as a function of stress state is as follows [13] : where 1 and 2 are constants that need to be calibrated by experiments, and are constants from the swift hardening law. Another fracture model was the BW model [13] . This function is based on three limiting cases: , , , which respectively correspond to the LP of -1, 0, and 1. From the early studies for these three limiting cases, the effect of ST on plastic strain can be written as Therefore, Eq. (12) can be written as [14] (12) It can be seen from Fig. 10, Fig. 11(b), Fig. 12(b), and Fig. 13(b) that ST and LP vary continuously during specimen deformation. To accurately calibrate the ductile fracture criterion, the integral mean values of ST and LP ̅ were calculated by Eq.
(13) and Eq. (14): The stress state and corresponding fracture strain used for calibration are shown in Table 3. Eq. (12) in MATLAB. The calibrated material constants are given in Table 3. The 3D fracture envelopes of the MMC and BW models are depicted in Fig. 14(a) and Fig. 14(b), respectively. It can be seen that the 3D fracture envelopes of both models agree well with the results of tension and torsion tests. When the ST of the two models is around -0.5, a large difference in fracture strain is observed due to the lack of valid experimental data for fitting. The fitting effects of different models were evaluated by 2 using Eq. (15): (15) where ̅ is the fracture strain from the experiment. The sum of the squared residuals, 2 , of the MMC and BW models are 0.8882 and 0.9901, respectively. It shows that the BW model fits better compared with the MMC model, may be related to the number of parameters in the model.  Since the LP hardly varies in axisymmetric tension and pure shear tests, the relationship between tension fracture strain or shear fracture strain and ST can be obtained when the LP is 1 or 0, respectively, in Fig. 15. It can be found that the fracture strain will decrease with the increase in ST under any deformation condition.
The shear fracture strain is smaller than the tension fracture strain as a whole. This indicates that cracks are more likely to occur during shear deformation. The experimental data points are close to the curve of the BW model.

Fig. 15
The relationship of fracture strain and ST.

Fracture prediction
The developed approaches were implemented into Abaqus 6.14 for FE simulation to verify the damage model further and observe damage distribution. Two developed approaches were applied to two selected tests (SRB, NR20), respectively.

Application of the MMC model
There are force-displacement responses from simulations compared with experiments in Fig. 16(a). The curve was obtained by FE simulation, and the scattered points show the experimental data. Better prediction in the weak stage was achieved by the MMC model compared with the results without considering damage effects in Fig. 11(a). However, compared with the experimental results, the fracture of the specimen in the simulation occurred earlier than that in the experiment. The crack of SRB initiated at the center of necking and propagated horizontally at first in Fig.   16(b). The crack of NR20 initiated in the middle of the thickness at the notch region and propagated horizontally.

Application of the Bai and Wierzbicki model
The results of the BW model were similar to those of the MMC model. The BW model was more accurate than the MMC model in predicting the fracture point on the force-displacement curve in Fig. 17(a). Through the analysis of the damage field shown in Fig. 17(b), it was found that the damaged area for the BW model was slightly smaller. The prediction ability of the two models for the damage of AA6063 was similar, although the accuracy of the BW model was slightly higher than that of the MMC model.

Conclusions
(1) The mechanisms of ductile fracture were investigated in various loading conditions. The initial ST of experiments varies from −0.33 at cylindrical compression tests to 0.67 at NR5 tension tests. The specimens with high ST tend to fail due to tension stress, while the specimens tend to fail along the direction of the maximum shear stress with the decrease in ST.
(2) The plasticity of AA6063 was probed through the SRB tension test. The force-displacement curve obtained from the experiment was consistent with the simulation results. (

Declarations
(1) This is an original article, and there is no multi contribution of one manuscript.
(2) The materials and data used in this paper are available, and other researchers can use them to repeat the experiments in this paper.
(3) All the authors have no objection to the author ranking. The cooperative units have no objection to the ranking of the units and there are no intellectual property disputes. If there is a change in the signature of the author and the unit in the process of revising the paper, the consent of all authors shall be required, and the certificate shall be issued by the first author or corresponding author.
(4) In case of any discrepancy with the statement, the author shall compensate the editorial department for the loss.