Modified GTN model for a broad range of stress states and application to ductile fracture

Accepted Manuscript Modified GTN model for a broad range of stress states and application to ductile fracture Wei Jiang, Yazhi Li, Jie Su PII:

S0997-7538(15)00171-0

DOI:

10.1016/j.euromechsol.2015.12.009

Reference:

EJMSOL 3263

To appear in:

European Journal of Mechanics / A Solids

Received Date: 28 February 2015 Revised Date:

22 August 2015

Accepted Date: 16 December 2015

Please cite this article as: Jiang, W., Li, Y., Su, J., Modified GTN model for a broad range of stress states and application to ductile fracture, European Journal of Mechanics / A Solids (2016), doi: 10.1016/ j.euromechsol.2015.12.009. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Modified GTN model for a broad range of stress states and application to ductile fracture Wei Jianga,*, Yazhi Lia, Jie Sua a

School of Aeronautics, Northwestern Polytechnical University, Shaanxi Xi’an 710072, China

RI PT

*Corresponding author. Tel.: +86 18292166256; E-mail address: [email protected] (W. Jiang).

Abstract:

Based on two ductile fracture mechanisms, a new modification to the Gurson–Tvergaard–Needleman (GTN) model is proposed for the simulation of ductile fracture behaviors under high, low and negative stress triaxiality

SC

loadings. In the modified GTN model, two distinctive damage parameters, respectively related to void growth mechanism and void shear mechanism, are introduced into yield function as internal variables of degradation process.

M AN U

The void volume fraction, similar to the original GTN model, is still adopted to capture the volumetric damage due to the void nucleation, growth and subsequent inner-necking, which plays the main role in tension. The additional shear damage, due to void rotation, distortion as well as “void sheet” effects, can accumulate in different stress states, even under negative stress triaxiality with the auxiliary of a new stress–state dependent function. The development of numerical procedures into finite element platform ABAQUS/Explicit allows the analysis of damage evolution, crack

TE D

initiation and growth. The validity of the new model is examined by comparing numerical results with the experimental data of specimens including axisymmetric tensile bar, flat grooved plate, torsion tube and compression cylinder. The results show that the modified GTN model performs well in predicting the overall load–

EP

displacement response and fracture path of various specimens under stress states covering a wide range of

AC C

triaxiality and Lode parameter combinations.

Keywords: ductile fracture; void volume fraction; shear damage

1

Introduction

At the microscopic level, one of the most important mechanisms capturing ductile fracture in metallic materials is the nucleation, growth and coalescence of micro voids. In structural materials, voids nucleate at inclusions and second-phase particles, by decohesion of the particle–matrix interface or by particle cracking. Void growth is subsequently driven by diffuse plastic deformation of the surrounding matrix. Void coalescence is the final stage of

ACCEPTED MANUSCRIPT ductile failure, during which plastic deformation localized at the microscale inside the inter-void ligament between neighboring voids [1–5]. Void dilation and elongation are two important void deformation modes that have been observed in the void growth phase, and this phase is highly influenced by stress states [6, 7]. Experimental studies by Barsoum and

RI PT

Faleskog [8] have shown that void dilation is significant at high stress triaxiality, whereas void elongation is predominant at low stress triaxiality. The dilating voids coalesce by reduction of the inter-void ligament, which necks down and lead to rupture, in the case that the void size increases significantly before localization. This is known as void internal–necking rupture mechanism and favored in the stress states of high triaxiality [5–9]. The

SC

second mode of void coalescence is shear coalescence, which is similar to shear banding but at the scale of the voids. At low stress triaxiality, void elongation, rotation and interaction with neighboring voids leads to shear

M AN U

localization in the inter-void ligaments causing ductile fracture. In particularly, when a second, smaller population of voids nucleate between the primary voids and link-up in shear band, this shear coalescence mechanism is called “void sheet” [8–12].

Many theoretical models have been developed to predict ductile fracture using the mechanism-based local approach concept. The early micromechanical models for ductile damage considered the growth of isolated voids in

TE D

a rigid perfectly plastic matrix [13, 14]. These studies revealed the dramatic effect of hydrostatic stress on void growth. Failure occurs when the void growth ration has reached a critical value which is assumed to be a material parameter. Thereafter, constitutive relations for porous plastic solids were developed by coupling the damage

EP

parameter (porosity) with stress components in yield functions [15, 16]. The most widely used porous material model for analyzing the ductile fracture is that of Gurson (1977) [15], which is the representative of the void

AC C

growth stage only. Based on micromechanical cell model studies, heuristic supplements have been made to account for void nucleation [17], void coalescence [18], void shape effect [19, 20], void size effect [21, 22], and anisotropy [23– 25].

One of the most referenced void–nucleation–growth–coalescence models is the Gurson–Tvergaard–Needleman (GTN) model [15, 17, 18]. Although it is well known to give accurate results in medium to high stress triaxiality, GTN model still suffers an apparent limitation that it is powerless for shear dominated loadings [26–35]. In the extreme case of pure shear, where the stress triaxiality is zero, the GTN model predicts no damage accumulation and gives no failure at all if the nucleation of voids is neglected. Nahshon and Hutchinson (2008) [26] introduced a phenomenological damage law into GTN model to account for the effective damage accumulation due to the shear

ACCEPTED MANUSCRIPT mechanism. A similar modification to GTN model for shear failure was pursued by Xue (2008) [27]. These extensions retain the isotropy of the GTN model by making use of the third invariant of stress to introduce shear dependence. Both of the two shear modifications have been continually adopted to simulate shear failure in literatures [28–33], but unfortunately, their disadvantages cannot be ignored. (1) Shear damage may introduce a

RI PT

mendacious contribute to the volumetric part of plastic strain. (2) Failures in axisymmetric compression states cannot be dealt with, although they are controlled by shear rupture mechanisms. In the recent study of Malcher et al. (2014) [34], an evolution law for shear damage, incorporating the damage nucleation, growth and post-coalescence effects, was proposed and combined with GTN model to capture the shear localization of plastic strain inside the

SC

inter-void ligament. Zhou et al. (2014) [35] proposed a new modification for shear effect by combining the damage mechanics concept of Lemaitre [36] with the GTN model. Both of these two studies improve the performance of

M AN U

GTN based model in predicting shear dominated failures.

In general, the coupled mechanism–based damage models have the ability to predict the correct fracture behavior under a specific range of combined stress triaxiality and Lode angle and are fairly accurate for loading conditions close to calibrating points. Therefore, the grand challenge is the development of a computational capability for predicting localization, crack formation and crack propagation under all stress states. The aim of this paper is to

TE D

develop a new modification and the relative calibrated procedures to the GTN model and make it possible for modeling ductile failures under a broad range of stress states. The organization of this paper is as follows: Section 2 explains the modified GTN model in detail after presenting the original GTN model and the shear extension by

EP

Nahshon and Hutchinson. Section 3 presents the FE analyses for a variety of laboratory specimens by the developed numerical procedures to assess its ability of predicting the fracture processes under various stress states.

2

AC C

Section 4 lists several conclusions drawn from the FE analysis.

Constitutive modeling of porous materials

2.1 Characterization of stress states For an isotropic material, stress states are characterized by the symmetric stress tensor (six components) or its eigenvalues (three principal stresses: σ1 , σ 2 , σ 3 ), while most of material models are formulated in terms of three invariants of the stress tensor σ , which can be defined by I1 = −3 p = 3σ m = σ : Ι = σ1 + σ 2 + σ 3

(1)

ACCEPTED MANUSCRIPT 1 1 1 2 2 2 J 2 = q 2 = s : s =  (σ 1 − σ 2 ) + (σ 2 − σ 3 ) + ( σ 3 − σ 1 )    3 2 6

(2)

J3 = det ( s) = (σ1 −σm )(σ2 −σm ) (σ3 −σm ) = s1s2s3

(3)

where I is the second order identity tensor; p is hydrostatic stress, which has a positive value in compression; σ m is

deviatoric stress tensor; s is the deviatoric stress tensor s = σ + pI

RI PT

mean stress, which has a positive value in tension; q is the equivalent stress; s1 , s2 , s3 are principal values of

(4)

Stress triaxiality ratio T, defined as the ratio of the mean stress to the equivalent stress, is widely used as the sole

σm q

=−

p q

(5)

M AN U

T =

SC

parameter to characterize the effect of triaxial stress states in ductile fracture.

However, experimental observations and numerical analyses indicate that stress triaxiality cannot fully describe ductile fracture envelope because multiple stress states with different principal stress values can result in the same stress triaxialty ratio [37–40]. For instance, experimental data by Wierzbicki et.al (2005) [41] showed that different critical strains in ductile fracture of 2024-T351 aluminum alloy under axisymmetric and plane stress states even

TE D

though the stress triaxiality ratio remains the same. Therefore, another parameter, e.g., the Lode parameter, which plays the role of a deviatoric state parameter,–should be introduced to distinguish different stress states having the same stress triaxiality ratio. In literatures, Lode parameter–has different representations [30, 42, 43]:–

27 J 3 3

, µ=

2σ 2 − σ 1 − σ 3 σ1 − σ 3

AC C

ξ = cos ( 3θ ) =

EP

 2  s1 − s3 1    1  s −s  π −   , θ L = θ − = tan −1   2 2 3 − 1   6    3  s2 − s3 2    3  s1 − s3

θ = cot −1 

2q

Where the range of θ is [0,

π

3

]; the range of θ L is [ −

(6)

(7) π π , ]; both ξ and µ are limited to [-1, 1]. 6

6

2.2 Gurson-Tvergaard-Needleman model The yield function of GTN model [15, 17, 18, 44, 45], which assumes isotropic hardening, is expressed by Φ ( σ, f ) =

q2

σm

2

 3q p  2 *2 + 2 q1 f * cosh  − 2  − 1 − q1 f = 0  2 σm 

(8)

ACCEPTED MANUSCRIPT where q is macroscopic equivalent stress, p is macroscopic hydrostatic pressure. Macroscopic stress σ is defined as the force per “current unit area,” comprised of voids and the solid (matrix) material. The current state is characterized by f, the damage parameter which can be interpreted as an effective void volume fraction, and σ m which is the current effective true stress governing flow of the matrix material. Tvergaard (1981, 1982) [44, 45]

better agreement with experimental data.

RI PT

introduced the constants q1 , q2 to account for void interaction effects due to multiple-void arrays and to give a

The f * function is adopted by Tvergaard and Needleman (1984) [18] to model the rapid loss of stress carrying

f ≤ fc fc < f < fF f ≥ fF

(9)

M AN U

f  f * − fc  f * =  fc + u ( f − fc ) f F − fc  f*  u

SC

capacity that accompanies void coalescence. The function is defined as

where f c is the critical void volume fraction at which voids begin coalesce; fu* = 1 q1 is the corresponding value of

f * at zero stress; f F denotes the void volume fraction at complete failure.

According to the hypothesis of generalized normality, the plastic flow rule is given by  3q p   f * q1 q2 ∂Φ &  3 =λ  2 s + sinh  − 2  I  σ ∂σ σm  2 σm    m

TE D

ε& p = λ&

(10)

where ε& p is the macroscopic plastic strain rate tensor; λ& is the plastic multiplier.

nucleation of new voids.

EP

The change in volume fraction of the voids is due partly to the growth of existing voids and partly to the

AC C

f& = f&growth + f&nucleation

(11)

Assuming that the matrix material is plastically incompressible, the void growth rate is a function of the plastic volume change and can be written as f&growth = (1 − f ) ε& p : I

(12)

Nucleation of voids can occur as a result of micro-cracking and/or decohesion of the particle-matrix interface. The plastic strain controlled nucleation is given by

f&nucleation = Aε&mp

(13)

ACCEPTED MANUSCRIPT The void nucleation intensity A is a function of ε mp , the equivalent plastic strain in the matrix material. It is assumed that void nucleation is taken to be a continuous process as suggested by Chu and Needleman (1980) [17], so that the coefficient A take the form

sN

  

2

  

(14)

RI PT

A=

 1  ε p −ε N exp  −  m  2  sN 2π 

fN

where f N is the total void volume fraction that can be nucleated; ε N is the mean equivalent plastic strain for void nucleation; and sN is the standard deviation of the distribution.

( )

SC

p The hardening of the matrix material is described through σ m = σ m ε m , ε mp is the equivalent plastic strain of

plastic work and the matrix plastic dissipation:

(1− f ) σmε&mp = σ : ε& p

M AN U

matrix. The evolution of ε mp is assumed to be governed by enforcing equality between the rates of macroscopic

(15)

2.3 Shear extension by Nahshon and Hutchinson (N-H)

TE D

Although the original GTN model can give accurate results in medium to high stress triaxialities, it is still unable to capture localization and fracture for low triaxialities in shear–dominated deformations. Nahshon and Hutchinson (2008) [26] introduced an additional term to quantify the void shearing effects at low triaxiality, which was directly

EP

added to the total void volume fraction in spite of not representing a physical accumulation of the porosity. The extended shear damage term that assumes linear dependence on the current porosity and the deviatoric part of plastic

AC C

strain rate can be expressed as s : ε& f&shear = kω f ω q

p

(16)

In the original work of Nahshon and Hutchinson (2008), Eq. (16) was proposed and calibrated in a pure shear condition that void nucleation was absence. Accordingly, the void volume fraction rate, f& , is given by the sum of the growth of voids as well as shear effects by s : ε& p f& = (1 − f ) ε& p : I + kω f ω q

(17)

where kω is a material parameter which is used to scale the magnitude of damage growth rate in pure shear loading;

ACCEPTED MANUSCRIPT ω , the function of the stress states, is used as a weight function, which is defined by

 27J  ω = 1 − ξ = 1 − cos (3θ ) = 1 −  33   2q  2

2

2

(18)

where, ξ is Lode parameter defined in Eq.(7). The value of ω changes from zero at generalized tension/compression

RI PT

(uniaxial tension/compression + hydrostatic states) to 1 at generalized shear (pure shear + hydrostatic state). Generally, the growth and internal necking of voids governs the rupture mechanism in tension dominated loading modes, while the internal shearing due to the effects of void elongation, rotation and “void-sheet” dominates for shear conditions [5–12]. Since the two failure mechanisms are predominant in different loading conditions, they should

SC

naturally be described by two separate damage parameters. This modification for incorporating shear effects by directly adding shear damage to the total void volume fraction cannot satisfy the assumption of the original GTN

M AN U

model that matrix material is plastically incompressible and the plastic volume change of the material is totally coming from void size changing. During this case, shear damage may introduce a mendacious contribute to the volumetric part of plastic strain. Furthermore, Nahshon and Hutchinson adopted a Lode dependent function, ω , to discriminate between axisymmetric and shear–dominated stress states. For stress states with negative triaxialities, the shortest distance of internal ligament can be reduced and void can be distortion and linking under compressive

TE D

loading without the increases of porosity. Unfortunately, the weight function, ω , has a

quite low value in

axisymmetric stress states with negative triaxiality, suggesting that the extended shear term work with limited effects. In addition, it can be seen from Eq. (17) that the second term may accumulate certain damage under

EP

negative triaxiality while the first term will give a reverse effect due to the negative plastic volumetric strain. Thus,

AC C

the generalized compression failure cannot be captured by the N-H version extended GTN model.

2.4 Modified constitutive relationship To overcome the above disadvantages, a modified GTN model is proposed, which includes two separate damage parameters. In this new model, an independent damage variable due to shear effects is introduced into yield function by incorporating the damage concept of Lemaitre. As a consequence, shear damage only affects the deviatoric deformation while the only cause of the plastic volume change is the porosity. A new function of triaxiality and Lode angle is adopted to generalize the shear damage accumulation obtained in a pure shear state to an arbitrary stress state.

ACCEPTED MANUSCRIPT 2.4.1 Yield function incorporating shear damage To characterize the losing of load carrying capacity due to shear effects, a new damage variable, Dshear , need to be introduced into yield function and couple with stress components. Following the principle of strain equivalence (Lemaitre, 1985, 2000) [36, 46], the only change in comparison with the original damage model is to replace the

RI PT

stress by the effective stress in the constitutive equations. The effective stress and its equivalent and hydrostatic components are defined as: σ% =

σ q p , q% = , p% = 1 − Dshear 1 − Dshear 1 − Dshear

q2

σ m (1 − Dshear ) 2

2

  3q  p 2  +  2q1 f cosh  − 2  − 1 − ( q1 f )  = 0    2 σ m (1 − Dshear ) 

(20)

M AN U

Φ=

SC

Therefore, the yield function is modified as

(19)

The new yield surface includes two separate damage parameters. The volumetric damage component is characterized by the void volume fraction, f, and the deviatoric damage component is described by shear damage variable, Dshear .

The evolution law for the plastic flow, assuming the hypothesis of generalized normality and associated flow rule,

TE D

is given by

 3q 2   q1q2 f ∂Φ &  3 p =λ  s + sinh − ε& p = λ&   I  2  2 σ (1 − D  σ 2 (1 − D  ∂σ σ m (1 − Dshear ) m shear )    shear )  m

(21)

EP

The evolution of matrix equivalent plastic strain, ε mp , should take the shear damage into account and can be represented as:

σ : ε& p (1 − f ) (1 − Dshear )

AC C

σ m ε&mp =

(2 2)

It is important to mention that the modified yield surface, described by Eq. (20), does not change the original GTN model under generalized tension where shear effects are not present. In a pure shear condition where q=0, if void nucleation is not invoked, Eq. (20) can be simplified as Φ=

q2

(1 − q1 f 0 ) (1 − Dshear ) 2

2

− σ m2 = 0

(23)

In the situation of zero initial porosity, Eq. (23) represent the Lemaitre model that damage is coupled with Mises

ACCEPTED MANUSCRIPT plasticity. In addition, it should also be noted that f = Dshear = 0 implies that the material is damage free, and the modified yield condition reduces to that of von Mises; f = 1 q1 or Dshear = 1 implies that the material is fully voided or at shear failure, and the material totally loses its load carrying capacity. Fig. 1 illustrates the yield surfaces for

M AN U

SC

RI PT

different levels of void volume fraction and shear damage in the p–q plane.

(b)

EP

TE D

(a)

(c)

AC C

Fig.1 Yield surfaces of the proposed model. (a) Yield surfaces of q1f=0.01,0.1,0.3,0.6,0.9 when Dshear is fixed as 0; (b) Yield surfaces of Dshear=0,0.2,0.5,0.7,0.95 when q1f is fixed as 0.01; (c) Yield surfaces of q1f=0.01,0.1,0.3 when Dshear changes as 0, 0.2, 0.5 and 0.7.

2.4.2 Damage evolution

In the current model, the damage evolution involves two aspects to capture the material failure: the accumulation of the void volume fraction and damage due to shear effects.

(1) Void volume fraction The accumulation of porosity, f, which is the same as the original GTN model, is composed of two parts, the

ACCEPTED MANUSCRIPT growth of existing voids and the nucleation of new voids. The void growth rate is still assumed as a result of plastic volume change. Since a phenomenological model will be used to describe the experimentally observed material degradation and softening due to shear effects, it is reasonable to exclude the porosity induced damage accumulation under a pure shear condition. Hence, a new void nucleation criteria is proposed here, which can be expressed as: f&nucleation = (1 − ωσ ) Aε&mp

RI PT

(24)

where, ωσ is a new stress–state dependent function and ωσ has a unit value in pure shear state. Void nucleation intensity A is still a function of ε mp and can be expressed as:   

2

  

p≤0

SC

  1  ε p −ε fN N exp  −  m  A =  sN 2π  2  sN   0 

p>0

(25)

M AN U

It should be mentioned that Eq. (25) implies voids nucleates only when the loading is tensile.

(2) Shear damage

The second damage evolution law is that induced by shear effects, which is defined by an independent scalar

s : ε& p D& shear = k ω Dshear ω σ q

TE D

variable, as

(2 6)

This modification follows the notion by Nahshon and Hutchinson (2008) [26] that the volume of voids undergoing

EP

shear may not increase, but the void distortion, rotation as well as smaller voids nucleating at internal shear bands contribute to material softening and constitute an effective increase in damage. The shear damage rate D& shear is

AC C

assumed to be linearly dependent on the current value of Dshear and the deviatoric part of plastic strain rate. In pure shear conditions, there is a relationship between stress tensor and its deviator: σ1 = −σ 3 ( σ1 > 0 ), σ 2 =0 , 2 σij =sij and q= 3σ1 ; there is ε&1p = −ε&3p ( ε&1p > 0 ), ε&2p = 0 and ε& p = ε&1p for strain rate components. Additionally, 3

the value of ωσ should be taken as unit. From the analysis above, the shear damage rate function can be reorganized as s : ε& p D& shear = k ω Dshear = k ω Dshear ε& p q

(27)

ACCEPTED MANUSCRIPT With the initial shear damage D0 , Dshear can be integrated and expressed explicitly by

(

Dshear = D0 exp k ωε p

)

(28)

For simplification, we consider a pure shear state in a full dense material (f0 = 0) that elasticity can be neglected n ( E → 0 , ε = ε p ) and it obeys a power-law hardening rule of σ m = σ 0 + K ( ε ) . σ 0 is initial yield strength; K and

RI PT

n are strain hardening parameters. The evolution of equivalent stress as a function of damage/strain, obtained from shear damage accumulation law (Eqs. (27) (28)) in combine with yield condition (Eq. (24)), could be explicitly expressed as:



n  1 Dshear   ln    (1 − Dshear ) D0    kω 

SC

 K = 1 + σ0  σ0  q



)n  (1 − D0 exp ( kω ε ) )

(29)

(30)

M AN U

 K = 1 + (ε σ0  σ0 q

A set of curves illustrating the relationship between normalized equivalent stress, shear damage and strain are

AC C

EP

TE D

plotted in Fig. 2.

Fig.2 The evolution of (a) equivalent stress normalized by initial yield strength; (b) shear damage as a function of strain, with σ 0 = 350MPa , K = 500MPa , n = 0.15 .

Since two separate scalar damage variables are considered, the present model has two distinct critical values: the critical void volume fraction, fc, and the critical shear damage value, Dc . Both two critical values should be determined experimentally: fc is obtained from a specimen subjected to tension dominant loading and Dc is obtained from a specimen subjected to a shear loading condition. In addition, similar to f * function, a new Dshear* function is adopted to simulate the post process of coalescence in micro-shear bands when the critical condition is reached. In

ACCEPTED MANUSCRIPT this stage, damage increases rapidly and the softening process accelerates so that material loses its load carrying capacity abruptly.

(31)

RI PT

Dshear*

 Dshear Dshear ≤ Dc  * D − Dc  =  Dc + u ( Dshear − Dc ) Dc < f < DF DF − Dc   D* Dshear ≥ DF  u

* where Du* = 1 is the corresponding value of Dshear at zero stress; DF denotes the shear damage at complete shear

failure; Dc is the critical shear damage which can be approximated by the following equation

)

2.4.3 Stress–state dependence function

M AN U

where ε f is the equivalent strain at the beginning of shear failure.

(32)

SC

(

Dc =& D0 exp kω ε f

The shear damage evolution, which is calibrated in a shear loading condition, need to be generalized to an arbitrary stress state. This is usually accomplished by introducing a Lode angle dependent function in literatures [26, 27]. The Lode angle dependent function proposed by Nahshon and Hutchinson (2008) [26] discriminates between

TE D

axisymmetric and shear–dominated stress states and expresses a quadratic relation with the normalized third invariant (Eq. (18)). An alternative Lode angle dependence function proposed by Xue (2008) [27] is defined by a linear

g (θ L ) = 1 −

6 θL

π

EP

expression of the normalized Lode angle, as:

(3 3)

AC C

where θ L is Lode angle defined in Eq.(6). As discussed in Section 2.3, these two types of Lode angle dependent functions give a zero value of g or ω in axisymmetric compression, preventing shear modification term from taking effect.

In this paper, a new stress–state dependent function is introduced as a weight function controlling shear effects in multiple stress states. The new function, ωσ , still has a unit value in pure shear and shear plus hydrostatic stress conditions. For the positive triaxiality, it is taken as the same form as the Lode angle dependent function of Nahshon and Hutchinson (2008) [26], therefore, ωσ still has a zero value in generalized tension. For a negative triaxiality, ωσ is dependent on the current triaxiality and Lode angle, giving a non-zero value in generalized

ACCEPTED MANUSCRIPT compression. The proposed stress state dependent function, ωσ , can be expressed as:

( (

 1−ξ 2   ωσ =  1 − ξ 2  1 − ξ 2 

) (1 − c ) + c ) 1 + ck T  − ck T 1

1

1

1

T < −k −k ≤ T ≤ 0

(34)

T >0

RI PT

where, c1 and k are constants which need to be calibrated by a group of experimental data with negative triaxiality.

M AN U

SC

Three types of stress–state dependent functions are illustrated in triaxial stress space by Fig.3.

(b)

AC C

EP

TE D

(a)

(c) Fig. 3 Stress–state dependent functions in triaxial stress space: (a) Lode angle dependent function proposed by Xue (2008) [27]; (b) Lode angle dependence function by Nahshon and Hutchinson (2008) [26]; (c) The proposed stress–state dependent function with c1=0.5 and k=0.5

3

Modeling of ductile fracture in 2024–T3 aluminum alloy The modified GTN damage model described in Section 2.4 is implemented in ABAQUS/Explicit via a user

ACCEPTED MANUSCRIPT defined material subroutine VUMAT. The numerical solution strategy adopted to perform simulations is following the work of Aravas (1987) [47] and Zhang (1995a, b) [48, 49] and illustrated in Appendix. 3.1 Material and Specimens A series of tests under various stress states were performed to verify the predictions of the overall elastic–plastic

RI PT

behaviors and the crack propagation paths by the present approach. Fig.4 illustrates the geometries and loading types of the specimens, including the smooth round bar (SRB), notched round bar (NRB), flat grooved plate (FGP) for tension, cylinder specimen for compression, and a thin–wall tube for torsion. SRB, NRB and FGP tensile tests provide medium to high triaxialities which change with notched (or grooved) radii (for SRB specimen, notched

SC

radius R=∞), while compressive tests give negative triaxialities. Lode angel θ in SRB, NRB, FGP and compression cylinder are 0, 0, π/6 and π/3, respectively. Pure shear condition was attained by testing a thin–wall

M AN U

tube under pure torsion load. The dimensions of the specimens are listed in Table 1. The materials of all the specimens are 2024 T3 aluminum alloy that is widely used in the aerospace industry. All the specimens are cut

from the same plate with a thickness 30mm. The tensile mechanical properties of 2024 T3 aluminum alloy are

AC C

EP

TE D

Young’s modulus, E = 71 GPa; Poisson ratio, ν = 0.33; initial yield strength, σ0 = 34.61MPa.

Fig.4 Geometries and loading types of the specimens

Table 1 Dimensions of the specimens (units: mm) Specimens

SRB NRB FGP Pure Torsion Tube Compression Cylinder

Gauge length 60 50 40 36 /

Diameter 9 9 / 12.5 15

Gauge section Width / / 40 / /

Thickness / / 2 1.2 /

Height

Notch radius

/ / / / 30,18.75,15,10

/ 12,6,3 12,4,2.4,1.6 / /

ACCEPTED MANUSCRIPT 3.2 Finite element models Four–node axisymmetric element with reduced integration (CAX4R) and eight–node solid element with reduced integration (C3D8R) were used in the finite element analyses. Due to the symmetry, half models were used to simulate FGP tension tests. Both 1/8 and full model were used for cylinder compression tests. Simulations of SRB,

RI PT

NRB tension tests as well as pure torsion tests were performed by full models. To properly capture the fracture surface, fine finite element mesh must be carefully designed to generate elements with nearly unit aspect ratio at failure in the rupture-critical locations. To improve computational efficiency, only the materials near the crack path were modeled using the modified GTN model for all the specimens except compression cylinders. Outside these

SC

regions, specimens were modeled with usual Mises plasticity. For compression cylinders, the compression platen was modeled by an analytical rigid body, and frictional surface contact was defined to model the interaction

M AN U

between the platen and the specimen. Since large slant fracture surfaces formed during failure in the axisymmetric compression tests, the modified constitutive relationship was applied to whole model to capture the correct shear

(b)

(c)

AC C

(a)

EP

TE D

fracture surface. Fig.5 illustrates the fine meshes adopted in the simulations.

(d) (e) Fig.5 Finite element meshes of the specimens: (a) SRB; (b) NRB–R6; (c) compression cylinder (H=15); (d) FGP–R12; (e) tension–torsion tube

ACCEPTED MANUSCRIPT 3.3 Parameter determinations The uniaxial true stress–strain relationship is approximated using the Ramberg–Osgood form fitted to the SRB tensile test data, which is of the form

where ε 0 =

σ0 E

n

(35)

RI PT

σ  ε +1 = α   ε0  σ0 

, α and n are fitted parameters.

In order to apply the present modified damage model, five groups of parameters should be first determined,

SC

including two constants of yield function (q1, q2); two hardening parameters of material (α, n); six parameters related to void volume fraction (f0, fc, fF, fN, εN and sN); four parameters related to shear damage (D0, Dc, DF and kω);

M AN U

and two parameter in the stress–state dependent function (c1, k). At first, the classic values ( q1 = 1.5 and q2 = 1.0 ) suggested by Tvergaard (1982) [45] were applied in current research. Then, two parameters, f0 and fN, can be obtained from the measurement of microstructure in the untested material. Chu and Needleman (1980) [17] suggested that fN is determined so that the total void volume nucleated is consistent with the volume fraction of second-phase particles. For the current material clusters of coarse Fe- and Si-rich intermetallic particles are

TE D

responsible for nucleation and the volume fraction of coarse particles was quantified as 2.0×10-2. In addition, the volume fraction of pores in the initial material was quantified as 0.4×10-2. So fN =0.02 and f0 =0.004 are used in the following simulations. Since void nucleation occurs at a very low level of plastic strain, εN and sN can be taken as

EP

0.2 and 0.1. At last, several tests are required to calibrate other parameters in the modified GTN model. A four–step procedure was followed: (1) data from a SRB tensile test was used to characterize the tensile stress–strain curve of

AC C

the undamaged material; (2) data from a NRB tensile test (R3) with a high stress triaxiality was used to identify the critical void volume fraction, fc, and void volume fraction at failure, fF; (3) a pure torsion test was used to calibrate shear damage parameters; and (4) the parameter c1 and k were calibrated using the data from a compression test (H=15mm) with negative triaxiality. All parameters characterizing the constitutive model of 2024–T3 aluminum alloy are summarized in table 2. Table 2 Calibrated parameters. Notations α n q1 q2

Parameters Strain hardening coefficient Strain hardening exponent Yield function constant Yield function constant

Value 0.5 6.1 1.5 1.0

Notations sN fN D0 Dc

Parameters Standard deviation on void nucleation Void volume fraction to be nucleated Initial shear damage Critical shear damage

Value 0.1 0.02 0.01 0.025

ACCEPTED MANUSCRIPT f0 fc fF εN

Initial void fraction Critical void volume fraction Void volume fraction at failure Mean nucleation strain

0.004 0.025 0.15 0.2

DF kω c1 k

Final shear damage at failure Shear coefficient Stress state dependence parameter Stress state dependence parameter

0.1 4 0.5 0.2

3.4 Results and discussions 3.4.1 Elastic-plastic response

RI PT

Fig. 6 compares the experimental load–displacement response of axisymmetric tensile specimens with the simulated results by the present approach. It can be seen that numerical results accord well with those of experiments. The void coalescence process is accompanied by an abrupt change in the slope of the local load–

SC

displacement curve. The effect of stress triaxiality on fracture strain under tensile loads is also obtained, just as what has widely been documented: higher triaxiality yields lower fracture strain.

M AN U

The simulated load–displacement curves of FGP specimens and the relevant experimental results are plotted in Fig.7. Both the stress traixiality and shear damage weight function (stress–state dependent function) have high values in the center of the FGP specimens, hence, void internal–necking rupture mechanism as well as void shearing rupture mechanism, simultaneously drive ductile failure initiation in these areas. The predicted load-displacement curves agree pretty well with the experiments.

TE D

Fig. 8 is the simulated torque–twist angle curve of pure torsion tube and the curve from experiment. Both are in accordance with each other. The degradation and failure process under pure shear load is well captured by the modified GTN model. It can be concluded that the modification of shear evolution, while phenomenological, is

EP

nevertheless formulated in accordance with the mechanism of softening in shear. The cylinder compression tests which can not be simmulated by the extended GTN model of Nahshon and

AC C

Hutchinson are also considered in this paper. Fig.9 compares the simulated load–displacement curves of compression cylinders with 4 different heights to the experimental results. Fig.9 shows that ductile failure under negative triaxiality loadings is treated well by the modified GTN model with the proposed stress–state dependent function.

As a summary, failure processes under torsion and compression loads that cannot be simulated by the original GTN model are well captured by the proposed model. With appropriate damage parameters, the modified damage model with an extra shear damage variable is able to give reasonable predictions for ductile metallic specimens covering a broad range of stress states.

RI PT

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

(a) (b) Fig.6 Elastic-plastic responses of axisymmetric tensile specimens: (a) SRB specimen; (b) NRB specimens.

AC C

EP

Fig.7 Elastic-plastic responses of FGP specimens.

Fig.8 Elastic-plastic responses of pure torsion tube

RI PT

ACCEPTED MANUSCRIPT

SC

Fig.9 Elastic-plastic responses of compression cylinders

3.4.2 Crack propagation path

M AN U

Fig. 10 (a)–(d) are the contour plots of stress triaxiality before crack initiation in axisymmetric tensile specimens. The highest triaxialities of SRB, NRB R12, R6 and R3, in sequence, are 0.74, 1.04, 1.29 and 1.37 at the center of the symmetric section. For SRB, NRB R12 and R6, the stress triaxialries at the free outer surface of symmetric section are very low, which are 0.21, 0.28 and 0.35, respectively. The triaxialities of NRB-R3, from the center to

TE D

the free outer surface along the symmetric section, are all at high levels ranging from 0.67 to 1.37. Void nucleation, growth and subsequent inner–necking result in failure initiation at the center of the symmetrical section subjected to the highest stress triaxiality.

Fig. 11 (a)–(d) illustrate the simulated crack growth paths of axisymmetric tensile specimens, from the center to

EP

the free outer surface through the minimum section. For specimens with large notch radii (smooth, R12, R6), the initial crack spreads in tensile mode that crack path is perpendicularly to the loading direction and the final

AC C

propagation is away from the mid-section on a conical surface. The transition from flat to slant fracture gives rise to a cup and cone appearance. For the specimen with smallest notch radius (R3), crack passes the minimum section in a zig-zag path but not form apparent shear lips. The simulated fracture surfaces are in accordance with the observed fracture morphology of NRB–R12 and R3 specimens shown in Fig. 12. In a word, cracks initiate at the center of the minimum section, where the stress triaxiality is the highest, and during the crack propagation, the transition of fracture mode from flat to slant will occur when stress triaxiality becomes quite low.

SC

RI PT

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

(a) (b) (c) (d) Fig. 10 Triaxiality distribution before fracture initiation in axisymmetric tensile specimens: (a) SRB; (b) NRB–R12; (c) NRB–R6; (d) NRB–R3.

AC C

(a) (b) (c) (d) Fig. 11 Simulated fracture surfaces in axisymmetric tensile specimens: (a) SRB; (b) NRB–R12; (c) NRB–R6; (d) NRB– R3;

–

–

(a) (b) Fig. 12 Experimental fracture surfaces in axisymmetric tensile specimens: (a) NRB–R12; (b) NRB–R3

Fig. 13 depicts the simulated results of FGP specimens, including final fracture surfaces by half models and

ACCEPTED MANUSCRIPT damage distributions in the mid-section before fracture occurs. Fracture initiates at the center of mid-section, where both two damages, shear damage as well as porousity, localizes in two diagonal shear bands. With damage accumulation, one of the shear bands dominates and leads to crack iniation. When comparing the maximum value of two damage in mid-section before fracture initiation, it is found that shear damages, in

RI PT

sequence, are 4.5, 3.4, 2 and 2 times as porousity for R12, R4, R2.4 and R1.6 specimens. Therefore, slant fracture prevails in all of FGP specimens and the crack is found to form at an approximately 45° angle to the

loading direction in R12 specimen and a smaller angle, 42°, 41° and 37°, in R4, R2.4 and R1.6 specimen. The fracture surface translates from slant to flat as the crack front approaches the end of the section due to the loss of

SC

plane strain stress state, leaving a small failure area of tensile mode. Details of the experimental failure surfaces of

TE D

M AN U

NGP specimens are shown in Fig. 14, which is in accordance with simulation results.

AC C

EP

(a)

(b)

M AN U

SC

(c)

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

(d) Fig. 13 Fracture surfaces and damage accumulation before fracture in flat grooved plates: (a) FGP–R12; (b) FGP–R4; (c) FGP–R2.4; (d) FGP–R1.6

Fig. 14 Experimental fracture surfaces in flat grooved plates.

In the thin–wall section of torsion tube, both triaxiality (T) and Lode parameter ( ξ ) remain nearly zero as plasticity evolves up to failure, suggestting that the thin–wall section undergoes a pure shear deformation during load proceeding. It is believed that shear fracture is driven by shear localization of plastic strain in void ligaments due to void rotation, distortion and the nucleation of second population voids inside inner shear band. Fracture initiates at the intersection of cylindrical surface and conic surface of the tube, where the euquivalent plastic strain

ACCEPTED MANUSCRIPT is the largest. Fig. 15 depicts the distribution of triaxiality (T), Lode parameter ( ξ ), equivalent plastic strain and shear damage before fracture initiation. Fig. 16 shows the calculated fracture surface in pure torsion specimen. Crack propagates from outside surface to inside of the tube wall and produces a flat fracture surface. Fig. 17 shows

RI PT

the experimental fracture surface, which is in accordance with simulation results.

(b)

M AN U

SC

(a)

EP

TE D

(c) (d) Fig.15 Contour plots before fracture initiation in pure torsion specimen: (a) Triaxialty; (b) Lode parameter; (c) Equivalent plastic strain; (d) shear damage.

AC C

Fig. 16 Simulated fracture surface of pure torsion tube.

Fig. 17 Experimental fracture surfaces of pure torsion tube

For compression cylinders, fracture initiates at the top or bottom of the free outer surface where both the stress– state dependent function and equivalent plastic strain have the highest values. It can be seen that the deformation

ACCEPTED MANUSCRIPT and damage are uniform in hoop direction at the beginning of the compression, but the deformation localizes at shear bands. Crack then propagates slantwise throughout the specimen, leading to a large shearing plane. Fig.18 depicts the simulation results of a compression cylinder with the height of 15 mm. Fig.18 (a) and (b) illustrate the contour plots of stress–state dependent function and equivalent plastic strain while (c) and (d) show the calculated

RI PT

fracture surfaces by 1/8 and full models. Fig. 19 plots the simulated fracture surfaces in compression cylinders with heights of 30mm, 18.75mm and 10mm. The experimental results of the compression cylinder with the height of 15

M AN U

SC

mm are given in Fig. 19. It can be seen that the simmulated fracture surface accord well with those of experiments.

(b)

EP

TE D

(a)

(c)

(d)

AC C

Fig. 18 Simulated results of compression cylinder with the height of 15 mm: (a) stress–state dependent function before fracture initiation; (b) equivalent plastic strain before fracture initiation; (c) fracture surface by a 1/8 model; (d) fracture surface by a full model.

(a)

(b)

(c)

RI PT

ACCEPTED MANUSCRIPT

M AN U

SC

Fig. 19 Simulated fracture surfaces in compression cylinders (a) H=10mm; (b) H=18.75mm; (c) H=30mm.

Fig. 20 Experimental fracture surfaces of compression cylinders with the height of 15 mm.

4

Conclusions

TE D

An extra damage variable related to void shear effects is introduced to Gurson–Tvergaard–Needleman (GTN) model, in combining with void volume fraction related to void growth and internal–necking mechanism, to capture the material degradation and failure process under various stress states. The implementation of the modified GTN

EP

model into ABAQUS/Explicit via user defined material subroutine (VUMAT) allows the analysis of damage evolution, crack initiation and growth as well as the overall elastic–plastic response in ductile metallic materials.

AC C

The accuracy of the new model is validated by comparing the results a series of numerical tests with experimental data. Several conclusions are summarized as follows: (1) The modification to GTN yield function in combine with two independent damage evolutions, while phenomenological, is nevertheless formulated in accordance with the mechanisms of softening in tension as well as in shear. (2) The additional shear damage term with the auxiliary of the new stress–state dependent function makes it possible for the computational approaches based on GTN model to be extended to failures under a broad range of stress states. (3) With carefully determined damage parameters, the modified GTN model has yielded reasonable predictions for

ACCEPTED MANUSCRIPT specimens with various geometries and loading types.

Appendix A

Numerical integration algorithms

The modified GTN damage model described in Section 2.4 is implemented in ABAQUS/Explicit via a user

RI PT

defined material subroutine VUMAT. In this section, the numerical solution strategy adopted to perform simulations is illustrated in detail. The procedure of integrating the rate-form constitutive equations includes two principal steps: an elastic prediction followed by a plastic correction where stresses and internal state variables are integrated using an implicit Euler method in combination with the return mapping algorithm [47–49]. Since a full

A.1 Elastic-plastic constitutive relationship

M AN U

SC

implicit integration method is used, all the equations should be constructed at the end of each time increment.

The constitutive relationship of the present damage model can be characterized by following equations: Strain decomposition: ε = εe + εp

Linear elasticity relationship: σ = Ce : ε e

(

)

Flow rule: ε& p = λ&

∂Φ ∂σ

TE D

1 2 3 Yield function: Φ p, q, H , H , H = 0

(A-2) (A-3)

(A-4)

(

α α p β Evolution laws of internal state variables: H& = h ε& , σ, H

EP

(A-1)

)

(A-5)

AC C

2   e where ε is the total strain; εe and εp are elastic and plastic strain components; Cijkl = 2Gδ ik δ jl +  K − G  δ ij δ kl is 3  

the fourth order isotropic elastic matrix, δ is Kronecker delta, G and K are the shear and bulk moduli; ε& p is plastic strain rate; λ& is the plastic multiplier; H1 = εmp , H 2 = f * and H 3 = Dshear* are three internal variables; α , β = 1, 2,3 .

A.2 Hydrostatic and deviatoric components The plastic strain rate is obtained by an associated flow rule and it can be decoupled into two scalar variables as ∂Φ 1 p = ε&v I + ε&dp n ε& p = λ& ∂σ 3

with

(A -6)

ACCEPTED MANUSCRIPT ε&vp = −λ&

∂Φ ∂Φ p , ε&d = λ& ∂p ∂q

(A -7)

where ε&vp and ε&dp are the volumetric and deviatoric parts of plastic strain rate induced by hydrostatic stress and equivalent stress, respectively; n=

3s is the vector in the deviatoric space normal to the yield surface; I is the second 2q

RI PT

order identity tensor. With the notation, n , stress tensor, σ , can also be rewritten as:

(A -8)

SC

2 σ = − pI + qn 3

Since we have Eqs. (A-6) and (A-8), the evolution laws of internal state variables can be reorganized as the functions of scalar variables. − pε&vp + qε&dp σ : ε& p = (1 − f ) (1 − Dshear ) σ m (1 − f ) (1 − Dshear ) σ m

M AN U

H& 1 = ε&mp =

(A-9)

H& 2 = f& = (1 − f ) ε& p : I + (1 − ωσ ) Aε&mp = (1 − f ) ε&vp + (1 − ωσ ) Aε&mp

(A-10)

s : ε& H& 3 = D& shear = kω Dshearωσ = kω Dshearωσε&dp q

(A-11)

A.3 Return mapping algorithm

TE D

p

Since linear isotropic elasticity is assumed before the material reaches the yield point, the stress tensor at the end

EP

of each increment, t + ∆t , is given by

AC C

σt +∆t = Ce : εet +∆t = σT − Ce : ∆εp

(A-11)

with σ T is the elastic trial stress, which can be expressed as

(

σ T = Ce : εet + ∆ε

)

(A-12)

where εet is the elastic strain at the beginning of time increment and ∆εε is the total strain increment. As the section A.2, the backward stress of Eq. (A-11) can be decomposed as

σt +∆t = σT − K∆εvpI − 2G∆εdpn

(A-13)

where ∆ε vp and ∆ε dp are the volumetric and deviatoric parts of plastic strain increments. It can be seen from the Eq.(A-13) that in the deviatoric stress space the backward stress returns to the yield surface along n , which implies

ACCEPTED MANUSCRIPT that n and sT are coaxial and n can be simply determined from the stress predictor σ T , as

n=

3 sT 2 qT

(A -14)

where sT and qT are determined from σ T . Based on the previous decomposition, the return mapping process can I

and n , respectively:

RI PT

also be decoupled into two parts along

pt +∆t = p T + K ∆ε vp , qt +∆t = qT − 3G∆ε dp

(A-15)

Eliminating ∆λ from flow rule (Eq. (A-7)) yields Eq. (A-16). By choosing ∆εvp and ∆ε dp as the primary unknowns,

 ∂Φ   ∂Φ  F1 = ∆ε vp  + ∆ε dp  =0    ∂q t +∆t  ∂p t +∆t

(

* F2 = Φ p, q, ε mp , f * , Dshear

)

t +∆t

=0

A.4 Numerical integration procedures

M AN U

incrementally by Newton–Raphson iterative method, as:

SC

Eq. (A-16) in combine with the yield function gives the basic nonlinear equations, which need to be solved

(A-16)

(A-17)

TE D

In VUMAT, the integration of the rate form constitutive equation is carried out incrementally at Gauss points. The solution of stress tensor and internal state variables are assumed to be known at the start of each increment (at t time) and need to be updated at the end of each strain increment (at t + ∆t time). The main steps of the implementation are

EP

summarized as below:

Step1: initiate the stress tensor, strain tensor, internal variables at t = 0 time:

AC C

Step 2: read the variables values of σ t , εt , ε m,p t , ft , Dshear, t at the start of an increment and the strain increment ∆εε ; Step 3: calculate the elastic predictor σTt +∆t , pTt +∆t , qTt +∆t :

σT = σt + Ce : ∆ε ,

1/ 2 1 3 p T = − σ T : I , q T =  s T : sT  ; 3 2 

Step4: calculate the value of yield function:

(

)

* F2 = Φ pT , qT , εm,p t , f *t , Dshear, t ,

If F2 < 0 , the current state is pure elastic and σt +∆t = σT , then go back to step 2.

ACCEPTED MANUSCRIPT If F2 ≥ 0 , the current state is elastic-plastic, then go to step 5 to calculate the plastic correction; Step 5: calculate the corrections cv and cd for ∆εvp and ∆ε dp which are required in solving Eqs. (A-16) and (A-17) by Newton–Raphson iteration. cv and cd are obtained through solving the linearized equations:  ∂F1     k +1  p     ( F )k   ∂∆ε d t +∆t  ( cv )t +∆t  1 t +∆t ,  = − k k k + 1     F ( )  ∂F2   ( cd )t +∆t   2 t +∆t    p    ∂∆ε d t +∆t  k

k

  t +∆t k

  t +∆t

where the upper index k is the iterative step counter.

= ∆ε vp

( ∆ε )

= ∆ε dp

p k +1 v t + ∆t

p k +1 d t + ∆t

(

)

(

)

+ ( cv )t +∆t , k +1

k

t +∆t

M AN U

( ∆ε )

SC

Step 6: update new guesses for ∆εvp and ∆ε dp :

+ ( cd )t +∆t ; k +1

k

t +∆t

RI PT

 ∂F 1   ∂∆ε vp    ∂F2 p   ∂∆ε v

(

Step 7: calculate the internal variables increments. It should be noted that ∆ε mp acquired by solving the following nonlinear equations in a local iteration:

)

=

(

− ptk++∆1t ∆ε vp

(

)(

k +1

(

t +∆t

k +1

)(

)

k +1

k +1

)

t +∆t k +1

EP

t +∆t

(

(

= kω ( Dshear )t +∆t ωσ ∆ε dp k +1

AC C

k +1





(

ptk++∆1t = p T t +∆t + K ∆ε vp

(

)

qtk++∆1t = q T t + ∆t − 3G ∆ε dp

k +1

t + ∆t

)

,

k +1

t + ∆t

;

Step 8: update the internal variables: p k +1 m t + ∆t

)

+ (1 − ωσ ) Atk++∆1t ∆ε mp

(σ m )tk++∆1t =σ m  (ε mp )t +∆t  ,

(ε )

(

+ qtk++∆1t ∆ε dp

1 − ft k++∆1t 1 − ( Dshear )t +∆t (σ m )t +∆t

∆ft k++∆1t = 1 − ft k++∆1t ∆ε vp

( ∆Dshear )tk++∆1t

)

TE D

(

k +1 ∆ε mp t +∆ t

( ) + ( ∆ε )

= ε mp

t

ftk+∆+1t = ft +∆ftk+∆+1t ,

p k +1 m t + ∆t

,

)

k +1

t +∆ t

,

,

)

k +1

t +∆t

,

)

k +1 t +∆t

k +1

, ∆ftk+∆+1t , ( ∆Dshear )t +∆t are

ACCEPTED MANUSCRIPT

( Dshear )tk++∆1t = ( Dshear )t + ( ∆Dshear )tk++∆1t ; Step 9: check the yield condition with current stresses and internal states variables. The global iterations terminate when the values of F1 –and F2 are less than a specified tolerance 10-7:

)

k +1

 ∂Φ  + ∆ε dp   t +∆t ∂q  t + ∆ t k +1

(

)

k +1

 ∂Φ    t +∆t ∂p   t +∆ t k +1

( ) ,( f ) ,( D )

F2 = Φ  ptk++∆1t , qtk++∆1t , ε mp 

k +1

* k +1

t +∆t

t +∆t

k +1 * shear t +∆t

RI PT

(

F1 = ∆ε vp

 , 

SC

If F1 > tolerence or F2 > tolerence , go back to step 5 and start a new iteration.

Else, the ∆εvp and ∆ε dp are converged solution to the basic nonlinear Eqs. (A-16) and (A-17)–and the

M AN U

iterating of current step time is finished. Output the stresses, strain and state variables by the following equations and then go back to step 2 for a next time increment.

k +1 pt +∆t = ptk++∆1t , qt +∆t = qtk++∆1t , ( ε mp )t +∆t = ( ε mp )t +∆t , ft +∆t = ftk+∆+1t , ( Dshear )t +∆t = ( Dshear )t +∆t , k +1

2 qt + ∆t n , 3

εt +∆t = εt + ∆ε

( ∆ε ) p

t + ∆t

=

TE D

σ t +∆t = − pt +∆t I +

1 ( ∆ε vp )t + ∆t I + ( ∆ε dp )t + ∆t n , 3

EP

εpt +∆t = εpt + ( ∆εp ) , t +∆t

AC C

ε e t + ∆ t = ε t + ∆t − ε p t + ∆t

Acknowledgments

We gratefully acknowledge Prof. M. Kikuchi (Department of Mechanical Engineering, Tokyo University of Science) for the support of all the experiments in this study.

Reference [1]. A. Needleman, A Continuum Model for Void Nucleation by Inclusion Debonding. J. Appl. Mech.54 (1987) 525–531. [2]. W.M. Garrison, N.R. Moody, Ductile Fracture. J. Phys. Chem. Solids, 48 (1987) 1035–1074. [3]. J. Koplik, A. Needleman, Void growth and coalescence in porous plastic solids. Int. J. Solids Struct. 24 (1988) 835–853. [4]. V. Tvergaard, Material failure by void growth to coalescence. Adv. Appl. Mech. 27 (1990) 83–151.

ACCEPTED MANUSCRIPT

RI PT

[5]. A.A. Benzerga, J.B. Leblond, Ductile fracture by void growth to coalescence. Adv. Appl. Mech. 44 (2011) 169–305. [6]. R. Kiran, K. Khandelwal, A triaxiality and Lode parameter dependent ductile fracture criterion. Engng. Fract. Mech. 128 (2014) 121–138. [7]. J. Besson, Continuum models of ductile fracture: a review. Int. J. Damage Mech. 19 (2010) 3–52. [8]. I. Barsoum, J. Faleskog, Rupture mechanisms in combined tension and shear—experiments. Int. J. Solids Struct. 44 (2007) 1768–1786. [9]. A. Weck, D.S. Wilkinson, H. Toda, E. Maire, Second and third visualization of ductile fracture. Adv. Eng. Mater. 8 (2006) 469 72.

AC C

EP

TE D

M AN U

SC

[10]. T. Pardoen, J.W. Hutchinson, An extended model for void growth and coalescence. J. Mech. Phys. Solids 48 (12) (2000) 2467–2512. [11]. V. Tvergaard, Behavior of voids in a shear field. Int. J. Fract. 158 (2009):41–49. [12]. V. Tvergaard, K Nielsen, Relations between a micro-mechanical model and a damage model for ductile failure in shear. J. Mech. Phys. Solids 58 (2010) 1243–1252. [13]. F.A. McClintock, A criterion for ductile fracture by the growth of holes. J. Appl. Mech. 35 (1968) 363–371. [14]. J.R. Rice, D.M. Tracey, On the enlargement of voids in triaxial stress fields. J. Mech. Phys. Solids 17 (1969), 201–217. [15]. A.L. Gurson, Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I– Yield Criteria and Flow Rules for Porous Ductile Media. J. Eng. Mater. Technol. 99 (1977) 2–15. [16]. G. Rousselier, Ductile fracture models and their potential in local approach of fracture. Nucl. Eng. Des. 105 (1987) 97–111. [17]. C.C. Chu, A. Needleman, Void nucleation effects in biaxial stretched sheets. J. Eng. Mater. Technol. 102 (3) (1980) 249–256. [18]. V. Tvergaard, A. Needleman, Analysis of the cup-cone fracture in a round tensile bar. Acta Metall. 32 (1) (1984) 157–169. [19]. M. Gologanu, J.B. Leblond, J. Devaux, Approximate models for ductile metals containing non-spherical voids — case of axisymmetric prolate ellipsoidal cavities. J. Mech. Phys. Solids 41 (11) (1993) 1723–1754. [20]. M. Gologanu, J.B. Leblond, J. Devaux, Approximate models for ductile metals containing non-spherical voids — case of axisymmetric oblate ellipsoidal cavities. J. Mech. Phys. Solids 116 (1994) 290–297. [21]. J. Wen, Y. Huang, K.C. Hwang, C. Liu, M. Li, The modified Gurson model accounting for the void size effect. Int. J. Plasticity 21 (2) (2005) 381–395. [22]. V. Monchiet, G. Bonnet, A Gurson-type model accounting for void size effects. Int. J. Solids Struct. 50 (2013) 320–327. [23]. A.A. Benzerga, J. Besson, Plastic potentials for anisotropic porous solids. Eur. J. Mech. – A/Solids 20 (3) (2001) 397–434. [24]. A.A. Benzerga, J. Besson, Anisotropic ductile fracture: Part II: theory. A. Pineau, Acta Mater. 52 (15) (2004) 4639–4650. [25]. Z.Y Chen, X.H. Dong, The GTN damage model based on Hill’48 anisotropic yield criterion and its application in sheet metal forming. Comp. Mater. Sci. 44 (2009) 1013 1021. [26]. K. Nahshon, J. Hutchinson, Modification of the Gurson Model for shear failure, Eur. J. Mech. – A/Solids 27 (1) (2008) 1–17. [27]. L. Xue, Constitutive modeling of void shearing effect in ductile fracture of porous materials Eng. Fract. Mech.

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

75 (11) (2008) 3343–3366. [28]. K.L. Nielsen, V. Tvergaard, Ductile shear failure or plug failure of spot welds modelled by modified Gurson model. Eng. Fract. Mech. 77 (7) (2010) 1031–1047. [29]. J. Jackiewicz, Use of a modified Gurson model approach for the simulation of ductile fracture by growth and coalescence of microvoids under low, medium and high stress triaxiality loadings Eng. Fract. Mech. 78 (2011) 487–502. [30]. Z.Y. Xue, J. Faleskog, J.W. Hutchinson, Tension–torsion fracture experiments – Part II: Simulations with the extended Gurson model and a ductile fracture criterion based on plastic strain. Int. J. Solids Struct. 50 (2013) 4258–4269. [31]. L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, Evaluation of shear mechanisms and influence of the calibration point on the numerical results of the GTN model. Int. J. Mech. Sci. 75(2013)407–422. [32]. M. Achouri, G. Germain, P.D. Santo, D. Saidane, Numerical integration of an advanced Gurson model for shear loading: Application to the blanking process. Comp. Mater. Sci. 72 (2013) 62–67. [33]. T.S. Cao, E. Maire, C. Verdu, C. Bobadilla, P. Lasne, P. Montmitonnet, P.O. Bouchard, Characterization of ductile damage for a high carbon steel using 3D X-ray micro-tomography and mechanical tests – Application to the identification of a shear modified GTN model. Comp. Mater. Sci. 84 (2014) 175–187. [34]. L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, An extended GTN model for ductile fracture under high and low stress triaxiality. Int. J. Plasticity 54 (2014) 193–228. [35]. J. Zhou, X.S. Gao, J.C. Sobotka, B.A. Webler, B.V. Cokeram, On the extension of the Gurson-type porous plasticity models for prediction of ductile fracture under shear-dominated conditionsInt. J. Solids Struct. 51 (2014) 3273–3291. [36]. J. Lemaitre, A continuous damage mechanics model for ductile fracture. J. Eng. Mater. Technol. 107 (1985) 83–89. [37]. Y. Bao, T. Wierzbicki, On fracture locus in the equivalent strain and stress triaxiality space. Int. J. Mech. Sci. 46 (2004) 81–98. [38]. X.S. Gao, T.T. Zhang, M. Hayden, C. Roe, Effects of the stress state on plasticity and ductile failure of an aluminum 5083 alloy. Int. J. Plasticity 25 (2009) 2366–2382. [39]. Y.J. Liu, Q. Sun, X.L. Fan, T. Suo, A stress-invariant based multi-parameters ductile progressive fracture model. Mater. Sci. Eng. A 576 (1) (2013) 337–345. [40]. X.L. Fan, Q. Sun, Y.J. Liu, A modified ductile fracture model incorporating synergistic effects of pressure and Lode angle. Int. J. Appl. Mech. 4 (2012) 1250022. [41]. T. Wierzbicki, Y.B. Bao, Y.W. Lee, Y.L. Bai, Calibration and evaluation of seven fracture models. Int. J. of Mech. Sci. 47 (2005) 719–743 [42]. L. Xue, Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading. Int. J. Solids Struct. 44 (2007) 5163–5181. [43]. S.M. Graham, T.T. Zhang, X.S. Gao, M. Hayden, Development of a combined tension–torsion experiment for calibration of ductile fracture models under conditions of low triaxiality. Int. J. Mech. Sci. 54 (2012) 172–181. [44]. V. Tvergaard, Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17 (1981) 389–407. [45]. V. Tvergaard, On localization in ductile materials containing spherical voids. Int. J. Fract. 18 (1982) 237–252. [46]. J. Lemaitre, R. Desmorat, M. Sauzay, Anisotropic damage law of evolution. Eur. J. Mech. A/Solids 19 (2000) 187–208. [47]. N. Aravas, On the numerical integration of a class of pressure-dependent plasticity models Int. J. Numer. Meth. Eng. 24 (1987) 1395–1416.

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[48]. Z.L. Zhang, Explicit consistent tangent moduli with a return mapping algorithm for pressure-dependent elastoplasticity models. Comp. Meth. Appl. Mech. Engng.121 (1995a) 1–4:29–44. [49]. Z.L. Zhang, On the accuracies of numerical integration algorithms for Gurson-based pressure-dependent elastoplastic constitutive models. Comp. Meth. Appl. Mech. Engng. 121 (1995b). 1–4:15–28.

ACCEPTED MANUSCRIPT

Highlights A modification to GTN model is proposed for the simulation of ductile fracture under high, low and negative stress triaxialities.

2.

Two distinctive damage parameters are introduced into yield function to capture the degradation process.

3.

A new function of triaxiality and Lode angle is adopted as shear damage weight function.

4.

The new model is validated by comparing numerical results with the experimental data of specimens including axisymmetric tensile bar, flat grooved plate, torsion tube and compression cylinder.

AC C

EP

TE D

M AN U

SC

RI PT

1.