Visualization of the damage evolution for Ti–3Al–2Mo–2Zr alloy during a uniaxial tensile process using a microvoids proliferation damage model


 Understanding the damage evolution of alloys during a plastic deformation process is significant to the structural design of components and accident prevention. In order to visualize the damage evolution in the plastic deformation of Ti–3Al–2Mo–2Zr alloy, a series of uniaxial tensile experiments for this alloy were carried out under the strain rates of 0.1–10 s−1 at room temperature, and the stress–strain curves were achieved. On the other hand, the finite element (FE) models of these uniaxial tensile processes were established. A microvoids proliferation model, Gurson–Tvergaard–Needleman (GTN) damage model, was implanted into the uniaxial tensile models, and the simulated stress–strain curves corresponding to different GTN parameter combinations were obtained. Based on the simulated and experimental stress–strain curves, the GTN parameters of this alloy were solved by response surface methodology (RSM). The solved GTN parameters suggest that higher strain rate can enhance the proliferation and coalescence of microvoids. Furthermore, the uniaxial tensile tests over different strain rates were simulated using the solved GTN parameters. Then, the damage processes were visualized and evaluated. The result shows that the degradation speed of this alloy is slow at the initial stage of the tensile deformation and then accelerates once the voids volume fraction reaches a critical value.


Introduction
Ti-3Al-2Mo-2Zr alloy is a near α titanium alloy with bits of β phase distributed on α matrix. It has wide application prospects in marine engineering, ship and chemical industry owing to its excellent combination properties, such as strong toughness, high fatigue resistance, middle strength and superior corrosion resistance [1,2]. This alloy is a prime candidate for Ti-6Al-4V alloy, mainly used in the skeleton parts and pressure vessels of underwater detectors [2]. When the detectors are subjected to violent shocks, these components will undergo destructive deformation to absorb the impact energy. Because of the high toughness and superior strength, this alloy is able to withstand huge deformation without breaking. In fact, the destructive deformation process is composed of three portions including elastic deformation, plastic deformation and damage degradation. Currently, the mechanical properties and the plastic deformation behaviors of this alloy have been reported [1,2]. However, little attention has been paid on the damage degradation behaviors of this alloy. It has been widely recognized that the damage degradation behavior of metals is of significance to the structural design and security evaluation of the components. Nevertheless, it is usually difficult to characterize the damage degradation process of metals because damage is essentially a nucleation and proliferation process of microvoids which cannot be detected directly [3]. Therefore, it is an important issue to find a suitable method to characterize the damage degradation behaviors of metals and to further understand the damage mechanism of this alloy.
Due to the fact that the essence of damage degradation is resulted from the nucleation and proliferation of microvoids, it is reasonable to characterize the degradation process of an alloy using the indicator of microvoids volume fraction. In order to achieve the evolution of this indicator during a certain deforming process, a reliable damage model is essential. Up till now, there have been a lot of damage models, such as coupled phenomenological models, uncoupled phenomenological models and micromechanical models [4]. The coupled or uncoupled phenomenological models have been increasingly applied in industry due to the advantage of computation speed. However, the coefficient solution in the phenomenological models is dependent on the specified experimental conditions, and so such models cannot intrinsically describe the damage degradation behavior of materials [4]. It also means that the phenomenological models have a relatively higher tolerance in the computation of damage evolution. Nevertheless, a micromechanical model can solve this confused issue [5]. On the basis of the micromechanical theory, Gurson-Tvergaard-Needleman (GTN) model was developed [6], and until now it has been widely applied to describe the damage behaviors of materials, such as in aluminum alloy [7], titanium alloy [8], shape memory alloy [9], high strength steel [7], etc. Once the GTN model was figured out and planted into the finite element (FE) software, the volume fraction of microvoids can be derived in real time during a deforming process. Generally speaking, the damage degradation evolution of materials can be characterized by the FE simulation coupled with GTN model. In GTN model, there are four vital parameters that are closely related to the evolution mechanism of microvoids. According to the intrinsic concept of GTN model, the evolution of microvoids directly affects the yield behavior of materials, meanwhile the latter affects the former backward. That is to say, there exists a mapping relationship between the macroscopic yield behavior and the four GTN parameters. How to describe the mapping relationship is the key to solve the GTN parameters.
Response surface methodology (RSM) is suitable to clarify the chaotic mapping relationship. RSM seeks to identify the interior relationships between response variables and control factors through concise polynomial equations [10]. Compared with conventional fitting methods like orthogonal experimental method and polynomial fitting, RSM can achieve better fitting effect based on less experimental results. On the other hand, compared to intellectual algorithms, the fitting result of RSM is more intuitive and the fitted model is more convenient to be implanted into FE models. In addition, unlike the traditional regression method whose fitting effect significantly relies on the complexity of the independent variables, RSM can handle the complex function relationship even though sometimes there is more than one independent variable and these variables influence each other. More importantly, the higher fitting precision can be achieved because the non-significant solution set will be deleted during the regression calculation of RSM [11]. Naturally, RSM will be adopted in this work to solve the GTN parameters of Ti-3Al-2Mo-2Zr alloy.
In this work, a series of uniaxial tensile experiments were first performed to obtain the real stress-strain data. Then, the FE models of the uniaxial tensile processes were established and the simulated stress-strain curves corresponding to different GTN parameter combinations were obtained. Based on the experimental stress-strain data and the simulated ones, the GTN parameters of Ti-3Al-2Mo-2Zr alloy were solved by RSM. As the uniaxial tensile FE models were simulated again using the solved GTN parameters, the damage degradation process of this alloy was characterized dynamically. Besides, the influence of strain rate on the damage degradation behavior was discussed in this work.
2 Methods and theoretical basis 2.1 Theoretical basis of GTN damage model GTN model, as one of the classic microscopic mechanical models, has been widely used to describe the effect of ductile damage on the plastic yielding of materials [12]. The ductile damage process that contains the microvoids evolution can be schematically described as shown in Figure 1. When the material is deformed under a tensile stress, microvoids existing in the material begin to selfproliferate and expand. As the deformation continues, some microvoids begin to merge with each other, forming microcracks. Then, macrocracks come into being with the propagation and coalescence of the microcracks, resulting in dramatic performance degradation.
In the theory of GTN model, the impact of microvoids on the plastic yielding potential Φ of materials can be expressed as equation (1) [13].
where σ eq is the macro von Mises equivalent stress; σ h is the macro hydrostatic stress; σ m is the equivalent stress of matrix material; q 1 , q 2 and q 3 are the calibration parameters of the model, and q 3 = q 1 2 . Here q 1 and q 2 are 1.5 and 1, respectively, for this alloy [14]; f * is the damage factor which can be expressed as the function of the voids volume fraction (VVF) as equation (2).
The VVF includes the contribution of newly nucleated microvoids and the growth of existing microvoids, and their relationships can be expressed as equation (3).
where df represents the total increment of VVF; df nucleation represents the increment of VVF that comes from newly nucleated microvoids; df growth represents the increment of VVF that comes from the growth of existing microvoids. ε d¯m pl is the microplastic strain increment of the matrix material; ε d kk p is the macroplastic strain increment. A is the nucleation coefficient of microvoids and it can be expressed as equation (4).
where f nucleation represents the VVF of the newly nucleated microvoids; ε p represents the equivalent plastic strain of matrix material; S N represents the standard deviation and ε N represents the equivalent plastic strain for the newly nucleated microvoids. According to refs [13,15,16], the values of ε N and S N can be taken as 0.3 and 0.1, respectively. Based on the theoretical description of GTN damage model, the control parameters of GTN model include q 1 , q 2 , q 3 , f 0 , f c , f f , f n , ε N and S N . Since some of the parameters have been obtained from references, the GTN parameters to be solved are f 0 , f c , f f , f n .

Theoretical basis of RSM
RSM, as a regression analysis method, was introduced by George E. P. Box and K. B. Wilson in 1951 [17]. It is dedicated to characterizing the quantified connection between the response variables and the independent variables. The main operation procedure of RSM is as follows. First, design a sequence of experiments to get the basic data for establishing response surface. Second, choose suitable regression equations to represent the response relationship. This relationship can be expressed as equation (5), where y represents the response variables and x i (i = 1, 2,…,n) represent the factors.
Finally, fit the regression equation using the obtained experimental data by least square method, working out the coefficients in the regression equation. Usually, for the sake of saving calculation time and avoiding over-fitting, linear or quadratic polynomial functions are widely favored for the regression equations. Quadratic polynomial functions are adopted in the present work. The function of the response surface can be expressed as follows: where R represents the response variables; b 0 ,b i ,b ii and b ij are the coefficients of response functions; X i and X j represent the independent variables.

Method for solving the GTN parameters based on RSM
One of the aims of this work is to figure out the values of GTN parameters. Hence, the GTN parameters can be regarded as the independent variables of the response surface. As introduced in the theory of GTN model, the GTN parameters determine the yielding behavior of alloys directly. Therefore, the indicators related to stress can be considered as the response variables. Since the stressstrain curve for a material gives the direct relationship between stress and deformation, the difference between experimental stress-strain curves and numerical ones can be defined as the response variables of RSM. The detailed expressions of the response values R are as follows: where ε E P and ε N P represent the experimental peak strain and numerical one, respectively; σ E P and σ N P represent the experimental peak stress and numerical one, respectively; ε E F and ε N F represent the experimental breaking strain and numerical one, respectively; σ E F and σ N F represent the experimental breaking stress and numerical one, respectively. These indicators are schematically shown in Figure 2. In Figure 2, the experimental stress-strain curve can be obtained from the tensile experiments, while the numerical stress-strain curve would be extrapolated from the FE simulation of the tensile processes. Thereupon, for each GTN parameter combination, there will be a specific numerical stress-strain curve and a set of corresponding response value R. The solution process of the GTN parameters can be schematically described as shown in Figure 3. Based on the tensile tests, the experimental curves can be obtained. Then, simulate the tensile process using different GTN parameter combinations; meanwhile, acquire the corresponding response variables based on the difference between the experimental curves and numerical ones. Construct the response surface using the input GTN parameters and the corresponding response variables. Calculate the accurate GTN parameters based on the established response surface.

Experimental procedures
Based on the theory of RSM, it is prerequisite to acquire the experimental stress-strain data for establishing the response surfaces. Consequently, a sequence of uniaxial tensile tests for Ti-3Al-2Mo-2Zr alloy were performed in the strain rate range of 0.1-10 s −1 at room temperature.
The experimental details are as follows. First, four standard tensile specimens were separated from a block billet. The as-received material of the billet is as-cast Ti-3Al-2Mo-2Zr alloy with the relative chemical compositions as shown in Table 1. The size of the block billet is 100 mm 3 × 100 mm 3 × 140 mm 3 , and its mass is 6.32 kg. Consequently, the mass density of this alloy was calculated as 4,510 kg·m −3 based on the density formula, where ρ is the mass density, m is mass and V represents the volume. The shape and size of the tensile specimens are shown in Figure 4. The uniaxial tensile experiments were carried out on a computer-controlled, servohydraulic WDW-100 tensile testing machine. Prior to the tensile experiments, the surfaces of the samples were grinded to avoid artificial microcracks. Subsequently, both ends of the samples were fixed on the clamps of the tensile machine. The experimental setup with a tensile specimen mounted in terms of grip setup is shown in Figure 5. Then, the specimens were stretched to break with a specific strain rate. During the tensile tests, the load and distance data were monitored and recorded in real time by a computer control data acquisition system.  From the load-distance data, the nominal stressstrain can be obtained by , where σ N represents the nominal stress; F is the applied load; S is the cross sectional area of the specimens; ε N is the nominal strain; d is the recorded distance; D 0 is the initial length of the specimens. Further, the nominal stressstrain data were converted into true stress-strain data based on the formulas: where σ T is the true stress and ε T is the true strain [18]. The true stress-strain curves for Ti-3Al-2Mo-2Zr alloy under varying strain rates are exhibited in Figure 6. In Figure 6, there is a linear stage on each stress-strain curve, i.e., the elastic stage of the material. Consequently, the Young's modulus of this alloy can be evaluated by linearly fitting the elastic stage of the stress-strain curves, and the result is 115 GPa. Besides, it can be seen from Figure 6 that the strain rate has little impact on the variety of peak stress, while it has obvious influence on the breaking stress, peak strain and breaking strain. With increase in strain rate, the values of both the peak strain and breaking strain decrease. These phenomena may be caused by the harder deformation coordination under higher strain rate. It is well-known that the multicrystal plastic deformation depends on dislocation glide and the rotation of grains. At higher strain rate, there is no sufficient time for the rotation of deforming grains, which can more easily lead to local stress concentration and nonuniform plastic deformation. Moreover, at higher strain rate, there is no sufficient time for the energy dissipation of material via plastic deformation. These reasons increase the fracture tendency of material over high strain rate. As for the peak stress, strain rate has no obvious impact on its variety because the peak stress of materials depends only on the inherent bearing capacity of the materials.
After the tensile experiments, the fracture morphologies of the tensile specimens were characterized on a VEGA3 TESCAN scanning electron microscope. The experimental accelerating voltage was set as 20 kV and the work distance was set as 26 mm.   Visualization of the damage evolution for Ti-3Al-2Mo-2Zr alloy  315 4 Results and discussion 4

.1 Numerical model of uniaxial tensile process
To calculate the response variables corresponding to various GTN parameters, it is necessary to carry out the numerical simulation of the uniaxial tensile processes. The simulation is carried out on the FE platform ABAQUS. For saving simulation time, the FE model of the uniaxial tensile process was simplified as an axisymmetric model. The simplified tensile model was meshed with the quadrilateral element (CAX4R), and 657 nodes and 548 elements were obtained. The simplified model and the boundary condition of the tensile simulation are shown in Figure 7. In order to ensure that the strain rate in the numerical simulation is consistent with that in the tensile experiments, the tensile speed exerted on the specimens should be calculated by , where l 0 represents the gauge length of the specimens, 25 mm here. GTN damage model is applied in the tensile simulation. The mechanical property parameters that are required for the tensile simulation of Ti-3Al-2Mo-2Zr alloy are presented in Table 2 [19,20]. According to the tensile experiments, it has been found that the deformation span is very short for the tensile process of this alloy. Therefore, in order to obtain more details on the damage evolution during the tensile process, the simulation is performed using the explicit analysis arithmetic with an automatic time step. The computational cost for the FE simulation is a personal computer with an eight-core processer. Four threads were applied during the computation of the uniaxial tensile FE model. The computational time lasts about four hours for each uniaxial tensile process.

Solution of the GTN parameters based on RSM
In order to construct the response surface of GTN parameters versus response variables, a series of simulation schemes need to be designed. Central composite design (CCD) method is applied in this work to design the simulation schemes. Based on the CCD method, a four-factors-three-levels simulation scheme was designed, as shown in Table 3. These factors correspond to the GTN parameters to be solved, i.e., the initial VVF ( f 0 ), critical VVF ( f c ), nucleation VVF ( f N ) and the final VVF ( f f ).
Taking the strain rate of 0.01 s −1 for instance, the detailed solution process of GTN parameters will be demonstrated. First, the uniaxial tensile process was simulated using the designed GTN parameter combination in Table 3 under the strain rate of 0.01 s −1 . The response values corresponding to different GTN parameters were calculated according to the previously introduced procedure, and the results were recorded as shown in Table 3.
Based on the factors and corresponding response values in Table 3, the final functions of the response surfaces were achieved by least squares method. The expressions are shown in equation (8).   Due to the fact that the response variables R in equation (8) represent the difference between experimental and simulated stress-strain curves, the minimum R-values correspond to the most possible GTN parameter values. In other words, the closer the input GTN parameter is to the actual values, the smaller the response variables R-values will be. Consequently, the GTN parameters will be solved by searching for a suitable parameter combination that makes the response variables reach the minimum almost simultaneously. Considering the significant difference among R 1 , R 2 , R 3 and R 4 , the response variables

Number
Factors Response values should be replaced by dimensionless variables ′ R 1 , ′ R 2 , ′ R 3 and ′ R 4 , as shown in equation (9). Repeating the above calculation procedures, the GTN damage parameters at different strain rates were figured out, and the results were recorded as shown in Table 4.
According to Table 4, the strain rate has significant influence on the GTN parameters, especially for the final VVF. With the increase in the strain rate, the value of the final VVF decreases. Whereas, the nucleation VVF and critical VVF increase with the increase in the strain rate. As for the initial VVF, there exists little discrepancy under different strain rates.
When the specimens deform under high strain rate, the slip speed of dislocations is rapid. The deformation will not be coordinated, leading to local stress concentration on the matrix. The high interior stress will stimulate the deformation and expansion of the microvoids [21,22]. In addition, as the deformation process continues, the microvoids will be more easy to connect under the  concentrated interior stress. When the necking occurs, the VVF reaches a critical value. The critical VVF is higher at high strain rate because the nucleation and expanding rate of voids is quicker at high strain rate. It is remarkable that the final VVF is smaller at higher strain rate, because there is no sufficient time for the nucleation and growth of the microvoids. For clearly displaying the impact mechanism of strain rate on the voids evolution, the scanning electron micrographs of the fracture morphologies for the tensile specimens are shown in Figure 8. There are obvious dimples in Figure 8 and it means that ductile damage is the dominant degradation mechanism for this alloy under the strain rate range of 0.1-10 s −1 . There are fine dimples in Figure 8a and b. In Figure 8c and d, the dimples are coarse and sparse. Therefore, a conclusion can be drawn that the higher strain rate contributes to the expansion and coalescence of the microvoids; meanwhile, the number of the microvoids will be less under high strain rate because of the insufficient deformation time.

Characterization of the damage degradation behavior of Ti-3Al-2Mo-2Zr alloy
In order to dynamically characterize the damage degradation process of Ti-3Al-2Mo-2Zr alloy, the uniaxial tensile processes were numerically simulated using the solved GTN parameters. Subsequently, the damage degradation behaviors of Ti-3Al-2Mo-2Zr alloy were discussed taking the strain rate of 0.01 s −1 as an example. Figure 9 shows the distribution and evolution of the VVF in the tensile process. According to the simulation results, it can be seen that the microvoids nucleate in the interior of the specimen at the beginning of the tensile process, and the microvoids begin to proliferate and propagate with the proceeding tensile deformation. When the true strain reaches 0.1 s −1 , necking occurs and the VVF increases to a critical level. With the ongoing tensile deformation, the VVF keeps increasing until the maximum VVF arrives at the final VVF (f f ). At this time, the specimen loses its overall loading capacity. The fracture position of the tensile specimen in the simulation is in good agreement with the experimental results.
In order to understand the damage process deeply, it is essential to analyze the stress evolution of the tensile process. Figure 10 shows the simulation result of equivalent Mises stress in the tensile process. It is obvious that the distribution of the equivalent stress is uniform at the initial stage of the tensile process. When the true strain reaches 0.1 s −1 , the equivalent stress shows a non-uniform distribution. As the deformation continues, the nonuniform stress distribution becomes more and more obvious, and higher stress concentrates in the middle and lower part of the specimen. It has been introduced in the theory of the GTN model that the stress concentration will promote the proliferation of microvoids. On the contrary, the increase in the VVF will make the loading capacity of materials decrease. Consequently, the mutual improvement between stress concentration and microvoids proliferation leads to the equivalent stress decrease in the later stage of the tensile deformation.
In order to reveal the relationship between damage degradation and plastic deformation, the simulation results of equivalent plastic strain are presented in Figure 11. It is apparent that the plastic deformation zone is very small, and it mainly occurs at the middle and lower part of the specimen, where VVF is higher. The evolution of equivalent plastic strain indicates that the strain is sensitive to the damage degradation of this alloy. Comparing Figure 11a-g with each other, it can be found that the strain distribution is relatively uniform at the initial stage of the tensile deformation. When the true strain reaches 0.1 s −1 , the plastic deformation area intends to contract inside because the VVF is higher at the center of the specimen; meanwhile, necking appears at the area where the equivalent plastic strain is high. Generally, deforming is more easy at the areas with higher VVF because the deformation resistance is lower. When the true strain reaches 0.15 s −1 , the maximum of plastic strain keeps basically stable. Further tensile deformation will make the specimen break off.
Not only the equivalent stress and plastic strain, but also the stress state is a key impact factor for the damage degradation behaviors of materials. Stress triaxiality is widely applied to represent the stress state of materials. High stress triaxiality values correspond to the tensile stress state, while low stress triaxiality values correspond to the compressive stress state. As for uniaxial tensile stress state, the stress triaxiality value is close to 0.3. Figure 12 exhibits the distribution and evolution of the stress triaxiality in the tensile process of Ti-3Al-2Mo-2Zr alloy. It is obvious that the stress triaxiality is 0.3 in most areas of the tensile specimen. The stress triaxiality at the fracture position is relatively high, which means the failure mechanism of this alloy is ductile damage. It is in agreement with the result that is observed from the EBSD images. More interestingly, the variety of stress triaxiality is complex at the fracture position. The center of the fracture position shows a higher stress triaxiality value, while the vicinities along the tensile axis show a lower stress triaxiality value. It indicates that the damage intents to initiate at the center of the specimen, and then expand along the direction perpendicular to the tensile axis. In addition, there is a visible change in the stress triaxiality distribution at the true strain of 0.1 s −1 . This suggests that the microvoids begin to proliferate and expand quickly at the true strain of 0.1 s −1 ; meanwhile, the VVF reaches a critical value.
In summary, the VVF is closely related to the stress, strain and stress state of the materials. The existence of microvoids will lead to the stress concentration, uneven plastic deformation and the change in local stress state. And the latter will contribute to the increase in VVF in return. This is the damage degradation process of materials during the uniaxial tensile deformation.
Since the evolution of VVF in the tensile process has been visualized by FE simulation, the degradation process of this alloy can be evaluated dynamically. Here the damage degradation degree is defined as the percentage of the current VVF to final VVF. Correspondingly, the damage degradation curves of Ti-3Al-2Mo-2Zr alloy at different strain rates were plotted as shown in Figure 13. According to Figure 13, it is apparent that all the degradation curves show a similar tendency, and the degradation speed is faster at higher strain rate. This point can be explained by the damage mechanism of this alloy. As discussed in the above Sections, the nucleation and proliferation of the microvoids is the dominant damage mechanism for this alloy. Therefore, in the early stage of the deformation, the progress of degradation process is slow because of the low nucleation speed of microvoids. In the middle and later stages of the tensile deformation, the degradation is faster because the nucleation and expansion of microvoids are promoted once the VVF reaches the critical value (f c ). Particularly, the promoting effect is enhanced under higher strain rates because the stress concentration and uneven plastic deformation is more severe at higher strain rate.

Conclusion
The damage degradation behaviors of Ti-3Al-2Mo-2Zr alloy have been studied by the FE simulation combined with uniaxial tensile experiments at the strain rate range of 0.1-10 s −1 . And the following conclusions have been drawn.
(1) For Ti-3Al-2Mo-2Zr alloy, the GTN parameters corresponding to different strain rates were worked out based on RSM.
(2) Strain rate has significant influence on the GTN parameters. Higher strain rate contributes to the growth and coalescence of microvoids. In addition, the final VVF is lower at the higher strain rate because of the insufficient deformation time.
(3) The evolution of microvoids volume fraction was visualized based on a FE platform and the solved GTN parameters. The damage degradation curves of Ti-3Al-2Mo-2Zr alloy were drawn out. The degradation speed of this alloy is slow at the beginning stage of the tensile process and then accelerates once the VVF reaches a critical value. Visualization of the damage evolution for Ti-3Al-2Mo-2Zr alloy  323