A robust approach for the determination of Gurson model parameters

Among the most promising models introduced in recent years, with which it is possible to obtain very useful results for a better understanding of the physical phenomena involved in the macroscopic mechanism of crack propagation, the one proposed by Gurson and Tvergaard links the propagation of a crack to the nucleation, growth and coalescence of micro-voids, which is likely to connect the micromechanical characteristics of the component under examination to crack initiation and propagation up to a macroscopic scale. It must be pointed out that, even if the statistical character of some of the many physical parameters involved in the said model has been put in evidence, no serious attempt has been made insofar to link the corresponding statistic to the experimental and macroscopic results, as for example crack initiation time, material toughness, residual strength of the cracked component (R-Curve), and so on. In this work, such an analysis was carried out in a twofold way: the former concerned the study of the influence exerted by each of the physical parameters on the material toughness, and the latter concerned the use of the Stochastic Design Improvement (SDI) technique to perform a “robust” numerical calibration of the model evaluating the nominal values of the physical and correction parameters, which fit a particular experimental result even in the presence of their “natural” variability.


INTRODUCTION
ecause of the increasing use of the aluminium alloys with high fracture resistance features in the field of the aerospace industry, the study of the techniques for the experimental determination of the fracture toughness and the research of the methodologies to transfer such results to the real structures have been subject of great interest for the scientific community in recent years.In fact, it is well known that the macroscopic parameters of the classical fracture mechanics theory, such as the stress intensity factor (K), the J integral [1], the CTOD (Crack Tip Opening Displacement), the CTOA (Crack Tip Opening Angle) [2], the extension of the plastic zone at the crack tip (r p ) [3][4][5][6][7], cannot be easily transferred from one geometry to another, since they strictly depend on it.Within this scenario, the parallel spreading of advanced techniques of numerical simulation has allowed the development of numerical models with which it is possible to obtain very useful results for a better understanding of the physical phenomena involved in the macroscopic mechanisms of propagation.Among the most promising models introduced in recent years the one proposed by Gurson-Tvergaard (GT) links the B propagation of a crack to the nucleation and growth of micro-voids in the material, and then it is able to connect the micromechanical characteristics of the component under examination to crack initiation and propagation up to a macroscopic scale.The three stages of nucleation, growth and coalescence of micro-voids are well established results of metallographic observation for polycrystalline metals at ductile failure.The simulation of these microstructural damage processes has been considered in various micromechanical and macromechanical model approaches in the literature.A macromechanical model can be exactly obtained by the statistical averaging of microscopic quantities in a homogenization process.The Gurson model [8], which derived a macroscopic yield function and an associated constitutive flow law for an ideally plastic matrix containing a certain volume fraction of spherical voids, is a well-known analytical approach to this problem.Empirical modifications of this approach have been proposed to improve the prediction at low void volume fractions [9] and to provide a better representation of final void coalescence [10].It must be pointed out that even if the statistical character of some of the physical parameters involved in this model has been put in evidence, no effective attempt has been made insofar to relate the corresponding statistic to the experimental results, as for example the R-Curve [1,11].In the first part of the work, after briefly recalling the ductile fracture mechanism and the classical approach to determine the parameters of the GT model, we pointed out the main aspects of the Stochastic Design Improvement (SDI) technique [12], implemented in the commercial code ST-ORM [13], which can be coupled with specialized FEM codes as LS-DYNA [14] or WARP 3D [15] for what concerns the deterministic structural part of the process; we then illustrated the use of the SDI technique to perform a robust numerical calibration of the parameters of the GT model.In the second part of the work, some numerical results regarding a simple structural component are shown, pointing out the influence of the statistical distribution of the physical parameters of the GT model on the toughness of the considered material.Material toughness was numerically evaluated both in terms of the critical value of the J-integral, by using the Equivalent Domain Integral method [16] implemented in the finite element code WARP 3D, and in terms of critical applied load, by using the explicit finite element code LS Dyna.In this case, as the load is monotonically increased by applying increments that are considered to be constant in time, the critical load value is univocally linked to a correspondingly crack propagation initiation time.For the validation of the numerical results, experimental data from literature are used [17][18].

DUCTILE FRACTURE MECHANISMS: OVERVIEW ON THE GURSON-TVERGAARD MODEL
s it is well known, a typical process of ductile fracture develops through the nucleation and the following growth and coalescence of cavities nucleated at crack tip.The best known nucleation mechanisms of such cavities are those linked to the presence of hard particles (inclusions or precipitates) in a ductile matrix (base material).In such cases, the brittle failure of a particle or its decohesion from the matrix, because of the interface breaking, leads to the nucleation of small voids, whose growth causes real cavities in the material.It is clear that the toughness of a material increases for smaller dimensions of such particles and for increasing mean distance and homogeneity of the distribution of such inclusions.From an analytical point of view, we can observe that a void nucleates when the elastic potential energy associated with the material surrounding the particles reaches a particular value, which can be determined by following one of the two approaches, which, even if rather different from each other, lead to similar or equivalent conclusions.The first one is based on the existence of a critical value of the strain in the direction of the applied load, and the second one on the existence of a critical value of the stress in the same direction [19]; once such critical values are overcome, the amount of elastic potential energy necessary for the voids nucleation is considered as reached.As those critical values for both stress and strain are usually low, the nucleation process develops rather easily; therefore, the ductile failure of metallic materials depends substantially on the way the cavities grow and coalesce, rather than on their nucleation.As a rule, a growing cavity is modelled as an ellipsoid of revolution, whose major axis (2b) is parallel to the direction of the load.As the strain increases, the cavity grows in volume while keeping its minor axis equal to the initial mean dimension R0 of the enclosed particle considered to be non-deformable, and b varying along with the deformation in the load direction, A according to the relation: As we can see, the bigger the initial radius of the particle the larger the cavity will be, for a given strain.
In analogy, if X 0 is the initial average distance between the particles, the current value of the average distance, X, can be written as: this equation states that the applied load leads to a reduction of the distance between the particles, which is added to the increasing effect of the cavity dimensions.
As a consequence of those assumptions, high-localised stresses are produced into the surrounding material, which can lead to a coalescence phenomenon of near cavities if a critical value of strain ( f ) is reached.This critical value is strictly related to the strain hardening coefficient of the material, n: high values of n mean that the strain peaks are transferred to less strained surrounding zones; consequently, the contraction of the material between the cavities becomes less localized, and the coalescence of the cavities delayed.
Once the critical strain is reached, the ratio b/X exhibits a limit value that can be expressed as follows: which, rearranging the terms, leads to an equation, which let us evaluate the critical strain: In the previous equation the critical strain, f  is implicitly showed as dependent on the initial volume fraction of the particles, f0.
Rice and Johonson [20], assuming that the crack tip is located at distance Xo from the nearest cavity with radius R 0 , formulated a bi-dimensional model for the localised strain and the coalescence process of cavities.
The crack tip opening displacement, at the conditions for which the coalescence between the growing cavity and the plastically deformed crack tip occurs, represents the CTOD i value at which the crack starts to growth.Intuitively, for a fixed space between the particles, the CTOD i is bigger if the volumetric ratio is smaller.Experimentally, it has been found a relation like CTODi/X0 = cost • εf and, according to some detailed finite element analyses, as well as experimental tests, that constant can be assumed to be function, g, of the square of the hardening coefficient of the material, n [21].Therefore, the following equation can be considered: According to this expression, the value of CTODi could be converted in an equivalent value of fracture toughness KIc, by means of the well-known relationship: where  represents the yielding stress of the material; from the Eqs.( 4) and ( 6) follows: The above equation states that in the limit state preceding the coalescence of voids and, therefore, the crack growth, the fracture toughness K Ic of the material depends on the power (-1/6) of the initial volumetric void fraction, f 0 , and on the basic material parameters (E,  , n) as well as on the initial dimension of the particles, R 0 .
Based on an equivalent approach, the GT model provides a constitutive relation for the increasing of the volumetric void fraction up to the starting of the coalescence phenomenon.The coalescence phenomenon is not accounted for in the model, as the basic assumption of the model theory is the homogeneity of the strain field, and therefore it does not allow reliable prediction in the case of localised strain phenomena.This limit of the Gurson model can be overcome if a critical value of the cavities volume fraction is introduced; the coalescence phenomenon starts when the actual volume fraction of cavities reaches this critical value, as explained below.
Fortunately, it was possible to find a strong dependency of this critical volume fraction, f c , from its initial value f 0 and therefore from the cavity radius, R0, and from the distance from the crack tip, X0.
The original Gurson model considers a rigid sphere of perfectly plastic material, encapsulating a spherical cavity under a homogenous strain.The derived model that can also consider material hardening, called the Gurson-Tvergaard model [10], can follow the growth of the volume fraction of existing cavities from f0 to fc and can be written as: where q 1 ,q 2 and q 3 are correction parameters introduced by Tvergaard [9] (in order to take into account the hardening behaviour of material), σ m is the average normal stress, σ is the equivalent Von Mises stress; with such model it is possible to predict the cavity growth ratio in the plastic field.
For real components, where many inclusions are to be found, the GT model considers that the volume fraction of voids increases over an increment of load because of both the continuing growth of existing voids and the nucleation of new voids; this is taken into account by increasing the void growth rate linked to the actual plastic strain flow at crack tip, (1 , by a quantity, ( ) Therefore, the introduced nucleation parameters are:  f N : particles volume fraction from which the nucleation of cavities can take place;   N and S N : mean value and standard deviation of the supposedly normal distribution of the critical strain value at which the nucleation takes place.
As it is possible to observe, the nucleation acceleration, driven by the parameter A, is limited to an internal strain condition in which the actual equivalent strain  is quite close to N  .Anyway, the nucleation phenomenology is of complex understanding and, even though some results from metallographic analyses about particle distribution and concentration are available, it is not sure if voids nucleate from all of them; this implies that the particles total volume fraction obtained by the metallographic analyses cannot be directly used as parameter in the GT model [22].Within this work, the values of the physical parameters of the GT model related to the nucleation mechanism were calibrated by adapting numerical and available experimental results by means of an iterative process that is described in the following pages.This crack growth mechanism leads to a numerical model characterised by cubic elements parallel to the crack plane; all those elements, which are equal in dimensions, exhibit material properties according to the GT model.Each element represents an elementary cell of material containing an initial cavity volume fraction, f o ; the cell dimension, D o , is comparable with the distance between the largest inclusions detected by microscopic analyses at the interior of the bulk material.According to the GT model, the volumetric fraction of the cavities increases with the stress-strain state; because of the induced softening of the surrounding material the critical value of cavity volume fraction, fc, can be reached and, in this case, the coalescence starts.From a numerical point of view, when the coalescence of the cavities starts, the cell elements reduce their ability to equilibrate the applied load, according to a given parameter,  , which can be expressed as follows: where D is the effective cell dimension, 0 D is the cell dimensions at the time at which f = f c and  is a numerical parameter.
From the considerations above, it follows that the GT model must be calibrated for every chosen material if it wants to be considered as a valid tool for fracture mechanics analysis within the non linear field.

CALIBRATION OF THE GT MODEL PARAMETERS.
he numerical calibration of the GT model parameters can be carried out with a two-step procedure [23][24].The first stage involves, on a micro-mechanical scale, the correction parameters, q1, q2 (q3=q1 2 [15]), the fixing of the physical parameters (, f c ) and of the geometric ones (D o , f o ) on the basis of experimental considerations.The second stage involves, on a macro-mechanical scale, the geometric parameters (Do, fo).For what concerns the parameters governing the nucleation mechanism of cavities, as already noted, they can be initially roughly determined by experimental data related to the initial distribution of the particles in a material sample, while their final refined values can be determined in the second stage of the calibration process.Beside these parameters, it is obviously necessary to know the material properties of the base material (E , n), or the relationshipfrom a standard tensile test.The validity and the reliability of the GT model are strongly dependent on the correct determination of all the aforesaid parameters and material properties, which can be summarized as follows: elastic-plastic material model parameters: n,    ; -Tvergaard correction parameters: q 1 , q 2 , q 3 ; nucleation process parameters: f N ,  N ,S N ; initial dimension of the elementary cell: D0; initial and critical cavities volume fraction: f0, fc; parameter .
In the present work, the first stage of the calibration was performed by comparing two different numerical models: the first one (reference model) was built using ANSYS code.The model (Fig. 1) was a cubic three-dimensional cell, with an initial characteristic dimension D o of the same order of magnitude of the distance between the largest inclusions (100 m), having the same material properties as the base material with a spherical cavity in the centre.A second FEM model (developed with WARP 3D v.14.2 code) consists of a unique 8-noded cubic finite element, whose material shows a constitutive law based upon the GT model, with parameters relative to the base material and to the cavity volume fraction equal to those adopted in the previous reference model.In the first stage of the calibration the stress-strain curve, obtained from the reference model under given boundary conditions, was compared with the curve obtained in the same conditions by the second model, on the basis of the maximum stress and of the work needed to deform the cell [23][24].All the combinations of q 1 , q 2 e q 3 that give the smallest difference of the said quantities in the second model, with respect to those obtained from the reference model, can be used as parameters of the GT model for the considered material.The second stage of the calibration process enabled to evaluate the parameters D0 and f0, as well as the nucleation parameters f N ,  N and S N , by fitting the numerical R-curves to the experimental one.
Anyway, metallographic analyses and numerical calculations suggest that an indicative value of D0 can be given by the Crack Tip Opening Displacement (CTODi) at the beginning of the propagation process [23][24], which depends on the average distance between the largest inclusions in the material.This value of the CTOD can be measured with standard static propagation tests carried out on a CT specimen [25].
T Figure 1: Elementary cell FE model.

THE SDI PROCESS
set of techniques was introduced in order to obtain a robust design, i.e. one whose behavior is rather insensitive to all variations of the main variables, or, what is the same thing, a design whose statistics are characterized by the smallest standard deviation, as a function of the statistics of input.This approach can be also linked to another very relevant question; the result of an experimental test carried on an assigned structure is the consequence of the particular and real values of all design variables, whose density functions are supposed to be known: when we try to correlate the test results to a numerical simulation procedure, we want, in effect to assess all other aspects stated, which are the values that the design variables had in the real structure tested in that particular experiment.Both these problems can be effectively dealt with by an SDI (Stochastic Design Improvement) [12] process, which is carried out by means of several MC series of trials (runs) as well as of the analysis of the intermediate results.In fact, input -i.e.design variables x -and output -i.e.target y -of a mechanical system can be connected by means of a functional relation of the type Which, in the largest part of the applications, cannot be defined analytically, but only ideally deduced because of its complex nature; in practice, it can be obtained by considering a sample xi and examining the response yi, which can be carried out by means of a simulation procedure and, first of all, by means of one of M-C techniques, as recalled above.
Considering a whole set of M-C samples, the output can be expressed by a linearized Taylor expansion centered about the mean values of the control variables, as where  i represents the vector of mean values of input/output variables, and where the gradient matrix G can be obtained numerically, carrying out a multivariate regression of y on the x sets obtained by M-C sampling.If y0 is the required target, we could find the new x0 values by inverting the relation above, i.e. by   as we are dealing with probabilities, the real target is the mean value of the output, which we compared with the mean value of the input, and considering that, as we shall illustrate below, the procedure will evolve by an iterative procedure, it can be stated that the relation above has to be modified as follows, considering the update between the k-th and the (k+1)-th step: The SDI technique is based on the assumption that the cloud of points corresponding to the results obtained from a set of MC trials can be moved toward a desired position in the N-dimensional space such as to give the desired result (target), and that the amplitude of the required displacement can be forecast through a close analysis of the points that are in the same cloud (Fig. 2): in effects, it is assumed that the shape and the size of the cloud don't change greatly if the displacement is small enough; it is therefore immediate to realize that an SDI process is composed by several sets of MC trials (runs) with intermediate estimates of the required displacement (Fig. 3).It is also clear that the assumption about the invariance of the cloud can be maintained just in order to carry out the multivariate regression which is needed to perform a new step -i.e. the evaluation of the G matrix -but that subsequently a new and correct evaluation of the cloud is needed; in order to save time, the same evaluation can be carried out every k steps, but of course, as k increases, the step amplitude has to be correspondently decreased.
It is also immediate that the displacement is obtained by changing the statistics of the design variables and, in particular, by changing their mean (nominal) values, as in the now available version of the method all distributions are assumed to be uniform, in order to avoid the gathering of results around the mode value.It is also pointed out that sometimes the process fails to accomplish its task because of the existing physical limits, but in any case SDI allows to quickly appreciate the feasibility of a specific design, therefore making easier its improvement.Of course, it may happen that other stochastic variables are present in the problem (the so called background variables): they can be characterized by any type of statistical distribution included in the code library, but they are not modified during the process.Therefore, the SDI process is quite different, for example, from the classical design optimization, where the designer tries to minimize a given objective function with no previous knowledge of the minimum value, at least in the step of the problem formulation.On the contrary, in the case of the SDI process, it is firstly stated what is the value that the objective function has to reach, i.e. its target value, according to a particular criterion that can be expressed in terms of maximum displacement, maximum stress, and so on The SDI process gives information about the possibility to reach the objective within the physical limits of the problem and determines which values the project variables must have in order to get it.In other words, the designer specifies the value that an assigned output variable has to reach and the SDI process determines those values of the project variables which ensure that the objective variable becomes equal, in the mean sense, to the target.Therefore, according to the requirements of the problem, the user defines a set of variables as control variables, which are then characterized from an uniform statistical distribution (natural variability) within which the procedure can let them vary, observing the corresponding physical (engineering) limits.In the case of a single output variable, the procedure evaluates the Euclidean or Mahalanobis distance of the objective variable from the target after each trial: where y i is the value of the objective variable obtained from the i-th iteration, y * is the target value and N is the number of trials per run.Then, it is possible to find, among the worked trials, that one for which the said distance gets the smallest value and, subsequently, the procedure redefines each project variable according to a new uniform distribution with a mean value equal to that used in such "best" trial.The limits of natural variability are accordingly moved of the same quantity of the mean in such way as to save the amplitude of the physical variability.If the target is defined by a set of output variables, the displacement toward the condition where each one has a desired (target) value is carried out considering the distance as expressed by: where k represents the generic output variable.If the variables are dimensionally different it is advisable to use a normalized expression of the Euclidean distance: where: But, in this case, it is of course essential to assign weight factors  k in order to define the relative importance of each variable.
Several variations of the basic procedures are available; for example, it is possible to define the target by means of a function which implies an equality or an inequality, too; in the latter case the distance is to be considered null if the inequality is satisfied.
Once the project variables have been redefined a new run is performed and the process restarts up to the completion of the assigned number of shots.It is possible to plan a criterion of arrest in such way as to make the analysis stop when the distance from the target reaches a given value.In the most cases, it is desirable to control the state of the analysis with a real-time monitoring, with the purpose to realize if a satisfactory condition has been obtained.

SENSITIVITY ANALYSIS
n this paragraph a sensitivity analysis of the nucleation parameters and of the main mechanical properties of the material on the material toughness of a thin plate made of aluminium alloy 2024 T3, with a through crack in the middle transversal section, was performed.
The values of the material properties used are:  = 324 MPa,  = 0.0045, n = 9.8; for what concerns the Tvergaard parameters they were set equal to q1 = 1.255, q2 = 1.20, q3 = 1.58, as obtained from the calibration of the model performed in [26]; the other necessary physical parameters values were derived from literature and put equal to f 0 = 0.09, D0 = 0.15 mm, fc = 0.2 and0.1; the first attempt value for the nucleation parameters were set equal to  N = 0.3, S N = 0.5, f N = 0.05, which represents the default values in the algorithms of the WARP 3D code.The dimensions of the plate were 800 mm in the longitudinal direction, 600 mm in the transversal direction and 1.6 mm in the thickness.
The static crack propagation test, from both a numerical and an experimental point of view, were developed by considering a uniform applied longitudinal displacement at both the ends of the plate.The experimental results are available in [18].At first, the experimental R-curve of the plate was compared with a numerical one (Fig. 4), obtained by using the WARP 3D code; the used FE model consisted of 33216 nodes and 23680 8-noded solid elements.The elements linked to the nodes belonging to the crack plane were modelled by considering a material model based on the GT criterion, as implemented in the used code.
As it is shown in the plot of the Fig. 4, in the range of the considered crack elongation, numerical results are in good agreement with the experimental ones: the two R-curves are quite parallel.Further developed numerical investigations show that the distance between the numerical and the experimental R-curves depends substantially on the values considered for the nucleation parameters in the GT model.

I
In order to clarify the influence of the nucleation parameters f N ,  N and S N on the material toughness at the beginning of the crack propagation (output variable), that is, on the start point of the R-curve, a sensitivity analysis was developed; results are reported as a pie-chart in Fig 5.
As it is possible to observe, the most meaningful influence on the output variable is due to the parameter fN; for the variables distribution the following realistic characteristics were assumed: f N  N(0.05; 0.01),  N  N(0.3; 0.06); S N N(0.5;0.1).As it is possible to observe from the plot of the standard deviation versus the MC iteration number, the reached value of the standard deviation of the crack initiation time is about 10.6, which implies a high probability of varying the critical remote stress of about ± 3% of its mean value, given the considered statistics of the nucleation parameters, i.e. the same order of magnitude of the difference between the numerical and the experimental results.
From this considerations, it is possible to conclude that the influence of the calibration of the nucleation parameters is quite low in the reference case, even in case their intervals of variation are rather large.After this, an investigation on the influence of the variability of the mechanical material properties on the same output variable was developed.The scale factor of the stress-strain curve along the stress axes was chosen as the input random parameter, with a normal distribution N (1;0.1).
In the following Figs.9-11 the cumulative density function, the mean and the standard deviation of the crack initiation time, related to the variability of the stress-strain curve, vs. the MC trial number, are respectively reported.In this case, the variability of the critical remote value can reach the 8% of its mean value, which is a rather appreciable difference to be taken into account when we try to fit the experimental results with numerical ones.

CONCLUSIONS
he use of the SDI approach to calibrate the physical and correction parameters of the GT model is an useful numerical tool that provides a numerical model able to well represent the real development of the considered phenomena even if the involved physical parameters are subjected to natural variability.Within this work we investigated on the feasibility to apply the SDI approach, and we used it to calibrate the physical nucleation parameters of the GT model.The results were satisfactory and the R-curve of the examined component, obtained by using the nominal value of the found nucleation parameters, is in good agreement with the experimental one.In the last paragraph, some sensitivity analyses showed that the natural scatter of the material properties and of the nucleation parameters around their "nominal" value do not influence considerably the toughness of the material in the considered reference case.

Figure 2 :
Figure 2: initial and final structural responses.

Figure 5 :
Figure 5: Sensitivity analysis results.The commercial code used to develop the sensitivity analyses is the St-Orm code, linked to the LS-Dyna code for what concerns the deterministic structural results.In this code the GT model can be related only to shell type elements under plane stress conditions.The cumulative density function (CDF), the mean and the standard deviation of the output variable obtained for an increasing number of MC trials, are respectively reported in the plots of Figs.6-8.The mean and the standard deviation of the crack initiation time show quite constant values after about eighty MC iterations.That reached mean value corresponds to a critical remote stress of 74 N/mm 2 , which differs from the experimental results of a few percent.As it is possible to observe from the plot of the standard deviation versus the MC iteration number, the reached value of the standard deviation of the crack initiation time is about 10.6, which implies a high probability of varying the critical remote stress of about ± 3% of its mean value, given the considered statistics of the nucleation parameters, i.e. the same order of magnitude of the difference between the numerical and the experimental results.From this considerations, it is possible to conclude that the influence of the calibration of the nucleation parameters is quite low in the reference case, even in case their intervals of variation are rather large.After this, an investigation on the influence of the variability of the mechanical material properties on the same output variable was developed.

Figure 6 :
Figure 6: CFD of the crack init.time [ms]

Figure 7 :
Figure 7: Mean of the crack init.time [ms] vs. MC trial.

Figure 8 :
Figure 8: Std dev. of the crack init time [ms] vs. MC trial.

Figure 10 :
Figure 10: Mean of the crack init.time [ms] vs. MC trial.