- Email: [email protected]
An extended GTN model for ductile fracture under high and low stress triaxiality
Contents lists available at ScienceDirect
International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas
An extended GTN model for ductile fracture under high and low stress triaxiality L. Malcher a,⇑, F.M. Andrade Pires b, J.M.A. César de Sá b a b
Department of Mechanical Engineering, Faculty of Technology, University of Brasília, Brasília, DF, Brazil IDMEC – Institute of Mechanical Engineering, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, Porto 4200-465, Portugal
a r t i c l e
i n f o
Article history: Received 8 April 2013 Received in final revised form 28 August 2013 Available online 12 September 2013 Keywords: Ductile fracture Volume void growth Shear mechanism Lode angle Stress triaxiality
a b s t r a c t This contribution provides an improvement on GTN model upon the prediction of fracture location within low level of stress triaxiality. In the proposition, two distinct damage parameters are introduced as internal variables of the degradation process and an effective damage is calculated as a sum of both contributions in the post-processed step. In the beginning, the volume void fraction, based on conservation mass law, is assumed as the first damage parameter, similar to Gurson’s original model. This volumetric damage contribution is able to capture spherical void growth, which plays the main role in tensile loading condition. The second damage parameter is proposed as a new shear mechanism, based on geometrical and phenomenological aspects and is also a function of the equivalent plastic strain, Lode angle and stress triaxiality. The shear damage parameter is formulated independent of the volume void fraction and requires a new nucleation of micro-defects mechanism to trigger the shear growth contribution, and hence is able to capture elongated (and rotation) void growth which is present in simple shear and combined shear/tensile or shear/compression loading conditions. Furthermore, the first and the second damage parameters are coupled in the yield function in order to affect the hydrostatic stress and deviatoric stress contributions, separately. In the first part of this paper, a review of Gurson’s model and its most famous version as GTN’s model is done. After that, the new contribution is presented and an implicit numerical integration algorithm is determined, based on the operator split methodology. The calibration strategy is discussed for determination of material parameters. Numerical tests are performed for a butterfly specimen using two types of materials (aluminum alloy 2024-T351 and steel 1045) under ranges of stress triaxiality between1=3 < g < 1=3 (shear/compression or shear/tensile). At the end, the behavior of internal variables is analyzed, such as: evolution of both damage parameters, evolution of the equivalent plastic strain, the reaction curve and the contour of the effective damage parameter. The results obtained are compared with experimental data and have shown that the present formulation performs well in the prediction of the fracture location and determination of the correct level of equivalent plastic strain at fracture under predominant shear loading condition. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction In many engineering applications, the characterization of material ductility can be done by a combination of equivalent plastic strain and stress triaxiality parameters (Bridgman, 1952; McClintock, 1968; Rice and Tracey, 1969; Johnson and Cook, ⇑ Corresponding author. Tel.: +55 6199944157. E-mail address: [email protected] (L. Malcher). 0749-6419/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijplas.2013.08.015
194
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
1985; Bao, 2003). McClintock (1968) and Rice and Tracey (1969) have also proposed an exponential expression, which is function of the equivalent plastic strain, in order to describe the behavior of ductile material in two dimensional spaces. Furthermore, the micro mechanical models for void growth together with their work have laid the foundation of modern ductile fracture mechanics (McClintock, 1968; Rice and Tracey, 1969; Hancock and Mackenzie, 1976). In the eighties, Johnson and Cook (1985) integrated the effects of the stress triaxiality, strain rate and temperature into a simple equation of a separable form. The large majority of engineering alloys contains several populations of inclusions corresponding to different length scales. Typically, it is possible to distinguish two main populations: one composed by primary inclusions, which are large particles embedded in the matrix, and one composed of secondary inclusions (or second phase particles), which can be 10–1000 orders of magnitude smaller. The phenomenon of ductile failure is usually induced by primary inclusions and second phase particles where micro-voids nucleate either by de-cohesion of inclusions (or second phase particles) from the surrounding matrix of by fracture of inclusions. The nucleated damage grows consistently with the applied stress state and the material degrades rapidly with the coalescence of multiple damage sites. Two ductile failure micro mechanisms have been identified in Fig. 1 (see Besson, 2010) and in both cases, failure by internal necking and failure by void sheeting, the coalescence mechanism is accelerated by second population particles. The ductile fracture phenomenon can be described, based on a micromechanical analysis of micro cavity growth, especially for the fracture computation within local approaches of fracture (see Pineau, 1981; Mudry, 1985; Rousselier’s, 1987; Besson et al., 2001) or based on the Continuum Damage Mechanics theory and a thermodynamic framework, either phenomenological or micromechanically based, as Lemaitre (1985) for damage caused by plastic flow, Chaboche (1984) and Murakami and Ohno (1981) for creep damage, Krajcˇinovic´ and Fonseka (1981) on micromechanical grounds. The formulations proposed by Lemaitre and Gurson are the most important coupled damage ductile models to describe the above two methodologies (see Chaboche et al., 2006). Since then, motivated by the limitations of these classical models, such as in prediction of the correct fracture location or in determination of the correct values of the internal variables at fracture, many researchers have proposed improvements in both methodologies, by introducing more effects in the constitutive formulation or in the damage evolution law like the pressure effect, temperature, Lode angle dependence, visco-plastic effects, crack closure effect, shear mechanisms, among others (Tvergaard and Needleman, 1984; Rousselier’s, 1980; Rousselier’s, 2001; Barsoum and Faleskog, 2007a; Barsoum and Faleskog, 2007b; Xue, 2007; Nahshon and Hutchinson, 2008; Nielsen and Tvergaard, 2008; Tvergaard, 2008; Lemaitre and Chaboche, 1990; Chaboche et al., 2006; Chaboche, 2003; Andrade Pires et al., 2003; Chaboche et al., 2006; Besson, 2010; Mirone et al., 2010; Li et al., 2011; Stoughton et al., 2011; Khan and Liu, 2012). In general, these classical coupled damage models have the ability to predict the correct fracture location under a specific range of combined stress triaxialities and Lode angle parameter (see Xue, 2007; Nahshon and Hutchinson, 2008; Teng and Wierzbicki, 2006; Malcher et al., 2012; Brunig et al, 2013) and are extremely accurate for loading conditions close to the calibration point (see Malcher et al., 2012). For example, within range of high levels of stress triaxialities, where the spherical void growth is the predominant mechanism, the models based on Gurson theory, like the Gurson–Tvergaard–Needleman model, have good performance in prediction of fracture location and parameters in fracture as equivalent plastic strain and displacement. However, under shear dominated loads, where failure is mainly driven by the shear localization of plastic strain of the inter-voids ligaments due to void rotation and distortion, the model does not perform well (see Engelen, 2005; Chaboche et al., 2006). Motivated by these short comings, in this contribution, a new extension to the GTN model is proposed in order to improve the ability to predict the correct fracture location and determinate the internal parameters at fracture. A new independent damage parameter is suggested to capture elongation of micro-defects and coupled to constitutive model to affect only the deviatoric stress part. A nucleation of general micro defects is introduced to trigger the shear mechanism and gives more accuracy to the model in prediction of ductile failure under mixed loading condition.
(a)
(b)
Fig. 1. Schematic representation of ductile failure micro mechanisms: (a) internal necking and (b) void sheeting (adapted from Besson (2010)).
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
195
2. Constitutive modelling 2.1. Gurson’s model The model proposed by Gurson (1977) is one of the first micromechanical based models for the description of ductile damage and fracture, which introduces a strong coupling between plastic strain and damage, in the presence of finite strains. The governing equations of the model were established by assuming a spherical cavity embedded in a cubic rigid-plastic matrix without hardening (see Fig. 2)a and use of the upper bound plasticity theorem. The degradation of the material is measured through the relation between the volume of void, V v oid , and the volume of a representative element, V RVE . This relation is denominated by the volume void fraction, f. Under tensile dominant loads, the macroscopic behaviour of the material can be represented by the reaction versus displacement curve depicted in Fig. 2b. In the elastic domain, the material is represented by stage (a), there is no appreciable change in the micro structure. Nevertheless, with the increase of the macroscopic load the nucleation of the micro voids is trigged due to existence of the localized plastic strain stage (b). In stage (c), the growth of micro voids is promoted by the high tensile hydrostatic stresses followed by coalescence of voids in stage (d). The evolution of the volume void fraction predicted by Gurson’s model follows as a direct consequence of the requirement for mass conservation of a rigid-plastic material assuming plastic incompressibility. Hence, the density of a representative volume element of a material with voids (Fig. 2a) can be determined by:
q ¼ qm
Vm : V RVE
ð1Þ
where q represents the density of the RVE; qm represents the density of the matrix of the material and V m is the volume of the matrix of the material. Thus, the relationship between the volume of the matrix of the material, V m , and the volume void fraction, f, can be established by:
Vm ¼ ð1 f Þ: V RVE
ð2Þ
Substituting Eq. (2) into Eq. (1), and by time differentiation of the result, the density rate of the representative volume element, q_ , can be expressed as a relation between the density rate of the material matrix, q_ m , and the volume void fraction rate, f_ , as follow:
q_ ¼ q_ m ð1 f Þ qm f_ :
ð3Þ
The matrix material is assumed to be plastically incompressible. In addition, the elastic volumetric strains are neglected by assumption. Therefore, the principle of mass conservation requires that q_ m ¼ 0. Thus, manipulating algebraically the above equations, the following expression can be obtained:
(a)
(b)
Fig. 2. (a) Representative volume element (RVE) with a spherical void; (b) schematic representation of the process of nucleation, growth, and coalescence of micro voids and the relationship with the macroscopic load (adapted from Gurson (1977) and Pineau and Pardoen (2003)).
196
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
q_ q_ f_ ¼ ¼ ð1 f Þ: qm
ð4Þ
q
The principle of mass conservation establishes that the volumetric strain rate is determined by:
q_ _ ¼ e ¼ e_ ev þ e_ pv ; q v
ð5Þ
where the elastic and plastic strain rate contributions are represented by e_ ev and e_ pv , respectively. In Gurson’s model the matrix material is assumed to be rigid-plastic, therefore, by disregarding the elastic contribution, Eq. (4) can be re-written as:
q_ f_ ¼ ¼ e_ pv ð1 f Þ:
ð6Þ
qm
The previous equation is the most significant contribution to the degradation of a porous material and expresses the evolution law for the volume void fraction. The original yield function derived by Gurson (1977) for a void-matrix aggregate is expressed by:
Uðr; k; f Þ ¼ J2 ðS Þ
1 3p 1 þ f 2 2f cosh r2y ; 3 2ry
ð7Þ
where J 2 represents the second invariant of the deviatoric stress tensor ðS Þ; p is the hydrostatic pressure, ry is the isotropic hardening rule, which can be defined as ry ¼ k r0 , where k represents the thermodynamical force associated to the isotropic hardening state variable and r0 is initial yield stress. According to the hypothesis of generalized normality, the plastic flow rule is given by
e_ p ¼ c_
@U 1 3p I; ¼ e_ pd þ e_ pv ¼ c_ S þ c_ f ry sinh 3 2ry @r
ð8Þ
where the plastic strain rate tensor, e_ p , involves two terms: the deviatoric, e_ pd , and volumetric plastic strains, e_ pv , and c_ represents the plastic multiplier. With the volumetric flow, e_ pv , constitutive equation (Eq. 8), it is possible to obtain the evolution law for the volume void fraction, f_ , after the substitution of e_ pv in Eq. (6):
_ y sinh f_ ¼ ð1 f Þe_ pv ¼ f f 2 cr
3p : 2ry
ð9Þ
2.2. Gurson–Tvergaard–Needleman (GTN)s model One of the shortcomings of the Gurson model is the fact that, whatever strain history the material might be subjected; no volume void fraction evolution will be predicted if the initial void ratio is zero. Therefore, in order to enhance the model, several mechanisms for damage nucleation have been proposed such that voids can nucleate depending on the strain history. One of the most well known nucleation laws was proposed by Chu and Needleman (1980) and later used by Tvergaard and Needleman (1984) in the GTN model. The damage evolution is represented by three simultaneous or successive mechanisms: nucleation, growth and coalescence of voids. The effective porosity, f , is determined by the following bilinear function:
f ¼
8
1 q1
fc
f < fc ðf fc Þ ; ðff fc Þ
f P fc
ð10Þ
where fc represents the critical volume void fraction and ff is the volume void fraction at fracture. The effective porosity, f , is obtained from both nucleation and growth mechanisms if the volume void fraction is less than the critical value, fc . The coalescence mechanism becomes active when the volume void fraction is higher than the critical value, fc . The volume void fraction rate, f_ , is given by the sum of the nucleation and growth mechanism as:
f_ ¼ f_ n þ f_ g :
ð11Þ
The nucleation mechanism can be driven either by plastic strain or hydrostatic pressure (see Tvergaard and Needleman, 1984). The definition of the nucleation mechanism based on the equivalent plastic strain is given by:
f_ n ¼
" 2 # fN 1 ep eN e_ p ; pffiffiffiffiffiffiffi exp 2 sN sN 2p
ð12Þ
where fN represents the volume fraction of all particles with potential for microvoid nucleation, eN and sN are the mean strain for void nucleation and its standard deviation. The variable ep represents the equivalent plastic strain and e_ p is the rate of the accumulated plastic strain. The evolution of the growth mechanism in the GTN model is given by the same expression as the
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
197
original Gurson model (see Eq. 6). The yield function of the GTN’s model, which assumes isotropic hardening and isotropic damage, is expressed by:
Uðr; k; f Þ ¼ J2 ðS Þ
1 q 3p 1 þ q3 f 2 2q1 f cosh 2 r2y ; 3 2ry
ð13Þ
where parameters q1 , q2 and q3 are introduced to bring the model predictions into closer agreement with full numerical analyses of a periodic array of voids. 2.3. Shear mechanism The original formulation of Gurson based models did not include shear effects, which excludes the possibility of predicting shear localization and fracture under conditions of low triaxiality, 0 6 g < 1=3. Under shear dominated loading conditions, the distortion of voids and inter-void linking promotes an effective increase in the material internal degradation and contributes to the material softening. Therefore, in order to improve Gurson based models predictive ability, under both zero and low levels of stress triaxialities, several researchers (Barsoum and Faleskog, 2007a; McVeigha et al., 2007; Xue, 2008; Nahshon and Hutchinson, 2008;Butcher et al., 2009; Lecarme et al., 2011;Stoughton et al., 2011) have suggested the introduction of shear effects. The formulation of shear mechanisms, which can be based on geometrical or phenomenological considerations, resulted in evolutions laws that include the influence of the third invariant of the deviatoric stress tensor, the plastic strain tensor and its rate. The shear damage mechanism proposed by Xue (2007) is based upon the solution of McClintock (1968) for the coalescence of holes in a shear band. Due to its geometrical and physical appeal, it will be revised here the shear damage mechanism proposed by Xue (2008) and also describe an extended version proposed by Butcher et al. (2009). The mechanism is based on geometrical considerations of a representative square cell, containing a circular void at the center, which is subjected to a simple shear strain (see Fig. 3). The length of the cell is equal to L and the radius of the central void is given by R. When the cell structure is loaded, the void rotates and elongates in the preferred direction. Due to the requirement of volume conservation of the cell structure, Xue (2008) assumes that the relative position of the void does not change with respect to the cell (see Fig. 3). As the shear strain increases, the distance between the free surface of the void and the boundary of the representative volume element decreases. Fig. 4 shows the cell structure in the initial configuration (a) and in the deformed configuration (a). The minimum distance between the free surface of the void and the boundary of the RVE, which is represented by the parameter a, can be expressed by the relation between the length of the cell and the radius of the void (see Fig. 3a):
a¼
L R: 2
ð14Þ
The application of a simple shear strain, c, leads to the appearance of a deformation angle, a, on the deformed configuration given by:
tan a ¼ c:
ð15Þ
The minimum distance at the deformed configuration can be related with the initial distance, a, and the deformation angle, a. In addition, it can also be related with the simple shear strain as:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 : 1 þ c2
a0 ¼ a cos a ¼ a
ð16Þ
(a)
(b)
Fig. 3. Void shear mechanism: (a) initial configuration; (b) deformed configuration (adapted from Xue (2007)).
198
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 4. Evolution of the Lode angle functions, g 0 , with regard to the normalized third invariant, n, proposed by Xue (2008) and by Nahshon and Hutchinson (2008).
An artificial strain can be defined (Xue, 2008), using the logarithmic definition of strain, which can be associated with the reduction of this minimum distance as:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a a
eart ¼ ln 0 ¼ ln 1 þ c2 :
ð17Þ
The fracture initiation in a shear band, according to McClintock (1968), can be defined by the boundary contact condition of the sheared void with the longitudinal direction of the shear band. For small volume void fraction, Xue (2008) expressed the failure macroscopic shear strain in the shear band as:
eshearband ¼
L : 2R
ð18Þ
Consequently, the damage associated with the shearing of the void, Dshear , is defined by the ratio of the artificial and the macroscopic shear strain in a shear band (Xue, 2008):
Dshear ¼
eart eshearband
¼
ln
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ c2 L 2R
:
ð19Þ
Xue (2008) performed a Taylor series expansion and simplified the expression of the artificial strain term to:
1 2
eart c2 :
ð20Þ
pffiffiffi The shear strain can be expressed as a function of the von Mises equivalent strain, c ¼ 3eeq . Therefore, for simple shear and for small volume void fractions, Eq. (19) can be approximated by: 1 2 c 3 2 ffiffiffiffiffiffiffiffi ¼ pffiffiffiffi f ð1=2Þ e2eq ; Dshear 1 p p p =f 2
ð21Þ
where f ¼ pR2 =L2 is the volume void fraction of the cell in a two dimensional problem. For the three dimensional case, with a spherical void of radius R at the center of the representative cell of length L, the void volume fraction is expressed by f ¼ 4pR3 =3L3 for the cell. A similar three dimensional relation can also be obtained:
Dshear ¼
ð1=3Þ 3 6 f ð1=3Þ e2eq : 2 p
ð22Þ
The evolution of shear damage can be represented in the rate form as:
D_ shear ¼ q4 f q5 eeq e_ eq ;
ð23Þ
where q4 and q5 are geometrical parameters that can be defined for two or three dimensional problems. For a two dimen ð1=3Þ sional problem, q4 ¼ p3ffiffipffi and q5 ¼ ð1=2Þ and for a three dimensional problem, q4 ¼ 32 p6 and q5 ¼ ð1=3Þ.
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
199
A modified shear damage expression was later derived by Butcher et al. (2009) that, contrary to Xue (2008), did not perform a Taylor series expansion of the artificial strain (Eq. 17) and expressed the failure strain with the logarithmic definition as:
Dshear ¼
eart eshearband
where the parameter
v ¼ Rx =Lx ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ c2 pffiffiffiffiffiffiffiffiffi ; ln 1=v
ln
¼
ð24Þ
v is the ligament size ratio defined for two or three dimensional problems, respectively, as:
4 k2 f p k1
12 ;
v ¼ Rx =Lx ¼
6 k2 f p k1
13 ð25Þ
:
The parameters k1 and k2 are the void and cell aspect ratios and defined as:
k1 ¼
Ry ; Rx
k2 ¼
Ly ; Lx
ð26Þ
where Ry and Rx represent the radii of the void in the direction y and x. The dimensions Ly and Lx are the length of the cell in the direction y and x. The ratios k1 and k2 are equal to one in Xue’s model (Xue, 2008), which implies that the rotation of the cell structure is proportional to elongation of the void. According to Butcher et al. (2009), these parameters can more generically be expressed as a function of the stress state and, as a result, the evolution of the ligament size ratio is related to the normal strains. Under the assumption of simple shear and small volume void fractions, the shear strain is related to the equivalent von pffiffiffi Mises strain as c ¼ 3eeq and the evolution of shear damage is given by:
D_ shear ¼
! 1 3eeq pffiffiffiffiffiffiffiffiffi e_ eq : ln 1=v 1 þ 3e2eq
ð27Þ
Butcher et al. (2009) have shown that the shear damage expression (Eq. 24) complies with McClintock criterion while Xue’s expression (Eq. 21) does not. In addition, it was emphasized that the simplifications proposed by Xue (2008) have a critical role on the shear damage criterion and evolution rule. 2.4. Lode angle function The shear damage evolutions, which were described for a simple shear loading condition in Section 2.3, need to be generalized for arbitrary stress states. This can be accomplished with the introduction of a Lode angle dependence function. The Lode angle, which is associated to the third invariant of the deviatoric stress tensor, is an essential parameter in the characterization of the effect of the stress state on ductile fracture (Kim et al., 2003; Kim et al., 2004; Bao and Wierzbicki, 2004; Gao et al., 2006; Barsoum and Faleskog, 2007a; Barsoum and Faleskog, 2007b; Bai and Wierzbicki, 2007; Mirone et al., 2010; Gao et al., 2009; Gao et al., 2011; Stoughton et al., 2011; Brunig et al., 2013). The Lode angle dependence function ranges between g 0 ¼0, for dominant tensile stress states, and g 0 ¼ 1, for shear dominant stress states. For intermediate values, ð0 < g 0 < 1Þ, there is a combined stress state and the function should define the relative magnitude of each stress condition. The Lode angle dependence function proposed by Xue (2008) is defined by a linear expression of the normalized Lode angle, as:
g 0 ¼ 1 h ;
ð28Þ
where g 0 represents the so-called Lode angle function and h is the normalized Lode angle, expressed by:
h ¼ 1 6jhj ¼ 1 2 arcosðnÞ;
p
p
ð29Þ
where h represents the Lode angle and n is the normalized third invariant, which is a ratio between the third invariant of the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=3 deviatoric stress tensor, r ¼ ½ð27=2Þ: det S , and the von Mises equivalent stress, q ¼ ð3=2ÞS : S :
n¼
3 r : q
ð30Þ
An alternative Lode angle dependence function as been proposed by Nahshon and Hutchinson (2008), which discriminates between uniaxial and biaxial tension and expresses a quadratic relation with the normalized third invariant:
g 0 ¼ 1 n2 :
ð31Þ
Expressions (28) and (31) can be used to activate the shear mechanisms, described in Section 2.3, whenever shear effects are present. Fig. 4 represents the shape of both functions with regard to the third invariant of the deviatoric stress tensor. The shear damage evolutions, expressed by Eq. (23) and (27), can be rephrased for arbitrary loading conditions as:
D_ shear ¼ g 0 q4 f q5 eeq e_ eq ;
ð32Þ
200
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
" D_ shear ¼ g 0
! # 1 3eeq pffiffiffiffiffiffiffiffiffi e_ eq : ln 1=v 1 þ 3e2eq
ð33Þ
3. Extended constitutive formulation Due to the limitation of Gurson based models, in the prediction of fracture onset under conditions of low stress triaxiality or capture Mode II and Mode III of crack initiation (see Socie and Marquis, 2000), several researchers (Barsoum and Faleskog, 2007a; McVeigha et al., 2007; Xue, 2008; Nahshon and Hutchinson, 2008) have proposed the introduction of shear effects (see Section 2.3) on the formulation. Although the results obtained with the modified GTN models (Xue, 2008;Nahshon and Hutchinson, 2008) have shown improvements in the prediction of damage, it has also been observed (Reis et al., 2010;Malcher et al., 2012), that both models have inherent limitations. In particular, the prediction of the location of fracture, the displacement to fracture and the equivalent plastic strain to fracture, for combined stress states, is not adequate. Therefore, in order to overcome these shortcomings, in this contribution, a new extended GTN model is proposed that incorporates a new nucleation law for second-phase particles, the yield surface is modified to include two distinct damage mechanisms (volumetric void growth and shear damage), a modified Lode angle dependence function is introduced and a new criterion for coalescence is proposed. It is important to observe that the proposed extension of GTN model is here based on the micromechanical description of shear mechanisms when there is elongation and rotation of the micro voids, as well as, empirical observation in order to introduce the effect of stress triaxiality and Lode angle dependence. 3.1. Nucleation mechanism The nucleation of voids associated with the GTN model, described by Eq. (12), was proposed by Chu and Needleman (1980). This specific form for the nucleation of primary voids is strain rate controlled and was introduced on a purely phenomenological basis. Nevertheless, as described in Section 2, engineering alloys loaded under shear conditions create localization bands due to the nucleation of secondary voids through a void sheeting mechanism. Although, secondary nucleation might be hard to detect in some materials, it will be introduced here a new independent nucleation mechanism that triggers void sheeting and localization. This nucleation mechanism is a result of second-phase particle deboning and cracking. Following the same approach of Chu and Needleman (1980), it will be considered a normal distribution of all secondphase particles with potential for nucleation as:
D_ n ¼
" 2 # DN 1 ep e0N e_ p pffiffiffiffiffiffiffi exp 2 s0N s0N 2p
ð34Þ
where DN represents the fraction of all second-phase particles with potential for nucleation, e0N and s0N are the mean strain for 0 0 second-phase nucleation and its standard deviation. This set of parameters DN ; eN ; sN needs to be calibrated under simple shear loading condition. The calibration of material parameters will be discussed latter in a specific section. The extended GTN model, proposed in this contribution, incorporates two independent nucleation mechanisms. The first one, which is the conventional nucleation mechanism of the GTN model (Eq. 12), triggers the evolution of the volume void fraction. The second, described by Eq. (34), triggers the evolution of the shear mechanism. The activation of these nucleation mechanisms under pure volumetric and shear conditions is relatively straightforward to establish. Nevertheless, under arbitrary stress states that may include combinations of tensile/shear or compressive/shear is not so easy to define. It is necessary to couple both mechanisms and also establish their relative magnitude. Here, it will be introduced the Lode angle function, g 0 , to combine both nucleation mechanisms. Therefore, Eq. (12) and Eq. (34) are re-defined as:
f_ n ¼ ð1 g 0 Þ
D_ n ¼ g 0
" 2 # fN 1 ep eN e_ p ; pffiffiffiffiffiffiffi exp 2 sN sN 2p
" 2 # DN 1 ep e0N e_ p : pffiffiffiffiffiffiffi exp 2 s0N s0N 2p
ð35Þ
ð36Þ
Under pure tensile loading conditions, the function g 0 is equal to zero and only primary nucleation of voids occurs (Eq. 35). For simple shear loading conditions, the function g 0 is equal to one and only secondary nucleation occurs (Eq. 36). For combined tensile/shear stress states, both mechanisms are active and the Lode angle function defines the relative importance of each component. Finally, if a combination of shear/compressive conditions is present there is no nucleation of primary voids and secondary nucleation takes place with the function g 0 defining the relative magnitude. 3.2. Incorporating shear effects As mentioned previously, several modified versions of the Gurson model, which include damage growth under low triaxiality straining for shear dominated stress states, have been proposed in the literature (Barsoum and Faleskog, 2007a;
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
201
McVeigha et al., 2007; Xue, 2008; Nahshon and Hutchinson, 2008). Nevertheless, the incorporation of shear effects on the Gurson model has been mainly accomplished through the introduction of additional terms on the evolution of the volume void fraction, f_ , as:
f_ ¼ f_ n þ f_ g þ f_ y :
ð37Þ
where the term f_ y does not represent a physical value of the porosity but ensures the detrimental effect of void distortion and inter-void linking, associated with low triaxiality, in the material. Therefore, the volume void fraction, f, in the modified versions of Gurson’s model, does not represent the plastic volume change of the material as in the original Gurson model. Alternatively, this scalar variable, f, measures the total accumulation of different types of damage in the material in an average sense. In contrast with this approach, in this work, it will be used two separate damage variables. The first one is the evolution of the volume void fraction employed in the GTN model, rewritten here with appropriate modifications, as:
f_ ¼ f_ n þ f_ g ¼ ð1 g 0 Þ
" 2 # fN 1 ep eN e_ p þ ð1 f Þe_ pv : pffiffiffiffiffiffiffi exp 2 SN SN 2 p
ð38Þ
The second variable is the evolution of damage due to shear effects, which is defined by an independent scalar variable, as:
D_ ¼ D_ n þ q6 D_ shear ¼ g 0
" 2 # DN 1 ep e0N e_ p þ q6 D_ shear ; pffiffiffiffiffiffiffi exp 2 s0N s0N 2p
ð39Þ
where D_ represents the evolution of the shear damage variable, D_ n represents its nucleation, which was introduced in Eq. (36), and D_ shear is the evolution of shear effects that can be defined based on geometrical considerations (see Eqs. 23 and 27) or phenomenological considerations (Nahshon and Hutchinson, 2008). The parameter q6 is a numerical constant, calibrated for each specific material, which defines the magnitude of the damage growth rate in shear. The extended GTN model proposed here has two scalar damage variables: a volumetric damage component characterized by the volume void fraction, f, and a deviatoric damage component described by shear damage, D. Each of these variables will be coupled with a specific component of the stress tensor: the hydrostatic pressure, p, will be related with the volume void fraction, f and the deviatoric component of the stress tensor, S, will be associated with the shear damage variable, D. The yield function of the model is therefore, defined by the following equation:
Uðr; k; f ; DÞ ¼
J 2 ðS Þ 1 q 3p 1 þ q3 f 2 2q1 f cosh 2 r2y : ð1 DÞ 3 2ry
ð40Þ
Following the principle of maximum dissipation, the yield function is taken as the dissipation potential of the model. Therefore, the evolution law for the plastic flow, assuming the hypothesis of generalized normality, is given by:
e_ p ¼ c_
@U S 1 q 3p I: ¼ e_ pd þ e_ pv ¼ c_ þ c_ q1 q2 f ry sinh 2 3 2ry @r ð1 DÞ
ð41Þ
where c_ represents the plastic multiplier, S represents the deviatoric component of the stress tensor and I is a second order identity tensor. From the previous equation, it is possible to observe that each component of the plastic strain rate tensor is affected by a different damage variable. In addition, only the deviatoric plastic strain rate, e_ pd , was altered, when compared with the Gurson model (Eq. 8), due to the introduction of a distinct shear damage variable. Remark 1. Due to the fact that we have proposed two distinct damage variables in the formulation of the extended GTN model, the evolution of the volumetric plastic strain, e_ pv , predicted by the model (see Eq. 41) will be different from previously proposed modifications of the Gurson model (Xue, 2008; Nahshon and Hutchinson, 2008). In these models, the volumetric plastic strain is coupled with an effective porosity (or damage variable), f, which includes both the volume void growth and shear damage. In addition, in these models, the deviatoric plastic strain rate, e_ pd , is not affected by the material degradation. The evolution law for the hardening variable, R, is determined by performing the derivative of the yield function with regard to the hardening force, k:
@U c_ 3q2 p 2 3q2 p þ ry 1 þ q3 f 2 2q1 f cosh ; R_ ¼ c_ ¼ q1 q2 fp sinh 2ry 3 2ry @k ð1 f DÞ
ð42Þ
where the term ð1 f DÞ is introduced to account for the softening effect on the material evolution law. The equivalent plastic strain rate, for the present model, can be determined from:
e_ p ¼
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ( 2 ) u2 2 p p S:S 1 3q2 p t e_ : e_ ¼ c_ þ ry q1 q2 f sinh 3 3 ð1 DÞ2 3 2ry
ð43Þ
202
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
It is important to mention that the extended GTN model, described by the previous set of equations, does not change the original model under stress states where shear effects are not present. The extension only modifies the predictions for stress states that include shear components, i.e. for problems with a Lode angle function, g 0 , different from zero. 3.3. Damage evolution The evolution law for the volume void fraction, f_ , has got two components (see Eq. 38): nucleation and growth. The most significant contribution to the evolution of spherical voids is the growth mechanism, f_ g , which depends on the evolution of the volumetric plastic strain, e_ pv , established on Eq. (41). Therefore, with the substitution of the rate of the volumetric plastic strain, e_ pv , on Eq. (38) it is obtained the evolution law for the volume void fraction, f_ , of the model:
f_ ¼ f_ n þ f_ g ¼ ð1 g 0 Þ
" 2 # fN 1 ep eN 1 3q2 p e_ p þ c_ f ð1 f Þq1 q2 ry sinh pffiffiffiffiffiffiffi exp I: 2 3 2r y SN SN 2 p
ð44Þ
Two possible evolution laws for the shear damage, based on geometrical considerations, were described by Eqs. (32) and (33). The evolution of damage in the material inevitably reduces the overall elastic properties. However, this effect is small when compared to the influence of damage on the plastic behavior. Therefore, the evolution of damage due to shear effects, employed in this work, will neglect the influence of damage on elasticity as is usually done in this type of model. The shear damage evolution laws are redefined as a function of both the accumulated plastic strain; ep , and the rate of the accumulated plastic strain, e_ p , instead of the equivalent strain, eeq , and equivalent strain rate, e_ eq :
D_ shear ¼ g 0 q4 Dq5 ep e_ p ; " D_ shear ¼ g 0
# 1 3ep e_ p ; pffiffiffiffiffiffiffiffiffi ln 1=v 1 þ 3ep2
ð45Þ
v ¼ q4 D
k2 k1
q5 :
ð46Þ
In addition, due to the introduction of two separate damage variables, the shear damage evolutions described by Eqs. (32) and (33), for the present model, do not depend on the current volume void fraction, f, but on the current value of the shear damage variable, D. With the substitution of the equivalent plastic strain rate, e_ p , established on Eq. (43), on the shear damage evolution law based on Xue’s work (Eq. 45), it can be obtained:
D_ shear
ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ( 2 ) u 2 S : S 1 3q p 2 : ¼ c_ g 0 q4 Dq5 ep t þ ry q1 q2 f sinh 3 ð1 DÞ2 3 2ry
ð47Þ
On the other hand, with the substitution of the equivalent plastic strain rate, e_ p , established on Eq. (43), on the shear damage evolution law based on Butcher’s work (Eq. 46), it can be obtained:
D_
shear
ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ( u 2 ) u2 1 3ep S : S 1 3q p 2 t : ¼ c_ g 0 pffiffiffiffiffiffiffiffiffi þ ry q1 q2 f sinh 3 ð1 DÞ2 3 2ry ln 1=v 1 þ 3ep2
ð48Þ
3.4. Modified lode angle dependence function The Lode angle dependence functions proposed by Xue (2008) and Nahshon and Hutchinson (2008), described in Section 2.4, were introduced to generalize a simple shear damage evolution for arbitrary stress states. They critically include the effect of the deviatoric stress tensor through the Lode angle. The relative magnitude of shear and volumetric effects under combined stress states is defined by these functions, whose behavior was presented in Fig. 4. Nevertheless, in this contribution, it is suggested the introduction of the stress triaxiality parameter, g, defined by
p 1=3trðrÞ g ¼ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; q
3=2S : S
ð49Þ
on the definition of the Lode angle dependence function, g 0 . This modification, which can be performed on either of the functions described in Section 2.4, has the following exponential form: 1
g 00 ðn; gÞ ¼ ½g 0 ðnÞjgjþk
ð50Þ
where g 00 is the modified function, g represents the stress triaxiality and k is a numerical constant that needs to be calibrated for each material. If either the Lode angle dependence function proposed Xue or Nahshon and Hutchinson, g 0 , is selected; the modified function is expressed by:
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
203
1 jgjþk 1 2 g 00 ðn; gÞ ¼ 1 h jgjþk ¼ 1 1 acosðnÞ ;
ð51Þ
1 g 00 ðn; gÞ ¼ 1 n2 jgjþk :
ð52Þ
p
Fig. 5 illustrates the behavior of the original functions on the space of fg 0 ; n; gg and Fig. 6 represents the behavior of the modified function on the same space. The influence of stress triaxiality will be dominant whenever the range of stress triaxiality is within the range of ½1=3; 1=3. In order to analyze the influence of the new exponential term on the Lode angle dependence function, let us restrict ourselves, for instance, to the function described by Eq. (51) that was proposed by Xue (2008). The constant k introduced in the exponent can be better understood through the analysis of Fig. 7a and b. Fig. 7a illustrates the shape of the modified Lode angle function, g 00 , when the stress triaxiality is zero, g ¼ 0, and for different values to the constantk. If the constant is equal to one, k ¼ 1, the value of the exponent will be equal to 1 (since the stress triaxiality is zero) and it is obtained the same evolution for the modified Lode angle function as the original one (see Fig. 5). For values of constant, k, higher than unity, the modified Lode angle function will be higher than the original function. If the value of the constant, k, is smaller than 1 but higher than zero, the modified Lode angle function will be lower than the original function. Fig. 7b illustrates the shape of the modified Lode angle function, g 00 , when the normalized third invariant is zero, n ¼ 0, and for different values to the constant k. If the constant, k, is equal to zero, it is possible to appreciate the influence of the stress triaxiality, g, on the modified Lode angle function. It can be observed that for high values of the parameter, k, the modified Lode angle function, g 00 , is less affected by both third invariant and stress triaxiality. Nevertheless, for low values of the constant, k, the impact of the value of the parameters is very significant. Therefore, it is recommended that the range of this value should be defined within the interval ½0:1; 1=3. Fig. 8 also presents the influence of the parameter k on the global shape of the function g 00 . In particular, Fig. 8 (a) represents the original Lode angle function (without the introduction of the dependence of stress triaxiality) proposed by Xue (2008) and the set of Figs. 8 (b)–(c)–(d) represent the global shape of the modified function for different values of the parameter k. When the value of k is low, it can be observed a stronger dependence on the stress triaxiality and when k is high, the modified function is less dependent on stress triaxiality, recovering in the limit the original Lode angle function. 3.5. Coalescence criterion The definition of a criterion for void coalescence is extremely important for the prediction of fracture onset and for the simulation of crack formation and propagation. The simplest criterion for void coalescence is to assume a constant critical value of the volume void fraction, fc . Once this value is reached, the mechanism of coalescence accelerates the rate of increase of the volume void fraction, f_ , which will lead to final failure. The void coalescence process can be simulated by the function, f , introduced by Tvergaard and Needleman (1984) described by Eq. (10). This criterion has been followed by several authors (Needleman and Tvergaard, 1987; Xia et al., 1995; Faleskog et al., 1998). Nevertheless, further research on the topic focused on whether the volume void fraction could be regarded as a constant for different loading conditions. Several researchers (Berzerga et al., 1999; Zhang et al., 2000; Pardoen et al., 2000; Kim et al., 2004) employing unit cell models have shown that the volume void fraction of voids, fc , is not sufficient to describe the ini-
Fig. 5. Three dimensional representation of the Lode angle functions, g 0 : (a) Xue’s model (b) Nahshon and Hutchinson’s model.
204
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 6. Three dimensional representation of the modified Lode angle dependence function, g 00 : (a) Xue’s model (b) Nahshon and Hutchinson’s model.
Fig. 7. Influence of the constant k on the behavior of the function g 00 (a) with regard to the third invariant fixing the stress triaxiality to zero (g ¼ 0), and (b) with regard to stress triaxiality fixing the third invanriant to zero (n ¼ 0).
tiation of fracture and established that it depends on several factors, such as the void shape and spacing, stress triaxiality, material hardening, etc. Recently, criterions based on either the effective strain or effective stress have been proposed to trigger the coalescence of voids. In particular, Gao et al. (2009) employed the Gologanu et al. (1995), with the modification proposed by Pardoen et al. (2000), to describe void growth and the macroscopic plastic response of cell elements containing non-spherical micro voids. They concluded that when the macroscopic effective strain of the element reaches a critical value, void coalescence occurs. On the other hand, Jackiewicz (2011) has proposed a coalescence criterion based on the assumption that a singular value of the effective stress triggers the coalescence of micro voids in materials. In the present model, it has used two separate scalar damage variables, namely the volume void fraction, f, and the shear damage variable, D. Therefore, it will have two distinct critical values: the critical volume void fraction, fc , which is the critical volume void fraction employed in the GTN original model, and the critical shear damage value, Dc , which is regarded as a material constant that needs to be obtained. In this case, each critical values will have to be determined under different conditions: the critical volume void fraction, fc , will be obtained from a specimen subjected to tensile dominant loading conditions (associated with high triaxiality) and the critical shear damage value, Dc , will be obtained from a specimen subjected to a simple shear loading condition (associated with low triaxiality). The coalescence criterion proposed here introduces an effective damage variable, Def , which is conveniently normalized, to combine both critical damage parameters (fc and Dc Þ. The determination of fracture onset is established whenever the effective damage variable, Def , reaches unity.
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
205
Fig. 8. Behavior of the modified Lode angle dependence function g 00 for different values of k: (a) without dependence, (b) k ¼ 0, (c) k ¼ 0:1, (d) k ¼ 0:4.
8D ; g 00 ¼ 1 > > > Dc < 1 þ Dfcc ffc þ DDc ; 0 < g 00 < 1 Def ¼ > > > :f ; g 00 ¼ 0 fc
ð53Þ
The coalescence criterion based on the effective damage, Def , has three possible cases: – Under generalized tension (g 00 ¼ 0Þ, where the spherical void growth drives the damage evolution, the effective damage is given by the ratio f =fc and fracture onset is predicted when the volume void fraction, f, reaches the critical volume fraction, fc ; – Under generalized shear (g 00 ¼ 1Þ, where the distortion of voids and inter-void linking promotes damage evolution, the effective damage variable is given by the ratio D=Dc and fracture is predicted when the shear damage variable, D, attains the critical shear damage, Dc ; – Under combined stress states, both the void growth under hydrostatic tension and shear localization compete with each other and the prediction of fracture onset is a combination of both contributions ðf =fc Þ and ðD=Dc Þ. In addition, an extra term is introduced, in the definition of the effective damage variable, in order to account the multiaxiality of the stress state (combined stress state), which can be a hard factor for the degradation of the materials for many researchers (see Dowling, 1999; Socie and Marquis, 2000). This term will play the role of an accelerator parameter in the prediction of fracture onset, and here is defined as ð1 þ fc =Dc Þ.
The effective damage variable is evaluated in the model in a post processed manner and does not affect the evolution of the two independent scalar damage variables. Critically, the criterion proposed allows fracture initiation at different values of the volume void fraction, f, and at different values of the shear damage variable, D. Furthermore, in the case of
206
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
axsymmetric uniaxial tension, this criterion (employed in conjunction with the extended GTN model proposed) recovers the original GTN model criterion with critical fc as a particular case, where g 00 ¼ 0 and the damage shear evolution is zero. In Box 1, the basic constitutive equations and evolution laws for the internal variables and damage are summarized: Box 1. GTN’s extended modified model including nucleation, growth and shear effects. (i) Elasto-plastic split of the strain tensor e ¼ ee þ ep (ii) Elastic law r ¼ D e : ee (iii) Yield function h i J2 2p Uðr; r; f ; DÞ ¼ 1D 13 1 þ q3 f 2 2q1 f cosh 3q r2y 2 ry (iv) Plastic flow and evolution equations for R, f and D e_ p ¼ c_ N R_ ¼ c_ @@kU p 2 fNffiffiffiffiffi e eN 1 e_ p þ ð1 f Þe_ pv f_ ¼ ð1 g 0 Þ S p exp 2 SN N 2p p 0 2 1 _ ¼ g DpNffiffiffiffiffi exp 1 e 0 eN e_ p þ ½g jgjþk q D_ shear D 0 s0 2p N
2
sN
0
6
where,
8 < q3 Dq4 ep :e_ p ; if Xue’s shear mechanism is choosen p _ shear ¼ D 3e _ p if Butcher’s shear mechanism is choosen : ln p1ffiffiffiffiffiffi p2 e ; 1þ3 e 1=v 1 h ; if Xue’s Lodeangle dependence function is choosen g0 ¼ 1 n2 ; if Nahshon’s Lode angle dependence function is choosen and, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e_ p ¼ 23 ðe_ p : e_ p Þe_ pv ¼ tr ðe_ p Þ (v) Loading/unloading criterion
c_ P 0, U 6 0, c_ U ¼ 0
4. Numerical integration algorithm In this section, the numerical solution strategy adopted to perform the numerical simulations is summarized. Algorithms based on operator split methodology are especially suitable for the numerical integration of the evolution problem and have been widely used in computational plasticity (see Simo and Hughes, 1998; De Souza Neto et al., 2008). This method, which is employed in the present development, consists of splitting the problem in two parts: an elastic predictor, where the problem is assumed to be elastic and, a plastic corrector, in which the system of residual equations comprising the elasticity law, plastic consistency and the rate equations is solved, taking the results of the elastic predictor stage as initial conditions. In the case of the yield condition has been violated, the plastic corrector stage is initiated and the Newton- Raphson procedure is used to solve the discretised system of equations. The Newton–Raphson procedure was chosen motivated by the quadratic rates of convergence achieved which results in return mapping procedures computationally efficient (see Simo and Hughes, 1998; De Souza Neto et al., 2008). The model was implemented in a quasi-static finite element framework based on the infinitesimal strain theory. The extension of the model to the finite strain range was done by adopting the well-established multiplicative hyperelasto-plastic framework (Peric0 et al., 1992; Eterovic and Bathe, 1990). Let us consider what happens to a typical Gauss point of the finite element mesh within pseudo-time interval [t n ; t nþ1 ]. Given the incremental strain, De, and the values of rn , epn , epn , Rn , fn and Dn at time t n , the numerical integration algorithm should obtain the updated values at the end of the interval, rnþ1 , epnþ1 , epnþ1 , Rnþ1 , fnþ1 and Dnþ1 , in a manner consistent with the constitutive equations of the model. 4.1. Elastic predictor The first step in the algorithm is the evaluation of the elastic trial state, where the increment is assumed purely elastic with no evolution of internal variables. The elastic trial strain and trial state variables are given by: trial trial eenþ1 ¼ een þ De; epnþ1 ¼ epn ; Rtrial nþ1 ¼ Rn
ð54Þ
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
207
Dtrial nþ1 ¼ Dn :
trial fnþ1 ¼ fn ;
The corresponding elastic trial stress tensor is computed e e trial rtrial nþ1 ¼ D : enþ1
ð55Þ
e
where D is the standard isotropic elasticity tensor. Equivalently, in terms of stress deviator and hydrostatic pressure, it can be written as: e trial Strial nþ1 ¼ 2Gednþ1 ;
e trial ptrial nþ1 ¼ K ev nþ1
ð56Þ
trial trial where eednþ1 is the deviatoric elastic trial strain tensor and eev nþ1 is the volumetric elastic trial strain. The material constants G and K represent the shear and bulk moduli. The trial yield stress is defined, in this case, as a function of the internal hardening variable at time t n :
rtrial ¼ ry ðRn Þ ¼ r0 þ kn ¼ r0 þ HRn y
ð57Þ
where kn is the thermodynamical force associated with the isotropic hardening internal variable, Rn . P The next step of the algorithm is to check whether rtrial nþ1 lies inside or outside of the trial yield surface. With variables R, e , f and D frozen at time tn , it is computed: trial
U
" !# 2 3q2 ptrial 1 nþ1 trial2 trial 1 þ q3 fnþ1 2q1 fnþ1 cosh ¼ rtrial y trial trial 3 2 r 1 Dnþ1 y J trial 2
ð58Þ
If Utrial 6 0, the process is indeed elastic within the interval and the elastic trial state coincides with the updated state at t nþ1 . In other words, there is no plastic flow evolution within the interval and the trial state is equal to real state,
ðÞnþ1 ¼ ðÞtrial nþ1
ð59Þ
trial
Otherwise, if U > 0, it is necessary to apply the plastic corrector or return mapping algorithm whose derivation is described in the following. 4.2. Return mapping algorithm Following a straightforward specialization of standard return mapping procedure for the present constitutive equations, leads to the numerical integration of the evolution equations for een , epn , Rn , fn and Dn having the trial state as the initial condition. The discretisation of the elastic strain tensor reads:
trial eenþ1 ¼ eenþ1 Dc
Snþ1 1 3q2 pnþ1 I : þ ry q1 q2 fnþ1 sinh 1 Dnþ1 3 2ry
ð60Þ
With the application of Hooke’s Law to the above expression, it is possible to determine the evolution of the stress tensor:
rnþ1 ¼ De : eenþ1 2GDc
Snþ1 1 3q2 pnþ1 I: DcK ry q1 q2 fnþ1 sinh 1 Dnþ1 3 2ry
ð61Þ
Eq. (61) can be split into a deviatoric and a hydrostatic contribution. The updating relation for each component is given by:
Snþ1 ¼ h
1þ
Strial nþ1
2G:Dc 1Dnþ1
i ;
ð62Þ
1 3q2 pnþ1 I pnþ1 I ¼ ptrial D c K r q q f sinh y nþ1 1 2 nþ1 3 2ry
ð63Þ
Furthermore, the discrete counterparts of the other variables of the problem read:
Rnþ1 ¼ Rtrial nþ1 þ
Dc 3q2 pnþ1 2 3q2 pnþ1 2 þ ry 1 þ q3 fnþ1 q1 q2 fnþ1 pnþ1 sinh 2q1 fnþ1 cosh 3 ð1 fnþ1 Dnþ1 Þ 2ry 2ry
trial fnþ1 ¼ fnþ1 þ 1 g 0nþ1
Dnþ1 ¼
Dtrial nþ1
" p 2 # fN 1 enþ1 eN 3q2 pnþ1 pffiffiffiffiffiffiffi exp Dep þ ð1 fnþ1 ÞDcry q1 q2 fnþ1 sinh 2 SN 2ry SN 2 p
" p 2 # DN 1 enþ1 e0N þ g 0nþ1 0 pffiffiffiffiffiffiffi exp Dep þ q6 DDshear 2 S0N SN 2 p
ð64Þ
208
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
where the incremental shear damage component, DDshear , can be either defined by Eq. (47) or (48). Due to the fact that the Lode angle dependence function proposed by Nahshon and Hutchinson (2008) is continuous and does not have singular points (see Fig. 4), it is more convenient for numerical implementation. Therefore, in present’s derivations, the Nahshon and Hutchinson function is adopted. Furthermore, the updating relation for the equivalent plastic strain can be obtained from:
epnþ1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 u 8 u > 2 > trial trial < = u Snþ1 : Snþ1 2 1 3q2 pnþ1 trial ¼ epnþ1 þ Dcu þ ry q1 q2 fnþ1 sinh : i2 t3 >h > 3 2r y : 1 þ 2G:Dc ; ð1 Dnþ1 Þ2
ð65Þ
1Dnþ1
The above equations must be complemented by the so-called consistency condition that guarantees that the stress state at the end of the plastic step lies on the updated yield surface:
Unþ1 ¼ h
1þ
J trial 2nþ1 i2
2GDc 1Dnþ1
ð1 Dnþ1 Þ
1 3q2 pnþ1 2 1 þ q3 fnþ1 2q1 fnþ1 cosh r2y : 2ry 3
ð66Þ
Since it is possible to express the deviatoric stress tensor, Snþ1 , as a function of the plastic multiplier, Dc, the shear damage, Dnþ1 , and the trial stress state, Strial nþ1 , (Eq. 62) it is possible to eliminate the deviatoric stress tensor from the initial system of equations, which is composed by ten equations in the three-dimensional case. The return mapping scheme can therefore be reduced to a set of only five coupled non-linear equations, which need to be solved for the unknowns pnþ1 , Rnþ1 , fnþ1 , Dnþ1 and Dc. After the solution of the system, all other variables need to be continently updated. The overall algorithm for numerical integration is summarized in Box 2. Box 2. Fully implicit Elastic predictor/Return mapping algorithm. (i) Evaluate the elastic trial state: Given the incremental strain De and the state variables at t n : trial trial eenþ1 ¼ een þ De; epnþ1 ¼ epn ; Rtrial nþ1 ¼ Rn trial fnþ1 ¼ fn ;
e trial Strial nþ1 ¼ 2Genþ1 ¼ ry Rtrial nþ1
Dtrial nþ1 ¼ Dn ;
e trial ptrial nþ1 ¼ K ev nþ1 ;
rtrial y
(ii) Check plastic admissibility: 2 3q ptrial J trial trial 2 trial nþ1 2q1 fnþ1 cosh 2r2 trial rtrial 6 0 THEN IF Utrial ¼ 2 trial 13 1 þ q3 fnþ1 y 1Dnþ1
y
ðÞtrial nþ1
(elastic step) and go to (v) Set ðÞnþ1 ¼ ELSE go to (iii) (iii) Return mapping (plastic step): Solve the system of equations for Dc, pnþ1 , fnþ1 , Rnþ1 and Dnþ1 8 h i 9 Jtrial > 2nþ1 3q2 pnþ1 1 2 > > > > 1 þ q f 2q f cosh r2y > h i 2 3 nþ1 1 nþ1 8 9 > > 3 2 ry > > 2GDc > > 0> ð1Dnþ1 Þ 1þ 1D > > > > > > > nþ1 > > > > > > > 0> = = < < 3q p 2 nþ1 pnþ1 ptrial þ D c K r q q f sinh y 1 2 nþ1 nþ1 2ry ¼ 0 > > > > trial > > > Df n Df g 0> > > > fnþ1 fnþ1 > > > > > ; > : > > > trial > > 0 > > Rnþ1 Rnþ1 DR > > > > ; : trial n shear Dnþ1 Dnþ1 DD q6 DD where, Df n ¼ 1 g 0nþ1 S
ep e 2 N 12 nþ1SN Dep 3q p Df g ¼ ð1 fnþ1 ÞDcry q1 q2 fnþ1 sinh 22rnþ1 y n h io 3q2 pnþ1 3q p Dc 2 þ 23 ry 1 þ q3 fnþ1 DR ¼ ð1fnþ1 Dnþ1 Þ q1 q2 fnþ1 pnþ1 sinh 2ry 2q1 fnþ1 cosh 22rnþ1 y fNffiffiffiffiffi p exp N 2p
ep
e0
ffi exp½ 12 ð nþ10 N Þ2 Dep DDn ¼ g 0nþ1 S0 DpNffiffiffiffiffi SN 2:p 8N ðjg 1 jþkÞ q5 > g nþ1 epnþ1 Dep ; < q4 Dnþ1 0nþ1 p shear DD ¼ ðjg 1 jþkÞ 1 3enþ1 > :Dep ; : g 0nþ1 nþ1 ln pffiffiffiffiffiffi 1=v 1þ3 epnþ1 2 g 0nþ1 ¼ 1 n2nþ1 8 1 > < 4 Dnþ1 k2 2 ; if 2D problem p k1 v¼ 1 > k2 3 : 6D p nþ1 k1 ; if 3D problem
if Xue’s mechanism is chosen if Butcher’s mechanism is chosen
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
(iv) Update the others state variables: 2 e nþ1
e
Dc4h
e trial nþ1
¼e
1þ
Strial nþ1 2GDc 1Dnþ1
i
ð1Dnþ1 Þ
1 3
þ ry q1 q2 fnþ1 sinh
3q2 pnþ1 2ry
209
3 I5
trial
S Snþ1 ¼ h nþ1 i 1þ
2GDc 1Dnþ1
rnþ1 ¼ Snþ1 þ pnþ1 I
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9ffi u 8 u > = h i2 > u < trial trial Snþ1 :Snþ1 3q p 2 Dep ¼ Dcu þ 13 ry q1 q2 fnþ1 sinh 22rnþ1 t3 >h 2GDc i2 y > 2 : 1þ 1D ; ð1Dnþ1 Þ nþ1
pnþ1
e
trial pnþ1
¼e
p
þ De
(v) Exit
The previous set of discrete equations needs to be solved for the unknowns pnþ1 , Rnþ1 , fnþ1 , Dnþ1 and Dc. The Newton– Raphson (NR) method will be used for solving the return mapping system of equations due to the asymptotic rate of quadratic convergence of the method. Let us rewrite the non-linear scalar system of residual equations, in the following form:
8 h i 9 J trial 3q p 2 2nþ1 >h 2q1 fnþ1 cosh 22rnþ1 r2y > 13 1 þ q3 fnþ1 i2 > 8 9 > > > y > > 2GDc > > r 1þ 1D 1Dnþ1 Þ > > > > ð D c > > > > nþ1 > > > > > > > > > > > > r > > > > p 3q p trial < =
> > > trial > > Df n Df g fnþ1 fnþ1 > > rR > > > > > > > > > > > > > > > > : ; > > > Rnþ1 Rtrial > > > rD nþ1 DR > > > > : ; trial n shear Dnþ1 Dnþ1 DD q6 DD
ð67Þ
To obtain a new guess for each variable of the problem, it is necessary to solve a linearized system of equations given by:
2 @r
Dc
6 @ Dc 6 @rp 6 @ Dc 6 6 @rf 6 @ Dc 6 6 @rR 6 @ Dc 4 @r D @ Dc
@r Dc @p
@r Dc @f
@r Dc @R
@r Dc @D
@r p @p
@r p @f
@r p @R
@r p @D
@r f @p
@r f @f
@r f @R
@r f @D
@r R @p
@r R @f
@r R @R
@r R @D
@r D @p
@r D @f
@r D @R
@r D @D
3j
3jþ1 3j 2 2 r Dc 7 dDc 7 6 7 7 6 7 6 dp 7 6 rp 7 7 6 7 7 6 7 6 7 6 7 :6 df 7 7 ¼ 6 r f 7 7 6 7 6 7 4 dR 7 5 4 rR 5 7 5 dD rD
ð68Þ
In the above system, the terms in the matrix are composed by the derivative of each residual equation (see Eq. 67) with regard to each variable of the problem (pnþ1 , Rnþ1 , fnþ1 , Dnþ1 and DcÞ at iteration j. The matrix is multiplied by a vector with the incremental values of each variable at iteration, ðj þ 1Þ. The vector on the right hand side represents the residual of each variable at iteration j. After solving the system for the unknowns and obtaining a new guess for each variable, the convergence needs to be checked. More details of the linearization procedure for the present model can be found in appendix. The overall algorithm for numerical integration is summarized in Box 3 in pseudo-code format. 4.3. Consistent tangent operator Under elastic loading condition, the tangent operator for this constitutive formulation is the standard linear elasticity tensor. Nevertheless, for the plastic step, the elasto-plastic tangent operator is obtained by the linearization procedure of the above system of residual equations. Hence, the first step for determining the operator is to differentiate the stress tensor update expression:
rnþ1 ¼ h
1þ
Strial nþ1
2GDc 1Dnþ1
i þ pnþ1 I
After some algebraic manipulation, the above equation can be rewritten by differentiate form as:
ð69Þ
210
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
drnþ1 ¼ h
2G
2GDc 1Dnþ1
1þ
trial i deednþ1
8 <
h : 1þ
8 92 < = 1 2G Dc trial trial i eednþ1 i eednþ1 dDc h dDnþ1 þ dpnþ1 I : 1 þ 2GDc ; ; 1 Dnþ1 ð1 Dnþ1 Þ2 92 =
2G
2GDc 1Dnþ1
1Dnþ1
ð70Þ The terms dDc, dDnþ1 and dpnþ1 can be obtained by the linearization procedure of the residual system of equations. Furthermore, the elasto-plastic operator can be determined as:
Dep ¼
drnþ1 trial deenþ1
ð71Þ
Box 3. The Newton–Raphson algorithm for solution of the return mapping system of equations. ð0Þ
ð0Þ
ð0Þ
ð0Þ
ð0Þ
ð0Þ
ð0Þ (1) Initialize iteration counter, k : ¼ 0, set initial guess for Dcð0Þ ¼ 0; pnþ1 ¼ pn ; fnþ1 ¼ fn ; Rnþ1 ¼ Rð0Þ n and Dnþ1 ¼ Dn corresponding residual:
8 h i 9 J trial > > 2nþ1 3q p 2 >h > 13 1 þ q3 fnþ1 2q1 fnþ1 cosh 22 rnþ1 r2y > i 2 9 > 8 > > y > > 2GDc > > r Dc > > > > 1þ 1D 1Dnþ1 Þ ð > > > > nþ1 > > > > > > > > > > > > r > > > > p 3q p = =
> > > trial n g > > > > fnþ1 Df Df f > > > > rR > > > > > > > nþ1 > > > > > ; > > : trial > > R R D R > > nþ1 rD nþ1 > > > > > > ; :D Dtrial DDn q DDshear nþ1
6
nþ1
(2) Perform Newton–Raphson iteration 2 @r
Dc
6 @ Dc 6 @rp 6 @ Dc 6 6 @rf 6 @ Dc 6 6 @rR 6 @ Dc 4 @r D @ Dc
@r Dc @p
@rDc @f
@r Dc @R
@r Dc @D
@r p @p
@rp @f
@r p @R
@r p @D
@r f @p
@rf @f
@r f @R
@r f @D
@r R @p
@rR @f
@r R @R
@r R @D
@r D @p
@rD @f
@r D @R
@r D @D
3j
2 2 3jþ1 3j r Dc 7 dD c 7 6 6 7 7 7 6 dp 7 6 rp 7 7 6 6 7 7 7 6 6 7 7 :6 df 7 ¼ 6 rf 7 7 7 6 6 7 7 7 4 dR 5 4 rR 5 7 5 rD dD
New guess for pnþ1 ; Rnþ1 ; fnþ1 ; Dnþ1 and Dc: ðjþ1Þ
Dcðjþ1Þ ¼ DcðjÞ þ dDcðjþ1Þ , ðjþ1Þ
ðjÞ
ðjþ1Þ fnþ1
ðjÞ fnþ1
ðjþ1Þ
Rnþ1 ¼ Rnþ1 þ dRnþ1 , ¼
þ
ðjÞ
ðjþ1Þ
pnþ1 ¼ pnþ1 þ dpnþ1 ðjþ1Þ
ðjÞ
ðjþ1Þ
Dnþ1 ¼ Dnþ1 þ dDnþ1
ðjþ1Þ dfnþ1
(3) Check for convergence
~ ¼h U
1þ
Jtrial 2nþ1 i2
2GDc 1Dnþ1
ð1 Dnþ1 Þ
1 3q2 pnþ1 2 2q1 fnþ1 cosh r2y 1 þ q3 fnþ1 2ry 3
~ 6 tol THEN RETURN to Box 2. IF U (4) GO TO ( 2Þ The residual system of equations for the plastic corrector algorithm differentiated at the converged state results in the following identity:
2 @r
Dc
6 @ Dc 6 @rp 6 @ Dc 6 6 @rf 6 @ Dc 6 6 @rR 6 @ Dc 4 @r D @ Dc
@r Dc @p
@r Dc @f
@r Dc @R
@rDc @D
@r p @p
@r p @f
@r p @R
@rp @D
@r f @p
@r f @f
@r f @R
@rf @D
@r R @p
@r R @f
@r R @R
@rR @D
@r D @p
@r D @f
@r D @R
@rD @D
3j
2
3
2
@r Dc @ ee trial
trial deednþ1
3
7 6 dnþ1 d Dc 7 7 6 @rp 7 6 trial 7 6 e trial deev nþ1 7 6 dpnþ1 7 7 7 6 @ ev nþ1 7 6 7 7 6 7 6 @r f e trial 7 ¼6 7 :6 dfnþ1 7 d e 7 7 6 e trial dnþ1 7 6 @ e 7 6 7 4 dRnþ1 7 7 5 6 dnþ1 7 7 6 0 5 5 4 dDnþ1 @r D e trial d e dnþ1 @ ee trial dnþ1
ð72Þ
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
211
In order to determine the terms, and, the inversion of the above system of derivatives is required and can be re-written as:
3
2
2
C 1;1 dDc 7 6 6 6 dpnþ1 7 6 C 2;1 7 6 6 6 dfnþ1 7 ¼ 6 C 3;1 7 6 6 7 6 6 4 dRnþ1 5 4 C 4;1
C 1;2
C 1;3
C 1;4
C 2;2
C 2;3
C 2;4
C 3;2
C 3;3
C 3;4
C 4;2
C 4;3
C 4;4
C 5;1
C 5;2
C 5;3
C 5;4
dDnþ1
3 2 @rDc e trial trial dednþ1 3 @ eednþ1 7 C 1;5 6 7 6 p trial 7 76 @rtrial deev nþ1 C 2;5 76 @ eev nþ1 7 7 76 6 @rf e trial 7 C 3;5 7 76 @ ee trial dednþ1 7 7 76 dnþ1 C 4;5 56 7 7 6 0 5 C 5;5 4 @rD e trial d e dnþ1 @ ee trial
ð73Þ
dnþ1
where, the matrix C is the inversion matrix of the differentiation residual equations in order to each variable of the problem. After some algebraic manipulation, the term, and can be written as: @r
@r
@r
c p f D dDc ¼ C 1;1 @ee Dtrial C 1;2 @ee trial C 1;3 @ ee trial C 1;5 @ e@re trial dnþ1
@r
c ddnþ1 ¼ C 5;1 @ee Dtrial dnþ1
@r
v nþ1 @r p C 5;2 @ee trial v nþ1
dnþ1
dnþ1
@r
f D C 5;3 @ee trial C 5;5 @e@re trial dnþ1
@r
@r
v nþ1
dnþ1
dnþ1
ð74Þ
c p f D dpnþ1 ¼ C 2;1 @ee Dtrial C 2;2 @ee trial C 2;3 @ee trial C 2;5 @e@re trial dnþ1
dnþ1
Substituting the above expressions into Eq. 70, the close form for the tangent operator can be determined. More details about how to obtain the derivative of each residual equation in function of the elastic strain tensor can be verified in the appendix. 5. Calibration strategy Regarding the determination of material parameters for the proposed constitutive model, two independent calibration points are required and a so-called multi-objective function is used. The first calibration point is taken from a specimen at high level of stress triaxiality like a smooth bar specimen under tensile loading condition. In this step, the hardening law, ry ðRÞ, for the undamaged model and the critical volume void fraction, fc , are determined as well as the set of parameters for nucleation of micro void mechanism ðfN ; SN ; eN Þ. The second calibration point can be taken from a specimen in simple shear loading condition, where the accelerator parameter, q6, and the critical shear damage, Dc, are determined as well as the set of parameters for the nucleation of micro defects mechanism ðDN ; S0N ; e0N Þ. In this case, the butterfly specimen is used under simple shear loading condition. Eq. 75 represents the multi-objective function used in parameters identification.
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # v # u N " Num u N " Num Exp u1 X f u 1 X g ðwÞ g Exp ð v Þ f i i i i t þ ; Sðv ; wÞ ¼ t N i¼1 N i¼1 fiExp g Exp i
ð75Þ
where Sðv ; wÞ is the multi-objective function, based on the least squares, fiNum and fiExp represent, respectively, numerical and experimental reaction forces for the calibration point in high level of stress triaxiality, g Num and g Exp and represent, respeci i tively, numerical and experimental reaction forces for the calibration point in low level of stress triaxiality, both calibration points determined for each incremental displacement of the process and N is the number of increments or the number of experimental points. The terms v and w represent vectors of material properties which are variables in the optimization process and can be written as:
v ¼ ðry ; fN ; SN ; eN Þ
and w ¼ ðDN ; S0N e0N ; q6 Þ:
ð76Þ
The search method used in the optimization process is based on gradient method and is expressed as:
v ðkþ1Þ ¼ v k þ ak
OSðv k ; wk Þ ; kOSðv k ; wk Þk
wðkþ1Þ ¼ wk þ ak
OSðv k ; wk Þ ; kOSðv k ; wk Þk
ð77Þ
where kOSðv k ; wk Þk represents the gradient of the multi-objective function. Based on the parameter identification method described previously, the material properties for an aluminum alloy 2024-T351 and a steel 1045 were obtained. Furthermore, regarding numerical tests, three different loading conditions were carried out: simple shear, shear/tensile and shear/compression. Assuming the limitation of the size of this paper, numerical tests are presented using only the effect of Xue shear mechanism. 5.1. Geometry and mesh definition Regarding the material properties for the first calibration point, a classical smooth bar specimen is used and the following dimensions were employed (see Fig. 9):
212
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 9. The geometry for the smooth bar specimen. Dimension in (mm). Taken from Teng (2008).
Fig. 10. Finite elements meshes for (a) aluminum alloy and (b) steel, regarding the gauge section.
In order to trigger necking, a dimensional reduction of 5% in the central diameter of the specimen is used. However, different gauges sections are taken regarding the experimental data (see Teng and Wierzbicki, 2006). For the aluminum alloy and steel, gauges sections of 25.4 mm and 20.6 mm are used, respectively. The standard eight-nodded axsymmetric quadrilateral element, with four Gauss integration points, is adopted. The initial meshes discretisation for both cases are illustrated in Fig. 10, where only one symmetric quarter of the problem, with the appropriate symmetric boundary conditions imposed to the relevant edges, is modeled. A total number of 1800 elements have been used in the discretization of both the smooth specimens, amounting to a total of 5581 nodes. For the second calibration point and all numerical tests that will be presented, a butterfly specimen is used. The specimen was initially designed by Bai (2008) and the geometry and general dimensions can be verified by Fig. 11. In this case, a three dimensional finite element mesh of 3392 twenty nodded elements, with nine Gauss integration points, is used amounting to 17,465 nodes (see Fig. 12).
5.2. First calibration point: smooth bar specimen (tensile loading condition) In the present section, the stress–strain curves, the parameters required for modeling micro void nucleation mechanism from the GTN model are calibrated by tensile tests in cylindrical smooth bars. Through experimental tests conducted on both materials (see Teng, 2008 and Bai, 2008), the reaction versus displacement curves are determined as well as the stress–strain curves for an elasto-plastic model of von Mises type. The inverse method is adopted in order to calibrate the material parameters for coupled damage model by forcing the numerical solution to be, as close as possible to the experimental results.
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
213
Fig. 11. The geometry for butterfly specimen. Dimension in (mm). Taken from Bai (2008).
Fig. 12. Finite elements mesh for butterfly specimen.
Figs. 13a and 14a show reaction curves for the model determined after the application of inverse method. A good agreement between the experimental and numerical results can be observed. Furthermore, the critical volume void fraction is also determined in the point where the model attains the displacement to fracture, experimentally observed (see Figs. 13b and 14b). The critical values obtained are fc ¼ 0:06 and fc ¼ 0:76, for aluminum alloy 2024-T351 and steel 1045, respectively. The results for the calibration procedure, in terms of stress–strain curve, can also be observed in Fig. 15, where both curves, for uncoupled and coupled damage models, were determined.
5.3. Second calibration point: butterfly specimen (simple shear loading condition) In this calibration point, the parameters related with the micro defects nucleation mechanism are determined as well as the critical value for the shear damage. Furthermore, details related to how to determine the accelerator parameter, q6, are described. The butterfly specimen is here used under simple shear loading condition and the displacement to fracture,
214
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 13. (a) Reaction versus displacement curve for GTN model and experimental results to aluminum alloy 2024-T351. (b) Critical volume void fraction parameter calibrated for the material.
Fig. 14. (a) Reaction versus displacement curve for GTN model and experimental results to steel 1045. (b) Critical volume void fraction parameter calibrated for the material.
experimental determined for both materials, was suggested by Bai (2008). An inverse method is also adopted, regarding the calibration of the parameters by forcing the numerical results to be as close as possible to the experimental data. The behavior of parameter q6 can be better understood by looking at Figs. 16 and 17, where the reaction versus displacement curve and the evolution of the shear damage parameter, assumes the definition of the shear mechanism by Xue, is observed for different values of q6. According to the value of the accelerator parameter, a different critical shear damage value, Dc, is also established. Fig. 18 represents the critical shear damage as a function of the value of q6. Tables 1 and 2 contain the material properties found for the proposed model using the optimization method described. Remark 2. In order to study the influence of the spatial discretization on the numerical results, several numerical simulations were performed using different mesh refinements. It was possible to conclude that, although the numerical results can be affected by the discretization, there is no strong dependence of the numerical results with the level of mesh refinement. The location of fracture onset was also not affected by the level of mesh refinement. This is mainly due to the fact that the level of damage attained by the specimens is still relatively low, for the applied displacement to fracture, and the overall softening effect is very small (see Figs. 13, 14, 16 and 17). Therefore, the meshes selected and used in this paper have
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
215
Fig. 15. Stress–strain curves determined for an uncoupled and coupled models: (a) for aluminum alloy 2024-T351 and (b) steel 1045.
Fig. 16. (a) Reaction versus displacement curve and (b) evolution of the shear damage parameter, for different values of q6 (steel 1045).
Fig. 17. (a) Reaction versus displacement curve and (b) evolution of the shear damage parameter, for different values of q6 (aluminum alloy 2024-T351).
216
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 18. Representation of the critical shear damage, Dc, as a function of the accelerator parameter, q6.
Table 1 Materials properties for aluminum alloy 2024-T351: (E = 72.400 MPa,
v = 0.33, q1 = 1.5, q2 = 1, q3 = 2.25 and k = 0.1).
fN
SN
eN
fc
DN
S0N
e0N
q6
Dc
S(u)
0.04
0.2
0.1
0.060
0.08
0.15
0.10
1.0
0.08
0.0012
Table 2 Materials properties for steel 1045: (E = 220.000 MPa,
v = 0.33, q1 = 1.5, q2 = 1, q3 = 2.25 and k = 0.1).
fN
SN
eN
fc
DN
S0N
e0N
q6
Dc
S(u)
0.05
0.2
0.1
0.076
0.10
0.15
0.10
1.0
0.16
0.0009
Table 3 Reference values for different loading scenarios of two materials. Angle (°)
0 10 5
Aluminum alloy 2024-T351
Steel 1045
uf
ef a
Fracture location
uf
ef a
Fracture location
0.70 0.50 1.00
0.22 0.26 0.22
Surface of the critical zone (node 2) Middle of the critical zone (node 1) Surface of the critical zone (node 2)
1.03 0.42 1.71
0.50 0.36 0.60
Surface of the critical zone (node 2) Middle of the critical zone (node 1) Surface of the critical zone (node 2)
a The values for the equivalent plastic strain were estimated by Bai (2008) using an equation which represents the three dimensional fracture locus. A mixed experimental and numerical procedure was applied.
given numerical results that do not change noticeably with the spatial discretization and the conclusions of the assessment are not perceptibly affected.
6. Numerical results For a consistent analysis of the proposed formulation at low level of stress triaxiality, some numerical tests are performed using the butterfly specimen and the implicit algorithm suggested in above sections. Three different loading conditions are taken as: simple shear, shear/tensile and shear/compression, and two materials as: aluminum alloy 2024-T351 and steel 1045. The performance of some parameters as equivalent plastic strain and displacement at fracture as well as the ability to predict the correct fracture location are evaluated. At the end, the numerical results determined by the new formulation are compared with results obtained by other shear mechanisms as Xue (2008) and Nahshon and Hutchinson (2008). 6.1. Evolution of equivalent plastic strain and damage parameters The experimental results obtained by Bai, 2008, which will be used as reference for comparison, are listed in Table 3. In particular, the displacement at fracture, uf , the expected value for the equivalent strain at fracture, ef , and the location of
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
217
Fig. 19. Evolution of damage parameter and evolution of equivalent plastic strain for 1045 steel: (a, b) simple shear, (c, d) shear/tensile and (e, f) shear/ compression loading conditions.
218
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 20. Evolution of damage parameter and evolution of equivalent plastic strain for 2024-T351 alloy: (a, b) simple shear, (c, d) shear/tensile and (e, f) shear/ compression loading conditions.
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
219
Fig. 21. Evolution of effective damage parameter for 1045 steel: (a) simple shear, (b) shear/tensile and (c) shear/compression loading conditions.
crack initiation are listed for each loading condition and material. The evolution of the internal variable are analyzed according the critical node for each experimental test (see Fig. 12 and Table 3). All numerical results obtained with the proposed model, which will be presented in this section, were conducted following the same approach. For each simulation, the butterfly specimen was subjected to a prescribed displacement that matches with the respective critical one experimentally observed that is listed in Table 3. The evolution of the damage variables and of the equivalent plastic strain for each loading condition and material are presented in Figs. 19–22, regarding the critical nodes (see Fig. 12 and Table 3). For a simple shear loading condition, it can be observed, in Figs. 19a and 20a, for the steel 1045 and the aluminum alloy, respectively, that there is an evolution of the shear damage parameter with the applied displacement and the volume void fraction does not grow and remains equal to zero. Hence, the introduction of the new damage variable allows the prediction of failure with the GTN original model and, in this case, plays the main damage role. Under a combined shear/tensile loading with an angle of 10°, it is possible to observe in Figs. 19c and 20c an evolution of both the shear damage variable and the volume void fraction variable. For this loading case, the prediction of crack initiation is established when the effective damage variable reaches unity (see Eq. 53). Due to the presence of a multi axial stress state, an additional factor is introduced by the term ð1 þ fc =Dc Þ in Eq. (53), to accelerate the prediction of fracture. In Figs. 19e and 20e, there is a combined shear/compression loading condition and, the shear damage parameter also plays a dominant role in the prediction of fracture. In this case, there is a crack closure effect and the degradation of material occurs due to the formation of shear bands, which can be captured by the proposed shear mechanism. The volume void fraction, in this case, is reduced to zero due to a negative hydrostatic pressure. In Fig. 19b, d and f, the evolution of equivalent plastic strain parameter is shown for different numerical parameter, q6, using the steel 1045. It can be observed that the numerical parameter does not have a strong influence over the evolution of
220
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 22. Evolution of effective damage parameter for 2024-T351 alloy: (a) simple shear, (b) shear/tensile and (c) shear/compression loading conditions.
Table 4 Numerical results for butterfly specimen using 1045 steel and different loading conditions. Angle (°)
Experimental data
Numerical results
uf
ep
q6
uf
ep
nav
hav
f
d
0
1.03
0.50
0.5 1.0 1.5
1.03 1.03 1.03
0.516 0.522 0.528
0.022 0.022 0.021
0.061 0.060 0.057
0.000 0.000 0.000
0.122 0.160 0.204
10
0.42
0.36
0.5 1.0 1.5
0.33 0.44 0.59
0.257 0.353 0.440
0.241 0.245 0.257
0.477 0.485 0.507
0.018 0.026 0.030
0.045 0.053 0.061
5
1.71
0.60
0.5 1.0 1.5
1.71 1.71 1.71
0.611 0.612 0.616
0.066 0.065 0.065
0.173 0.173 0.173
0.000 0.000 0.000
0.100 0.126 0.153
this internal variable. However, it is expected that the parameter q6 affects the evolution of the shear damage variable and effective damage, which can be seem in Fig. 21 that presents the evolution of the effective damage parameter for different values of q6, for both combined loading conditions and simple shear, for the 1045 steel. According to the value of the numerical parameter, the failure condition is established. From the analysis of Fig. 18, for q6 ¼ 0:5, q6 ¼ 1:0 and q6 ¼ 0:5, the failure
221
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228 Table 5 Numerical results for butterfly specimen using aluminum 2024-T351 alloy and different loading conditions. Angle (°)
Experimental data
Numerical results
uf
ep
q6
uf
ep
nav
hav
f
d
0
0.70
0.22
1.0 2.0 3.0 4.0
0.70 0.70 0.70 0.70
0.292 0.298 0.305 0.318
0.018 0.017 0.017 0.017
0.048 0.048 0.047 0.046
0.000 0.000 0.000 0.000
0.084 0.107 0.137 0.179
10
0.50
0.26
1.0 2.0 3.0 4.0
0.55 0.63 0.75 0.75
0.230 0.271 0.336 0.337
0.250 0.254 0.257 0.264
0.486 0.492 0.494 0.502
0.013 0.017 0.021 0.021
0.032 0.039 0.051 0.056
5
1.00
0.22
1.0 2.0 3.0 4.0
1.00 0.98 0.95 0.93
0.414 0.424 0.432 0.455
0.066 0.065 0.064 0.063
0.176 0.173 0.169 0.165
0.000 0.000 0.000 0.000
0.084 0.110 0.140 0.190
Fig. 23. Effective damage contour for butterfly specimen using 1045 steel.
condition is met when the shear damage variable is equal to Dc ¼ 0:12, Dc ¼ 0:16 and Dc ¼ 0:20, respectively. Based on the results presented in Fig. 21, it can be observed that when q6 ¼ 1:0, for the 1045 steel, the constitutive formulation predicts a crack formation closer to experimental evidence. For the aluminum alloy 2024-T351, the evolution of equivalent plastic strain is represented by Fig. 20b, d and e and a similar behavior is also observed, since the numerical constant, q6, has a small impact in the evolution . For this material, it is possible to conduct a similar analysis for the value of the constant q6, though the examination of Figs. 18 and 22, and conclude that q6 ¼ 1:0 also predicts failure closer to experimental evidence than the other values (see Fig. 22). The set of all results can also be analyzed in Table 4 and Table 5 which list all the numerical results and expected values, experimentally observed, for the equivalent plastic strain and displacement at fracture. Based on numerical results presented, it can be concluded that the new formulation has the ability to predict the correct moment to crack formation by appropriately calibrating the numerical constants and parameters of the model. Both the equivalent plastic strain and the displacement, calculated by present formulation, are in close agreement with the experimental data for both loading conditions and materials applied (see Table 4 and Table 5).
6.2. Prediction of the correct fracture location Another important feature to be analyzed, in order to validate the proposed constitutive formulation, is the ability to predict the correct fracture location. Based on experimental tests performed by Bai (2008), using the butterfly specimen (see Table 3), it can be concluded that under a simple shear loading condition, the micro crack is initiated in the surface of the
222
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Fig. 24. Effective damage contour for butterfly specimen using aluminum 2024-T351 alloy.
Fig. 25. Damage parameter contour for butterfly specimen using 1045 steel. (a) Nahshon and Hutchinson shear mechanism, (b) Xue shear mechanism and (c) new proposition. Section AA in the critical zone.
critical zone (node 2 of Fig. 12). However, when combined shear/tensile loading condition is applied, the crack is formed in the middle of the thickness (node 1 of Fig. 12) and grows toward the surface of the critical zone. Under a combined shear/ compression loading (5°), the surface of the critical zone is also the location of crack formation (node 2 of Fig. 12). Figs. 23 and 24 present the contour of effective damage for the steel 1045 and the aluminum alloy 2024-T351, respectively, at fracture. It is possible to conclude that the proposed damage formulation has the ability to predict the correct fracture location in all loading conditions. The predictive ability of the new model, in terms of fracture location, can also be compared against two recent extensions of the GTN model. The accuracy of the three models in the prediction of the fracture location is evaluated here for specifically a combined shear/tensile loading of 10° for the steel 1045. Fig. 25a illustrates the contour of the damage parameter for Nahshon and Hutchinson shear mechanism (Nahshon and Hutchinson, 2008), Fig. 25b for Xue shear mechanism (Xue, 2008) and Fig. 25c for the proposed model. It can be observed that only the new model predicts the initiation of failure in agreement with experimental evidence (node 1 – middle of thickness). The prediction of the GTN model extended with Xue’s shear
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
223
Fig. 26. Projection of three dimensional fracture loci for steel 1045.
Fig. 27. Projection of three dimensional fracture loci for aluminum 2024-T351 alloy.
mechanism is in complete disagreement with experimental evidence and the prediction of the GTN model extended with Nahshon and Hutchison shear mechanism is somewhat spread along the critical section, which may suggest an inaccuracy to the model. 7. Representation in three dimensional fracture locus In order to qualitatively judge the constitutive model, the numerical results, listed in Tables 4 and 5, are represented, within the projection of the so-called three dimensional fracture locus, on the space of stress triaxiality versus equivalent plastic strain. The fracture locus was originally proposed by Bao and Wierzbicki (2004), for the aluminum alloy 2024T351, and by Bai (2008), for the 1045 steel. In the Bao’s proposition (2004), the fracture locus is only dependent on the stress triaxiality and the equivalent plastic strain. However, Bai (2008) proposed a surface represented by the interpolation of the equivalent plastic strain, stress triaxiality and Lode angle. The fracture locus for both steel 1045 and 2024-T351 aluminum alloy were calibrated by classical specimens in different loading condition (see Bao and Wierzbicki, 2004; Bai, 2008). Fig. 26 represents the projection of 3D fracture locus on the space of equivalent plastic strain versus stress triaxiality. It can be observed that the numerical results, obtained by the proposed model, have a good agreement with the calibrated curve for different levels of the Lode angle. In the same figure, the numerical results obtained by the Nahshon and Hutchinson and Xue formulations are also plotted, using only simple shear and combined shear/tensile loading conditions. In both cases, these formulations do not present an uniform behavior in the prediction of the expected equivalent plastic strain at fracture. In Fig. 27, the results for the aluminum 2024-T351 alloy are reported. It can be observed that simple shear and combined shear/tensile loading conditions have a good agreement with the expected values (less than 20 percent of difference). However, regarding the combined shear/compression loading condition, the numerical equivalent plastic strain is in disagreement with the result reported by Bai (2008). The numerical prediction of the proposed formulation can be obtained for values of q6 less than 1. In the presented simulations, values of q6 were limited in 1, regarding the complexity of the calibration procedure. 8. Conclusion In this contribution, a new formulation was proposed to improve the original GTN model, regarding the ability to predict ductile fracture under a low level of stress triaxiality. Firstly, a new shear mechanism was proposed that is a function of the equivalent plastic strain, stress triaxiality and Lode angle. This mechanism can capture the elongation and rotation of microdefects, when shear loading condition is present. Furthermore, a new micro-defects nucleation mechanism was proposed
224
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
which is responsible for triggering the evolution of the shear damage parameter, since the new mechanism is independent on the volume void fraction. Then, the new damage parameter was coupled in GTN constitutive formulation in such a way that only affects the deviatoric stress contribution. Thus, the new model has two independent damage parameters: first one affecting only the hydrostatic stress component and the other affecting the deviatoric stress component. Numerical tests were conducted, with an implicit integration algorithm, in order to evaluate the formulation ability to predict the crack formation. A butterfly specimen was employed and two different materials: the steel 1045 and the aluminum alloy 2024-T351 were used. In all loading conditions, the model behaves well, either in the determination of the correct level of equivalent plastic strain and displacement at fracture, or in prediction of the location for crack formation. The introduction of two damage parameters affecting separate components of the stress tensor critically affects the evolution of internal variables and allows more accurate values at the time of crack formation. Furthermore, the introduction of a new micro-defects nucleation mechanism facilitates the calibration model and thus an improved performance for a wide range of stress triaxiality. The introduction of the stress triaxiality dependence in the evolution of shear damage parameter also enhanced the prediction of the fracture location under combined loading conditions, since this parameter influences the behavior of material under low stress triaxiality. An effective damage variable is determined in post-processed step as a function of both volume void fraction and shear damage parameter. A penalization factor is introduced in order accelerate the damage evolution due to the presence of multi axial loading conditions. In spite of the best performance of this formulation when compared to the models available in the literature, the introduction of more parameters that need to be calibrated requires special attention. In particular, two calibration points are required to fully define the model. A calibration point at high triaxiality, which was already required in GTN original model, and now a new point at low triaxiality, to obtain the parameters that govern the new shear damage evolution law. In summary, the new model was formulated in order to perform well in all loading conditions and for different materials. From the results presented, it is possible to conclude that the objective was achieved for the cases tested. Acknowledgments Supported by Portuguese Science and Technology Foundation (FCT), under scholarship no. SFRH/BD/45456/2008 and under Grant no. PTDC/EME-TME/71325/2006. Appendix A Derivative of each residual equation in relation of each variable of the problem. (a) Derivatives for rDc
@r Dc @ Dc @r Dc @p @r Dc @f @r Dc @R
@r Dc @D
trial trial 8G3 eednþ1 : eednþ1 h i3 ; ð1 Dnþ1 Þ2 1 þ 2G1Dc Dnþ1 3q2 pnþ1 ; ¼ q1 q2 fnþ1 ry sinh 2ry 2 3q2 pnþ1 ¼ r2y q3 fnþ1 q1 cosh ; 3 2ry 2 3q2 pnþ1 2 ¼ rHy 1 þ q3 fnþ1 2q1 fnþ1 cosh 3 2ry 3q2 pnþ1 ; q1 q2 fnþ1 pnþ1 H sinh 2ry 8 9 > > < = 1 4GDc ¼ J trial ; h i h i 2nþ1 2 3 > > 3 2GDc :ð1 Dnþ1 Þ2 1 þ 2GDc ; ð1 D Þ 1 þ nþ1 1Dnþ1 1Dnþ1
ða:1Þ
@rp 3q2 pnþ1 ; ¼ K ry q1 q2 sinh @ Dc 2ry
(b) Derivatives for rp
@r p 3 3q2 pnþ1 ; ¼ 1 þ DcKq1 ðq2 Þ2 fnþ1 cosh 2 @p 2ry @r p 3q2 pnþ1 ; ¼ DcK ry q1 q2 sinh @f 2r y @r p 3q2 pnþ1 3q p 3q2 pnþ1 2 nþ1 cosh ; ¼ DcK ry q1 q2 fnþ1 H sinh @R 2ry 2ry 2r y @r p ¼ 0; @D
ða:2Þ
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
225
(c) Derivatives for r f
" p 2 # p @ Df n fN 1 enþ1 eN Dep enþ1 eN @ Dep trial ; ¼ khnþ1 k pffiffiffiffiffiffiffi exp 1 2 @ðjÞ SN SN SN @ðjÞ SN 2 p
ða:3Þ
and,
@ Df g 3q2 pnþ1 ; ¼ ð1 fnþ1 Þry q1 q2 fnþ1 sinh @ Dc 2ry @ Df g 3Dc 3q2 pnþ1 ; ¼ ð1 fnþ1 Þ q1 ðq2 Þ2 fnþ1 cosh @p 2 2r y @ Df g 3q2 pnþ1 3q2 pnþ1 þ ð1 fnþ1 ÞDcry q1 q2 fnþ1 sinh ; ¼ Dcry q1 q2 fnþ1 sinh @f 2ry 2r y @ Df g 3q2 pnþ1 3q p 3q2 pnþ1 þ 2 nþ1 cosh ; ¼ ð1 fnþ1 Þ DcHq1 q2 fnþ1 sinh @R 2ry 2ry 2ry
ða:4Þ
@ Df g ¼ 0: @D (d) Derivatives for rR
@r R 1 3q2 pnþ1 2 3q2 pnþ1 2 þ 1 þ q3 fnþ1 ¼ q1 q2 fnþ1 sinh 2q1 fnþ1 cosh ry ; 3 @ Dc ð1 fnþ1 Dnþ1 Þ 2ry 2ry @r R Dy 3q2 pnþ1 3q2 pnþ1 3q2 pnþ1 fnþ1 pnþ1 cosh q1 q2 fnþ1 sinh ; ¼ @p ð1 fnþ1 Dnþ1 Þ 2ry 2ry 2ry @r R Dy 3q2 pnþ1 2 3q2 pnþ1 2 þ ¼ q q f sinh r 1 þ q f q f cosh nþ1 y nþ1 1 2 3 1 nþ1 2 3 @f 2ry 2r y ð1 fnþ1 Dnþ1 Þ Dy 3q2 pnþ1 2 3q2 pnþ1 þ ry 2q3 fnþ1 2q1 cosh ; q q fnþ1 sinh 3 ð1 fnþ1 Dnþ1 Þ 1 2 2ry 2ry ( @r R Dc 3q ðq Þ2 3q2 pnþ1 2H 3q2 pnþ1 2 þ 1 þ q3 fnþ1 ¼1 2q1 fnþ1 cosh 1 22 fnþ1 ðpnþ1 Þ2 H cosh 3 @R ð1 fnþ1 Dnþ1 Þ 2ðry Þ 2ry 2ry H 3q2 pnþ1 ; þ2q1 q2 fnþ1 pnþ1 sinh ry 2ry @r R Dy 3q2 pnþ1 2 3q2 pnþ1 2 þ 1 þ q ¼ q q f p sinh f 2q f cosh ry : nþ1 nþ1 1 2 nþ1 3 nþ1 1 2 3 @d 2ry 2ry ð1 fnþ1 Dnþ1 Þ ða:5Þ (e) Derivatives for r D
" p 2 # p @ DDn DN 1 enþ1 eN Dep enþ1 eN @ Dep p ffiffiffiffiffiffi ffi 1 exp k ; ¼ k1 htrial nþ1 2 @ðjÞ SN SN SN @ðjÞ SN 2p
ða:6Þ
and
@ DDg @g0 p p @ Dep e De þ q4 Dnþ1 q5 g 0 ðep þ Dep Þ ¼ q4 Dq5 ; nþ1 @ Dc @ Dc @ Dc @ DDg @g0 p p @ Dep e De þ q4 Dnþ1 q5 g 0 ðep þ Dep Þ ¼ q4 Dq5 ; nþ1 @p @p @p @ DDg @g0 p p @ DDg ¼ q4 Dq5 ðe De Þ ; nþ1 @p @f @f @ DDg @g0 p p @ DDg ¼ q4 Dq5 ðe De Þ ; nþ1 @p @R @R @ DDg @g0 p p @ Dep ðq5 1Þ e De þ q4 Dq5 p p ¼ q4 q5 Dnþ1 g 0 ep Dep þ q4 Dq5 : nþ1 nþ1 g 0 ðe þ De Þ @d @D @D
ða:7Þ
226
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
(f) Derivatives for Dep in relation of each variable of the problem trial S trial @ Dep pffiffiffi D 4G nþ1 : S nþ1 pffiffiffi ¼ a Dc ; 3 @ Dc c 3 a ½1 þ ð 2GDc Þ ð1 Dnþ1 Þ3 1Dnþ1 @ Dep 1 3q2 pnþ1 3q2 pnþ1 2 cosh ; ¼ Dc pffiffiffi ðq1 Þ2 ðq2 Þ3 ðfnþ1 Þ ry sinh @p 2ry 2ry 3 a ( )
2 @ Dep 2 3q2 pnþ1 ¼ Dc pffiffiffi ðq1 Þ2 ðq2 Þ3 fnþ1 ðry Þ2 sinh ; @f 2ry 9 a @ Dep 2 3q2 pnþ1 3q2 pnþ1 3q p 3q2 pnþ1 2 sinh 2 nþ1 cosh ; ¼ Dc pffiffiffi ðq1 Þ2 ðq2 Þ3 ðfnþ1 Þ ry H sinh @f 2ry 2ry 2ry 2ry 9 a 8 9 > > < = @ Dep 1 2 4GDc trial ¼ Dc pffiffiffi S trial : S : h i nþ1 nþ1 2 3 2GDc > > @p 3 a : ð1 D4nþ1 Þ½1 þ ð1Dnþ1 Þ ; ð1 þ Dnþ1 Þ3 1 þ ð 2GDc Þ
ða:8Þ
1Dnþ1
(g) Derivatives for
g 00 nþ1
in relation of each variable of the problem
" # gnþ1 pnþ1 @g 0 q6 2G trial ¼ g 0 lnð1 khnþ1 kÞ ; ð1 Dnþ1 Þ @ Dc ðkgnþ1 k þ q6 Þ2 kgnþ1 k qtrial nþ1 ( ) gnþ1 1 @g 0 q6 2G trial ; ¼ g 0 lnð1 khnþ1 kÞ 1þ 1 Dnþ1 @ Dc ðkgnþ1 k þ q6 Þ2 kgnþ1 k qtrial nþ1 " # gnþ1 pnþ1 @g 0 q6 2G : ¼ g 0 lnð1 khtrial nþ1 kÞ 2 kg 2 trial @ Dc k q ðkgnþ1 k þ q6 Þ nþ1 nþ1 ð1 Dnþ1 Þ
ða:9Þ
(h) Derivatives of each residual equation in relation of the elastic trial strain
@r Dc @ eed trail nþ1
( ¼
)2
2G
2GDc Þ ½1 þ ð1D nþ1
1 ee trail : Id ; 1 Dnþ1 d nþ1
@p ¼ k; @ eev trail nþ1 @r f @ eed trail nþ1
" p 2 # @khtrial fN 1 enþ1 eN nþ1 k pffiffiffiffiffiffiffi exp ¼ Dep ; 2 SN @ eed trail nþ1 SN 2p
khtrial nþ1 k @r d @ eed trail nþ1
" p 2 # p fN 1 enþ1 eN Dep enþ1 eN @ Dep pffiffiffiffiffiffiffi exp ; 1 2 SN SN SN @ eed trail SN 2 p nþ1
¼
p @khtrial dN 1 enþ1 eN nþ1 k p ffiffiffiffiffiffi ffi 2 @ Dep ; exp 2 SN @ eed trail nþ1 SN 2p
k1 htrail nþ1 k " q4 dnþ1
q5
ða:10Þ
" p 2 # p dN 1 enþ1 eN Dep enþ1 eN @ Dep pffiffiffiffiffiffiffi exp 1 ; 2 SN SN SN @ eed trail SN 2 p nþ1 @g 0
@ eed trail nþ1
p
p
p
e De þ g 0 ð e
# @ Dep þ De Þ e trail : @ ed nþ1 p
(i) Others derivatives
92 8 = @ Dep 2Dc < 2G h i p ffiffiffi ¼ ee trial : Id ; trial @ eednþ1 3 a : 1 þ 2GDc ; dnþ1
ða:11Þ
1Dnþ1
0 1 trial trial h @khnþ1 k 2 B C nnþ1 ¼ nþ1 @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA e trial ; trial 2 @ eednþ1 @ ednþ1 khtrial nþ1 k p 1n nþ1
ða:12Þ
227
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
2
3
trial trial e trial 1 9 deteednþ1 deteednþ1 nnþ1 27 6 d e trial 7 ¼ e 4 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 5 ednþ1 5 : I ; q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi dnþ1 3 trial 2 2 @ eednþ1 3 e trial 3 e trial e trial e trial e : e e : e dnþ1 dnþ1 2 dnþ1 2 dnþ1
8 < htrial k @k q6 q7 nþ1 htrial k htrial k ¼ 1 k þ g ln 1 k 0 nþ1 nþ1 trial trial : kg kgnþ1 k þ q6 @ eednþ1 @ eednþ1 @g 0
ða:13Þ
9
= gnþ1 3pnþ1 ð2GÞ3 trial i eednþ1 : Id : 2 trial 3 h kgnþ1 k 2 q 1 þ 2GDc ; nþ1 k þ q6 nþ1 q6
1Dnþ1
ða:14Þ
References Andrade Pires, F.M., César de Sá, J.M.A., Costa Sousa, L., Natal Jorge, R.M., 2003. Numerical modelling of ductile plastic damage in bulk metal forming. Int. J. Mech. Sci. 45, 273–294. Bai, Y., Wierzbicki, T., 2007. A new model of metal plasticity and fracture with pressure and Lode dependence. Int. J. Plast. 24, 1071–1096. Bai, Y., 2008. Effect of Loading History on Necking and Fracture. Ph.D Thesis, Massachusetts Institute of Technology. Bao, Y., 2003. Predition of Ductile Crack Formation in Uncrecked Bodies. Ph.D Thesis, Massachusetts Institute of Technology. Bao, Y., Wierzbicki, T., 2004. On fracture locus in the equivalent strain and stress triaxiality space. Int. J. Mech. Sci. 46 (81), 81–98. Barsoum, I., Faleskog, J., 2007a. Rupture in combined tension and shear: experiments. Int. J. Solids Struct. 44, 1768–1786. Barsoum, I., Faleskog, J., 2007b. Rupture in combined tension and shear: micromechanics. Int. J. Solids Struct. 44, 5481–5498. Besson, J., Steglich, D., Brocks, W., 2001. Modeling of crack growth in round bars and plane strain specimens. Int. J. Solids Struct. 38 (46–47), 8259–8284. Besson, J., 2010. Continuum models of ductile fracture: a review. Int. J. Damage Mech. 19, 3–52. Bridgman, P., 1952. Studies in Large Plastic and Fracture. McGraw-Hill Book Company, London. Brunig, M., Gerke, S., Hagenbrock, V., 2013. Micro-mechanical studies on the effect of the stress triaxiality and the Lode parameter on ductile damage. Int. J. Plast. 50, 49–65. Butcher, C., Chen, Z., Bardelcik, A., Worswick, M., 2009. Damage-based finite-element modeling of tube hydroforming. Int. J. Fract. 155, 55–65. Chaboche, J.L., 1984. Anisotropic Creep Damage in the Framework of Continuum Damage Mechanics. Nucl. Eng. Des. 79, 309–319. Chaboche, J.L., 2003. Damage mechanics. In: Milne, I., Ritchie, R.O., B. Karihaloo (Eds.), Comprehensive Structural Integrity, vol. 2. Elsevier-Pergamon, pp. 213–284. Chaboche, J.L., Boudifa, M., Saanouni, K., 2006. A CDM approach of ductile damage with plastic compressibility. Int. J. Fract. 137, 51–75. De Souza Neto, E.A., Peri’c, Owen, D.R.J., 2008. Computational methods for plasticity: theory and applications. John Wiley & Sons Ltd. Dowling, N., 1999. Mechanical Behavior of Materials. Eng. Methods for Deformation, Fracture, and Fatigfue, second ed. Prentice Hall. Engelen, Roy A.B., 2005. Plasticity-induced Damage in Metals/nonlocal modelling at finite strains. PhD Thesis, Technische Universiteit Eindhoven, Eindhoven. Eterovic, A.L., Bathe, K.-J., 1990. A hyperelastic based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int. J. Numer. Meth. Eng. 30, 1099–1114. Faleskog, J., Gao, X., Shih, C.F., 1998. Cell Model for Nonlinear Fracture Analysis I, Micromechanics Calibration. J. Fract. 89, 355, 73. Gao, X., Kim, J., 2006. Modeling of ductile fracture: significance of void coalescence. Int. J. Solids Struct. 43, 6277–6293. Gao, X., Zhang, G., Roe, C., 2009. A study on the effect of the stress state on ductile fracture. Int. J. Damage Mech. 19, 75–94. Gao, X., Zhang, T., Zhou, J., Graham, S.M., Hayden, M., Roe, C., 2011. 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. 27 (2), 217–231. Gologanu, M., Leblond, J.B., Perrin, G., Devaux, J., 1995. Recent extensions of Gurson’s model for porous ductile metals. In: Suquet, P. (Ed.), Continuum Micromechanics. Springer-Verlag. Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth – Part I. Yield criteria and flow rules for porous ductile media. J. Eng. Mater. Tech. 99, 2–15. Hancock, J.W., Mackenzie, A.C., 1976. On the mechanisms of ductile failure in high-strength steels subjected to multi-axial stress-states. J. Mech. Phys. Solids 24, 147–160. Jackiewicz, J., 2011. 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, 487–502. Johnson, G.R., Cook, W.H., 1985. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng. Fract. Mech. 21 (1), 31–48. Khan, A.S., Liu, H., 2012. A new approach for ductile fracture prediction on Al 2024–T351 alloy. Int. J. Plast. 35, 1–12. Krajcˇinovic´, D., Fonseka, G.U., 1981. The Continuous Damage Theory of Brittle Materials – Part 1: General Theory. J. Appl. Mech. 48, 809–815. Kim, J., Gao, X., Srivatsan, T.S., 2003. Modeling of crack growth in ductile solids: a three-dimensional analysis. Int. J. Solids Struct. 40, 7357–7374. Kim, J., Gao, X., Srivatsan, T.S., 2004. Modeling of void growth in ductile solids: effects of stress triaxiality and initial porosity. Eng. Frac. Mech. 71, 379–400. Lecarme, L., Tekoglu, C., Pardoen, T., 2011. Void growth and coalescence in ductile solids with stage III and stage IV strain hardening. Int. J. Plast. 27, 1203– 1223. Lemaitre, J., 1985. A continuous damage mechanics model for ductile fracture. J. Eng. Mater. Technol. Trans. ASME, vol. 107, pp. 83–89. Lemaitre, J., Chaboche, J.L., 1990. Mechanics of Solid Materials. Cambridge Univ, Press. Li, H., Fu, M.W., Lu, J., Yang, H., 2011. Ductile fracture: experiments and computations. Int. J. Plast. 27, 147–180. Malcher, L., Andrade Pires, F.M., César de Sá, J.M.A., 2012. An assessment of isotropic damage constitutive models under high and low stress triaxialities. Int. J. Plast. 30-31, 81–115. McClintock, F.A., 1968. A Criterion for Ductile Fracture by the Growth of Holes. J. Appl. Mech. 35, 363–371. McVeigha, C., Vernereya, F., Liu, W.K., Moran, B., Olson, G., 2007. An Interactive Micro-Void Shear Localization Mechanism in High Strength Steels. J. Mech. Phys. Solids 55, 225–244. Mirone, G., Corallo, D., 2010. A local viewpoint for evaluating the influence of stress triaxiality and Lode angle on ductile failure and hardening. Int. J. Plast. 26, 348–371. Mudry, F., 1985. Methodology and application of local criteria for prediction of ductile tearing. D.Reidel Publishing Co. Murakami, S., Ohno, N., 1981. A continuum theory of creep and creep damage. In: Ponter, A.R.S. (Ed.), Proceedings of the IUTAM Symposium on Creep in Structures, Leicester, Berlin, pp. 422–443. Nahshon, K., Hutchinson, J., 2008. Modification of the Gurson model for shear failure. Eur. J. Mech. A/Solids 27, 1–17. Needleman, A., Tvergaard, V., 1987. An Analysis of Ductile Rupture Modes at a Crack Tip. J. Mech. Phys. Solids 35, 151, 83. Pardoen, T., Hutchinson, J.W., 2000. An extended model for void growth and coalescence. J. Mech. Phys. Solids 48, 2467–2512. Peric, D., Owen, D.R.J., Honnor, M.E., 1992. A Model for Finite Strain Elasto-Plasticity Based on Logarithmic Strains: Computational Issues. Comput. Meth. Appl. Mech Eng. 94, 35–61.
228
L. Malcher et al. / International Journal of Plasticity 54 (2014) 193–228
Pineau, A. (1981). Review of fracture mechanisms and local approaches to predicting crack resistance in low strength steels. In: François, D. et al. (Ed.), Advances in Fracture Researches. New-York, Pergamon Press, ICF5, Cannes. Pineau, A., Pardoen, T., 2003. Damage mechanics. In: Milne, I., Ritchie, R.O., B. Karihaloo (Eds.), Comprehensive Structural Integrity, vol. 2. ElsevierPergamon, pp. 686–783. Reis, F.J.P., Malcher, L., Andrade Pires, F.M., César de Sá, J.M.A., 2010. A modified GTN model for the prediction of ductile fracture at low stress triaxialities. Int. J. Struct. Integrity 1 (4). Rice, J.R., Tracey, D.M., 1969. On the ductile enlargement of voids in triaxial stress fields. J. Mech. Phys. Solids 17, 201–217. Rousselier’s, G., 1980. Finite deformation constitutive relations including ductile fracture damage. In: Nemat-Nasser (ed.) Three-Dimensional Constitutive Relations and Ductile Fracture. North-Holland Publ. Comp. 1981, 331–355. Rousselier’s, G., 1987. Ductile fracture models and their potential in local approach of fracture. Nucl. Eng. Des. 105, 97–111. Rousselier’s, G., 2001. The Rousselier model for porous metal plasticity and ductile fracture. In: Lemaitre, J. (Ed.), Handbook of Materials Behavior Models, vol. 2. Academic Press, New York, pp. 436–445 (Chapter 6.6). Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Springer-Verlag, New York. Socie, D.F., Marquis, G.B., 2000. Multiaxial Fatigue. SAE International, Warrendale. Stoughton, T.B., Yoon, J.W., 2011. A new approach for failure criterion for sheet metals. Int. J. Plast. 27, 440–459. Teng, X., Wierzbicki, T., 2006. Evaluation of six fracture models in high velocity perforation. Eng. Fract. Mech. 73 (12), 1653–1678. Teng, X., 2008. Numerical prediction of slant fracture with continuum damage mechanics. Eng. Fract. Mech. 75, 2020–2041. Tvergaard, V., Needleman, A., 1984. Analysis of the cup-cone fracture in a round tensile bar. Acta Met. 32, 157–169. Xue, L., 2007. Ductile Fracture Modeling – Theory, Experimental Investigation and Numerical Verification, Ph.D Thesis, Massachusetts Inst. of Technology. Xue, L., 2008. Constitutive modeling of void shearing effect in ductile fracture of porous materials. Eng. Fract. Mech. 75, 3343–3366. Xia, L., Shih, C.F., Hutchinson, J.W., 1995. Computational Approach to Ductile Crack Growth under Large Scale Yielding Conditions. J. Mech. Phys. Solids 43, 389–413. Zhang, K.S., Bai, J.B., Francois, D., 2000. Numerical analysis of the influence of the lode parameteron void growth. Int. J. Solids Struct. 38 (32–33), 5847–5856.