Numerical investigation and experimental validation of physically based advanced GTN model for DP steels

Materials Science & Engineering A 569 (2013) 1–12

Contents lists available at SciVerse ScienceDirect

Materials Science & Engineering A journal homepage: www.elsevier.com/locate/msea

Numerical investigation and experimental validation of physically based advanced GTN model for DP steels Joseph Fansi a,b,c,n, Tudor Balan b, Xavier Lemoine b,c, Eric Maire d, Caroline Landron d, Olivier Bouaziz c,f, Mohamed Ben Bettaieb e, Anne Marie Habraken a a

University of Lie ge, Departement ArGEnCo, Division MS2F, Chemin des Chevreuils 1, Lie ge 4000, Belgium Arts et Me´tiers ParisTech, LEM3, UMR CNRS 7239, 4 rue A. Fresnel, 57078 Metz cedex 03, France c ArcelorMittal R&D Global Maizie res S.A., voie Romaine, Maizie res-Le s-Metz 57238, France d INSA de Lyon, MATEIS CNRS UMR5510, 7 Avenue Jean Capelle, Villeurbanne 69621, France e Ensicaen, 6 Boulevard du Mare´chal Juin, 14050 CAEN Cedex 4, France f Ecole des Mines de Paris, Centre des Mate´riaux, CNRS UMR 7633, BP 87, Evry Cedex 91003, France b

a r t i c l e i n f o

abstract

Article history: Received 23 November 2012 Received in revised form 11 January 2013 Accepted 12 January 2013 Available online 23 January 2013

This numerical investigation of an advanced Gurson–Tvergaard–Needleman (GTN) model is an extension of the original work of Ben Bettaiebet al. (2011 [18]). The model has been implemented as a user-defined material model subroutine (VUMAT) in the Abaqus/explicit FE code. The current damage model extends the previous version by integrating the three damage mechanisms: nucleation, growth and coalescence of voids. Physically based void nucleation and growth laws are considered, including an effect of the kinematic hardening. These new contributions are based and validated on experimental results provided by high-resolution X-ray absorption tomography measurements. The current damage model is applied to predict the damage evolution and the stress state in a tensile notched specimen experiment. & 2013 Elsevier B.V. All rights reserved.

Keywords: Damage GTN model Nucleation Coalescence DP steel Tomography

1. Introduction Many models for material degradation, commonly called damage, have been proposed in the past. The damage modeling is usually divided in two groups: the phenomenological and the micromechanical approaches. The phenomenological approaches initiated by Kachanov [1] coupled the elastoplasticity theory with damage in the framework of the so-called continuum damage mechanics. This approach assumes the existence of a classical true stress tensor s computed from macroscopic loading and macroscopic area measurements, and an effective stress tensor s effec theoretically closer to the actual average microscopic stress state existing between defects. The effective stress tensor s effec can be related to the true stress s by an effective stress operator M , depending on a damage parameter D which characterizes the state of damage of the material.

s ef f ec ¼ M ðDÞ : s

ð1Þ

The operator M (D) takes into account the area of the microvoids and microcracks, stress concentrations due to microcracks n Corresponding author at: University of Lie ge, Departement ArGEnCo, Division MS2F, Chemin des Chevreuils 1, Liege 4000, Belgium. E-mail addresses: [email protected], [email protected] (J. Fansi).

0921-5093/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.msea.2013.01.019

and the interactions between neighboring defects. The best known contribution on this type of approach is the work of Lemaitre [2]. In his basic model, the isotropic damage depends on a scalar variable D, neglecting the microcracks orientation. He defines the effective stress by the following relation:

s ef f ec ¼

1 s 1D

ð2Þ

Cordebois et al. [3] modified the isotropic damage parameter D (Eq. (1)) to a second order tensor D to introduce the anisotropic damage. An elastic energy equivalence hypothesis is introduced and effective values for stress and strain tensors are defined. In the micromechanical approach, initially introduced by Gurson [4] for ductile fracture, the yield function in Eq. (3) is modified by the presence of voids in the material; when the porosity f starts to increase, material softens and loses its capability to carry loads:   s2eqv 3sm 2 F p ¼ 2 þ 2f cosh  ð3Þ 1f ¼ 0 2sy sy where

 seqv designates the macroscopic von Mises equivalent Cauchy stress.

2

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

 sy is the yield stress of the fully dense matrix, initially equal to s0 (initial yield stress).  sm ¼ ð1=3Þtrðs Þ corresponds to the macroscopic mean stress.

 f defines the void volume fraction or porosity. In its original form, the Gurson model represents only the void growth stage. It is based on a representation of the voided material which consists of a hollow sphere surrounded by the matrix. The latter is assumed to be rigid-plastic and isotropic. It is described by a standard von Mises yield criterion and an associated flow rule. In order to minimize the differences between the model and finite element cell results, Tvergaard et al. [5] introduced new damage parameters (q1, q2, q3) into the expression of the yield surface Eq. (4). Later, Tvergaard et al. [6,7] introduced the nucleation and the coalescence stages. These extensions of the initial Gurson model led to the so-called GTN model (standing for Gurson–Tvergaard– Needleman) Fp ¼





s2eqv 3q2 sm 2 þ 2q1 f cosh 1q3 f ¼ 0 2sy s2y

ð4Þ

with q1 ¼ 1:5, q2 ¼ 1, q3 ¼ q1 2 . One limitation of the GTN model is to consider the voids as spherical. In order to overcome this limitation, Gologanu et al. [8] extended the GTN model by taking into account the void shape effects. Another limitation of the GTN model is to consider the matrix isotropic in the yield function. Observations have proved that plastic anisotropy of the matrix surrounding the voids can influence the effective stress function and damage evolution. As a consequence the isotropic equivalent stress is often replaced by an anisotropic function adopting the Hill 48 criterion. Benzerga et al. [9] introduced two material parameters in the yield function to reflect the plastic anisotropy. An additional significant extension is the hardening behavior. The initial Gurson model describes a perfect plastic behavior. Leblond et al. [10] extended in an accurate theoretical way the model to the case of a matrix with isotropic hardening. Their model introduces two new parameters in the yield function representing the plastic hardenable hollow sphere. These parameters are depending on the deviatoric and hydrostatic parts of the macroscopic plastic strain. Ne gre et al. [11] used the Hollomon law to describe the yield stress which is a function of equivalent plastic strain. Mear and Hutchison [12] extended the Gurson model to the kinematic hardening case. Other contributions use a yield surface function of plastic strain rate and allowed the introduction of the viscoplasticity in the model. The improvement of microstructure observation method enhanced the void nucleation law from the initial GTN model. Recently Helbert et al. [13] and Bouaziz et al. [14] with help of X-ray tomography method accurately modeled the voids nucleation stage. Pardoen [15] improved the void growth modeling and the Thomason coalescence criterion [16]. Due to its micromechanical roots and to the explicit use of the void volume fraction as a state variable, the GTN model has been chosen in the current research to introduce results from recent experimental X-ray tomography measurements. The extended GTN damage model presented in this paper was developed a few years ago by Ben Bettaieb et al. [17–18], who focused on dual-phase steel behavior and introduced the kinematic hardening and the plastic anisotropy of the matrix in the Gurson yield function. Their second contribution was to integrate experimental nucleation/growth observations in the GTN model. Through this approach, they also introduced a more accurate way to compute the density of voids and the evolution of the mean radius. In that way, the usual GTN damage parameters (q1, q2, q3) and the porosity become variables explicitly depending on triaxiality and plastic strain. Validation of Ben Bettaieb’s model was performed on experimental results for tensile specimens with square cross-section and very large notch radius.

Recently, further X-ray tomography measurements have been carried out and investigated by Landron et al. [19] on in-situ tensile notched axisymmetric specimens of DP steels. The experiments revealed a strong dependency between the density of voids, the backstress, and the triaxiality for these grades. Motivated by these new experimental observations and the industrial needs, an extension of the Gurson–Tvergaard–Needleman–Ben–Bettaieb (GTNB) model has been developed. The main goal of this paper is to correlate the experimental results on a notched specimen and the model predictions. To reach this goal, the GTNB model has been adapted as ‘‘User-defined Material model subroutine’’ (VUMAT) in the Abaqus/ explicit FE code [20–21]. The model has been enriched by adding a coalescence model and a recent void nucleation law integrating the backstress variable [19]. These enhancements have been done based on X-ray tomography observations and measurements. In order to accurately correlate the finite element predictions with the experiments, a post-processing method has been developed to compare consistent results between experiments and simulations. The numerical void volume fraction definition is the same as the one used in the test related to the number of cavities and their size. After this introduction, the second section of this paper introduces the model, which is an improved version of the GTNB model. The tensile experimental tests on notched samples are presented in Section 3. The fourth section describes the identification of the model parameters for the selected DP steel. The fifth section presents the finite element (FE) simulation of the tensile test. Preceding the conclusions, the last section is dedicated to the comparison between the measurements performed by the highresolution X-ray absorption tomography method and the FE predictions.

2. Improved GTNB model The total macroscopic strain e is divided into elastic and plastic contributions e e and e p :

e ¼ ee þ ep

ð5Þ

Hooke’s elasticity law simply states:

s ¼ C e : ee

ð6Þ

where s is the macroscopic (matrix and cavity) Cauchy stress tensor and C e is the elastic stiffness tensor (chosen isotropic and linear). The anisotropic yield function [18] reads:   s 2eqv 3q s m n n2 F p ¼ 2 þ 2q1 f cosh  2 ð7Þ 1q3 f r 0

ksy

sy

where 



 s ¼ s X is the shifted stress tensor defined as the difference 

  

between the Cauchy stress tensor s and the backstress tensor X of the macroscopic medium (matrix and void). s eqv is the anisotropic equivalent shifted stress (with respect to the quadratic Hill criterion). s eqv is computed by qffiffiffiffiffiffiffiffiffiffiffiffiffiffi s eqv 1=2ðs : H : s Þ where H is the Hill matrix of anisotropy coefficients defined as a function of the Lankford coefficients r0 , r45 , r90 (see Appendix). sy defines the yield stress of the dense matrix, it is a function describing the matrix hardening. s m corresponds to the macroscopic mean shifted stress equal to ð1=3Þtrðs Þ. n f is a function of the porosity (or void volume fraction) f ; its definition will be given in Section 2.2. f is the void volume

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

fraction defined as the ratio between the volume of voids V v and the total volume V m þ V v , f¼



Vv V m þV v

ð8Þ

where V m is the matrix volume. Finally, q1 , q2 and q3 ¼ q1 2 are three material parameters introduced by Tvergaard and Needleman [6]. In most cases, these parameters are kept constant and take the following values: q1 ¼ 1:5 and q2 ¼ 1. In the current model, the values of these parameters evolve during deformation to account for the plastic incompressibility of the matrix in the context of physically based porosity evolution laws [17,18]; more details are given in Section 2.1. The parameter k is reflecting the influence of the plastic anisotropy. Derived by Benzerga and Besson [9] from a micromechanical analysis, it is a function of the Lankford coefficients r0 , r45 , r90 . For isotropic materials, k is equal to 2 (see Appendix).

The plastic flow rule is defined by the normality relationship as associated plasticity was chosen: 8 < l_ ¼ 0 if F p o 0 @F p e_ p ¼ l_ ð9Þ @s : l_ Z 0 if F p ¼ 0 where l_ is the plastic multiplier. Consequently, the material is modeled by the quadratic Hill anisotropic yield function coupled with the mixed hardening law of Ben Bettaieb et al. [17]. The isotropic hardening model is defined by the well known Swift law:   sy ¼ K e0 þ epm n ð10Þ where K, n and e0 are material parameters and epm represents the equivalent plastic strain in the dense matrix. The kinematic hardening is described by a variant of the Armstrong–Frederick saturating model [13], adapted to damaged materials:    n n ð11Þ X ¼ 1q1 f X n ; X_ ¼ C x Ssat e_ p X n epeqv where C x and Ssat are material parameters, e_ p is the macroscopic plastic strain rate and epeqv is its equivalent value defined as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi e_ peqv ¼ 2 e_ p : H 1 : e_ p ð12Þ

3

where A is a material constant and eN is the critical strain value for which nucleation is supposed to start. The parameter eN0 designates the value of this critical strain for pure shear loading. In its current version, the modified GTNB model includes a new nucleation law recently proposed by Landron et al. [23]. This law defines the rate of the numerical void density N by: dN Bs  s N ¼ 1þT ð15Þ sc sX N0 depeqv where B and N0 are material constants and sc is the critical shear stress value that the Martensite/Ferrite interface can support without breaking. The quantities s and X represent the macroscopic stress and backstress, in the context of the uniaxial tensile loading that served to derive the law. In the current work, a straightforward extension of this one-dimensional law is adopted by replacing s and sX by the quadratic Hill equivalent stress seqv ðs Þ and shifted equivalent stress s eqv ðs X Þ, respectively. Finally, the triaxiality is defined as T ¼ sm =seqv . After the nucleation process, void growth takes place. It is important to underline that in the GTNB model the voids are spherical when nucleated and remain at this form during the growth stage. This strong hypothesis allows calculating a mean void radius R for the whole void population. The evolution of the mean void radius R is defined by the Rice and Tracey model [24] modified by Bouaziz [25], Maire [26], and Huang [27] to take into account nucleation and different void sizes:     dR sm 1=4 3 sm 1 dN ¼ a exp ðRR0 Þ ð16Þ R H 2 seqv N depeqv seqv depeqv The second term in Eq. (16) is the reduction of the average radius of the cavities due to nucleation, as compared to the initial Rice and Tracey model [24]. Indeed, it is easy to check that Eq. (16) reduces to the classical Rice and Tracey model when the nucleation rate dN is equal to zero. In Eq. (16) aH is a material constant introduced by Huang [27] in order to fit the model to experimental values. R0 is the initial mean radius of cavities, at nucleation. To complete the model, the porosity f must be related to the numeric void density and to the radius of the voids (supposed identical and spherical), starting from its definition (8) in terms of matrix and void volumes. Indeed, the total void volume is related to the mean void radius and the number of voids per volume unit V v ¼ ð4=3ÞNpR3

ð17Þ

1

is the pseudo-inverse of Hill’s anisotropy matrix (its where H expression is provided in the Appendix). According to Arndt et al. [22], Eq. (11) respects the initial Gurson model approach when the kinematic hardening is used. Finally, the work equivalence principle is used:

s : e_ p ¼ ð1f Þsy e_ pm

ð13Þ

Combining Eq. (8) with Eq. (17), one obtains: f¼

Nð4=3ÞpR3 V m þNð4=3ÞpR3

On the other hand, the plastic incompressibility of the matrix keeps V m constant during the loading: V_ m ¼ 0

2.1. Void nucleation and growth modeling The original GTNB model [18] introduced a void nucleation law based on the metallurgical laws proposed by Helbert [13] for a biphasic titanium and further enhanced by Bouaziz et al. [14] for dual-phase steels. In this approach, the numerical void density N (number of voids per unit volume) is related to the triaxiality T and to the macroscopic equivalent plastic strain epeqv by the following relationship ! ! p p N¼A

eeqv eeqv exp ; eN eN

eN ¼ eN0 expðT Þ

ð14Þ

ð18Þ

ð19Þ

It is noteworthy that in the classical Gurson model, this equation governs the void growth stage. An originality of the approach adopted in this work is that both nucleation and growth are governed by physically based equations fitted to the experimentally revealed evolution of the porosity, without invoking plastic incompressibility. Consequently, Eq. (19) was used by Ben Bettaieb et al. [18] to update the value of the q-parameters in the yield function, thus giving them a more physical significance. The details of the equations enforcing the incompressibility condition (19) and its numerical implementation can be found in Refs. [17,18]. Only one parameter, namely q2 , can be determined with this condition. The two other parameters are determined as q1 ¼ 1:5q2 and q3 ¼ q21 according to Tvergaard and Needleman [7].

4

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

current work, the constitutive model has been implemented in Abaqus/Explicit through a VUMAT user subroutine [29].

2.2. Coalescence modeling Tvergaard and Needleman [5–7] used the Brown and Embury [28] original contribution to extend the Gurson model for void growth to account for the coalescence phenomenon and final material failure. The yield function presented in Eq. (7) is only valid for void volume fraction below a critical value f c and it is n modified when f Zf c through the function f ðf Þ expressed by 8 8 f > > f of c n < f f f c  < fcþ f c r f of f f f c f > > : f n f f f c : f Zf n

f

ð20Þ

f

where n

ff ¼

1 q1

ð21Þ

and f f is the ultimate value of void volume fraction when the material completely fails. The discrete, incremental constitutive equations are integrated with an implicit integration scheme for all of the state variables, except for the porosity, which is treated explicitly in order to facilitate the implementation and comparison of various nucleation, growth and coalescence laws. The complete set of equations is solved by an iterative Newton–Raphson method [18]. In the

Fig. 1. Optical micrograph of the DP steel’s microstructure [19,23]. Ferrite appears in light gray and Martensite in dark gray.

3. Experimental tests A detailed description of the experiments has been already published in Ref. [23]. We just remind here briefly what is necessary, to understand the experimental quantities that were compared with the model. A rolled 3-mm thick DP sheet steel with 11% of martensite is used for the experimental validation. This material was thoroughly experimentally investigated in Refs. [19,23]; its bi-phased microstructure (with Ferrite in light gray and Martensite in dark gray) can be observed in Fig. 1. In-situ tensile tests have been performed on cylindrical notched specimens by using high-resolution X-ray absorption tomography. The specimen has been mounted on a rotation stage of the tomograph and the applied displacement speed was between 1 and 5 mm/s. Fig. 2(b) shows the notched axisymmetric specimen used for the experiments. Several quantities were measured during the test and used for the validation of the model introduced in Section 2. Most measurements focus on the central cross-section of the specimen, using the classical post-processing techniques initiated by Bridgman [30] for the analysis of tensile test on notched specimens. The total tensile load provided by a load cell is divided by the smallest cross-section area and is reported as the corresponding average axial stress. The number and average radius of the voids are measured in a spatial volume located at the center of the specimen, i.e., at the intersection of the symmetry axis with the plane of the weakest cross-section. The location of this particular material point is automatically updated during the test. In the mean time, the current radius of this minimal cross-section is recorded as well as the local profile of the neck—from which the local neck radius is determined [19]. Due to the small dimensions of the sample and of the testing machine, the grip displacement cannot be reliably used as the driving parameter for the numerical simulation and the comparison to the experiments. Accordingly, the previously enumerated quantities were extracted from the numerical simulation in a manner similar to the experimental treatment, in order to allow a consistent confrontation of the numerical and experimental results.

Fig. 2. Experimental test setting. (a) In-situ X-ray experimental device with 1-mm notched sample [19,23] and (b) the specimen design.

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

The values of the minimal cross-section radius (rsection) and the notch radius (rnotch) (Fig. 3(a)) are extracted from the outer shape of the specimen. These parameters are necessary to calculate the axial stress, the Bridgman approximation for the triaxiality [30] and the axial strain at each loading step. The average stress calculation uses the measurement of the minimal cross-section radius of the specimen (rsection) and the measured force (F) given by the sensor during the tensile test:

saxial ¼

F

with rini and rsection are the initial and the current radii of the minimum cross-section, respectively (see Fig. 3(a)). Specific postprocessing of the FE results was applied in order to use the same definitions of triaxiality, axial stress and strain during the entire test for consistent comparisons. The quantification of void nucleation during the tensile test has been possible by counting the number of cavities in a fixed spatial volume located at the center of the specimen of dimensions 0.3  0.3  0.3 mm3 (Fig. 3(b)) at each step of deformation:

ð22Þ

pr2section

N¼

The Bridgman formula for the triaxiality T was rewritten in terms of the measured rsection and rnotch [19,30]. Its expression is   1 pffiffiffi r T ¼ þ 2ln 1þ section ð23Þ 3 2r notch The average axial strain is defined as an average value over the entire minimal cross-section

eaxial ¼ ln

5

r 2ini 2 r section

ð24Þ

Number of cavities 0:027mm3

ð25Þ

The total volume of the detected cavities is used to calculate the experimental porosity f at each step of deformation, as well as the mean void radius R corresponding to the assumption of identical spherical voids. Landron et al. [31] underlined that the number of detected cavities depends on the resolution of the measurement, as well as the mean void radius R, while the measured porosity (which is the quantity that actually enters the elasto-plastic model) is almost independent of this resolution.

Fig. 3. (a) The geometry parameters required for axial stress and triaxiality computation (using Bridgman’s method). (b) Studied spatial volume for the porosity measurement at the center of the specimen, of dimensions 0.3  0.3  0.3 mm3.

Table 1 Material parameters for the developed model, corresponding to the considered DP steel. Elasticity E (GPa) 210

Anisotropy (Lankford coefficients)

n 0.35

r0 0.85

Isotropic hardening, Eq. (10)

r45 1.04

r90 0.94

Kinematic hardening, Eq. (11)

K (MPa)

n

e0

Cx

Ssat (MPa)

891

0.245

0.02

92.04

58.02

Nucleation and growth parameters Eqs. (14)–(16) aH B (mm  3) R0 (mm)

sc (MPa)

N0 (mm  3)

A (mm  3)

eN0

0.0021

0.8

1100

1300

4500

Coalescence parameters Eq. (20) and initial q-values fc ff f0

0.55

q1(t ¼0)

q2 (t ¼ 0)

q3 (t¼ 0)

2  10  5

1.5

1

2.25

0.001

4500

0.15

6

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

4. Identification of the model for the studied DP steel Table 1 summarizes the whole set of material parameters required to define the GTNB model. The first and the second categories of parameters correspond to the elasticity and the plasticity modeling. They have been identified by using uniaxial tensile tests, as well as monotonic and Bauschinger shear tests. Three tensile tests were performed in the rolling, diagonal and transverse directions, in order to determine the anisotropy parameters. The shear and reverse shear tests were performed in the rolling direction. The reverse shear tests served for the identification of the kinematic hardening parameters, while the monotonic shear test allowed identifying the isotropic hardening parameters over a wider strain range as compared to the tensile test. Fig. 4 compares the predicted and the experimental stress– strain curve for the tensile test in the rolling direction. The third set of parameters concerns the void nucleation and growth. As already mentioned, B, R0 and N0 are obtained by fitting Eq. (16) to the experimentally measured evolution of the number of cavities for smooth and notched samples [19]. The value of aH has been identified during

800

Axial stress [MPa]

700 600 500

the quantification of the cavities growth step [19]. An average value of the critical shear stress sc has been obtained by calculating the value of the stress at the Martensite/Ferrite interface at the experimentally observed nucleation strain (eN ) for a smooth specimen (eN ¼ 0.17) and notched specimen (eN ¼0.02) [23]. The nucleation strain represents the moment when the void density started to increase. The initial porosity (f0) was directly measured by X-ray tomography. Concerning the coalescence parameters (fourth part of Table 1). The critical value of void volume fraction (f c ) was measured by X-ray tomography, when the first coalescence phenomena were observed [19,31]. The ultimate value of void volume fraction (f f ) and the initial values of the Tvergaard– Needleman damage parameters (q1 , q2 , q3 ) were estimated from the literature [6,9]. In order to perform axisymmetric simulations, the Lankford coefficients defined for the Hill 48 criterion [18] have been set equal to 1 to obtain the von Mises hypothesis, which is a reasonable approximation for the relatively isotropic DP steels, such as the one studied here. As already mentioned, Landron et al. [31] have noticed a strong influence of the X-ray diffraction resolution on the number and average radius of the measured voids, and accordingly on the corresponding material parameters B and N0. However, the ratio (B/N0) appears to be almost insensitive to the resolution of the experimental equipment or to the sample geometry used for the experiment. Accordingly, we have considered this ratio as a single material parameter equal to 3.46 for all of the simulations.

400 Experiment Simulation

300

5. Finite element simulation of the tensile test

200 100 0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

0.18

Axial strain Fig. 4. Axial stress–strain curves comparison between the GTNB model prediction and the experiment on a uniaxial tensile test.

The in-situ tensile test presented in Section 3 has been simulated with the improved GTNB model (Section 2). The mesh generation was constrained by the specific post-processing procedures developed in order to extract from the numerical simulation the physical quantities corresponding to those obtained from the in-situ testing. The need for accurate average values across the weakest cross-section, within a fixed space volume at the root zone of the sample (Fig. 3(b)), and at the central zone of the neck,

Fig. 5. Finite element model of the tensile test on the notched axisymmetric sample. (a) Boundary conditions and velocity field applied; (b) mesh design zoomed on the sub-volume measurement area.

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

imposed a refined mesh at each of these locations. Thus an adapted mesh design has been built up as shown in Fig. 5. A very fine element size of 10 mm  10 mm is adopted in the neighborhood of the central point of the sample. Thanks to this choice, a large number of integration points were available in the subvolume of interest, for accurate average values to be obtained. A relatively small element size is maintained across the entire cross-section corresponding to the symmetry plane of the sample, as well as at the central part of the neck. The nodes of the three central-most elements in the neck region were used to determine the notch radius [30], r notch (see Fig. 3(a)). Homogeneous velocity boundary conditions were applied on the top surface, along with usual symmetry boundary conditions in axial and radial directions. The finite element type chosen in Abaqus/explicit [29] software is a reduced integration linear axisymmetric element with four nodes and one integration point (CAX4R). The mesh design shown in Fig. 5 will be called ‘adapted mesh’ hereafter. Fig. 6 shows the isovalues of axial stress component, triaxiality factor (T), numerical density of voids per mm3 (N), mean void radius (R) and porosity (f). These isovalues are shown at the moment when the maximum plastic strain is close to unity in the minimal cross-section. The damage variables and the axial stress are maximum at the center of the necking section due to the high concentration of the plastic deformation and the triaxiality.

5.1. Post-processing of FE values Fig. 6 underlines a strong heterogeneity in the cross-sections of the sample. In order to compare the experimental results with the numerical predictions at different levels of the axial strain in the minimum cross-section, average values of specific scalar fields

7

of interest were extracted from the heterogeneous field distributions. The weighted average method [32] was used to extract such average values over the central sub-volume and over the crosssection. Due to the choice of axisymmetric element configuration the central sub-volume is cylindrical rather than hexahedral (Fig. 3(b)), with its half-height and radius equal to 0.15 mm. The loss of accuracy due to this approximation can be considered negligible as the experimental measurements [19] show that the void density is sufficiently homogeneous within the measurement domain. Thus, its exact shape and dimensions would not influence the result. The average values of f, N, R and in general for any scalar field a are calculated with this expression Pnelement

Si ai

1 aav ¼ Pi ¼nelement i¼1

Si

ð26Þ

where ai is the value of the scalar field a in the finite element i in the domain. Si represents either the area of the cross-section that corresponds to the finite element i, or for the volume average calculations Si designates the element volume. Surface averages are performed for the calculations of the axial stress and strain, while volume averages are used for the porosity-related quantities. It is noteworthy that the set of finite elements that contribute to the volume average is evolving in time, since the average is computed over a fixed spatial domain—not a material domain. The area of the necking section is calculated by measuring the radius of the minimal cross-section (rsection) during the testing. In order to compare the triaxiality to the experiments, Eq. (24) is also applied to the numerical results. For this, the notch radius rnotch from FE results is required at different strain levels. Three nodes on the notch radius that are the closest to the symmetry

Fig. 6. Simulation results isovalues when maximum plastic strain is near unity: axial strain, triaxiality, axial stress, numerical void density, porosity and mean void radius.

8

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

plane are chosen and used to fit a circle whose center (a, b¼0) is retrieved by the circle equation (Eq. (27)), qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 r notch ¼ ðxaÞ2 þ ðybÞ ð27Þ where rnotch (in plane X-Y) evolves during the tensile test, and (x,y) are the coordinates of the three nodes. Using rsection and rnotch, the average values (saxial , eaxial , T) can be computed in a similar way as in the experimental results. 5.2. Mesh size influence In order to verify the mesh convergence for the simulation of the notched tensile test of interest, six different element sizes (see Fig. 7) have been compared. The mesh size element influence has been conducted by comparing the distribution and the evolution of the local axial strain, porosity (f), density of voids (N), and mean void radius (R). Fig. 8 shows the distribution of the axial strain and of the void radius as a function of the distance from the center (x/rsection ¼0) to the notched side (x/rsection ¼1) of the current necking section. Fig. 9 shows the distribution of the porosity at two different loading steps. Except for the coarsest mesh (Fig. 7(f)), almost no influence of the mesh size has been detected during the test simulations. The local axial strain distribution (Fig. 8(a)) is quite

homogeneous along the cross-section for all deformation steps. A small heterogeneity of the mean void radius is observed (Fig. 8(b)), with the maximum value located close to the symmetry axis. One can note that the maximum is not exactly located at the center of the specimen but at about 1/10 from the center, probably due to the numerical treatment of the axial symmetry. The analysis of the porosity (f) in Fig. 9 points out a mesh sensitivity after 0.5 of axial strain. This deformation level corresponds to the beginning of coalescence according to Landron et al. [19,23,31]. The graphic at 0.75 of deformation shows that the values of f decrease when the mesh size increases. The heterogeneity distribution becomes stronger when the deformation increases. The maximum value is located close to the center (x/rsection ¼0) for low axial strain, but moved at around 0.1 from the symmetry axis when the axial strain reached 0.75. In summary, a strong heterogeneity of the distribution of the damage variables (porosity f, density of voids N) through the sample cross-section was observed. The element size had no influence on the axial strain and equivalent void radius distributions, except for the 167-mm element size which is clearly too coarse to provide accurate results. The damage state variables such as the porosity (f) and the density of voids seem to be more mesh sensitive. As stated earlier, in the current study the element size choice has been restricted by the experimental conditions to densities which verify by far these numerical considerations.

Fig. 7. Element sizes used during the mesh size influence sensitivity on notched tensile simulation: (a) 15 mm, (b) 25 mm, (c) 35 mm, (d) 45 mm, (e) 50 mm, (f) 167 mm.

0.80 0.0024

Axial strain

0.60 0.50 0.40

Adapted 15µm 25µm 35µm

0.30

45µm

Mean voids radius (R) [mm]

0.70 0.0022 0.0020

Adapted 15µm 25µm

0.0018 35µm 0.0016

45µm

0.0014

50µm

0.10

0.0012

167µm

0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.0010 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.20

50µm 167µm

X/Rsection

X/Rsection

Fig. 8. Element mesh size influence and heterogeneity at an average axial strain of 0.75: (a) on local axial strain (eaxial ), (b) on mean void radius (R).

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

0.025

0.025

Adapted 0.020

15µm

0.020 167µm Porosity f

Porosity f

25µm 0.015

35µm 45µm

0.010

0.015

Adapted 15µm

0.010

50µm

25µm 35µm

167µm

0.005

9

0.005

45µm 50µm

0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 X/Rsection

0.000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 X/Rsection

Fig. 9. Element mesh size influence and heterogeneity on porosity (f): (a) at an average axial strain of 0.5, (b) at an average axial strain of 0.75.

1400

σAxial in MPa

1200 1000 800 600

Experiment

400

Complete GTNB (geometrical Method)

200

Complete GTNB (weighted average Method)

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

εAxial Fig. 10. Axial stress–strain curves comparison between the GTNB model prediction and the experiment. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

6. Comparison of simulation results and experiments Now that the measurements and the output quantities extraction from the finite element analysis have been presented, this section is focused on comparing the experimental and the simulation results. The average stress–strain response and triaxiality are first analyzed, preceding the confrontation of the experimental and calculated damage evolutions. 6.1. Stress–strain response and triaxiality Fig. 10 presents the comparison between the experimental and numerical average stress values. The blue curve is the experimental response from the in-situ tensile test, and the red curves are the numerical results. The dashed curve shows the stress value computed with Eq. (22) in a similar way as in the experimental results. The red solid curve is the stress response using the weighted average method in Eq. (26). The comparison underlines that the numerical stress values lie below those obtained by the experiment. The experimental points seem to fluctuate slightly under 0.3 of axial strain and higher stresses are measured as compared to the standard tensile test shown in Fig. 4. The reduced size of the sample and possible geometrical inaccuracies due to its particular shape, along with the use of a miniature insitu machine may explain the differences between the two stress measurements—besides the stress heterogeneity within the

notched specimen. The two numerical averaging methods also show a little stress difference, except for the final stages where the differences become very large. This observation shows the importance of taking into account the element volume weight in the sub-volume quantification and to consider a sufficiently large number of elements for the average calculation. Fig. 11(a) compares the evolution of rnotch obtained by measurement and the finite element simulation. The notch radius from experiment and finite element simulation are in good agreement, a maximum of 11% difference can be detected when the axial strain reaches 50%. This difference could be due to the material behavior model, but also to the approximate description of the notch radius with only three nodes of the finite element mesh. It is also important to recall that the experimental values are averaged over several tested specimens, where differences of the same order were observed between the different experiments [19]. The evolution of the triaxiality approximated with Bridgman’s formula (Eq. (23)) is shown in Fig. 11(c). Computed with the current GTNB model, it is above the experimental average triaxiality of several specimens. The gap between the experimental and the simulation results is linked with the differences in the rnotch measurements and predictions. The radius imperfection of the specimen during the cutting process cannot be reproduced in the axisymetric (2D) simulation and some discrepancies happened within real rnotch. These preliminary results, concerning standard mechanical quantities, put into perspective the accuracy that should be expected between the predicted and experimental results in the context of such large strain heterogeneous tests confronted to advanced constitutive models. 6.2. Porosity, void density and mean radius Figs. 12–14 show the evolution of the damage-related variables of the model with the axial strain—which plays the role of the loading parameter for all the comparisons to the experimental results. Simulations are performed with both nucleation equations presented in Section 2, i.e. Eq. (14) and Eq. (15). Fig. 12 shows the experimental and predicted evolutions of the numerical density of voids. The two predictions lay relatively close to each other, and both approximate well the experimental curve. It is noteworthy that the parameter identification for Eq. (14) was performed using smooth tensile samples [18], while both smooth and notched samples were used for the parameter identification of Eq. (15) [23]. Consequently, further work would be required to

10

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

1.2

Rnotch in mm

1

Fitted_From_experimental measurement Fitted_From_finite_Element

0.8 0.6 0.4 0.2 0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

εAxial

Fig. 11. Comparison of rnotch (a) and rsection (b) evolution obtained from measurement and finite element simulation during the uniaxial tensile loading. (c) Triaxiality comparison, according to Bridgman’s approximation, between the current GTNB model prediction and the test.

Fig. 12. Comparisons between predictions of simulation with Landron, Eq. (15), and Bouaziz, Eq. (14), nucleation laws introduced in the GTNB model and the test measurement X-Ray tomography observations [23]. The number of voids per unit volume is calculated using the weighted average method Eq. (26).

investigate the respective accuracy of the two nucleation equations. After an axial strain of approximately 0.5, the number of voids drastically increases, since the coalescence phenomenon catalyzes the nucleation of new generations of small voids. When coalescence takes place and leads to micro-cracks interconnecting the larger voids, the assumption of spherical voids becomes more and more inadequate—as indicated by the sample’s micrograph at 0.83 strain in Fig. 12. Nonetheless, the experimental numerical

Fig. 13. Comparison between predictions of simulation with Landron, Eq. (15), and Bouaziz, Eq. (14), during the notched tensile simulation, evolution of mean void radius.

density of voids is measured up to such very large strain values and can be used for the validation of the models. The mean void radius is calculated with Eq. (16). Given that the second term in this equation is related to nucleation, the two predictions of the mean void radius are not identical, as shown in Fig. 13. The prediction using Eq. (15) gives slightly higher mean void radius values. Note that with respect to all of the other variables related to damage, the mean void radius remains almost constant during the deformation for both nucleation laws, in agreement with the experiments. When coalescence develops at

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

11

spherical voids; consequently, a simple phenomenological description is adopted to describe the phenomenon. In addition to the average values used for the confrontation to experiments, the numerical simulation illustrated the heterogeneity of most state variables, with the maximum values of, e.g., porosity and plastic strain, located near the center of the specimen. This heterogeneity, which increases during the loading, exhibited little mesh sensitivity prior to the development of coalescence. Consequently, the adopted approach appears as a good candidate to describe void nucleation and growth with physicallyinspired state variables and evolution laws. The accurate description of coalescence requires further improvements and remains under investigation, as well as its coupling with a fracture criterion, in the light of the recent observations of fracture strains and microstructure of void populations in DP steels. The experimental validation will be further extended to other sample geometries in flat sheet steels, as well as industrial applications. Fig. 14. Comparison between predictions of simulation with Landron, Eq. (15), and Bouaziz, Eq. (14), during the notched tensile simulation, evolution of porosity.

large strains, however, the experimental void radius increases. The calculated ones evolve with the same very small growth rate. This is consistent with the phenomenological description of the coalescence in the present model, and indicates that the physical meaning of the quantity R is lost during the coalescence step. This observation should be also corroborated with the fact that the approximation of spherical voids of identical size becomes more and more questionable as coalescence appears. Eventually, the combined effect of the numerical void density and of the mean void radius provides the prediction of the void volume fraction shown in Fig. 14. In the range of strains preceding coalescence, the void volume fraction evolves from the initial value of 2  10  4, up to values of 2  10  3 and more. The two predictions have similar accuracies although with different trends. Both seem able to correctly describe the experiments. The parameter identification has a non-negligible influence on the final predictions, and future investigations will be devoted to better understand its effect.

Acknowledgments The authors acknowledge ArcelorMittal Research S.A., the ANRT–Agence Nationale de la Recherche et de la Technologie (France), and the Interuniversity Attraction Poles Program-Belgian State–Belgian Science Policy P7 INTEMATE for their support. As Research Director of FRS-FNRS, AM Habraken thanks this fund for its financial support.

Appendix Hill matrix (H ) 2 G þ H H 6 H HþF 6 6 6 G F H ¼6 6 0 0 6 6 4 0 0

G

0

0

0

F F þG

0 0

0 0

0

2N

0

0

0

2L

7 7 7 7 7 0 7 7 7 0 5

0

0

0

0

F¼

7. Conclusions The GTNB model, extension of the classical Gurson–Tvergaard– Needleman model, has been completed with a new void nucleation law developed by Landron et al. [23] and was implemented in the commercial finite element code Abaqus. This physically based law incorporates the effect of the kinematic hardening on nucleation, along with an improved evolution of the mean radius of the voids, inferred from X-ray diffraction tomography observations in dual-phase steels. The model was extended to include the coalescence stage. The damage model was applied for the simulation of a tensile notch test performed on DP steel and was validated with respect to direct porosity measurements obtained by X-ray diffraction tomography. An adapted mesh was designed to ensure an accurate extraction of average values over the same volumes/ areas of observation as in the actual experiments. The predicted porosity evolution is well validated up to a strain of 0.5. Furthermore, each of the two ingredients of the porosity evolution – number of voids and their mean radius – are in good agreement with the experimental evolutions, thus confirming the importance of this physically inspired description. For larger strain levels, the apparition of the coalescence weakens the physical meaning of these quantities and of the hypothesis of

0

3

0 0

2M

2r 0 2 2r 0 ;G¼ ;H¼ ; L ¼ M ¼ N ¼ ðF þ GÞðr 45 þ 0:5Þ r 90 ð1þ r 0 Þ ð1 þ r 0 Þ ð1 þ r 0 Þ

Pseudo-invert hill matrix "

½B11 

0

0

½B22 

2

#

2=N

6 ½B22  ¼ 4 0 0 2 HþGþ4 F 1 6 H2 F2 G ½B11  ¼ 4 9 ðH G þ H F þ G FÞ G2 H2 F H 1 ¼

;

0 2=L 0

0

3

0 7 5

2=M

3

H2 F2 G

G2 H2 F

Hþ4 Gþ F

F2 H2 G 7 5

F2 H2 G

F þGþ4 H

K calculation (from Benzerga and Besson [9]) The coefficient k introduced into the cosh of the yield function. s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     h1 þ h2 þ h3 1 1 1 þ 0:8 k¼ 1:6 þ þ h4 h5 h6 h1 h2 þh2 h3 þ h1 h3     2 r 0 r 90 2r 0 2 3ðr 0 r 90 1Þ ; h2 ¼ h1 1 3 1þ r 0 r 0 r 90 2r 0 2      3r 0 ðr 90 1Þ 3ðr 90 þ 1Þ ; h4 ¼ h1 0:5 h3 ¼ h1 1 r 0 r 90 2r 90 2 r 0 r 90 2r 90 2       3r 0 ðr90 þ 1Þ ð2r 45 þ1Þðr 90 þ 1Þ ; h6 ¼ h1 0:5 h5 ¼ h1 0:5 r 0 r 90 2r 0 2 r0 r 90 2r 0 2 h1 ¼ 

The coefficient k is equal to 2 when the plasticity of the matrix is isotropic.

12

J. Fansi et al. / Materials Science & Engineering A 569 (2013) 1–12

References [1] L. Kachanov, Izv. Akad. Nauk SSR otd Tech. Nauk 8 (1958) 26–31. [2] J. Lemaitre, A course on damage mechanics, Springer-Verlag, Berlin, 1992. [3] J.P. Cordebois, Crite res d’instabilite´ plastiques et endommagement ductile en grandes de´formations. Ph.D. Thesis, University Pierre et Marie Curie, France, 1983. [4] A.L. Gurson, J. Eng. Mater. Technol. 99 (1977) 2–15. [5] V. Tvergaard, Int. J. Fracture 17 (1981) 389–407. [6] V. Tvergaard, Int. J. Fract 18 (1982) 237–252. [7] V. Tvergaard, A. Needleman, Acta Metall. 32 (1984) 157–169. [8] M. Gologanu, J.B. Leblond, J. Devaux, J. Mech. Phys. Solids 41 (11) (1993) 1723–1754. [9] A.A. Benzerga, J. Besson, Eur. J. Mech. A: Solids 20 (2001) 397–434. [10] J.B. Leblond, G. Perrin, J. Devaux, Eur. J. Mech. A: Solids 14 (1995) 499–527. [11] P. Ne gre, D. Steglich, W. Brocks, M. Koc- ak, Comput. Mater. Sci. 28 (3–4) (2003) 723–731. [12] M. Mear, J.W. Hutchinson, Mechanics of Materials 4 (1985) 395–407. [13] A.L. Helbert, X. Feaugas, M. Clavel, Acta Metall. 46 (1998) 939–951. [14] O. Bouaziz, E. Maire, M. Giton, J. Lamarre, Y. Salingue, M. Dimechiele, Rev. Me´t. 105 (2008) 102–107. [15] T. Pardoen, Comput. Struct. 84 (2006) 1641–1650. [16] P.F. Thomason, Ductile Fracture of Metals, Pergamon Press, Oxford, 1990. [17] M. Ben Bettaieb, X. Lemoine, L. Duchˆene, A.-M. Habraken, Int. J. Numer. Methods Eng. 85 (8) (2012) 1049–1072.

[18] M. Ben Bettaieb, X. Lemoine, O. Bouaziz, A.-M. Habraken, L. Duchˆene, Mech. Mater. 43 (2011) 139–156. [19] C. Landron, E. Maire, O. Bouaziz, J. Adrien, L. Lecarme, A. Bareggi, Acta Mater. 59 (20) (2011) 7564–7573. [20] J. Fansi, M. Ben Bettaieb, T. Balan, X. Lemoine, A.-M. Habraken, Key Eng. Mater. (2012) 77–80. [21] J. Fansi, A.-M. Habraken, T. Balan, X. Lemoine, C. Landron, E. Maire, O. Bouaziz, M. Ben Bettaieb, Prooceedings of the 11th Biennial Conference on Engineering Systems Design and Analysis, Transactions of the ASME, 2012. [22] S. Arndt, B. Svendsen, D. Klingbeil, Tech. Mech. 17 (1997) 323–332. [23] C. Landron, O. Bouaziz, E. Maire, J. Adrien, Scr. Mater. 63 (2010) 973–976. [24] J.R. Rice, D.M. Tracey, J. Mech. Phys. Solids 17 (1969) 201–217. [25] S. Allain, O. Bouaziz, Mater. Sci. Eng. A 496 (2008) 329–336. [26] E. Maire, O. Bouaziz, M. Dimechiele, C. Verdu, Acta Mater. 56 (2008) 4954–4964. [27] Y. Huang, Trans. ASME 58 (1991) 1084–1086. [28] L.M. Brown, J.D. Embury, Proceedings of the 3rd International Conference on Strength of Metals and Alloys, 1973 pp. 164–169. [29] Hibbitt-Karlsson & Sorensen, Abaqus User’s Manual, Version 6.11 (2011). [30] P.W. Bridgman, Rev. Mod. Phys. 17 (1945) 3–17. [31] C. Landron, E. Maire, J. Adrien, O. Bouaziz, M. Di Michiel, P. Cloetens, H. Suhonen, Nuclear Instruments & Methods in Physics Research Section B 284 (2012) 15–18. [32] Y. Ling, AMP J. Technol., Harrisburg (Pennsylvania) 5 (1996) 37–48.