Numerical evaluation of the ductile fracture for AA6101-T4 and AISI 4340 alloys using the Lemaitre and Gurson models

This paper provides a numerical assessment of the ductile fracture of the AA6101 and AISI 4340 alloys using the Lemaitre and Gurson models.The simulations were carried out to verify the evolution of the accumulated plastic strain at fracture, assuming the high range of the triaxiality ratio, applying monotonic tensile loads on smooth and notched cylindrical specimens. The numerical strategy was structured for the Lemaitre and Gurson models, from establishing a non-linear system of equations using an implicit integration algorithm and solving the non-linear equation system using the Newton-Raphson method. When assessing the numerical and experimental results, it was observed that the cumulative plastic strain at the fracture decreases with the increasing levels of the triaxiality ratio for both alloys and models. On the one hand, the Lemaitre model was more optimistic than the experimental results. On the other hand, Gurson’s model proved to be more conservative in its prediction of ductile fracture. Regarding determining the fracture onset, in general, both models showed good predictive capacity. However, the numerical results for the aluminum alloy presented more in agreement with experimental data than for the AISI 4340 alloy. In the end, assuming the determination of the level of displacement at fracture, the Gurson’s model has the best performance. In this sense, for a high level of triaxiality ratio, the Gurson porous material can be recommended to describe the mechanical behavior of the material at fracture.


INTRODUCTION
Using lighter and more reliable metallic materials in different engineering areas has become necessary in recent decades. Consequently, the search for solutions that tend to optimize mechanical components is intrinsically related to the structuring of constitutive models capable of describing, with a certain precision, the moment, region, and exact location of the onset of ductile fracture of the constituent materials. The search for solutions implies that the need to design devices or components with knowledge of the state of maximum stress or maximum level of degradation must be identified for a correct design and desired operational safety. This paper considers two alloys with different ductility levels, which are widely used in numerous areas of engineering. Due to its low ductility, the AA6101-T4 aluminum alloy tends to fracture already in the forming process. Studies related to the fracture of this material have become relevant for applications in the design of lighter structures and mechanical components [1]. Furthermore, the annealed AISI 4340 alloy is considered a high-strength material, widely used in the offshore industry. It presents high ductility levels, a characteristic that makes the description of the mechanical behavior of this material highly relevant [2].
In this sense, analyzes regarding the mechanical behavior of materials have been carried out since the last century and date back to the studies by von Mises [3]. With the work of Kachanov [4], Damage Mechanics emerged, which studies the behavior of materials, admitting that they have damage or defects in their microstructure. When subjected to efforts exceeding their elastic limit, they will result in permanent deformations, also called plastic strains. In general, Damage Mechanics has as its primary objective to study and develop mathematical models that can predict the useful life of a given material or component [5]. However, it can be analyzed under two different approaches: i. One based on the continuous damage mechanics, where the Lemaitre model [5] is the main representative. In this approach, the damage depends on the accumulated plastic strain, and the evolution law is based on the thermodynamics of solids. Then, the so-called critical damage criterion is assumed, and a crack in the mesoscopic scale starts when the damage value reaches a critical level. In this sense, several models have been proposed [6][7][8][9][10][11][12][13][14][15][16][17][18][19]; ii. The second one is based on micromechanics of defects, where Gurson's model is the main contribution [20]. The approach assumes the material's porosity as the damage variable and is one of the most important reference in the area. Later, other formulations were published to improve the original proposal.
Given the importance of Lemaitre's approaches [5] and Gurson's [20] models, this contribution proposes to numerically study its predictive capacity regarding the determination of the level of displacement at fracture and the fracture onset. For this, two materials, such as aluminum AA6101-T4 and annealed AISI 4340 alloys, are used to manufacture cylindrical specimens in a high of triaxiality ratio level.

Mathematical Modeling and Numerical Strategy
This section presents the mathematical formulations for Lemaitre's and Gurson's models and their corresponding numerical strategies for integrating the flow equations and obtaining the return algorithms. Both damage models are built based on an isotropic damage variable, considering isotropic hardening for the studied alloys.

Lemaitre's Model
The mathematical structuring of the constitutive models, and the description of the variables involved in the study, are of fundamental importance for understanding the research approach. Initially, the mathematical aspects of the model based on the Continuous Damage Mechanics proposed by Lemaitre are presented [5], and later the mathematical aspects of the elastoplastic model published by Gao [40].
In the construction of the Lemaitre damage mathematical model [5], there is a need to define the effective stress according to equation (1).
where, σ ef is the effective stress of Lemaitre [5], σ is the Cauchy stress tensor and D is the isotropic damage variable. Thus, substituting equation (1) in the generalized Hooke's Law, the modified Hooke's Law (equation 2) is obtained.
The expression above (equation 3) can be written as follows: Where, D e is the elasticity material tensor and ε e is the elastic strain tensor. The yield function (equation 4) for the Lemaitre model [5] can be written as: Where the term J 2 (S) is the second invariant of the deviatoric stress tensor, σ y0 represents the initial yield stress of the material, ε p represents the accumulated plastic strain, and H I is the isotropic hardening modulus, which is a function of the accumulated plastic strain, being represented by H I ε p ( ) .
The flow vector (see equation 5), calculated based on the associative plasticity can be written as: Where N represents the flow vector and q is the von Mises equivalent stress. Thus, the plastic flow (equation 6) rule is then determined: Where ! ε p represents the plastic strain rate tensor and ! γ is the so-called plastic multiplier rate. The rate of evolution of the accumulated plastic strain, which assumes the role of an internal variable of isotropic hardening, is calculated according to the Prandt-Reuss, equation (7).
Where ! ε p represents the rate of evolution of the accumulated plastic strain. The law of evolution of the damage variable is given by equation (8).
Where ! D is the rate of evolution of the isotropic damage variable, Y represents the energy released due to the damage, S is the damage denominator, and s is the damage exponent. Terms (S, s) are material parameters and need to be calibrated. Thus, Y can be mathematically calculated by equation (9).
Where p is the hydrostatic pressure, G and K represent, respectively, the shear and bulk modulus. Thus, Figure 1 contains the mathematical equations that define the Lemaitre model [5] with isotropic hardening and isotropic damage.

Return map algorithm for the Lemaitre model
Lemaitre's model [5] is integrated through an implicit Euler strategy [41], where the operator split methodology is used, assuming that the problem can be split into two parts: a) trial state, in which the prescribed strain increment is considered entirely elastic, and b) plastic corrector, in which the hypothesis of full elastic strain increment is not valid. Moreover, in the plastic corrector, a system of non-linear equations is solved, using the Newton-Raphson method. The model is then added to an academic finite element tool called Hyplas, using the large strain transformation algorithm proposed by De Sousa Neto [42]. The system of nonlinear equations to be solved is: where, R σ n+1 , R D n+1 and R Δγ are residual equation, defined by the evolution law of the internal variables. The solution of the nonlinear system of equations shown in equation (10) is performed using the Newton-Raphson method, presented in Figure 2.

Gurson's mathematical model
Gurson's model [20] is based on the micromechanics of defects, where the porosity is taken as the damage variable. The yield function by Gurson is a function of the second invariant of the deviatoric stress tensor, porosity, hydrostatic stress, and hardening curve. Besides that, the evolution of e trial + Δγ D e : N n+1 the porosity is a function of the volumetric part of the plastic strain tensor, which represents the growth of the defect mechanisms. Among them, the Gurson model stands out [20], and other phenomenological extensions are proposed in the literature [21,43,44,45]. For the construction of Gurson's mathematical model [20], Hooke's Law in its original and generalized form is considered ( see equation 11).
Where is the Cauchy stress tensor, De is the elasticity material tensor, and ε e is the elastic strain tensor. The yield function for the Gurson model is defined according to equation 12: Where the term J 2 (S) is the second invariant of the deviatoric stress tensor, represents the material's yield stress or hardening rule, p is the hydrostatic pressure, and f is the volume void fraction or porosity.
The flow vector, calculated based on the associative plasticity can be written as shown in equation 13: N d and N v are the deviatoric and volumetric contributions of the flow vector, S is the deviatoric stress tensor, and I is the second-order identity tensor. Thus, the plastic flow rule is then determined by equation 14: where ! ε d p and ! ε v p represent the deviatoric and volumetric parts of the plastic strain rate tensor, respectively, the rate of accumulated plastic strain, which assumes the role of an internal variable of isotropic hardening, is also calculated according to the Prandt-Reuss equation (see equation 15).
Where ! ε p is the rate of accumulated plastic strain. The evolution law of porosity or void volume fraction, according to Gurson [20], is defined in equation 16:

Return algorithm for the Gurson model
In the construction of the return algorithm, the following system of non-linear equations must be solved (equation 17), having as variables: σ n+1 ,ε n+1 p , f n+1 , and Δγ Where, R σ n+1 , R ε n+1 p , R f n+1 and R Δγ are residual equations based on the evolution equations of the internal variables. In solving the non-linear system of equations, the Newton-Raphson method is applied; Figure 4 presents the Newton-Raphson strategy.

METHODOLOGY OF TESTS
In this part, it is presented the methodology of the experimental tests and the description of the specimens in the high level of triaxiality ratio.

Material description
AISI 4340 and AA6101-T4 alloy test specimens were manufactured to evaluate the Gurson and Lemaitre models in a high triaxiality ratio, considering smooth and notched cylindrical bars. A total of 7 specimens were manufactured: 4 for annealed AISI 4340 alloy and 3 for AA6101-T4 alloy. For the AA6101-T4 alloy, the experimental data were taken from Malcher [46], assuming a smooth cylindrical specimen and two notched cylindrical specimens with a notch radius of 6 and 10 mm. Besides that, for the AISI 4340 alloy, the experimental data were taken from Morales [2], assuming a smooth cylindrical specimen and three cylindrical notched specimens with a notch radii of 4, 6, and 10 mm. All specimens were submitted to a tensile test under strain control with a strain gauge of 25 mm of length gauge. The tests were performed regarding an MTS 810 of 100 kN of maximum force.

Geometry of specimens
The triaxiality ratio is conceptually defined as the ratio between the hydrostatic pressure, p, and the equivalent von Mises stress, q, represented in equation (18).
Bridgman [45], related the triaxiality ratio η, with the constructive parameters used to make the specimens: notch radius R, and radius of the central section of the specimen , according to equation (19).
Based on Bridgman's equation [45], cylindrical specimens were prepared, designed to allow an initial triaxiality ratio in the critical region of 0.33(smooth), 0.50, 0.60, and 0.70 (notched), with a cylindrical of the entire cross-section. The smooth and notched cylindrical specimens were machined with a total length of 128 mm and 120 mm, respectively, and a useful length of 40 mm. In Figure 5, the representation of the critical region of the cylindrical specimens is schematically shown.
In Figure 6. The technical drawings of the manufactured test specimens are shown.

NUMERICAL AND EXPERIMENTAL RESULTS
This section shows the numerical and experimental results for the specimens manufactured for the two studied alloys. Initially, the characteristics  of the discretization of the geometries and the corresponding calibration to obtain the parameters of each material are presented. The force reaction curves are shown, with a comparison between the experimental and numerical values, contrasted with the damage evolution. Finally, the evolution of the accumulated plastic strain is also described.

DISCRETIZATION AND CALIBRATION OF MATERIAL PARAMETERS
The finite element meshes were defined from axisymmetric modeling and structured with quadrilateral elements of eight nodes using reduced integration. Figure 7 shows the finite element meshes used in the numerical simulations, each mesh formed by 2146 nodes and 675 elements, with refinement in the central part to allow greater discretization of the results in this region of interest.
For both materials, the smooth cylindrical specimens, under tensile load, resulting in the experimental reaction and optimized curve, through the multivariable search method based on the gradient proposed by Machado and Malcher [46] and Machado [47], for calibration of material properties. In this case, the Kleinermann and Ponthot [48] curve parameters were also calibrated, equation (20), used to represent the behavior of the isotropic hardening curve of materials through four parameters: σ y0 , σ ∞ , ξ, δ. Figure 8 and Figure 9 show the graphic result of the parametric identification process, where it can be seen that the numerically calculated reaction force versus displacement curve approximates the experimentally determined curve.  Figure 10 shows the hardening curves determined according to Lemaitre [5] and Gurson [20] models. A difference between the curves is observed based on the nature of the damage evolution equations proposed by both models. For Lemaitre [5], the damage is a function of the accumulated plastic deformation and the energy     [20], the damage is associated with the ratio between the number of voids and the representative total volume, deduced from the micromechanics of a spherical defect.

Reaction curves and damage evolution
This item presents the numerical results for the notched specimens with R = 10 mm and R = 6 mm manufactured with the AA6101-T4 alloy and the results of the notched specimens manufactured with annealed AISI 4340 alloy with R = 10 mm, R = 6 mm, and R = 4 mm.
It is noteworthy that in a region of high triaxiality η ≥ 1 3 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟, the increase in specimen notch severity and the radius reduction result in an increase in the triaxiality ratio. The initial triaxiality ratio values are shown in Table 2.
Based on the parameters of the materials obtained previously, for the AA6101-T4 alloy, the notched specimens were simulated, subjected to tensile force and the damage evolution laws proposed by Lemaitre [5] and Gurson [20] to assess the predictive capacity of the models.   Figure 11 shows the reaction curve, damage evolution at the critical node, and damage contour for the cylindrical notched specimen with R = 10 mm (AA6101-T4). The simulations were carried out until the models reached the critical damage values. Thus, the experimentally observed displacements are compared with the numerically calculated ones, the location of the maximum damage value (critical node) along the boundary of the finite element mesh.
It can be seen that the Lemaitre model [5] predicted the displacement required for fracture onset, more significant than the experimentally measured value, which shows that the model had an optimistic behavior when the ductile fracture began. Gurson's model [20] predicted the numerically calculated displacement for fracture onset, lower than the experimentally measured value, which shows that the model had a conservative predictive capacity for fracture onset. As for the location of the critical node, that is, the node with the maximum damage, both models indicated that the beginning of the fracture occurs in the center of the specimen, which is in agreement with the experimental observations. Figure 12 shows the behavior for the notched specimen with R = 6 mm (AA6101-T4). Again, the Gurson model [20] presented a behavior closer to the real one, given the comparison between the numerically calculated and the experimentally observed displacements. In this case, also, both models could indicate the center of the specimen as the starting point of the fracture, compatible with the experimental observations.
It is also possible to highlight the difference between the reaction force levels calculated by the models and the levels observed experimentally. Both models calculated reaction force levels smaller than the experimental levels. This behavior is typical of models based on the von Mises yield function, Figure 10. Isotropic hardening curves for AA6101-T4 and AISI 4340 alloys, regarding Lemaitre and Gurson models.  where the triaxiality ratio or the hydrostatic pressure effect is neglected in determining the material's plastic flow law.
For the annealed AISI 4340 alloy, the notched specimens with radii of 10, 6, and 4 mm were simulated. For the annealed AISI 4340 alloy, the notched specimens with 10, 6, and 4 mm radii were simulated. Figure 13 shows the results for the specimen with R = 10 mm. Better performance of Gurson's model [20] concerning the prediction of the correct displacement for the fracture can again be observed. Lemaitre's model [5] was extremely optimistic, calculating a displacement for the fracture of about 30% greater than that observed through the experimental data. The location of the critical node, for both models, approached the starting point of the ductile fracture, experimentally verified. Figure 14 presents the results for the notched specimens with R = 6 mm. The predictive ability of the Gurson model [20] was again more significant than the predictive capacity of the Lemaitre model [5].
It is important to highlight that for both simulations, the calculated reaction force levels were very close to those experimentally measured, which differs from the result obtained for the AA6101-T4 alloy. This fact can be explained by the level of ductility of the materials, where the models based on Mises have better behavior regarding the level of strength in less ductile materials. Observing the level of accumulated plastic strain at fracture calculated by Lemaitre and Gurson's models, it can be seen that the AA6101-T4 alloy has a higher level of ductility than the AISI 4340 alloy (see Table 1).
For the notched specimen with R = 4 mm (Figure 15), it is observed that the Lemaitre model [5] behaves very optimistically, and the Gurson model [20] behaves conservatively in determining the displacement necessary for the beginning of the fracture. Gurson [20] reached the critical damage for a displacement at fracture less than the experimental observation, and Lemaitre [5] had the opposite behavior.

Evolution of accumulated plastic strain
This section shows the numerical results obtained for the accumulated plastic strain at fracture. The fracture initiation is set for the same displacement when the critical damage is reached. Figure 16, Figure 17, and Figure 18 show the behavior of the accumulated plastic strain at fracture, for the smooth and notched specimens, according to Lemaitre [5] and Gurson [20], for AA6101-T4 alloy. It can be observed that both models reach levels close to accumulated plastic strain at fracture, however, at different displacements, which shows that the rate of evolution of plastic strain for Gurson [20] is always higher than for Lemaitre. Figure 19, Figure 20, Figure 21, and Figure 22 also show the behavior of the accumulated plastic strain at fracture, for the smooth and notched specimens, according to Lemaitre and Gurson, now for AISI 4340 alloy. It can be observed that both models reach levels close to accumulated plastic strain at fracture, however, at different displacements, which shows that the rate of evolution of plastic strain for Gurson [20] is always higher than for Lemaitre.
Remarkably, the results were obtained for different displacement levels, which may indicate that even for materials with lower ductility, the evolution rate of plastic strain evolution, according to Gurson [20], is always greater than for materials with lower      Figure 19. Evolution of the accumulated plastic strain for AISI4340, regarding the smooth specimen.
ductility Lemaitre [5]. Figure 23 shows the behavior of the accumulated plastic strain at fracture, which can be a function of the triaxiality ratio. In this sense, for η ≥ 1/3, the accumulated plastic strain at fracture decreases with the increase of the triaxiality ratio. This behavior has been observed for both AISI 4340 and AA6101-T4 alloys. Nevertheless, due to ductility, the AA6101-T4 reaches levels of accumulated plastic strain higher than the AISI 4340. Table 3 presents the displacement at fracture determined experimentally and numerically by Lemaitre and Gurson.

CONCLUSIONS
It can be observed that both Lemaitre and Gurson damage models are able to capture ductile failure of AA6101-T4 and AISI 4340 alloys in high levels of triaxiality ratio. However, based on the level of displacement at fracture, the Gurson model had the best performance. The numerical results by Gurson were less than 18% of the experimental displacement observations. The performance of the Lemaitre model was extremally optimist, assuming values around 64% far from the experimental Figure 20. Evolution of the accumulated plastic strain for AISI4340, regarding the notched specimen with R=10 mm.  data. Regarding identifying the correct fracture locus, both models had a satisfactory performance, showing the center of the specimens as the site of fracture initiation, according to experimental observations. Furthermore, analyzing the localization of the maximum level of the accumulated plastic strain for the AA6101-T4, both models perform well and behave similarly. However, for the AISI 4340, the Gurson model had the best performance, remarking on the influence of the level of ductility in the predictive performance of the models. In this case, for Lemaitre, the accumulated plastic strain was spraided around the region in the center of the specimens; meanwhile, for Gurson, the variable was concentrated at the geometrical center of the specimens. Figure 23. Accumulated plastic strain at fracture as a function of the initial triaxiality ratio, for both, AA6101-T4 and AISI4340 alloys.