Comparison of Modified Mohr–Coulomb Model and Bai–Wierzbicki Model for Constructing 3D Ductile Fracture Envelope of AA6063

Ductile fracture of metal often occurs in the plastic forming process of parts. The establishment of ductile fracture criterion can effectively guide the selection of process parameters and avoid ductile fracture of parts during machining. The 3D ductile fracture envelope of AA6063-T6 was developed to predict and prevent its fracture. Smooth round bar tension tests were performed to characterize the flow stress, and a series of experiments were conducted to characterize the ductile fracture firstly, such as notched round bar 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-Wierzbicki (BW) model to develop the 3D ductile fracture envelope. Finally, two ductile damage models were proposed based on the 3D fracture envelope of AA6063. Through the comparison of the two models, it was found that BW model had better fitting effect, and the sum of squares of residual error of BW model was 0.9901. The two models had relatively large errors in predicting the fracture strain of SRB tensile test and torsion test, but both of the predicting error of both two models were within the acceptable range of 15%. In the process of finite element simulation, the evolution process of ductile fracture can be well simulated by the two models. However, BW model can predict the location of fracture more accurately than MMC model.


Introduction
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. Cold working is generally performed at room temperature and leads to an improvement in the mechanical properties and a decrease in plasticity due to work hardening. In fact, cold working can easily result in ductile damage, meaning it is critical to investigate the damage forming mechanism at room temperature [1,2]. It is the main means to study the damage forming mechanism of metal by establishing mathematical model [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 ductile damage is caused by the nucleation, growth, and aggregation of micro-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. Crockroft and Latham deemed that the maximum tension stress is the main factor in the failure of materials, and the Crockroft-Latham ductile failure criterion was established based on this [7]. Chen et al. [8] believed that the equivalent stress is also one of the factors that leads to material failure, and based on the work of Crockroft and Latham, a normalized Crockroft-Latham ductile fracture criterion was thus established. Meanwhile, in accordance with various applications, different ductile failure criterions have been established. For example, the Brozzo model [9] and the Oyane model [10] were used to predict sheet metal forming and round bar drawing. While the failure criterion is simple, it cannot predict the complex deformation path and large plastic deformation.
In terms of the second approach, Gurson first introduced the void volume fraction f into the Mises yield criterion and established the Gurson plastic potential [11]. Following this, Tvergaard et al. [12] established the modified Gurson plastic potential by considering the effects of the nonuniform stress field around the voids on the softening behavior of the material. Needlemen et al. [13] held a view that the nucleation of the void conforms to the normal distribution controlled by strain, and the nucleation of the void will lead to a change in the void volume fraction f. Through the research of above authors, the so-called GTN constitutive model which takes into account the nucleation, growth and convergence of the voids was established. Following further research on the void evolution mechanism, a number of more accurate GTN constitutive models have been developed. The advantage of this model is that it is more accurate for the damage prediction under a relatively low ST. However, 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. According to the continuum theory, Lemaitre coupled the equivalent stress σ with the damage factor D to describe the dissipation potential in the deformation process, and established a coupled phenomenological damage model [14]. Later, it was found that the application condition of Lemaitre's damage model was limited by the range of ST. Considering this, Xue established the relationship between the damage factor D and ST and LP, and revised the Xue phenomenological damage model [15]. He et al. coupled the damage variable D with dislocation density, grain size and recrystallization fraction, and established the damage constitutive equation of 52100 bearing steelballs [1].
Given that the coupled phenomenological damage model requires numerous experiments to calibrate the parameters, large number of researchers have established the uncoupled phenomenological damage model. According to the uncoupled phenomenological damage model, there is a linear relationship between the damage parameter D and plastic strain. The expression of damage parameter D and equivalent plastic fracture strain ε f is as follows: where η is ST, and θ is LP. Therefore, it is essential to construct a 3D fracture envelope in the space ε f , η, θ for uncoupled damage. Bai-Wierzbicki first constructed the 3D fracture envelope in the space ε f , η, θ , which defines the fracture strain as a function of ST and LP [16]. The Mohr-Coulomb model was formed based on the maximum shear stress to predict shear damage. Bai-Wierzbicki used the relationship between principal stress and ST and LP to transform the Mohr-Coulomb model into the MMC model [17]. Working on the basis of the MMC model, Lou et al. established the concept of ST threshold and proposed that when the ST is lower than the threshold, ductile fracture will not occur [18]. Elsewhere, Mohr et al. transformed the Mohr-Coulomb model into the space of ST, LP and equivalent plastic strain ε p devise the mixed stress/strain version and established the KHSP fracture criterion [19]. The phenomenological damage model parameters 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 [20]. The fracture of materials with positive ST can be obtained by using smooth round bar (SRB) tension test, notched round bar (NR) tension test, and plate tension test [21][22][23]. The fracture of materials with nearly zero ST can be obtained by the torsion and notched tube tension tests [24,25]. The fracture of materials with negative ST can be obtained by a cylindrical compression test.
The phenological damage models of aluminum alloys have 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 [26]. Many researchers used the phenomenological damage model to the cold deformation of aluminum alloys. Mirnia et al.  [28]. However, a 3D fracture envelope has not yet been constructed for AA6063 alloy.
In this paper, the fracture characteristics of AA6063 under different deformation conditions were studied, and two uncoupled phenomenological damage models based on MMC and BW ductile fracture criterion was established to describe the effect of ST and LP on plastic damage. By comparing the two types of phenological damage models, the prediction ability of the different ductile fracture criterion for AA6063 were assessed. To obtain a wide range of ST variation, tension, torsion, and compression tests were conducted to determine the relationships between the fracture strain and ST with 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 forcedisplacement 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.

Experimental Process
To observe the fracture of the specimen under as many different stress states as possible, various fracture tests were designed in the initial stress state space, as shown in Figure 1. The tests included tension, torsion, and compression tests. The plasticity and fracture behavior of the material under high ST were obtained by the SRB tension tests. The design of the SRB tension test piece was based on GB/T 228-2002, as shown in Figure 2(a). To   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. Figure 2 shows the shape and dimension of the sample of the room-temperature fracture test. The smallest cross-sectional diameter of the three types of NR was 5 mm, and R in Eq. (2) were 20 mm, 10 mm, and 5 mm, respectively. The dimension and shape of these specimens are shown in Figure 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. Through the torsion test, the specimen reached a pure shear state. The shape and dimension of the torsion specimen are shown in Figure 2(b). MTS809 axial torsional test system was used in this experiment. The test speed was 0.1 rad/s, and the gauge length was 50 mm. Figure 2(c) shows the compression test specimens. To study the effect of height-to-diameter ratio on the plasticity and ductile fracture of materials, two ratios were considered (1.5 and 1). Through the compression test of the cylinder, the damage evolution and repair under negative ST were studied. To reduce the friction between the cylindrical specimen and experimental equipment, the lubricant was used before the experiment. To maintain the quasi-static condition, a compression speed of 0.015 mm/s was used.

Experimental Results
Macroscopic fracture surfaces were investigated to study the fracture mechanism of specimens. Figure 3 displays the macroscopic fracture phenomena of different tension specimens. The macroscopic fracture surface of these specimens with high initial ST mainly presents the "cup cone" shape. 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 cup-cone fracture is caused by the dominant effect of the tension stress in the center region of the fractured surface [29]. Figure 4 shows the force-displacement relationships of SRB, NR5, NR10, and NR20. The softening stage of the material can be seen from the image. The SRB can be regarded as a notched bar with an infinite notch radius. It can be seen from Figure 4 that the fracture elongation increases with the increase in notch radius, while the peak load increases with the decrease in notch radius. The load of the curves decreases obviously after necking. Most experimental results show that the ST increases evidently when necking occurs [23,24]. The increase in ST will affect the plasticity of materials. This is because of the fact that the regions with high triaxiality tend to have the characteristics of small plastic deformation and large volume deformation, which will lead to the release of more elastic strain energy at the points with high ST. When more elastic strain energy is released, the stress will be concentrated at this point, which will hinder the metal flow and will lead to a decrease in metal plasticity [30]. Therefore,  as Figure 4 shows, the NR5 with the largest initial ST was the first to enter the softening stage and to fracture, while the slowest to enter this stage and to fracture was the SRB with the minimum initial ST. Figure 5 shows the force-displacement relationship of specimens after cylindrical compression tests with different height-to-diameter ratios. It can be seen from Figure 5 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 force-displacement curve continuously increases. Figure 6(a) shows the torque-twist angle curve obtained from the torsion test under a rate 0.1 rad/s, and Figure 6(b) shows the morphology of the fracture surface of the specimen. It can be observed that the fracture surface is very smooth without any macroscopic cracks due to the shear stress acting along the cross section in the pure shear state. Compared with the tension fracture test with high ST, the specimens fracture along the direction was perpendicular to the maximum tension stress, while in the torsion fracture test with low ST, the specimens fracture along the direction of maximum shear stress.

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 Figure 7 was calculated according   to the SRB force-displacement curve in Figure 4. According to the stress-strain curve in Figure 7, Young's modulus of elasticity was 14611 MPa. The true stress-strain curve in Figure 7 was calculated until the ultimate tension strength, thereby excluding the effects of damage when fitting the hardening model in Section 4. In the true stress-strain relationship, the reduction of the cross section of the specimen under tension is considered. As such, the load on the unit cross section of the specimen increased, which led to the true stress-strain curve being higher than the engineering stress-strain curve. The commercial software Abaqus 6.14 was used to simulate the SRB tension test. The geometric model is shown in Figure 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 (Figure 8(b)) by comparing the force-displacement curves from the experimental and simulation results until they matched, as shown in Figure 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 Figure 8(b). Additionally, due to the needs of the ductile fracture model in Section 5, the Swift hardening law was used to fit the experimental data as following [31]: where K is a material constant, ε 0 is a prestrain-like material constant, and m 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 Figure 8(b). The fracture experiments were simulated 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 [32,33], the critical element of specimen was usually selected on the element with the largest equivalent plastic strain before fracture, and the selection of critical  elements is shown in Figure 9. Since the real mesh size is very small, the grid division in Figure 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 Figure 8(a). The force-displacement curve obtained from the experiment is consistent with the simulation results. Additionally, 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 Figure 10, where LP remains 1 while ST increases. This indicates that the critical element is continuously elongated along the axial direction [23].
The force-displacement curves of NR from simulations compared with the corresponding experiments are shown in Figure 11(a). The simulation results of the softening stage are not consistent with the experimental data because the plastic model considers only the hardening law. The location of critical elements for obtaining stress state variables was at the axes in the notched region where the crack occurred, as shown in Figure 9(b-d).
The variation of ST and LP with equivalent plastic strain is shown in Figure 11(b). It can be seen that the growth rates of the ST of specimens with different notched radii are similar.
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 Figure 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 Figure 9(e)-(f ). The variation of ST and LP with equivalent plastic strain is shown in Figure 12(b).
The torque-rad curve obtained from torsion test is shown in Figure 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 Figure 9(g). The variation of ST and LP with plastic strain is shown in Figure 13(b). It was found that ST and LP approach zero, which leads to the failure of the specimen by shear stress.

Figure 9
The finite element mesh with critical elements: a SRB, b NR20, c NR10, d NR5, e cylinder with a height-to-diameter ratio of 1.5, f cylinder with a height-to-diameter ratio of 1.0, and g torsion test

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: where σ m and σ are the mean stress and equivalent stress. Another important parameter is LP, which is related to the third stress invariant as follows: where J 2 and J 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 MMC criterion was derived from the Mohr-Coulomb criterion. The expression of the fracture strain as a function of stress state is as follows [23]: where c 1 and c 2 are constants that need to be calibrated by experiments, K and m are constants from the Swift Therefore, Eq. (12) can be written as [19]: where D 1 , D 2 , D 3 , D 4 , D 5 and D 6 are constants, which need to be calibrated by experiments. It can be seen from Figures 10, 11(b), 12(b), and 13(b) that ST and LP vary continuously during specimen deformation (Table 3). To accurately calibrate the ductile fracture criterion, the integral mean values of ST (η) and LP ( θ ) were calculated by Eq. (13) and Eq. 14: GA optimization technique was used to calculate the parameters of Eqs. (8) and (12) in MATLAB. The calculated model parameters are shown in Table 4. The 3D fracture envelope of the MMC and BW models are depicted in Figure 14(a) and (b), respectively. It can be seen that the 3D fracture envelope of both models agree well with the results of tension and torsion tests.
The fitting effects of different models were evaluated by R 2 using Eq. 15: To compare the prediction ability of the two models, the prediction error of the two models were calculated according to Eq. 16: Figure 15 shows the comparison of the prediction errors between the two models, with the results indicating that, on the whole, the prediction ability of the two models is similar. While the two models returned several errors in the SRB tensile test and torsion tests, there were fewer errors with the other ductile fracture tests. This may because of the serious necking of the specimen in the SRB test, which led to the inaccurate measurement of the fracture strain. Secondly, the ductile fracture strain under shear conditions is less which results in the inaccurate prediction ability of the models for specimens with (16)   low and medium ST. The prediction ability of the two models was similar for the fracture strain of the NR tensile test. However, the BW ductile fracture criterion was superior to that of the MMC in predicting the fracture strain in the SRB tensile test and the torsion tests.
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 Figure 16. It can be found that the fracture strain will decrease with the increase in ST under any deformation condition. According to Bai-Wierzbicki [16,17], when the ST is low η < −0.33 , the material is in a negative hydrostatic stress state. In this stress state, the nucleation and growth of the voids will be restrained, which makes the occurrence of ductile fracture very unlikely. However, at a higher ST η > 0.4 , the material is in a high-ST stress state, wherein the fracture of the material is caused by the nucleation, growth, and propagation of the voids. In this stress state, the material is highly prone to form the ductile fracture. Therefore, under any deformation condition, the fracture strain will decrease with the increase in ST. 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.

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 simulate two selected tests (SRB, NR20), respectively.

Application of the MMC Model
There are force-displacement responses from simulations compared with experiments in Figure 17(a). The solid 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 Figure 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 Figure 17(b). The crack of NR20 initiated in the middle of the thickness at the notched region and propagated horizontally.

Application of the BW 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 Figure 18(a). Through the analysis of the damage field shown in Figure 18(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 failure due to tension stress, while the specimens tend to failure along the direction of the maximum shear stress with the decrease in ST.
(2) The plastic model of AA6063 was established via a SRB tensile test. By combining the ductile fracture test with the FE simulation, the ductile fracture criterion of AA6063 were established for both the BW model and the MMC model. (3) The 3D fracture envelope was represented by two models: the MMC model and the BW model. Compared with the experimental results, the two models can describe the softening phenomenon of materials due to the initiation and propagation of damage. The fitting effect and the prediction ability of the BW model were slightly better than those of the MMC model. While the fracture location predicted by the two models was different than the actual situation, the accuracy of the BW model was higher than that of the MMC model.