A consistent use of the Gurson-Tvergaard-Needleman damage model for the R-curve calculation

The scope of the present work is to point out a consistent simulation procedure for the quasi-static fracture processes, starting from the micro-structural characteristics of the material. To this aim, a local nineparameters Gurson-Tvergaard-Needleman (GTN) damage law has been used. The damage parameters depend on the micro-structural characteristics and must be calculated, measured or opportunely tuned. This can be done, as proposed by the author, by using an opportunely tuned GTN model for the representative volume element simulations, in order to enrich the original damage model by considering also the defect size distribution. Once determined all the material parameters, an MT fracture test has been simulated by a FE code, to calculate the R-curve in an aeronautical Al-based alloy. The simulation procedure produced results in a very good agreement with the experimental data.


INTRODUCTION
luminium alloys are widely used in aircraft structures, transportation industries and civil engineering due to their relative lightness and versatility.In particular, the fuselage skin, which thickness ranges from 1 to 3 mm, is often made of aluminium alloys capable of fully ductile behaviour.Thus, in order to better take advantage from the material metallurgical characteristics, the damage tolerance design concept is often used in the aircraft structures sizing.Following this approach, the presence of macroscopic cracks during the ordinary service life can be tolerated, as long as they are under control.For this reason, in this field the study of both the static and the dynamic response of cracked panels in ductile regime, in terms of residual strength it assumes a relevance that can hardly be overestimated [1].The socalled R-curve, that represents the critical stress intensity factor versus the crack length for a given structure, is a widely used tool to design by following the above criteria.Unfortunately the R-curve behaviour depends both from material and geometrical properties and then it changes as the test set-up changes [2,3].A reliable calculation model of the R-curve is, for this reasons, very important for the design purposes, in order to drastically decrease the amount of the very expensive physical tests.From a micro-mechanical point of view, ductile failure is typically characterized by three coupled stages: nucleation, growth and coalescence of voids, which are induced in the metal alloy matrix by the presence of a variety of microscopic defects.In order to consider the effects of these voids evolution on the stress-carrying capability of a mechanical continuum during simulations, damage mechanics models are often used.Many models for ductile fracture growth have been proposed in the past [4 -11] and, for their relative simplicity and efficiency the Gurson-Tvergaard-Needleman (GTN) model and many of its variants are proposed by several authors and also included in some public FE codes, like WARP3D [12].These kinds of models for ductile fracture include the micro-mechanical effects of void nucleation, growth, and coalescence of micro-voids in the constitutive equation used to describe the mechanical continuum.The coalescence mechanism, in particular, induces a strain softening at the large-scale material response.Therefore, the A material becomes susceptible of bifurcations, instabilities and localized deformation modes.This is because the governing equations become hyperbolic in some unknown area of the continuum and make ill-posed the boundary value problem.As a consequence, numerical solutions chronically suffer from mesh dependency [13 -15].A popular and well proven in literature way to avoid this problem is the use some nonlocal formulations of the damaged material model.This is done by replacing the local variables denoting stress and strain fields, as well as the damage variables and the other internal ones, with nonlocal variables, which are space-weighted averages of the local ones.The averaging volume is related to a characteristic length l c depending from the failure evolution at the micro-structural scale.The nonlocal formulations are capable to regularize ill-posed problems associated with damage, but require a very large increasing of the computational amount with respect to the local damage formulations [16 -19].From the computational point of view, the regularization is strictly related to the characteristic length; in fact, with l c a constraint to the material law is introduced in addition to the large scale material behavior [20].Also the local models can be regularized, and the easiest way to do this in a FE analysis is introducing the characteristic length not at the material model level, but more roughly at the element level.In other words, the element dimension will be strictly related to the characteristic length lc.In many works the latter computational strategy is adopted [21 -24].They typically present many free parameters that should be calibrated and fixed in a physically consistent way in order to make the calculation procedure effective.In the present work, a complete procedure to calculate the R-curve starting from a local form of the GTN model is presented and a numerical application on a widely used aeronautical alloy, 2024A-T351, is reported, using material and test data provided by industrial research offices.Differently from the purely phenomenological approach that is usually employed to calculate all the model parameters, a procedure to consistently calculate some of these parameters on the basis of objective micro-structural data has been performed.In order to enrich the original damage model also the defect size distribution has been considered.With such methodology, the arbitrariness of the parameters tuning is strongly reduced, and a double result is attained: the number of the FE run needed to simulate the R-curve behaviour drastically decreases, and the transferability of the whole methodology to different geometries and boundary conditions is improved.Once determined all the material parameters, to calculate the R-curve an MT fracture test has been simulated by the free FE code WARP3D.

THE GTN MODEL FOR DUCTILE FRACTURE
n a homogeneous ductile metal model, the total strain amount usually leaves unchanged the volume, because the plastic part of the strain is dominant with respect to the elastic one.On the other hand, in a material model containing voids which matrix phase is ductile, the volume can globally change, due to the local plastic flow arising around the voids boundary and, consequently, the void volume change.Then, the material model response to an imposed global strain will be a stress-strain curve presenting a softening curve, nevertheless the matrix material model have hardening behaviour.In such model, the voids can grow until the global load carry capability becomes negligible.The Gurson-Tvergaard (GT) model is able to explain the local strength decreasing during the fracture process of ductile metals in the intermediate phase between the nucleation and the coalescence of voids.Standing this characteristic, in the void growth model the number of voids is kept constant.Nucleation and coalescence are taken into account before the homogenization of the material model, applying some opportune corrections, originally proposed by Needleman et al., directly to the stress-strain cell response.Updated in this way the GT model, it is usually named GTN model.

Gurson-Tvergaard model
The homogenization technique is based on the stress-strain characterization of a Representative Volume Element (RVE), that can be roughly defined as the minimum material volume containing all the micro-structural information of a heterogeneous material, related to the specific problem under investigation (see, e.g., [25 -27]).In the Gurson-Tvergaard void growth model, the RVE of the material is considered as a cubic volume with a single void, yet existing in the virgin material.The initial void volume fraction f0 should be chosen as the 'equivalent' volume fraction corresponding to the physical distribution of voids inside the RVE (Fig. 1).The resulting homogenized material model was defined by modifying the analytical solution of the single void cell (Fig. 1) performed by Gurson and limited to a rigid-plastic material model of the matrix.The corrections take into account for the hardening of the matrix material and for the presence of void cell array instead of a single void cell.They are represented by the coefficients q 1 , q 2 , q 3 in the critical surface definition, originally introduced by Tvergaard.The homogenized constitutive law is defined by (tensors are indicated with bold characters):   1 : Where: Eq. ( 1) defines the plastic surface; Eq. ( 2) is the plastic flow rule; Eq. ( 3) is the void growth rate definition; Eq. ( 4) imposes the equivalence between micro and macro-mechanical plastic work; Eq. ( 5) is the global stress-strain relationship; Eq. ( 6) is the plastic hardening power law for the matrix material. 0  0 : yield equivalent stress and strain; q1, q2, q3 : Tvergaard correction coefficients.

Voids nucleation model
The nucleation of the first void in an array of cells is implicit in the above model.With the increasing load, further voids can nucleate in the RVE, and this fact is indirectly taken into account in the homogenized material law by increasing the void growth rate defined in Eq. ( 3).
In the present study the following relation is used: where: The corrected void growth rate becomes: The expressions (7-9) relate the voids growth rate due to the nucleation process to the internal equivalent strain rate d .
This acceleration is driven by the parameter A and it is limited to an internal strain condition in which  is quite close to a certain value N, representing the matrix strain value where the nucleation takes place.
Three new parameters are so introduced in the constitutive law (1-6): fN : magnitude of the volume fraction rate increasing; N : main value of the internal strain for which the nucleation takes place SN : standard deviation of the Gaussian distribution.
The void nucleation law completes the analytical description of the GTN model.

Cell dimension
It is well known that the undefined equilibrium problem in a material domain which constitutive law has softening branches loses ellipticity near the critical point and then it is hill-posed [27].
If the finite element method is used to search a solution, disregarding the above problem, the resulting discrete equation can be solved, but the solution is chronically mesh dependent and the convergence can't be reached.This occurs because near the softening condition the strain field would become discontinuous and, correspondently, in the discrete problem the strain tends to concentrate in a small zone, which size depends on the elements size.Smaller the element size, smaller the zone affected by the strain concentration.This behaviour is commonly named strain localization.
If the softening law derives from a physical model of the material, the incongruence can be explained noting that, at the micro-structural scale, the material domain can't be longer be considered homogeneous, so the softening at the large scale material law can be attributed to the micro-structural geometric changes.Apparently, the consequence of this fact is that we can't use local homogenization techniques if the global material law presents a softening behaviour.In practice, using a local FE model, we can avoid the mesh dependence by keeping fix the element dimensions in the fracture zone, regardless with the macroscopic geometry of the domain.In particular, the dimension D 0 to be mandatory controlled is the cell length in the direction perpendicular to the crack plane.The correct value of the element dimension depends on the scale of the local fracture process related to the particular material under consideration.In this way, the element dimension becomes a further parameter of the material constitutive law.In other words, we couldn't use a local homogenized law because the analytical equilibrium problem is hill-posed, but the error that cannot be eliminated in the FE solution of this problem can be driven in order to reach a numerically correct result.

Voids coalescence
In the constitutive law (1-6) the internal strain state is represented by a unique variable  .This simplification can be accepted only in the early phases of the void growth, in which the internal strain (and stress) distribution doesn't vary too much.In the following load phase the matrix material model becomes unrealistic.In fact, in this phase the fracture advances between two consecutive voids (coalescence) and then the RVE model with the single equivalent void itself becomes unrealistic.The voids coalescence, and the consequent crack advancing, is modelled in the FE code by eliminating the corresponding element (named killed element) when a critical volume fraction value fc is reached.Many studies have demonstrated that the critical value should be much less than 1, which is the value rigorously provided by the theory.For many ductile materials the critical volume fraction can be set-up between 0.15  0.25, and a phenomenological correlation furnishes the following empirical law [21,22]: Along the crack plane, the element traction forces are still present in this condition (f = fc) and they are gradually decreased until zero is reached, using a multiplicative coefficient , related to the cell dimension D0 and depending from a further parameter , to be set.
In Eq. ( 11) D and Dc indicates respectively the actual and the critical value of D0, averaged over the cell volume.
CELL CALIBRATION he material model described above needs the calculation, the measurement or the phenomenological tuning of (at least) 13 parameters:  The internal plastic-hardening parameters N,  0 ,  0 .
 The nucleation parameters fN , N , SN .
 The cell height D0.  The initial and critical void volume fractions f 0 , f c .
 The force release parameter .

Critical void volume fraction and force release parameter
These parameters are related to the coalescence phase.For the present calibration they have been chosen as mean values from the ones present in literature [21,22].
f c = 0.2; Notwithstanding the above position, the release forces parameter  contributes to form the energy released in the element extinction process, so it can't be rigorously considered a merely calculation parameter.The coefficient defined in (11) depends on the component of the strain orthogonal to the crack plane.In fact, we can write: where  rel is the stress value present during the force release process (coalescence).Thus, the corresponding energy density w released in the coalescence phase can be evaluated, for small , c: where c is the stress corresponding to the critical strain c (reached when f = fc).
Then, the parameter  is related to the energy released during the coalescence process.

Cell height and initial void volume fraction
Both the cell height D0 and the initial void volume fraction can be estimated considering as a starting point the defect size distribution.As shown above, the cell dimension, and in particular the cell height, is a fundamental parameter because, when the strain localization occurs, the cell height coincides with the localized strip of material.So, from this point, the most of the energy released is proportional to this parameter.On the other side, in the present model the energy released is also influenced by the coalescence parameters and by the shape of the stress-strain curve (that depends on the nucleation and Tvergaard parameters).Then, the correct cell height should be chosen taking into account the influence of all the above parameters in the fracture process.
If we consider the cell as a representative volume element (as it is usually done) it should contain a sufficiently representative voids distribution.Of course, a greater RVE contains more micro-structural information, but it is less useful in the present model because the strain would localize in a strip with fixed height.In conclusion, the cell height can't be precisely defined without a critical consideration of all the constitutive law parameters influence.
As an approximated evaluation, we can consider the bigger defects above the last 20% of the total volume fraction as it were the driving defects of the localization process.In this way, with reference to the Tab.1, an RVE, which dimension is coincident with the cell height, will contain one defect with 10+ m diameter and the appropriate number of minor defects as deduced from the size distribution (Tab.2).The cell height and the initial void volume fraction can be finally estimated:

Category
The D value is in good agreement with the literature values that range between 80÷200 m.
This procedure can appear quite arbitrary.In fact, if the categories were chosen in a different way (for instance the last category could be 8+ or 12+), the resulting D and f0 parameters were different.On the other side, this can be avoided if we consider for all the distribution, as stated above, the bigger defects above the last 20% of the total volume fraction; in this case, is arbitrary the 20% value.
As discussed above, a certain degree of arbitrariness is un-eliminable because the cell height (and consequently the initial volume fraction) has the double role of dimension of the RVE and dimension of the strain localization zone.Then, once established the cell height on the basis of a reasonable RVE choice, the effect of the strain localization should be taken into account by tuning both the nucleation parameters (that can be considered responsible for the minor defects growth) and the coalescence parameters f c , .
As an example, in the pictures below (Fig. 2) is reported the equivalent strain distribution inside two different twodimensional RVEs with the same initial volume fraction (see Fig. 1 for the unstrained configurations).The first, is the simplified Gurson RVE with a single void; the second, has a void distribution compatible with the Tab. 2 distribution.It can be deduced that the localization zone dimension is larger in the first that in the second RVE.This fact influences the global cell response (Fig. 3).A more general discussion on this topic can is presented in [28].Supposing that the second RVE represents a sufficiently accurate evaluation of the material behaviour, the pre-defined cell height seems too large in relation to the 'true' localization zone, as it results comparing the different entity of energy involved in the load process.Once established a compromise value D = 100 m < 130 m, the residual difference should be taken into account by tuning, as written above, the nucleation and the coalescence parameters.

Tvergaard correction coefficients
The correction coefficients q 1 , q 2 , q 3 can be calculated from a single void RVE model.The global cell behaviour has been calculated for two different load conditions: 1) Imposed global strain x > 0 and free surface condition for y and z direction.
In the following pictures 4-5 are presented the cell mesh and the resulting strain distribution in the x direction for the two considered load cases.It can be noted that in the first case the x-direction strain is about uniform, but in the second there is a consistent strain concentration on the void boundary.This agrees with the constitutive law (1-6), that attributes the void volume fraction growth to the mean stress (Eq. 3) that in the first case is much lower than in the second one.
The global RVE responses have to be compared with the eq.(1-6) stress-strain curves resulting from an equivalent load process.This comparison allows to determine the Tvergaard correction coefficients.In the present work, the analytically derived relation is adopted [23], that: In the subsequent refinement of the model calibration the parameter q3 could be left free.The two remaining free parameters q 1 , q 2 are determined by imposing on the homogenized law the condition that the characteristic values: are equal to that resulting from the FEM calculation.The first load case (uniaxial strain) is quite insensitive to the parameter variation, due to the very slow void growth rate present in this condition, so the calibration has been done with reference to the second load case.Finally, the resulting values are: In Fig. 6, the comparison between the FE model of the cell and the homogenized law (1-6) stress-strain curve is reported for the two load cases.

Nucleation parameters
In the present work, the nucleation parameters have been considered as correction parameters, not directly related to the microstructure.Then, the parameters fN, N , SN have been used to fit the experimental-numerical results in the residual strength curve, and their values are reported in the following section.In Fig. 7 the coalescence law (11) has been included into the constitutive law, with the imposed values f c = 0.2 for the starting point, and  = 0.1 for the curve slope, as established before.It is shown that the energy related to the coalescence process is a non negligible part of the total cell energy.

R-CURVE DETERMINATION
he R-curve calculation procedure has been tested on an M(T) specimen model, which experimental results have been published in [2].The M(T) specimen used for the residual strength results is made by a 1.28 mm thick 2024 sheet material.The T panel width is H = 500 mm and the initial half crack length is a 0 = 49.78mm.Testing and data evaluation were done according to ASTM specification E561-86 for R-curve determination.The residual strength tests were done under displacement control.

Mesh generation
The FEM mesh has been generated within the ANSYS code pre-processor module.All the elements are eight nodes brick, in order to be allowable for the subsequent calculations with the WARP3D [12] solver.The specimen has been modelled taking into account the three symmetry conditions.Along the crack plane, all the elements are cubic with length D0 = 106.667m (see the calibration section).Along the thickness t, six elements have been generated, so that t/2 = 0.64 mm = 6D 0 .In Fig. 8 some particulars of the mesh generation are reported.
It should be noted that, as shown in the last picture, there the row of elements before the initial crack front is not present.Further, the elements that will be killed during the crack advance simulation are indicated with a different colour, that indicates two different materials.In fact, for computational reasons, the killable elements have to be separated from the non-killable ones; otherwise, the solver can try to kill the above elements in contrast with the crack plane definition.This contrast will produce a program error (the limitation is related only to the force-release kill procedure).The elements immediately before the initial crack front have been deleted in order to avoid an undesired 'collaboration' to the initial crack strength.The as generated mesh, and the boundary conditions have been 'translated' into the corresponding WARP3D commands by using a Fortran routine.

Results: Load vs. crack length curve
From the experimental test data, the following material parameters have been chosen for the numerical analysis: The other constitutive law parameters values, calculated in the above section, are: After the solution has been calculated, by using a Fortran routine, the WARP3D output file can be read and interpreted.For every step, the total reaction value Ry in the traction direction and the total number of killed elements Nk are stored.Then, with simple data manipulations, the gross stress S and the half crack length a are pointed out: The resulting curve S(a) has to be compared with the experimental curve in order to perform the tuning of the nucleation parameters.

Nucleation parameters tuning
From the afore cited experimental test data, the following experimental stress vs. crack length curve results (Fig. 9).
By varying the nucleation parameters N, SN, fN, different curves have been calculated .The parameter fN gets down the curve when it is increased; S N has the opposite effect, with a little influence on the shape of the curve; the increasing of N can slightly move the curve maximum value versus the increasing a.

Final results and R-curve
The ASTM specification E561-86 for R-curve determination imposes that the curves have to be expressed in terms of equivalent stress intensity factor and equivalent crack length, in order to take into account for the plastic deformation around the crack tip.For the M(T) central crack specimen, the equivalent SIF is defined by:   eq eq eq eq eq a K sec S a a S a H The equivalent crack length has two alternative definitions: the one, based on the Crack Opening Displacement in the traction direction of the centre of the crack; the other, which has been adopted in the present work, based on the plastic radius evaluation.The plasticity contribution to the equivalent crack length is done by: 2 1 2 Then, the equivalent crack length is defined by: 0 ; Δ eq p eq eq a a r a a a     The plastic radius definition (13) contains the equivalent SIF and, vice-versa, the expression (12) contains the equivalent crack length.It is then necessary to calculate iteratively the two expressions ( 12), (13)(14) up to the convergence is reached.The current plastic radius value is then done, in an iterative formulation by: The convergence is reached when the (i+1)-th value of the plastic radius differs from the i-th value for a small enough quantity.By using the S(a) values calculated before, and after the above descript transformation (S, a)  (K eq , a eq ), we can compare the calculated and experimental R-curves, as reported in Fig. 10.As a result, the calculated R-curve is in very good agreement with the experimental test.In the following pictures, some results from the FE solution are reported.The WARP3D solution data stored during the calculation has been 'translated' in ANSYS commands with a Fortran routine, in order to use its postprocessor module.In Fig. 11 the calculated stress distribution in the traction direction (y) for the half-crack length increasing a = 25 mm is reported for the whole model and near the crack front.It is possible to appreciate the thorough-the-thickness distribution (remember that symmetry conditions have been imposed).In the following Fig. 12, the crack front shape is pointed out.It results from the elements extinction procedure described above.In Fig. 13, the enhanced deformed configuration is reported, for a = 25 mm.In the first picture, the deformation around the crack tip is pointed out and it is possible to appreciate the strain localization effect.In the last picture, the residual strain near the initial crack length is shown.Considering that the crack surface is stress free, the deformation is totally plastic.Further, the difference between the residual strain on the symmetry surface (the visible surface) and the external one is shown, were for the external surface (near the plane stress condition) the plastic strain is much higher.

Further
flow stress (internal variable of the model);  : current matrix equivalent strain (internal variable); f : voids volume fraction (internal variable); N : hardening coefficient;

Figure 4 :
Figure 4: FE model of the single void cell.

Figure 5 :
Figure 5: Strain (x) distribution for the two load cases.

Figure 6 :
Figure 6: Global cell response compared with eq.(1-6) for the two load cases.

Figure 7 :
Figure 7: Global cell response compared with eq.(1) and the coalescence law.

Figure 10 :
Figure 10: Comparison between calculated and experimental R-curve.

Figure 11 :
Figure 11: Calculated stress distribution in the traction direction (y) for a = 25 mm.

Figure 12 :
Figure 12: Crack front shape resulting from elements extinction.

Figure 13 :
Figure 13: Deformed configuration near the crack tip and the initial crack length.

T351 All Three Directions
T 2024A- The best nucleation parameter values are: