A finite strain continuum damage model for simulating ductile fracture in steels

Engineering Fracture Mechanics 116 (2014) 172–189

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

A finite strain continuum damage model for simulating ductile fracture in steels Kapil Khandelwal a,⇑, Sherif El-Tawil b a b

Dept. of Civil & Env. Engrg. & Earth Sci., U. of Notre Dame, Notre Dame, IN 46556, United States Dept. of Civil and Env. Eng., U. of Michigan, Ann Arbor, MI 48109, United States

a r t i c l e

i n f o

Article history: Received 10 June 2013 Received in revised form 19 November 2013 Accepted 17 December 2013

Keywords: Ductile fracture Damage mechanics Plasticity Stress triaxiality A36 steel Finite element analysis

a b s t r a c t A new coupled plasticity-damage model is proposed within a finite deformation framework for modeling ductile fracture in ASTM A36 structural steels. Damage mechanics principles of effective stress and strain equivalence are employed to formulate a new constitutive model for simulating damage due to the physical processes associated with microvoid nucleation, growth and coalescence. A scalar damage variable is used to track the micro-structural changes that occur during the ductile fracture process. The model is calibrated and validated by comparing its response to the results obtained from experimental testing of four symmetrically and asymmetrically notched ASTM A36 steel specimens. The proposed model is shown to successfully model failure due to ductile fracture under stress states typically observed in structural engineering applications. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Ductile fracture in steel is a multistep process resulting from microvoid nucleation, growth and coalescence of microvoids in a plastically deforming material [1]. Microvoids typically nucleate at inclusions either by decohesion/debonding of the inclusion matrix interface or by fracture of the inclusion itself [2]. Void nucleation is followed by a void growth stage where microvoids grow under the applied stress state, and the final phase of ductile fracture occurs when adjacent microvoids coalesce together into a crack. The microvoid growth stage is highly influenced by the state of stress in the material, measured in terms of stress triaxiality (T := rh/rm) and defined as the ratio of the hydrostatic stress (rh) and Mises stress (rm). At high stress triaxiality (T > 1.0), void growth is usually volumetric and microvoid coalescence is due to internal necking of the intervoid ligaments. However, at low stress triaxiality (T < 1.0), microvoids may also elongate and final fracture initiation is by shear localization between the microvoids [3–5]. For instance, high triaxiality is observed in connection regions and shear links [6], while reduced beam sections are under low triaxiality [7]. This paper is concerned with constitutive modeling of damage attributed to microvoid growth in structural steels. The ductile fracture process can be modeled using either uncoupled or coupled models. The constitutive relationship inherent to uncoupled models is assumed to be unaffected by material damage, such as that caused by nucleation and growth of microvoids and usually employ J2 plasticity as a constitutive model. These models assume that fracture initiation occurs upon the onset of microvoid coalescence, detected by means of a fracture criterion in the post processing step [8–11]. Coupled models can be broadly classified into micromechanics and damage-mechanics based models. The former are derived by explicitly modeling the microstructure through unit cell simulations [12]. The latter simulate microvoid evolution ⇑ Corresponding author. E-mail addresses: [email protected] (K. Khandelwal), [email protected] (S. El-Tawil). 0013-7944/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2013.12.009

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

173

Nomenclature a0, at, a2, a3 constitutive material parameters ad, ap damage and plastic internal variables be elastic Finger tensor adc critical damage parameter Cp plastic Cauchy tensor Dint internal energy dissipation rate en ; ec nucleation and coalescence strains De smoothening parameter F, Fe, Fp total, elastic and plastic deformation gradient c plastic consistency parameter J Jacobian determinant k, l bulk and shear modulus ^ a ; a ¼ 1; 2; 3 Eigen values and Eigen vectors of elastic Finger tensor kea ; n W, We, Wp total, elastic and plastic free energy functions /p plastic potential P2 ; P1 pullback and pushforward operators 2 rh, rm hydrostatic and Mises stress T stress triaxiality s; s Kirchhoff and effective Kirchhoff stress tensor sv ol ; sdev volumetric and deviatoric component of Kirchhoff stress tensor T 1T 2 threshold functions  g ; g ~ ; a; s; p  model implementation parameters ee ; b; K p ; rmax ; rY ; a saturation hardening parameters fp hardening function

by introducing damage internal variables into the constitutive relationship [13,14]. Irrespective of these differences, the constitutive equations in the case of coupled models have internal variables and evolution equations to simulate damage due to microvoid nucleation, growth and coalescence during the ductile fracture process. The Gurson model [12] and its extensions [15–17] belong to the category of micromechanical models. Gurson [12] first proposed an approximate yield criterion and flow rules for rigid plastic porous material using an upper bound variational approach. Subsequently, Tvergaard [18] performed bifurcation studies on unit cells and modified the Gurson yield criterion to account for the void interaction effects. Tvergaard and Needleman [19] further modified the Gurson model to account for the rapid loss of load carrying capacity of the material at the onset of microvoid coalescence. This form of Gurson model, also referred to as Gurson–Tvergaard–Needleman (GTN) model, is popularly used to model ductile fracture in metals [20]. The relative success of the GTN model is primarily due to the incorporation of a sufficient number of adjustable parameters to enable curve fitting of the desired experimental results. However, the GTN model has nine material parameters and these parameters are difficult to calibrate [21]; improper choice of model parameters can easily lead to convergence issues [22]. Furthermore, its implementation is computationally expensive as four coupled nonlinear equations are needed to be solved at each integration point during finite element simulation, and without careful implementation to ameliorate numerical difficulties, convergence may be difficult to obtain [22]. The Gurson model has also been recently extended to account for the Lode parameter [23] and to account for damage due to shear [24], with a corresponding increase in the number of model parameters that are again difficult to calibrate. GTN model parameters are dependent on the state of stress within the material [25]. Therefore, parameters calibrated to a particular experiment cannot be used to simulate fracture in situations which differ largely from the original experiment used for model calibration. Also, in the GTN model, the void nucleation strain is independent of hydrostatic stress. However, void nucleation models proposed by Argon et al. [26] and Goods and Brown [27], and metallurgical studies [28–30] have shown that void nucleation strain is sensitive to hydrostatic stress and that the nucleation strain decreases with increase in mean stress. Furthermore, the chosen parameters are sometimes found to be inconsistent with metallurgical results and may have no physical significance [31]. With the growing interest in the structural engineering community to model failure of steel structural systems, computationally efficient models are needed to simulate ductile fracture in structural steels [32]. In this paper a new continuum damage mechanics-based plasticity model is proposed within a finite deformation framework for modeling the micromechanical process of ductile fracture in ASTM A36 structural steels. Damage mechanics principles of effective stress [33] and strain equivalence [34] are employed to formulate a constitutive model for simulation of damage due to microvoids nucleation, growth and coalescence. The microscale response of the damage model arises from the addition of internal state variables known as damage variables that represent the change in microstructure of the material undergoing a deformation process.

174

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

In the proposed model, a scalar damage variable is used to represent the damage that arises due to micro-structural evolution during the ductile fracture process. In particular, the three stages of ductile fracture initiation: microvoid nucleation, growth and coalescence are modeled by an appropriate evolution function for the damage variable. There is a distinct computational advantage over the GTN model as only two nonlinear equations need to be solved at each integration point as opposed to four nonlinear equations in the GTN model. The void nucleation strain is taken to be a function of stress triaxiality, and the model is calibrated and validated by comparing its response to the results obtained from experimental testing of four symmetrically and asymmetrically notched steel bar specimens. The paper is organized as follows: Section 2 presents the kinematics of finite deformation, and the kinetics of damage and plasticity is described in Section 3. In Section 4 numerical implementation of the proposed model is described. Section 5 presents prediction capabilities of the proposed model and comparison with experimental results. Finally, important conclusions of this study are summarized in Section 6. 2. Kinematics Multiplicative decomposition of the deformation gradient tensor (F) that maps the infinitesimal material neighborhood from the initial configuration ðB0 Þ to the current configuration ðBt Þ is employed (Fig. 1). In this multiplicative plasticity model for finite deformation, an intermediate stress-free configuration ðBst Þ, referred to as structural space, is introduced [35]. Using this multiplicative decomposition, the deformation gradient is decomposed into elastic and plastic parts as follows:

F ¼ Fe  Fp

ð1Þ p

In case of metal plasticity, from the micromechanical point of view, F is an internal variable related to the amount of plastic flow (dislocation movement) associated with the underlying crystalline structure. However, from the phenomenological standpoint the intermediate configuration defines the local, stress free unloaded configuration. The elastic deformation gradient (Fe) represents the elastic distortion in an infinitesimal neighborhood of a material point due to the stretch and rotation of the microstructure. Using the intermediate configuration the following tensors can be defined using the pull back ðP2 Þ and push forward ðP1 2 Þ operations [35] (Fig. 1): def

e

C p1 ¼ F p1  F pT ¼ P2 ½b  e def

b ¼ Fe  F

eT

p1 ¼ P1  2 ½C

ð2Þ

where be is the elastic Finger tensor that characterizes the elastic deformation with respect to stress free intermediate conp1 figuration, P2 is the contravariant pullback, P1 is the inverse of plastic Cauchy 2 is the contravariant pushforward, and C tensor that characterizes the plastic deformation at a material point. Full kinematic description of this multiplicative framework can be found in Ref. [35] (chapter 91). 3. Kinetics The concept of effective stress [33] and the strain equivalence principle [13] are used to formulate the proposed model. Damage due to microvoids nucleation, growth and coalescence is modeled using a scalar damage internal variable ad. Then, according to the effective stress concept, the actual stress resisted by the undamaged or intact material surface is expressed as follows:

s def ¼

s

ð3Þ

1  ad

 is the effective Kirchhoff stress tensor and s is the Kirchhoff stress tensor. Note that the Kirchhoff stress tensor is where s related to the Cauchy stress tensor (r) as follows: s = Jr, where J = det F. The strain equivalence principle [34], which states that ‘‘any strain constitutive equation for a damaged material may be derived in the same way as for the virgin material except that

Fig. 1. Kinematics of finite deformation: multiplicative decomposition of deformation gradient.

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

175

the usual stress is replaced by the effective stress’’, is then used to formulate material constitutive law. Three superscripts – e, p and d – are used in the following discussion, where e represent variables associated with elastic process, p represent variables associated with plastic processes, and d represent variables associated with damage processes. For an isotropic material, the functional form of the Helmholtz free energy is given by [36,37]:

W ¼ Wðbe ; ad ; ap Þ

ð4Þ

where ad is an internal damage variable described above and ap is a strain-like internal variable that describe the state of material induced by dislocation pileups and other micro-defects. In this study the following functional form of the Helmholtz free energy is employed:

W ¼ ð1  ad ÞWe ðbe Þ þ Wp ðap Þ

ð5Þ

A similar free energy decomposition was proposed by Mahnken [38]; however, the damage model in that work was formulated in the intermediate configuration. Steinmann [37] proposed an alternate form of Helmholtz free energy where the damage is only considered to be a function of the deviatoric stress; however, such a formulation is inadequate to model damage due to microvoid growth under high triaxiality. With the free energy potential given by Eq. (5), the dissipation inequality, in terms of energy dissipation rate, is expressed as:

Dint ¼



s  2ð1  ad Þbe 

@ We e @b



  e  1 e1 e @W : d þ 2ð1  ad Þb  e :  b  F  C_ p1  F T þ We a_ d þ fp a_ p P 0 2 @b

ð6Þ

Using the Coleman–Noll argument [39], the Kirchhoff stress tensor is given by:

s ¼ 2ð1  ad Þbe 

@ We e @b

ð7Þ

The internal mechanical energy dissipation rate can now be expressed as follows:

Dint ¼

  e  1 e1 e @W :  b  F  C_ p1  F T þ We a_ d þ fp a_ p P 0 2ð1  ad Þb e 2 @b

ð8Þ

Thus, for any kinetically admissible deformation process, the inequality given in Eq. (8) should be satisfied. For describing the evolution of plastic internal variable (ap) a potential approach is adopted, wherein a plastic potential function is introduced in terms of the effective Kirchhoff stress as follows:

; fp Þ ¼ ks dev k  /p ðs

rffiffiffi 2 ðry þ fp Þ 3

ð9Þ

dev is the deviatoric component of the effective Kirwhere ry is the yield strength, fp = fp(ap) is the hardening variable and s . Thus, Eq. (9) represents the von Mises yield criterion in terms of the effective Kirchhoff stress tensor; chhoff stress tensor s ; fp Þ is a convex function of s . Evolution of the plastic internal variables is then expressed using associative flow also, /p ðs rules [36]:

a_ p ¼ c

rffiffiffi @/p 2 ¼ c 3 @fp

ð10Þ

1 e1 @/p c @/p  b  F  C_ p1  F T ¼ c ¼  2 @s 1  ad @ s

ð11Þ

where c > 0 is a scalar consistency parameter. The plasticity part of the model is completed by specifying the Kuhn–Tucker complementarity conditions:

c > 0; /p < 0 and c/p ¼ 0

ð12Þ

and the plastic consistency condition:

c/_ p ¼ 0

ð13Þ d

A typical approach for modeling the evolution of the damage internal variable (a ) is to introduce an additional potential in terms of the variable conjugate to ad. Such an approach was suggested by Lemaitre [13] for small deformation plasticitydamage models, and has also been used by Steinmann [37] and Mahnken [38] for large deformation elasto-plastic damage models. In particular, a quadratic damage potential is introduced in these approaches. Such a model trivially satisfies the dissipation inequality but does not explicitly consider damage due to the deviatoric and hydrostatic components of stress. In contrast, the damage evolution law proposed in this paper explicitly considers both the effect of deviatoric and volumetric components of stress in a phenomenological sense. The evolution law proposed in this study is given in Eq. (14):

176

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

a_ d ¼ a0 T 1 ðap Þa_ p þ a1 TT 1 ðap Þa_ p þ a2 expða3 ap ÞT 2 ðap Þa_ p def ksv ol k T¼ ksdev k

ð14Þ

where svol and sdev are the volumetric and deviatoric components of the Kirchhoff stress tensor s and a0, a1, a2, and a3 are material parameters. The first term in the evolution of ad (Eq. (14)) represents the damage attributed to the elongation of voids under applied shear, and a0 can be considered as a shear parameter. The second term represents the damage due to void growth under high triaxiality, and a1 can be considered as a triaxiality parameter. In addition, a threshold function, T 1 ðap Þ, is introduced so that damage is activated only when the nucleation strain (en) is exceeded. This quantity is taken equal to the plastic strain at which the microvoids nucleate. Thus, the first and second terms of the damage evolution equation model the material damage due to elongation and volumetric growth of microvoids, respectively, after microvoid nucleation. The last term represents material damage due to microvoid coalescence, and a2 and a3 can be considered as void coalescence parameters. An exponential function is chosen to model the rapid disintegration of material that occurs at this stage. This is done through the threshold function, T 2 ðap Þ; which is activated only when a coalescence strain (ec) is attained. The plastic strain limit ec represents the plastic strain value at which the microvoids start to coalesce due to localization in the intervoid matrix. Also, in the model implementation void coalescence strain ec is assumed as a material constant, while nucleation strain en is assumed a material-specific parameter that varies with triaxiality. The functional form of threshold functions is chosen as Hermitian polynomials to facilitate numerical implementation because of its smoothness. The functional forms of the threshold functions are given as follows [38]:

8 T 1 ðap ; en ; DeÞ ¼ 0 if ap < en > > < h i 2 p p T 1 ðap ; en ; DeÞ ¼ ða Dee2n Þ 3  2ðaDe en Þ if en < ap < en þ De > > : and T 1 ðap ; en ; DeÞ ¼ 1 if ap P en þ De

ð15Þ

8 T ðap ; ec ; DeÞ ¼ 0 if ap < ecl > > < 2 h i 2 p p T 2 ðap ; ec ; DeÞ ¼ ða Dee2c Þ 3  2ðaDe ec Þ if ec < ap < ec þ De > > : and T 2 ðap ; ec ; DeÞ ¼ 1 if ap P ec þ De

ð16Þ

In Eqs. (15) and (16) the parameter, De, represents the range over which the damage activation process is smoothened out and can be regarded as a material parameter. 4. Numerical implementation The above described hyperelastic–plastic-damage model is implemented in the explicit finite element program LS-DYNA [40]. A quadratic logarithmic model for free energy, We, in terms of the Eigen values (kea , a = 1, 2, 3) of the elastic Finger tensor be is assumed with the following functional form: e

1 2

2

2

2

2

We ðb Þ ¼ kðln ke1 þ ln ke2 þ ln ke3 Þ þ l½ðln ke1 Þ þ ðln ke2 Þ þ ðln ke3 Þ  

l 3

ðln J e Þ

2

ð17Þ

where k is the bulk modulus and l is the shear modulus of the material. This free energy based on logarithmic strain meae 1=2 sures satisfies the growth conditions, i.e. We ? 1 as J e ¼ ke1 ke2 ke3 ¼ ðdet b Þ ! 0 and We ? 1 as Je ? 1. Moreover, We is not a convex function of Je and this functional form cannot be used for problems with very large elastic strains [36]. Although plastic strains can be arbitrary large, this model provides an excellent approximation for moderately large elastic strains, as discussed later on. To advance the solution within an incremental scheme in a finite element framework, the flow rule and the evolution equations for the hardening and damage variables (Eqs. 10, 11, and 14) have to be integrated over a finite time step Dt = tn+1  tn; such a computation is typically carried out at an integration point. The known values at time tn include the Cauchy stress p d tensor (rn), deformation gradient Fn and internal variables ðC p1 n ; an ; an Þ. Within the user defined material subroutine, the known values at time tn+1 include the deformation gradient (Fn+1). The objective of the integration scheme is to compute p d the Cauchy stress tensor (rn+1) and internal variables ðC p1 nþ1 ; anþ1 ; anþ1 Þ at time tn+1. In the present work, integration of the damage and the plastic internal evolution equations is carried via the implicit Euler backward method [41] and the integration algorithm closely follows the implementation suggested by Simo [36] for finite deformation plasticity. The following additional terms are introduced to facilitate implementation of the algorithm:

2

3 2 3 ee1 ln ke1 6 e7 6 7 ¼ 4 e2 5 ¼ 4 ln ke2 5; ee3 ln ke3

e def

e

2

3 1 b 6 7 b ¼ 4 b2 5; 3 b  def

2

3 g 1 6 g ¼ 4 g 2 7 5; g 3 def

2

3

2 3

1 g~ 1 g 6~ 7 6 7  g ¼  ¼ 4 g2 5; 1 ¼ 4 1 5 kgk 1 g~ 3 def

ð18Þ

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

  1 def a ¼ k1  1 þ 2l I  1  1 3

177

ð19Þ

def b ¼ aee

s ¼

3 X ^a  n ^aÞ ba ðn

ð20Þ

a¼1

1 3

  ð1T bÞ1  g def ¼b s ¼

ð21Þ

3 X

ga ðn^ a  n^ a Þ

a¼1

¼ p

3 X ^a  n ^aÞ kð1T ee Þðn

ð22Þ

a¼1

^ a g are the eigen pairs of be and I is a 3  3 identity matrix. The plastic potential can now be expressed in terms of where fkea ; n  as follows: deviatoric component of stress g p

k  ; fp Þ ¼ kg / ðs

rffiffiffi 2 ½rY  fp ðap Þ 3

ð23Þ

The integration of the plasticity models is usually carried out in two steps [41]. In the first step known as the ‘‘elastic’’ or ‘‘trial’’ step, all inelasticity is frozen. If the yield condition is not exceeded then the trial step gives the correct state. However, if the yield condition is exceeded the algorithm proceeds to the ‘‘plastic’’ step, where further computations are carried out assuming that plastic flow occurs and the consistency condition is enforced in this step. The two steps of this algorithm are briefly described next. 4.1. Trial step: freeze plastic flow tr

tr

tr

etr

p d For the case of no plastic flow the trial values: ctrnþ1 , C p1 nþ1 , anþ1 , anþ1 , and bnþ1 are given as follows:

ctrnþ1 ¼ 0

ð24Þ

tr

p1 C p1 nþ1 ¼ C n tr

apnþ1 ¼ apn tr

adnþ1 ¼ adn tr

etr

T bnþ1 ¼ F nþ1  C p1 nþ1  F nþ1

Let

n

ð25Þ

J nþ1 ¼ detðF nþ1 Þ tr ^ tra keanþ1 ; n nþ1

o

2 6

ð26Þ

be the eigen pairs of

tr ln ke1nþ1

3

etr bnþ1 ,

then using Eqs. (18)–(22), the following relationships can be obtained:

7

7 eenþ1 ¼ 6 6 ln ke2nþ1 7 tr

4

tr

tr ln ke3nþ1

ð27Þ

5

  tr ¼ aeetr and g tr  1 1T b tr 1  trnþ1 ¼ b b nþ1 nþ1 nþ1 nþ1 3 rffiffiffi  tr i tr 2h  trnþ1 k  /pnþ1 ¼ kg rY  fp apnþ1 3

ð28Þ

ð29Þ

tr

Now if /pnþ1 6 0, this implies that the trial step is admissible and the following elastic updates are carried out:

 trnþ1 ¼ p

3    X tr  tra ^ tra  n k 1T eenþ1 n nþ1 nþ1 a¼1

strnþ1 ¼

3 X



gtranþ1 n^ tranþ1  n^ tranþ1

a¼1



ð30Þ

178

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

strnþ1 ¼ strnþ1 þ p trnþ1

ð31Þ tr

ptr

apnþ1 ¼ anþ1 and adnþ1 ¼ adnþ1

ð32Þ

tr

p1 C p1 nþ1 ¼ C nþ1

snþ1

ð33Þ

  tr trnþ1 ¼ 1  adnþ1 s

rnþ1 ¼

ð34Þ

1 snþ1 J nþ1

ð35Þ

tr

If /pnþ1 > 0, this implies that the trial step is not admissible and the algorithm proceeds to the second step, i.e., the plastic step. 4.2. Plastic flow In this step, the integration of flow rules is carried out using the exponential integrator and the consistency condition is enforced as described next. 4.2.1. Integration of flow rule Integration of flow rule using an exponential integrator yields:

C p1 nþ1 e

bnþ1

! 2Dcnþ1 ¼  exp  nnþ1  F nþ1  C p1 n 1  adnþ1 ! 2Dcnþ1 etr ¼ exp  nnþ1  bnþ1 1  adnþ1 F 1 nþ1

etr

ð36Þ e

s where n ¼ ksk , Dcnþ1  Dtcnþ1 and bnþ1 ¼ F nþ1  C p1  F Tnþ1 . Now using the following spectral forms for nn+1 and bnþ1 : n

nnþ1 ¼

3 X

g~ anþ1 ðn^ anþ1  n^ anþ1 Þ

a¼1 e bnþ1

¼

2 ^ anþ1 keanþ1 ðn

ð37Þ

^ anþ1 Þ n

and combining with Eq. (36) gives: etr bnþ1

" 3 X ¼ exp a¼1

etr

Also, bnþ1 tr

! # 2Dcnþ1 e2 ^ anþ1  n ^ anþ1 Þ ~ ga k ðn 1  adnþ1 nþ1 anþ1

ð38Þ

 P  2tr ^ tra  n  tra Þ , and the uniqueness of the spectral decomposition yields: ¼ 3a¼1 keanþ1 ðn nþ1 nþ1

2 keanþ1

! 2Dcnþ1 2 ke ¼ exp g~ a 1  adnþ1 nþ1 anþ1

^ anþ1 ¼ n ^ tra n nþ1

ð39Þ

a ¼ 1; 2; 3 Taking the logarithm of both sides of Eq. (39) gives: tr

eenþ1 ¼ eenþ1 

Dcnþ1 g~ nþ1 1  adnþ1

ð40Þ

Combining Eqs. (20), (21), (28) and (40) gives:

nþ1 ¼ b tr  2lDcnþ1 ag ~ nþ1 b nþ1 1  adnþ1 ! 2lDcnþ1 ~ trnþ1 kg  trnþ1 k g~ nþ1 kg nþ1 k þ ¼g 1  adnþ1

ð41Þ

Eq. (41) implies that:

 nþ1 k ¼ kg  trnþ1 k  kg

g~ nþ1 ¼ g~ trnþ1

2lDcnþ1 1  adnþ1 ð42Þ

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

179

4.2.2. Enforcing the consistency condition Enforcing the consistency condition at time t = tn+1 gives: rffiffiffi

 nþ1 k  /pnþ1 ¼ kg

 

2 ry þ fp apnþ1 ¼ 0 3

ð43Þ

Combining Eqs. (42) and (43), the consistency condition can be written as follows:

/pnþ1

 trnþ1 k  ¼ kg

2lDcnþ1  1  adnþ1

rffiffiffi  

2 ry þ fp apnþ1 ¼ 0 3

ð44Þ

4.2.3. Integration of damage and plastic internal variables Integrating Eq. (10) using the Euler backward method gives: p nþ1

a

rffiffiffi 2 ¼a þ Dc 3 nþ1 p n

ð45Þ

Integrating Eq. (14) using the Euler backward method gives:

adnþ1 ¼ adn þ where T nþ1

rffiffiffi    

2 Dc ða0 þ a1 T nþ1 ÞT 1 ðapnþ1 Þ þ a2 exp a3 apnþ1 T 2 apnþ1 3 nþ1  nþ1 k  nþ1 k kp kp ¼ ¼  nþ1 k kg kg  trnþ1 k  2lDcdnþ1

ð46Þ

1anþ1

4.2.4. Solution by Newton–Raphson method Eqs. (44) and (46) are nonlinear in Dcn+1 and adnþ1 and these equations are solved together with the Eq. (45) to obtain Dcn+1 and adnþ1 . For solving these equations, Eq. (44) is written as follows:

 trnþ1 k  F 1 ¼ kg

2lDcnþ1  1  adnþ1

rffiffiffi  

2 ry þ fp apnþ1 ¼ 0 3

ð47Þ

Similarly, Eq. (46) is written as follows:

F2 ¼ a

d nþ1

rffiffiffi  

2 a  Dc ða0 þ a1 T nþ1 ÞT 1 ðapnþ1 Þ þ a2 exp a3 apnþ1 T 2 ðapnþ1 Þ ¼ 0 3 nþ1 d n

ð48Þ

Thus, Eqs. (44) and (46) can be considered as two nonlinear simultaneous equations in two variables adnþ1 and Dcn+1, and can ð0Þ ð0Þ be solved by the Newton–Raphson method. Starting with the initial values of Dcnþ1 ¼ 0 and adnþ1 ¼ adn , the following iterations are carried out to obtain the solution:

2  3 2 @F 1 ðkþ1Þ d Dcnþ1 6 7 4 @ Dcnþ1 4  ðkþ1Þ  5 ¼ @F 2 d adnþ1 @ Dcnþ1

@F 1 @ adnþ1 @F 2 @ adnþ1

31 5 

"

ðkÞ

ðkÞ

ðkÞ

ðkÞ

Dcnþ1 ;adnþ1

ðkÞ

F 1 ðadnþ1 ; Dcnþ1 Þ

 F 2 ðadnþ1 ; DcðkÞ nþ1 Þ

# ð49Þ

The incremental updates are obtained as follows: ðkþ1Þ ðkÞ ðkþ1Þ Dcnþ1 ¼ Dcnþ1 þ dðDcnþ1 Þ ðkþ1Þ

ðkÞ

ðkþ1Þ

adnþ1 ¼ adnþ1 þ dðadnþ1 Þ

ð50Þ

The derivatives entering in Eq. (49) are given as follows:

@F 1 2l 2 dfp ¼ þ @ Dcnþ1 1  adnþ1 3 dDcnþ1

ð51Þ

2lDcnþ1 @F 1 ¼  2 @ adnþ1 1  ad

ð52Þ

nþ1

rffiffiffi

@F 2 2 ða0 þ a1 T nþ1 ÞT 1 ðapnþ1 Þ þ a2 expða3 apnþ1 ÞT 2 ðapnþ1 Þ ¼ 3 @ Dcnþ1 " # rffiffiffi rffiffiffi  p  @T nþ1   @T 2 2 @T 1 2 p p p a2 a3 expða3 anþ1 ÞT 2 ðanþ1 Þ  Dc ða0 þ a1 T nþ1 Þ þ a1 T 1 anþ1 þ a2 exp a3 anþ1 þ 3 nþ1 3 @ Dcnþ1 @ Dcnþ1 @ Dcnþ1 ð53Þ @F 2 ¼1 @ adnþ1

rffiffiffi   @T nþ1 2 Dc a1 T 1 apnþ1 3 nþ1 @ adnþ1

ð54Þ

180

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

The derivatives in these equations ative

dfp dDcnþ1

@T nþ1 , @T 1 , @T 1 , @T 2 @ adnþ1 @ Dcnþ1 @ adnþ1 @ Dcnþ1

and

@T 2 @ adnþ1

can be calculated from Eqs. (15), (16) and (46); the deriv-

can be evaluated after the hardening function is specified. After solving these equations, the relevant updates are

obtained using Eqs. (39)–(42), (45) and (46).

5. Model results 5.1. Single element response To investigate the effect of various model parameters, a parametric study is carried out using a single element of size 6.35 mm (0.25 in.), with one point integration. The model set up is shown in Fig. 2. The bulk modulus and shear modulus are taken as follows: k = 166625.6 MPa, and l = 76904.1 MPa. A saturated hardening model is assumed as follows:

fp ðap Þ ¼ K p ap þ ðrmax  rY Þð1  expðaap ÞÞ

ð55Þ

A critical value of the damage parameter, adc , is introduced and the element is removed (is assumed to have failed) when this critical damage parameter is reached. For this parametric study adc is taken as 0.90 so that the element is removed when the stress state is close to zero. The various cases considered for this parametric study are given in Table 1 and the corresponding results are shown in Figs. 3–5. These figures show that a wide range of material responses can be simulated and that the model has the flexibility to address micro-structural changes due to void growth and coalescence. For example, as shown in Fig. 5, the two stages of material response can be modeled: the first descending branch of the curve corresponds to damage due to microvoid growth processes, while the second descending branch corresponds to the rapid disintegration of the material during the microvoid coalescence stage. Clearly, the rate of damage evolution can be controlled by appropriate selection of the model parameters.

5.2. Validation with A36 steels 5.2.1. Experimental test To calibrate and validate the proposed damage model, tests on symmetric and asymmetric notched specimens are carried out. The experimental study involved tension tests on ASTM A36 steel coupons that are notched to a desired geometry as shown in Figs. 6 and 7. The chemical composition of the tested steel by weight percentage is as follows: C = 0.16%, P = 0.017%, S = 0.009% and Si = 0.19%. The specimens are extracted from flat hot rolled bars of dimension 45.7  10.2 mm. The specimen shown in Fig. 6 has a symmetric notch while the specimen in Fig. 7 has an asymmetric notch. These two different notch patterns are chosen to promote different stress conditions in the material. In particular, specimens with symmetrical notches are subjected to higher triaxiality as compared to specimens with asymmetrical notches, which are subjected to lower triaxiality. Thus, these specimens will experience different modes of damage, and are chosen to calibrate the proposed damage model to represent a wide range of triaxiality conditions that can be experienced in structural applications. Two different thicknesses of 6.35 mm and 12.7 mm are used for each of the test specimens. The specimens are tested to failure using a servo-hydraulic load frame under displacement control (Fig. 8) and the data of interest from the tension test is the force versus displacement response of the specimens. An extensometer of 38.1 mm gage length is used for measuring displacement and a loading rate 0.76 mm/min is used in the tests. To facilitate further discussions, symmetrically notched specimens are designated as SN-X, where X represents the thickness of the specimen (X = 1 implies a thickness of 6.35 mm, whereas X = 2 implies a thickness of 12.7 mm). Asymmetrically notch specimens are designated as ASN-X. The yield strength of the steel is 289.5 MPa, and the plastic hardening curve is obtained from a uniaxial test and is shown in Fig. 9.

6.35 mm

6.35 mm

Prescribed Displacement

6.35 mm Fig. 2. Model setup for parametric study.

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

181

Table 1 Parameters considered in the parametric studies. Analysis case

Model parameters

Case 1

a1 = a2 = a3 = 0, ry = 344.7 MPa, rmax = 448.1 MPa, a = 40, Kp = 200 MPa en = 0.2, De = 0.05, adc ¼ 0:9 a0 = 0, 0.25, 0.5, 1, 2.5, 5, 10 a0 = a2 = a3 = 0, ry = 344.7 MPa, rmax = 448.1 MPa, a = 40, Kp = 200 MPa

Case 2

en ¼ 0:2; De ¼ 0:05; adc ¼ 0:9 a1 = 0, 0.25, 0.5, 1, 2.5, 5, 10

ry = 344.7 MPa, rmax = 448.1 MPa, a = 40, Kp = 200 MPa en ¼ 0:2; De ¼ 0:05; adc ¼ 0:9

Case 3

a0 = a1 = 0.5, a2 = 1.0 a3 = 1, 5, 10, 15

700

a0 = 0.0 a0 = 0.25

Stress (MPa)

600 500 400

a0 = 0.5

300 200

a0 = 2.5

100

Element Failure

a0 =10.0

0 0.0

2.0

a0 = 5.0

4.0

6.0

8.0

Displacement (mm) Fig. 3. Parametric study case (a) – effect of parameter a0.

700

a1 = 0.0 a1 = 0.25 a1 = 0.5 a1 =1.0

Stress (Mpa)

600 500 400 300

a1 = 2.5

200

a1 = 5.0

100

a1 = 10.0

0 0.0

2.0

4.0

Element Failure 6.0

8.0

Displacement (mm) Fig. 4. Parametric study case (b) – effect of parameter a1.

5.2.2. Model calibration Calibration of the developed constitutive model is carried out using finite element analysis of the notched specimens used in the experimental part of this study. Failure of steel members results in softening response and is characterized by formation of localization zones or bands associated with high deformations where most of the energy dissipation occurs. The size of such localization zones cannot be resolved by classical continuum mechanics theories due to the lack of an intrinsic length scale in these theories; see for instance [42] and references therein. This leads to finite element results that are mesh dependent. Indeed, the finite element results fail to converge no matter how small an element size is used. In this work, this pathology is addressed by assuming that the material parameters depends on the mesh size and the constitutive model is calibrated for finite elements of size 0.635 mm. Although not presented in this work, a parametric study can be carried out to established mesh dependence of critical model parameters [22]. Model parameters are obtained from finite element

182

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

700

a0 = a1 a2 = 0

Stress (Mpa)

600 500 400 300

a3 = 1.0

200 100

a3 = 15

0 0.0

2.0

a3 = 10

a3 = 5

4.0

6.0

8.0

Displacement (mm) Fig. 5. Parametric study case (c) – effect of parameter a3.

Fig. 6. Symmetric notch (SN) geometry.

Fig. 7. Asymmetric notch (ASN) geometry.

simulations of specimen SN-1, wherein multiple simulations are carried out to adjust the model parameters in order to match the experimental results. As the model parameters are assumed to represent the constitutive behavior under different stress states, the calibrated parameters should be able to model other stress states. Finite element models of the notched bar employ selectively reduced fully integrated eight node isoparametric solid elements [40]. Fig. 10 shows the finite element

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

183

Extensometer with gage length of 38.1 mm

Fig. 8. Experimental test set-up.

mesh in the notch region of SN-1 and ASN-2 specimens. It should be noted that the parameter calibration is carried out based on only SN-1 specimen without any consideration of the performance of the model for the other specimens. The void nucleation strain, en, is taken as a function of stress triaxiality. Fig. 11 shows the bilinear relationship that is used to model void nucleation strain as a function of stress triaxiality. This relationship is modeled after the results obtained by Bandstra et al. [43], where it was shown that the effective strain to failure first decreases rapidly with triaxiality and then at higher triaxiality the decrease is at a slower rate. Fracture is simulated by removing affected elements when the damage variable, ad, reaches a critical value of adc ¼ 0:90. Specimen SN-1 is used for calibration of model parameters: a0, a1, a2, a3, en and ec. Table 2 shows the material parameters, which give a good match as compared to the experimental results obtained for SN1. As shown in Fig. 12(a), with the above model parameters, a good match is obtained between experimental and finite element results for this specimen. Based on the parametric studies during the model calibration process, it is observed that the model parameters are fairly independent of each other. Indeed, only a0 and a1 parameters are inter-dependent as the first two terms in the damage in the damage evolution Eq. (14) contributes simultaneously to damage. The void coalescence parameters, a3 and a4, have little influence on a0 and a1 parameters. 5.2.3. Model validation Model validation is carried out by using the model parameters obtained in Section 5.2.2 to simulate the response for other test cases i.e. for specimens SN-2, ASN-1 and ASN-2. Fig. 12(b)–(d) shows the comparisons between finite element and experimental results. A good match between the finite element and experimental results is obtained for all the cases. These validation studies show that the model parameters can be considered as constitutive properties, and therefore, can be used to represent material response under a variety of stress states. The results from J2-plasticity model in which damage is not modeled are also included in Fig. 12. These results clearly show that for asymmetric notches (ASN-1 and ASN-2) there is no softening because of geometric nonlinearity and softening can be mainly attributed to material damage. However, for symmetric notches (SN-1 and SN-2) there is some softening due to geometric nonlinearity; but in this case also most of the softening is due to material damage. Plots of stress triaxiality at various locations in the specimens are shown in Fig. 13. The average triaxiality on the fracture plane (as shown in Figs. 14 and 15), together with the triaxiality value at the center and the edge of the fractured plane are plotted in these figures. It can be observed that stress triaxiality varies from an average value of about 1.5 for symmetric

1400

Stress (MPa)

1200 1000 800 600 400 200 0 0.0

0.5

1.0

1.5

2.0

Effective Plastic Strain Fig. 9. Hardening curve for ASTM A36 steel.

2.5

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

Fig. 10. Finite element mesh in the notch region: (a) SN-1 specimen and (b) ASN-2 specimen.

0.35 0.30

Nucleation Strain

184

0.25 0.20 0.15 0.10 0.05 0.00 0.0

1.0

2.0

3.0

4.0

Stress Triaxiality Fig. 11. Variation of nucleation strain with stress triaxiality.

Table 2 Model parameters for ASTM A36 steels for mesh size of 0.635 mm. Model parameters Bulk modulus Shear modulus Hardening-curve Yield stress, ry Shear parameter, a0 Triaxiality parameter, a1 Void coalescence parameter, a2 Void coalescence parameter, a3 Void coalescence strain, ec Void nucleation strain, en

166625.6 MPa 76904.1 MPa Fig. 9 289.5 MPa 0.175 0.325 1.75 2.5 0.5 Fig. 11

185

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

(b) 60

(a) 30

50

Load (KN)

25

Load (KN)

20 15 Experiment

10

40 30 Experiment 20 Damage Model

Damage Model 5

10

J2-Plasticity

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

J2-Plasticity

0 0.0

3.5

0.5

1.0

Displacement (mm)

20

40

Load (KN)

(d) 50

Load (KN)

(c) 25

15 Experiment

10

Damage Model 5

1.5

2.0

2.5

3.0

3.5

Displacement (mm)

30 Experiment

20

Damage Model 10

J2-Plasticity

J2-Plasticity

0

0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

0.0

5.0

0.5

1.0

Displacement (mm)

1.5

2.0

2.5

3.0

3.5

4.0

4.5

Displacement (mm)

Fig. 12. Comparison of damage model and experimental results: (a) SN-1 specimen, (b) SN-2 specimen, (c) ASN-1 specimen and (d) ASN-2 specimen.

(b) 4.0 3.5

3.0

Center

2.5

Stress Triaxiality

Stress Triaxiality

(a) 3.5 Average

2.0 1.5 1.0

Edge

0.5

Center

3.0

Average

2.5 2.0 1.5 1.0

Edge

0.5 0.0

0.0 0.0

1.0

2.0

3.0

0.0

1.0

Displacement (mm)

(c) 2.5

3.0

4.0

(d) 3.5

2.0

Stress Triaxiality

Stress Triaxiality

2.0

Displacement (mm)

Average

1.5

Center 1.0 0.5

Edge 0.0

3.0

Center

2.5 2.0

Average

1.5 1.0

Edge

0.5 0.0

0.0

1.0

2.0

3.0

Displacement (mm)

4.0

0.0

1.0

2.0

3.0

4.0

Displacement (mm)

Fig. 13. Variation of stress triaxiality: (a) SN-1 specimen, (b) SN-2 specimen, (c) ASN-1 specimen and (d) ASN-2 specimen.

5.0

186

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

Plane of Fracture

(a)

(b) Fig. 14. Fracture in symmetric notches: (a) undeformed specimen and (b) fracture mode.

Plane of Fracture

(a)

(b) Fig. 15. Fracture in asymmetric notches: (a) undeformed specimen and (b) fracture mode.

Location of Fracture initiation

Fracture initiation

(a)

Fracture initiation

(b) Fig. 16. Fracture initiation in symmetric notches: (a) SN-1 specimen and (b) SN-2 specimen.

187

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

Location of Fracture initiation

Location of Fracture initiation

Fracture initiation

(a)

Fracture initiation

(b) Fig. 17. Fracture initiation in asymmetric notches: (a) ASN-1 specimen and (b) ASN-2 specimen.

(a) 0.8 0.6

: Edge

0.5

: Edge

: Center

0.7

: Center

History Variables

History Variables

(b) 0.8

: Center

0.7

0.4 0.3 0.2 0.1

: Center

0.6

: Edge

0.5

: Edge

0.4 0.3 0.2 0.1

0.0 0.0

1.0

2.0

0.0 0.0

3.0

1.0

Displacement (mm) : Center

(d) 0.8

0.8

: Center

0.7

0.7

: Edge

0.6

: Edge

History Variables

(c) 0.9 History Variables

2.0

3.0

Displacement (mm)

0.5 0.4 0.3 0.2

: Center : Center

0.6

: Edge

0.5

: Edge

0.4 0.3 0.2 0.1

0.1 0.0 0.0

1.0

2.0

3.0

Displacement (mm)

4.0

0.0 0.0

1.0

2.0

3.0

4.0

Displacement (mm)

Fig. 18. Evolution of history variables at selected locations: (a) SN-1 specimen, (b) SN-2 specimen, (c) ASN-1 specimen and (d) ASN-2 specimen.

188

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

notches to an average of about 0.9 for asymmetric notches. A large variation of triaxiality within the section can also be observed. The locations of fracture initiation for various cases are shown in Figs. 16 and 17. Failure initiation occurs at the edge of specimens near the notch where triaxiality is high during the initial loading stage leading to smaller values of void nucleation strain in these regions. The location of fracture initiation is consistent with the results reported by [8,10,44], which showed fracture initiates at locations of highest triaxiality. Thus, the proposed damage evolution model is able to successfully represent the ductile fracture process under various states of stress. Finally, Fig. 18 shows the evolution of plastic and damage internal variables at selected locations – edge and center – of the specimens. Clearly, the evolution of damage variables matches the fracture locations shown in Figs. 16 and 17. 6. Concluding remarks The study presented in this paper was concerned with the development of a damage mechanics based constitutive model that can be efficiently used to simulate ductile fracture in structural steels. Specifically, a damage model based on the concept of effective stress and strain equivalence principle in a multiplicative plasticity framework is proposed. The damage evolution law is formulated to take into account the changes that arise due to micro-structural evolution during the ductile fracture process. In particular, damage due to the three stages of ductile fracture initiation – microvoid nucleation, growth and coalescence – are modeled. There are six tunable parameters – a0, a1, a2, a3, en and ec – in the proposed damage model as opposed to nine tunable model parameters in the GTN model. The parameter en controls damage initiation, a0 and a1 control damage under shear and high triaxiality, and finally the parameters a2, a3 and ec control failure due to rapid microvoid coalescence. Success of the proposed model is due to the proper tuning of these parameters; however, the calibration of these parameters is relatively easy as the parameters are independent of each other. In contrast, the GTN model parameters can be interrelated as shown in a recent study by Kiran and Khandelwal [25], and therefore the estimation of GTN model parameters is nontrivial. Furthermore, as only two nonlinear equations are solved as opposed to four equations in the GTN model, the proposed model is computationally efficient. A time stepping algorithm for numerical integration of the proposed damage model is presented and implemented in the explicit finite element code LS-DYNA. Parametric studies with varying model parameters showed that the proposed model can account for various stages of fracture. Experimental testing of four symmetrically and asymmetrically notched ASTM A36 steel bars is carried out, and results from the tests are used to calibrate and validate the proposed model. Finally, a set of model parameters are proposed for ASTM A36 steel and the simulation results from the finite element analysis showed that the proposed constitutive model is able to successfully predict the failure due to ductile fracture under varied stress-states. For the limited set of validation studies the model performed adequately; however, more experimental data and corresponding analyses are needed to cover the range of stress not included in the present study. It should also be noted that the presented constitutive model parameters are mesh dependent, and the model has to be recalibrated if a different material or mesh size is employed. Although the proposed model is used for simulation of fracture in ASTM A36 structural steels, the model can be calibrated for simulating ductile fracture in other structural steels, including ASTM A572 and ASTM A992 steels, which are also commonly used in the steel construction industry. Acknowledgements The presented work is supported in part by the US National Science Foundation through Grant CMS-0928547. Any opinions, findings, conclusions, and recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the sponsors. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

Wilsdorf HGF. Void initiation, growth, and coalescence in ductile fracture of metals. J Electron Mater 1975;4:791–809. Anderson TL. Fracture mechanics: fundamentals and applications. 3rd ed. Boca Raton, Florida: CRC Press, Taylor and Francis; 2005. Rice JR, Tracey DM. On ductile enlargement of voids in triaxial stress fields. J Mech Phys Solids 1969;17:201–17. Barsoum I, Faleskog J. Rupture mechanisms in combined tension and shear—micromechanics. Int J Solids Struct 2007;44:5481–98. Danas K, Ponte Castaneda P. Influence of the Lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials. Int J Solids Struct 2012;49:1325–42. Chao S-H, Khandelwal K, El-Tawil S. Ductile web fracture initiation in steel shear links. J Struct Engng 2006;132:1192–200. Kanvinde AM. Micromechanical simulation of earthquake induced fracture in steel structures, PhD thesis. Department of Civil & Environmental Engineering, Stanford University, Stanford, CA; 2004. Kiran R, Khandelwal K. A micromechanical model for ductile fracture prediction in ASTM A992 steels. Engng Fract Mech 2013;102:101–17. Kiran R, Khandelwal K. Experimental studies and models for ductile fracture in ASTM A992 steels at high triaxiality. J Struct Engng 2013. http:// dx.doi.org/10.1061/(ASCE)ST.1943-541X.0000828. Toribio J. A fracture criterion for high-strength steel notched bars. Engng Fract Mech 1997;57:391–404. Kiran R, Khandelwal K. Fast-to-compute weakly coupled ductile fracture model for structural steels. J Struct Engng 2014. http://dx.doi.org/10.1061/ (ASCE)ST.1943-541X.0001025. Gurson AL. Continuum theory of ductile rupture by void nucleation and growth. I. Yield criteria and flow rules for porous ductile media. Trans ASME Ser H, J Engng Mater Technol 1977;99:2–15. Lemaitre J. A continuous damage mechanics model for ductile fracture. Trans ASME J Engng Mater Technol 1985;107:83–9. Rousselier G. Ductile fracture models and their potential in local approach of fracture. Nucl Engng Des 1987;105:97–111.

K. Khandelwal, S. El-Tawil / Engineering Fracture Mechanics 116 (2014) 172–189

189

[15] Perrin G, Leblond JB. Analytical study of a hollow sphere made of plastic porous material and subjected to hydrostatic tension – application to some problems in ductile fracture of metals. Int J Plast 1990;6:677–99. [16] Gologanu M, Leblond JB, Devaux J. Approximate models for ductile metals containing nonspherical voids-case of axisymmetric oblate ellipsoidal cavities. Trans ASME J Engng Mater Technol 1994;116:290–7. [17] Gologanu M, Leblond J-B, Devaux J. Approximate models for ductile metals containing non-spherical voids. Case of axisymmetric prolate ellipsoidal cavities. J Mech Phys Solids 1993;41:1723–54. [18] Tvergaard V. On localization in ductile materials containing spherical voids. Int J Fract 1982;18:237–52. [19] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32:157–69. [20] Rakin M, Gubeljak N, Dobrojevi M, Sedmak A. Modelling of ductile fracture initiation in strength mismatched welded joint. Engng Fract Mech 2008;75:3499–510. [21] Xue Z, Pontin MG, Zok FW, Hutchinson JW. Calibration procedures for a computational model of ductile fracture. Engng Fract Mech 2010;77:492–509. [22] Mahnken R. Aspects on the finite-element implementation of the Gurson model including parameter identification. Int J Plast 1999;15:1111–37. [23] Nahshon K, Hutchinson JW. Modification of the Gurson model for shear failure. Eur J Mech A/Solids 2008;27:1–17. [24] Gao X, Tingting Z, Jun Z, Graham SM, Hayden M, Roe C. On stress-state dependent plasticity modeling: significance of the hydrostatic stress, the third invariant of stress deviator and the non-associated flow rule. Int J Plast 2011;27:217–31. [25] Kiran R, Khandelwal K. Gurson model parameters for ductile fracture simulation in ASTM A992 steels. Fatigue Fract Engng Mater Struct 2013. http:// dx.doi.org/10.1111/ffe.12097 [in press]. [26] Argon AS, Im J, Safoglu R. Cavity formation from inclusions in ductile fracture. MTA 1975;6:825–37. [27] Goods SH, Brown LM. Overview No. 1: the nucleation of cavities by plastic deformation. Acta Metall 1979;27:1–15. [28] Beremin FM. Cavity formation from inclusions in ductile fracture of A508 steel. MTA 1981;12:723–31. [29] Le Roy G, Embury JD, Edwards G, Ashby MF. A model of ductile fracture based on the nucleation and growth of voids. Acta Metall 1981;29:1509–22. [30] Thomason PF. Ductile fracture of metals. 1st ed. Oxford, UK: Pergamon Press; 1990. [31] Thomason PF. A view on ductile-fracture modeling. Fatigue Fract Engng Mater Struct 1998;21:1105–22. [32] Khandelwal K, El-Tawil S. Collapse behavior of steel special moment resisting frame connections. J Struct Engng 2007;133:646–55. [33] Kachanov L. Rupture time under creep conditions. Int J Fract 1999;97:11–8. [34] Lemaitre J. Evaluation of dissipation and damage in metals submitted to dynamic loading. In: Proceedings of ICM-1, Kyoto Japan; 1971. [35] Gurtin ME, Fried E, Anand L. The mechanics and thermodynamics of continua. Cambridge University Press; 2010. [36] Simo JC. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput Methods Appl Mech Engng 1992;99:61–112. [37] Steinmann P, Miehe C, Stein E. Comparison of different finite deformation inelastic damage models within multiplicative elastoplasticity for ductile materials. Comput Mech 1994;13:458–74. [38] Mahnken R. A comprehensive study of a multiplicative elastoplasticity model coupled to damage including parameter identification. Comput Struct 2000;74:179–200. [39] Coleman BD, Noll W. The thermodynamics of elastic materials with heat conduction and viscosity. Arch Ration Mech Anal 1963;13:167–78. [40] Hallquist J. LS-DYNA theory manual. Livermore Software Technology Corporation, 7374 Las Positas Road, Livermore, CA 94551; 2012. [41] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [42] Belytschko T, Moran B, Liu WK. Nonlinear finite element analysis for continua and structures. West Sussex, England: Wiley; 2000. [43] Bandstra JP, Goto DM, Koss DA. Ductile failure as a result of a void-sheet instability: experiment and computational modeling. Mater Sci Engng A 1998;249:46–54. [44] Toribio J, Ayaso FJ. Macro- and micro-scopic approach to fracture of high strength steel notched bars. In: Nishida S-I, editor. Macro and microscopic approach to fracture. Boston: WIT Press; 2004.