Relevance of Incorporating Cavity Shape Change in Modelling the Ductile Failure of Metals

The purpose of this paper is to assess the relevance of considering the cavities shape change in the context of physically based modelling of the ductile rupture inmetals. Two thermomechanical models have been used in this study: the GTNmodel, developed by Gurson, Tvergaard, and Needleman for spherical cavities keeping their shape unchanged during loading, and the GLD model, proposed by Gologanu, Leblond, and Devaux for ellipsoidal cavities that can change their shape. The GTN and GLD models have been extended to take into account material thermal heating due to plastic dissipation. These two constitutive laws have been implemented into the commercial finite element code Abaqus/Explicit in order to simulate the necking of a round bar and the failure of a sheet deep drawing. The results showed the importance of incorporating the shape effects of the cavities for a correct description of the material failure.


Introduction
The ductile fracture of the materials is associated with the degradation of the mechanical properties under the effect of three physical mechanisms: nucleation of the cavities, their growth under appropriate loading, and the coalescence of the neighbouring voids at an advanced deformation stage.Several nucleation criteria have been proposed in the literature.Some of these criterions have been formulated by exploiting the energetic criterion of germination, initially enunciated by Gurland and Plateau [1,2].Needleman and Rice proposed a nucleation criterion depending on the inclusions rate using a phenomenological approach [3].Other nucleation criteria have experimental origins [4].The modelling of void growth has undergone several important contributions.Gurson [5] proposed a kinematically based limit analysis approach of a hollow sphere with an isotropic rigid-plastic matrix and developed a macroscopic yield criterion depending on the porosity.Tvergaard and Needleman extend this model to incorporate cavities interaction and to describe cavity coalescence in a phenomenological way [6,7].The new model is known as GTN model.Recently Nahshon and Hutchinson [8] proposed an interesting modification of the porosity evolution law to account material damage under conditions of low or vanishing triaxiality.With the same aim Bai and Wierzbick [9], Xue [10], Xue and Wierzbicki [11], Danas and P. Ponte Castañeda [12], and Zhou et al. [13] studied the dependency of ductile fracture on the Lode angle in metals by incorporating shear mechanisms effects in the porosity evolution law.Although these new expressions do not agree with the classical porosity law derived from homogenization and having a physical interpretation, these extensions offer good results and improve them substantially at low triaxialities; see, for instance, [14][15][16][17][18].All these models have been proposed by considering the matrix behaviour governed by von Mises-like model.carried out a series of unit cell calculations for initially spherical voids embedded now in an orthotropic matrix.Cazacu et al. [22] assumed the matrix obeying Tresca's criterion when performing limitanalysis of the hollow sphere.Gologanu et al. [23,24] extend the Gurson analysis to introduce cavity shape effects in the socalled GLD model.They proposed a model incorporating the porosity and a new internal variable: the shape parameter.The model has been assessed by Pardoen and Hutchinson [25] in comparison to unit cell prediction with different spheroidal cavities in the case of a hardenable matrix.Ould Ouali and Aberkane [26][27][28] studied namely the microstructure evolution during rolling of an A1050 aluminium sheet and introduced thermomechanical coupling in the GLD model.Madou and Leblond's [29,30] extended the GLD model to general ellipsoidal voids having three distinct axes.
The numerical simulation of the mechanical behaviour of metals undergoing forming operations has been a great contribution to the optimization of manufacturing processes.Indeed, the study of the material failure during such operations must prevent the arising of the internal or surface cracks as well as their propagations before causing the total rupture of the structure.Given the high number of parameters controlling the deformation mechanisms under such solicitations the recourse to modelling is necessary.In the present work, the modelling includes the effects of the microstructure evolution in the material on the overall response of the part being formed.However, it is generally observed in material forming processes low triaxiality which can be an obstacle to use the classical physically based models.The purpose of the study is to assess the relevance of taking into account the evolution of cavities shape parameter as proposed in the GLD model to correctly predict the ductile failure of metals.Benzerga et al. [32] showed that the consideration of cavity shape change is currently the most appropriate physical description for modelling ductile fracture at low triaxiality.Two applications have been chosen: necking of a round bar under tensile loading and the deep drawing of an aluminium sheet.Comparison of the numerical and experimental results shows that neglecting cavity shape change leads to poor numerical predictions.

Mathematical Modelling
2.1.Voids Growth 2.1.1.Gurson and Tvergaard (GT) Model.Gurson [5] developed a mechanical model for porous materials by limit analysis of a spherical Representative Volume (RVE) containing a spherical confocal void (Figure 1).Tvergaard [7] proposed an extension of the Gurson plasticity criterion by introducing three coefficients  1 ,  2 , and  3 =  2  1 to account for the interaction between cavities.It also replaced the yield stress  0 with the flow stress  in order to incorporate the isotropic work hardening of the matrix.The plastic potential is in the following form: where ∑  is the macroscopic von Mises equivalent stress, ∑  the macroscopic mean stress, and  the porosity of the material.The matrix hardening was introduced using the equivalence between the macroscopic and the microscopic plastic dissipation as proposed by Gurson [5]: where Ė  is the macroscopic plastic strain rate tensor and ε  the cumulated plastic strain rate of the matrix.The evolution law of the porosity can be approximated using the condition of plastic incompressibility of the matrix.In the case where nucleation of new voids is not taken into account, this is given by ḟ 2.1.2.The GLD Model.Gologanu et al. [23,24] extended the Gurson's analysis to a prolate or oblate ellipsoidal RVE containing a confocal cavity and subjected to some asymmetric loadings (Figure 2).The authors assumed that the voids can change their shape under the effect of loading.Such modelling requires the introduction of a new internal variable: the shape parameter  = ln  = ln(/).
The expression of the plastic potential is given by the following [21,25,26]: where ‖.‖ is the von Mises norm and X a second-order tensor defined by X = (1/3)(−  ⊗   −   ⊗   + 2  ⊗   ).∑ ℎ is the hydrostatic stress weighted by coefficients depending on the cavity geometry.It governs the voids growth [19]. and g are the coefficients of the GLD related to porosity and shape parameter.The coefficient q is the equivalent of the Tvergaard coefficient in the GT model [33]: where The evolution law of the shape parameter  is obtained by making some analysis and calibrations with cells calculations [23,24]: where ℎ 1 and ℎ 2 are two parameters depending on , , and triaxiality Τ = ∑  / ∑  .

Coalescence Stage.
Prior to coalescence, the GT or GLD model is used to describe the material behaviour [3,31].These models cannot capture any realistic void volume fraction during the final coalescence stage.So, the description of the material capacity loss needs to use a coalescence model.In this study, the function  * proposed by Tvergaard and Needleman is adopted [7]. * allow the reproduction of the rapid void volumes growth at the material coalescence beginning from some critical value   of : where  = (  −   )/(  −   ) describe the rapid increase of the porosity at the coalescence.  and   are the porosities at the onset of coalescence and at final failure, respectively.In the model, the critical porosity can be determined using the Thomason's criterion.This criterion has been formulated with the assumption that the onset of coalescence occurs when the limit plastic charge is reached at the material ligament between voids (see Figure 3).The coalescence criterion is therefore written in the case of the GTN model, as follows: Equation ( 8) can be reformulated in accordance with the dimensionless variables  = / = ((3/2) − ) 1/3 and  = / =  0 exp((3/2)  ) used in the GLD model.We write the following [19,28]:

Thermal Heating due to Plastic Dissipation.
The GTN and GLD models have been extended to take into account the thermal heating due to plastic dissipation.This extension is formulated using the hypothesis of a rigid-plastic behaviour of the matrix according to Gurson assumption.The thermoplastic coupling is introduced by enforcing the dependence of the flow stress , characterizing the isotropic matrix hardening, on equivalent plastic strain   and the temperature T [27,28].The following hardening law is chosen: where  0 is the initial yield stress and  hardening exponent. is a material parameter.
In the absence of external heat sources, the temperature evolution law is solved with the general heat equation under adiabatic conditions (neglecting the conductivity term): where  is the material density,  is the specific heat capacity, and   = (1 −  − /  ) is the expression of the Taylor-Quinney coefficient for porous material.
2.4.The Thermoelasticity Law.For the sake of consistency, the following thermoelastic law is introduced [27,28]: where E ℎ = ( −   )I is the strain tensor due to thermal expansion and Σ some objective derivative of the stress tensor.I is the second-order unit tensor.

Numerical Implementation.
The GTN and GLD models have been implemented into the finite element code Abaqus/Explicit using the Vectorized User MATerial (VUMat) subroutine and used to carry out the numerical simulations.The "elastic predictor-plastic correction" algorithm proposed by Aravas [34] is adopted in the local resolution of the constitutive laws GTN and GLD.The details of the implementation of these two behaviour laws, as well as the validation in comparison with the predictions of cell calculations using an explicit technique, are given in [26,27,33].
Aravas's algorithm is based on the decomposition of the stress tensor into its hydrostatic   and deviator parts   and the plastic strain increment tensor Δ  into its volumetric and deviator parts: Δ p and Δ  are used by Aravas as primary variables.So, the system is reduced to two non-inear equations: Φ (, , ) = 0 After the resolution of ( 15) and ( 16), the hydrostatic stress Σ  and the equivalent stress Σ  are updated: Then the state variables are actualized according to the following: where ℏ E = /  and ℏ T = / are, respectively, the plastic and thermal modulus.Details of the implementation procedure are done in the flowchart of Figure 4 and can be found in [19,20,24].
It is important to notice the GLD growth model implementation requires specific treatment due to the "limits" shapes that the cavities can take during the loading process: penny-shaped crack ( → 0, S→-∞), "sandwich" structure ( ≡ finite, S →-∞), spherical ( ≡ finite, S=0) or cylindrical ( ≡ finite, S → + ∞) shapes, or a needle ( → 0, S → + ∞).Specific numerical treatment of the particular expressions of the plastic potential Φ according to each case should be carried out.

Results and Discussion
In this section, we discuss the results obtained with the GTN and GLD models extended to take account the thermal heating due to plastic dissipation and to predict the onset of coalescence using Thomason's criterion.In order to highlight the effect of cavity shape change, we make use of the same  material parameters for both the GTN and GLD models with the exception of those related to the state variable describing the cavity shape change: the cavity shape parameter is not described in the GTN model.Two applications were chosen: tensile of a round specimen and deep drawing of a sheet.

Necking of Cylindrical Bar.
The GTN and GLD models are used to simulate the rupture of a round bar subjected to tensile loading.It has been found experimentally that the failure of the bar occurs according to the so-called "cup cone" mode [7,35].Indeed, Thomason observed that the crack appears in the centre zone of the specimen and then spreads horizontally in the plane perpendicular to the loading axis.Near of the outer surface, the crack deviates with an angle of 45 ∘ from the propagation direction forming a cup cone fracture surface.The height/radius ratio of the cylindrical bar is taken equal to  0 / 0 = 4. Three hundred and twenty (320) axisymmetric elements with reduced integration (CAXR) are used to mesh the quarter of the specimen.The initial and deformed meshes, as well as the loading conditions, are shown in Figure 5.The constituent material of the specimen is a standard steel with a rigidity / 0 = 500 and Poisson's ratio ] = 0.28.The hardening exponent  = 0.1.
The material is assumed to be initially dense ( 0 = 0).We use the same coalescence and nucleation parameters as those used by Tvergaard and Needleman when studying necking of the cylindrical bar [7].These parameters are done in Table 1.
The evolution of the nominal stress   = / 2 0  0 versus the longitudinal nominal strain   = / 0 and the radial nominal strain (reduction of radius)   = Δ/ 0 calculated with the two models GTN and GLD are represented in Figure 6.During the test, the nominal stress   increases to a maximum that corresponds to the onset of necking, then decreases [7].This diminution of the nominal stress is accentuated at the end of the test because of the coalescence of the cavities until the final rupture.A difference appears in nominal stress levels reached by the GTN and the GLD models.This difference can be explained by the effect of the cavity shape change on the resistance of the materials observed, namely, in cell calculations [19,[24][25][26]33].These authors have related these results to the morphological anisotropy generated by the cavities shape that the GLD model allows to take into account.
In the case of a simulation with the GLD model, the maximum force is equal to   max = 1.344 2 0  0 whereas that predicted by the GTN model is   max = 1.311 2 0  0 .The two maxima are obtained for the same radial nominal strain  neck  = 0.0847 and for close longitudinal strains of about  neck  = 0.04.These main results are reported in Table 2.

𝜀 fail
and  fail  are, respectively, the longitudinal and the radial nominal strain of the bar at failure.A qualitative explanation of this difference was given by Aravas and Ponte Castañeda concerning this test using the model proposed by Aravas and Ponte Castañeda [36].In the case of the GLD model, the justification can be done regarding the evolution of the microstructure (the porosity  and the shape parameter ) at the neck.Figure 7 represents the evolution of the nominal von Mises stress ∑  / 0 , the porosity , and the shape parameter  versus the nominal radial strain in three elements (A), (B), and (C) of the bar located at the neck (see Figure 7).It is well-known that the effective nominal stress depends on the net area of the section (without voids) perpendicular to the load at the neck that supports the load in the case of tensile loading.This area decreases when the total area of the section of the bar decreases (necking) and/or when the porosity in the plane carrying the load increases.As previously mentioned, the section of the bar for which the difference between the maximum force predicted by the GLD and GTN models is observed for the same radius reduction ΔR/R 0 = 0.0847.Consequently, the difference in the intensity of the loading force can only be justified by the second point: the surface occupied by the voids at the neck of the bar.We will study this point in more detail by considering the material microstructure evolution of the bar (Figures 7(b) and 7(c)).
(i) Figure 7(b) shows that the porosity predicted by the GTN and GLD models are the same for the three elements located at the neck of the bar.These curves coincide almost at least until the nominal radial strain   = 0.0847 which corresponds to the onset of necking.(ii) In Figure 7(c), we can see that the voids tend to lengthen in the direction of loading, in the case of the GTN model, since the shape parameter S increases.
For the GTN model, the cavities keep the spherical shape.It follows that, for the same porosity, the area occupied by the voids in the section at the neck is smaller in the case of the GLD model: the cavities shape have changed from spherical (initial shape  = 0) to elongated voids ( > 0).
So, for the same volume fraction of voids, the surface of the voids in the section of the neck supporting the load is smaller in the case of the GLD model compared to the GTN model.Indeed, during the simulation with the GTN model the increase of the porosity is done by the increase of the diameter of the void.So, the effective surface is smaller than that predicted by the GLD model.Knowing that the effective stress is directly related to the applied force, the difference between the effective surfaces explains the difference between the maximum forces achieved.

Deep
Drawing of a Sheet.The deep drawing of a thin sheet of aluminium alloy is studied.We are interested in the Swift deep drawing test which is applied to an aluminium sheet having a thickness of 1 mm, until failure [31].The dimensions of the sheet and the tools are shown in Figure 8. Due to geometry and applied load symmetries, the test is modelled under axisymmetric conditions.The mesh of the sheet is composed of 640 axisymmetric elements with reduced integration CAXR with 04 layers of elements on the thickness.The friction between the sheet and the tools has been determined experimentally by the CETIM/Senlis (France) and is equal to   = 0.17 [31].The clamping force applied by the blanking tool is   = 5 .
The material properties of the aluminium alloy sheet are Young's modulus  = 70, Poisson's ratio ] = 0.33, density  = 2700/ 3 , and yield strength  0 = 114.The hardening exponent is  = 0.146.In the absence of experimental data on the microstructure of the aluminum alloy, we have adopted the initial porosity:  0 = 0.1%.This level of initial porosity is generally observed in such material.For choosing the initial value of the shape parameter, a parametric study has been carried out using initial oblate, spherical, and prolate cavities.It is found that a value of  0 = 1.38 gives good results in this case.
Figure 9 shows a comparison between the experimental results obtained by the CETIM/Senlis with the numerical predictions of the GTN and GLD models.The curves describe the evolution of the deep drawing force as a function of the punch displacement.Good agreement is found between experimental and numerical results in the case of the GLD model.However, a difference appears during the simulation with the GTN model.
Table 3 resumes the main values of the forcedisplacement curve during deep drawing of the aluminium sheet: maximal force and corresponding displacement.In  the estimation of the error, the experimental results are taken as reference.The error made during the simulation with the GTN model (13.93% and 10.6%) is more important than that made using the GLD model (1% and 0.87%).An explanation of these differences can be found regarding the orthotropic behaviour of the sheet experimentally identified by the CETIM/Senlis [31].In the GTN model, the matrix behaviour is assumed to be isotropic.In the GLD model, a morphological anisotropy is introduced by the change of the ellipsoidal cavities shape.
Figure 10 shows comparison of the stamp shape obtained experimentally and that described numerically with the GLD model.We observe that the two geometries are practically identical and that the rupture appears in the same region of the sheet.

Conclusions
Two physically based constitutive laws, the GTN and the GLD models, have been implemented in the ABAQUS calculation via the VUMAT subroutine.Both models have been extended to account for thermal heating due to plastic dissipation under adiabatic conditions.The theoretical and numerical aspects of these formulations are detailed.These models have been used for the simulation of a round bar under tensile loading and deep drawing of aluminium sheet.From the comparison of the predictions of these models, we can state that incorporating cavity shape changes is an important step in the correct modelling of the ductile failure of metals As a perspective for this work, we plan to carry out experimental tests and unit cell calculations in order to verify the pertinence of the shape parameter evolution law under different stress states.Especially in the case of thermomechanical loading, another interesting extension is to introduce the orthotropic material behaviour existing, namely, in thin structures.Finally, the incorporation of shear mechanisms will also significantly improve the numerical predictions of ductile fracture arising in some particular engineering problems.

2 MathematicalFigure 1 :
Figure 1: Geometry of the RVE corresponding to the Gurson model.

Figure 4 :
Figure 4: Flowchart of the implementation of the GLD model.

Figure 7 :
Figure 7: Evolution of the stress state and the microstructure at the neck.(a) Equivalent nominal stress ∑  / 0 .(b) The porosity .(c) The shape parameter .

Figure 8 :Figure 9 :
Figure 8: Schematization of the deep drawing test: initial mesh and boundary conditions.

FailureFigure 10 :
Figure 10: Distribution of the porosity in the sheet at failure and experiment [31].

Table 1 :
Material parameters used with the GTN and GLD models.

Table 2 :
Representative values of the stress-strain curves corresponding to a cylindrical bar under tensile load.

Table 3 :
Representative points of the force-displacement curve during deep drawing of the sheet.