Numerical prediction of slant fracture with continuum damage mechanics

Available online at www.sciencedirect.com

Engineering Fracture Mechanics 75 (2008) 2020–2041 www.elsevier.com/locate/engfracmech

Numerical prediction of slant fracture with continuum damage mechanics X. Teng * Impact and Crashworthiness Laboratory, Room 5-218, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 15 December 2006; received in revised form 18 September 2007; accepted 3 November 2007 Available online 12 November 2007

Abstract Ductile specimens always exhibit an inclined fracture surface with an angle relative to the loading axis. This paper reports a numerical study on the cup–cone fracture mode in round bar tensile tests and the slant fracture in plane-strain specimens based on continuum damage mechanics. A combined implicit–explicit numerical scheme is first developed within ABAQUS through user defined material subroutines, in which the implicit solver: Standard, and the explicit solver: Explicit, are sequentially used to predict one single damage/fracture process. It is demonstrated that this numerical approach is able to significantly reduce computational cost for the simulation of fracture tests under quasi-static or low-rate loading. Comparison with various tensile tests on 2024-T351 aluminum alloy is made showing good correlations in terms of the load–displacement response and the fracture patterns. However, some differences exist in the prediction of the critical displacement to fracture. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Tensile tests; Cup–cone fracture; Damage evolution rule; Ductile fracture; ABAQUS

1. Introduction Ductile fracture is governed by nucleation, growth, and coalescence of microscopic voids and cracks, accompanied with large plastic deformation. In this deterioration process, a material element gradually loses its load-carrying capability until complete fracture. A number of numerical procedures have been developed to model ductile crack formation and growth in structural components. These numerical methods can be categorized into three types: (1) sudden failure approach, (2) porous plasticity models, and (3) continuum damage mechanics (CDM) [1]. In the first approach, a fracture criterion is uncoupled with a plasticity constitutive model, and damage accumulation is separated from the calculation of stresses and strains. An element is assumed to abruptly fail as soon as its damage indicator reaches a critical value. Under isothermal condition, the removal of an element usually occurs at the point when its stress reaches a peak value. This approach has been widely adopted in leading commercial finite element codes such as ABAQUS and LS-DYNA, because of *

Tel.: +1 617 253 2104; fax: +1 617 253 8125. E-mail address: [email protected]

0013-7944/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2007.11.001

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2021

Nomenclature A initial yield stress B hardening coefficient C fourth order tensor of material elastic constants cD correction to the damage indicator cp correction to the effective plastic strain p cr correction to the internal variable D damage indicator Dcr critical damage indicator E Young’s modulus of a virgin material f yield function G shear modulus Geff effective shear modulus Idev deviatoric fourth order unit tensor I second order unit tensor J Jacobian consistent tangent matrix J3 third invariant of a stress tensor K bulk modulus Keff effective bulk modulus N plastic flow direction tensor n hardening exponent p effective plastic strain pf critical effective plastic strain at fracture q von Mises equivalent stress R(r) isotropic hardening function Rd residual function about the damage indicator Rf residual function about the yield stress f Rr residual function about the internal variable r rf critical internal variable at fracture r internal variable S0 material constant s deviatoric stress tensor s1, s2, s3 three principal, deviatoric stress components Y damage strain energy release rate e total strain tensor ee elastic strain tensor ep plastic strain tensor k plastic multiplier r stress tensor rh hydrostatic stress n dimensionless third invariant of a stress tensor ð~Þ effective variable considering the reduction of a resisting area D( ) increment of a variable ð_Þ differentiation with respect to ‘‘time”

easy calibration and simple formulation. However, the sudden failure approach often fails to capture sophisticated fracture patterns such as slant rupture in quasi-static tests, since softening due to damage is neglected. The second approach, represented by Gurson’s [2] and Rousselier’s [3] porous plasticity model, takes into account effects of void growth on the yield surface. Both models are capable of simulating various fracture

2022

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

patterns, e.g. the cup–cone fracture mode in round bar tensile tests [4,5]. However, both models are suffering from some deficiencies. For instance, a total of nine material parameters in Gurson–Tvergaard–Needleman model have to be calibrated from tests besides the true stress–strain curve for a fully dense matrix material. This would be a formidable task for industrial applications, although empirical values for some parameters are available. Continuum damage mechanics introduces a damage indicator as an internal variable to describe the progressive degradation of mechanical properties of materials. In contrast to micromechanics-based porous plasticity models, continuum damage mechanics is formulated basically in a phenomenological way, and is relatively simple. This approach will be exploited in the present paper to study ductile fracture. A number of formulations based on continuum damage mechanics have been proposed for various applications in the literature. Some of early contributions in this area were collected in three books by Kachanov [6], Lemaitre [7], as well as Lemaitre and Chaboche [8]. Here, we attempt to summarize the recent progress on the application of continuum damage mechanics to ductile fracture. A comprehensive survey on various formulation and applications of continuum damage mechanics is beyond the subject of this study. Lemaitre [9] first formulated a continuum damage mechanics framework for ductile fracture based on the thermodynamics of irreversible processes. Basic features of continuum damage mechanics were summarized by Chaboche [10,11]. Wang [12] proposed a new damage evolution rule for ductile fracture, in which damage accumulates in a nonlinear way with respect to plastic strains. Bonora [13] identified three possible trends of damage evolution for various metals. Hambli [14] introduced Lemaitre’s coupled damage model into ABAQUS/Standard through a user defined material subroutine and simulated sheet-metal cutting processes. Børvik et al. [15] replaced Lemaitre’s evolution law of damage with the one suggested by Johnson and Cook [16]. The coupled material model was applied to simulate the perforation process of a target plate under projectile impact within LS-DYNA [17]. Andrade Pires et al. [18] improved the original Lemaitre’s damage model by introducing ‘‘crack closure effects”. The applicability of the model was demonstrated through the calculation of damage accumulation in three typical metal forming processes using an in-house explicit finite element code. Recently, Chaboche et al. [19] proposed a modified continuum damage model by considering the plastic compressibility in the context of void growth. In these studies, attention was directed to the modification, improvement, or generalization of continuum damage models for the prediction of ductile fracture. Either the calibration of material properties or the verification of the applicability of continuum damage models to a variety of problems was omitted. Both topics are important to convince the industry of practical values of continuum damage mechanics. By contrast, the present paper focuses on the calibration of damage parameters, the implementation into commercial finite element codes, and the simulations of a number of ductile fracture tests. Based on these results, one would be able to critically evaluate the applicability of continuum damage mechanics to ductile fracture problems. To predict plastic deformation all the way to crack formation and propagation, explicit numerical procedures are often used to solve governing equations with finite element codes. Because of the removal of failed elements and the coupling of damage with stress response, it would be extremely difficult for implicit numerical procedures to obtain convergence solutions. Explicit numerical procedures are essentially developed for transient, dynamic problems. It would be expensive if an explicit code is used to simulate the damage and fracture process of a specimen under quasi-static loading. To reduce computational cost, two techniques are commonly used: mass scaling by artificially increasing mass density and/or time scaling by artificially increasing an applied velocity. This inevitably introduces dynamic effects in the form of inertial force and kinetic energy into quasi-static problems. By contrast, the present paper for the first time proposes an implicit–explicit numerical procedure within ABAQUS for low-rate loading problems. The implicit code ABAQUS/Standard is first applied to model a plastic deformation process until fracture initiates or convergence becomes very slow. The explicit code ABAQUS/Explicit will then take over the calculation to simulate the remainder of the process characterized by crack formation and growth. In such a way, one would be able to substantially reduce computational cost and at the same time ensure the accuracy of solutions. To demonstrate the applicability of the newly implemented damage model and the implicit–explicit numerical scheme, a series of fracture tests on 2024-T351 aluminum alloy performed by Bao [20] are studied numerically. The experiments considered in the present paper include tensile tests on one smooth and two notched round bars, and a transverse plane-strain specimen. Although these tests are routinely conducted in labs to calibrate plasticity and fracture properties, the numerical modeling of the fracture processes and patterns still

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2023

poses a challenge. One of the outstanding difficulties is to predict realistic fracture modes such as the cup–cone fracture pattern in round bar tensile tests and the slant fracture mode of transverse plane-strain specimens under tension. Conventional finite element analysis always predicts a flat fracture surface perpendicular to the loading direction. The present paper first develops a numerical algorithm based on the formulation originally proposed by Lemaitre [9]. This numerical algorithm is implemented into ABAQUS/Standard and at the same time ABAQUS/Explicit by means of user defined material subroutines. In the second part, a true stress–strain curve and damage parameters are calibrated for 2024-T351 aluminum alloy from round bar tensile tests. With the implemented subroutines and the obtained material data, an attempt is made to numerically predict the cup–cone fracture mode of round bars under tension and the slant fracture mode of a transverse plane-strain specimen under tension. The applicability of the present formulation is examined by comparing the load–displacement response and the fracture patterns to test results. Meanwhile, the efficiency of the new numerical scheme is demonstrated through the simulations. 2. Basic formulation of continuum ductile damage mechanics The present paper considers a hypoelastic–plastic material model, in which the total strain tensor e can be decomposed into an elastic part ee and a plastic part ep: e ¼ ee þ ep :

ð1Þ

From the hypothesis of strain equivalence [9], the stress tensor can be written by r ¼ ð1  DÞC : ee ;

ð2Þ

where r is the stress tensor; the damage indicator D is introduced to represent the reduction of the resisting area of an material element; and C is the fourth order tensor of material elastic constants, given by C ¼ 2GIdev þ KI  I;

ð3Þ dev

where G and K are the shear and the bulk modulus; I is the deviatoric fourth order unit tensor and I is the second order unit tensor. Microvoids and microcracks are assumed to evenly distribute in all directions and thus the damage indicator D can be expressed by a scalar quantity. Physically more-sound anisotropic damage tensors were also proposed in the literature, e.g. see [21,22]. However, this leads to a more complicated formulation and calibration procedure than a scalar damage. Eq. (2) indicates that the reduction of the effective resisting area due to damage is equivalent to the degradation of elastic modulus. The damage indicator can be thus expressed in another form: e E ; ð4Þ E e are the virgin and the damaged Young’s modulus, respectively. If the variation of Young’s where E and E modulus during a loading process is measurable, damage accumulation can be quantified. Also from the hypothesis of strain equivalence, the isotropic yield function can be modified to account for damage q f ¼  RðrÞ 6 0; ð5Þ 1D D¼1

where R(r) is the isotropic hardening function; r is the internal variable; and q is the von Mises equivalent stress defined by rffiffiffiffiffiffiffiffiffiffiffi 3 s : s; ð6Þ q¼ 2 where s is the deviatoric stress tensor. Note, that ~q ¼ q=ð1  DÞ is the equivalent value of the effective stress tensor. Eq. (5) indicates that the material yield surface shrinks as damage grows. The normality rule gives of N e_ p ¼ k_ ¼ k_ ; or 1D

ð7aÞ

2024

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

r_ ¼ k_

of _ ¼ k; oR

ð7bÞ

where k is the plastic multiplier; and N is the plastic flow direction tensor given by N¼

of 3 ~s 3 s ¼ ¼ : q 2q o~ r 2~

ð8Þ

Taking the inner product of both sides of Eq. (7a) yields _ k_ ¼ r_ ¼ ð1  DÞp; where p_ is the rate of the accumulated plastic strain, defined by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 p p p_ ¼ e_ : e_ : 3

ð9Þ

ð10Þ

Note, that p > r is always valid. Without damage, the internal variable r reduces to the accumulated plastic strain p. In some formulations, the internal variable was taken to be identical to the accumulated plastic strain, e.g. see [18]. Lemaitre [9] proposed that the damage indicator evolves with respect to the accumulated plastic strain: Y _ D_ ¼ p; S0

ð11Þ

where S0 is the material constant and has the dimension of energy per unit volume; and Y is the damage strain energy release rate, given by ~2 1 q Y ¼ ee : E : ee ¼ Rv ; 2 2E

ð12Þ

where  2 2 rh ; Rv ¼ ð1 þ mÞ þ 3ð1  2mÞ 3 q

ð13Þ

and rh is the hydrostatic(mean) stress. Y is conceptually similar to stress intensity factors (SIF) and J-integral in elastic fracture mechanics, all of which are associated with the energy release rate. It should be pointed out that this fracture criterion does not discern the difference in damage accumulation between tension and compression, since the quadratic form of the stress triaxiality is defined. Lemaitre’s continuum damage model requires the calibration of two material parameters: S0 and the critical damage at the point of fracture Dcr. By comparison, nine material constants need to be identified from tests in the Gurson–Tvergaard–Needleman model. While the material parameters for a number of metals have been determined, e.g. see [8], the applicability of Lemaitre’s damage model to a variety of practical problems remains to be verified. This completes the description of Lemaitre’s coupled damage-plasticity constitutive model. The formulation of this model within the framework of the thermodynamic laws was presented by Lemaitre and Chaboche [9]. In summary, damage introduces two effects into material constitutive relationships: the degradation of Young’s modulus and the shrinkage of the yield surface. It is realized that various formulations for continuum damage mechanics have been developed, which more or less improve the original model of Lemaitre. However, the present paper does not intend to propose a new constitutive equation but to critically evaluate the application of Lemaitre’s damage model to ductile fracture. It is expected that some conclusions obtained can be extended to similar damage models. 3. Implicit integration algorithm This section develops a numerical integration algorithm for the above coupled damage-plasticity model. The algorithm is based on the elastic predictor, return mapping scheme, which is widely used in computational

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2025

plasticity. At each ‘‘time” increment, the stress tensor and material state variables such as the accumulated plastic strain and the damage indicator will be updated for a given elastic/plastic strain tensor increment De. For the conciseness, the subscript n + 1 will be omitted in the following derivation and all the variables without the subscript will be calculated at tn+1. Assume first that the increment is purely elastic within the time interval [tn, tn+1]. The trial effective stress tensor ~strial can thus be written as ~strial ¼ 2Gðejeln þ DeÞ;

ð14Þ

el

where ejn is the elastic strain tensor at tn; De is the deviatoric part of the strain tensor increment; and G is the shear modulus. If the equivalent value of this trial effective stress tensor, ~qtrial , satisfies the yield condition f 6 0, the assumption is valid and then all the material state variables related to plastic deformation do not evolve. Otherwise, one has to introduce a ‘‘plastic” correction to the strain increment. The solutions for the stress tensor and the material state variables have to be found at tn+1. In the presence of plastic deformation in the current increment, the effective stress tensor is given by ~s ¼ ~strial  2GDpN;

ð15Þ

where Dp is the increment of the accumulated plastic strain; and the notation D() indicates an increment D() = ()n+1  ()n. Note, that the plastic flow direction vector N is of the following useful property: N¼

3 s 3 ~s 3 ~strial ¼ ¼ : q 2~ 2q 2~ qtrial

With the aid of the above expression, Eq. (15) can be transformed into   2 ~ q þ 2GDp N ¼ ~strial : 3

ð16Þ

ð17Þ

Taking the inner product of both sides with themselves yields ~ q þ 3GDp ¼ ~ qtrial :

ð18Þ

The yield condition, Eq. (5), requires f ~ q  RðrÞ ¼ ~ qtrial  3GDp  RðrÞ ¼ 0:

ð19Þ

At the same time, Euler backward time discretization of Eqs. (9) and (11) gives the following equations: Dr ¼ ð1  DÞDp;

ð20aÞ

DD ¼ yð~ q; rh ÞDp;

ð20bÞ

where yð~ q; rh Þ ¼

Y ð~ q; rh Þ : S0

ð21Þ

The problem is reduced to find the solutions for the increments: Dr, Dp, and DD from Eqs. (19) and (20a). Newton’s iteration method is applied in this study. In each iteration, we have dR cr þ 3Gcp ¼ 0; dr Rr þ cr  ð1  DÞcp ¼ 0; oy dR Dpcr ¼ 0; Rd þ cD  ycp  o~ q dr Rf þ

ð22aÞ ð22bÞ ð22cÞ

where the notation c() = ()i+1  ()i denotes a ‘‘correction” to the solution at the (i + 1)th iteration; and Rf, Rd, and Rr are three residual functions, defined, respectively, by

2026

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Rf ¼ RðrÞ þ 3GDp  ~ qtrial ; Rr ¼ Dr  ð1  DÞDp;

ð23aÞ ð23bÞ

Rd ¼ DD  yDp:

ð23cÞ

From Eq. (23a), one is able to obtain the explicit expressions for the ‘‘corrections” for each variable: cp ¼

Rf þ dR R dr r ; 3G þ ð1  DÞ dR dr

cr ¼ ð1  DÞcp  Rr ; oy dR Dpcr  Rd : cD ¼ ycp þ o~ q dr

ð24aÞ ð24bÞ ð24cÞ

Hence, expensive matrix inversions are not needed. In the beginning of the iteration process, an elastic increment is assumed so that all the ‘‘plastic” increments and the residual functions except Rf vanish. The initial corrections are given by cp ¼

Rf ; 3G þ ð1  DÞ dR dr

ð25aÞ

cr ¼ ð1  DÞcp ;

ð25bÞ

cD ¼ ycp :

ð25cÞ

The above integration algorithm is summarized in Table 1. Although only one scalar damage variable is added, the formulation of the coupled damage-plasticity model is much more complicated than conventional isotropic elastic-plasticity models, in which only one nonlinear equation need to be solved. To efficiently integrate the coupled constitutive model, many algorithms were proposed in the literature. For example, Doghri [23] developed a numerical scheme for a coupled damage model with nonlinear isotropic and kinematic hardening, in which a nonlinear system of 14 equations have to be solved. de Souza Neto [24] successfully reduced the problem to only one nonlinear equation for the case with isotropic hardening. The present algorithm is similar to them but is relatively simple and stable for a variety of problems. ABAQUS/Standard uses Newton’s iteration method to solve nonlinear global equilibrium equations. This requires a material subroutine to provide a consistent tangent stiffness matrix (consistent Jacobian matrix) at each integration point. An accurate definition of the matrix is required for rapid convergence. Consistent tangent stiffness matrices are determined when the implicit numerical integration of the local constitutive equations has converged. For nonlinear isotropic and kinematic hardening materials coupled with damage, Doghri [25] derived an exact expression for the consistent tangent stiffness matrix. A simplified matrix was obtained by de Souza Neto [24] for coupled damage-isotropic hardening plasticity materials. Even so, the expression is Table 1 Summary of the integration algorithm Elastic predictor. Compute ~strial from Eq. (14) and ~qtrial from Eq. (6) if the yield condition is satisfied, i.e. Eq. (5) is valid, then Update the solutions: ~s ¼ strial , r = rn, p = pn, D = Dn, and exit else Compute the initial corrections c() by assuming the elastic increment from Eq. (25a) Loop over plastic iterations: Update variables Dp, Dr, and DD Calculate the residual functions from Eq. (23a) if tolerance criteria are satisfied on the residuals then Exit the loop else Compute corrections c(): cp from Eq. (24a); cr from Eq. (24b); and cD from Eq. (24c) end if end if

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2027

quite complex. In this study, an approximate expression for the consistent tangent matrix is developed by following the derivation given in ABAQUS user’s manual [26] for the conventional isotropic plasticity model. This approximate matrix takes into account the degradation of the material constants. The expression is given by "  # eff 4 dR J ¼ 2Geff Idev þ K eff I  I þ  3Geff N  N ð26Þ 9 dr where Geff and Keff are the effective shear and bulk moduli, defined by Geff ¼ Gð1  DÞ

~ q ; ~ qtrial

ð27Þ

and K eff ¼ Kð1  DÞ; and (dR/dr)eff is the effective hardening modulus, defined by  eff 3ð1  DÞG dR dR dr ¼ : dr 3ð1  DÞG þ dR dr

ð28Þ

ð29Þ

If the present increment is purely elastic, the exact consistent tangent stiffness matrix is given by J ¼ ð1  Dn ÞC:

ð30Þ

While the present expression for the elastic–plastic consistent tangent matrix is approximate, it is able to give convergent solutions after a limited number of iteration. 4. Numerical implementation Either an implicit finite element code or an explicit code can be used to model a damage/fracture process. It seems that no one used both types of codes to solve one single problem in the literature. Implicit finite element codes solve a system of equilibrium equations at each time/time-like increment in an iterative way. Programs always tend to take a large increment as long as solutions converge and thus these implicit codes are efficient to simulate a large plastic deformation process. However, if a problem involves strong nonlinearity and/or complex contact, it would be difficult for implicit codes to obtain convergence results. In the present case, the quick loss of load-carrying capability of materials at the point of fracture initiation would lead to a premature end of calculation. Meanwhile, implicit finite element codes usually prohibit to completely remove failed elements because of numerical stability. In the literature, it is often assumed that failed elements, whose accumulated damage has reached or exceed a critical value, are of ceratin residual stiffness and still participate in the rest of calculation, e.g. [27]. This inevitably introduces numerical errors. Therefore, implicit finite element codes have some limitations in the calculation of crack formation and propagation. As opposed to implicit codes, explicit finite element programs integrate governing equations in an explicit way. While explicit programs are developed essentially for transient, dynamic problems, they are also powerful to determine solutions for strong nonlinear problems and large-scale models. Free from any convergence problems, explicit numerical schemes are able to remove failed elements to represent crack propagation. However, explicit numerical procedures are conditionally stable. This requires that time increments be sufficiently small and thus computational cost would be prohibitive for quasi-static problems. To speed up calculation, large elements have to be defined, and mass scaling and/or time scaling is often applied. These techniques inevitably introduce dynamic effects represented by inertial force and kinetic energy. Therefore, explicit finite element codes are not efficient for the simulation of plastic deformation of a specimen under low-rate loading. The present paper proposes a combined implicit–explicit numerical scheme within ABAQUS. In contrast to conventional numerical procedures in the literature, the implicit solver: Standard, and the explicit solver: Explicit, are sequentially applied to simulate one single damage/fracture process in this scheme. Basically,

2028

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Standard is first used to calculate plastic deformation and damage accumulation up to fracture initiation. Then, Explicit continues the simulation until the complete failure of specimens. To fulfill this task, the ‘‘importing” technique has to be applied to transfer a deformed mesh and its associated material state variables such as plastic strains from Standard into Explicit. At the same time, two user defined material subroutines, one for Standard (UMAT) and the other for Explicit (VUMAT), have to be developed to calculate damage accumulation. Such a numerical procedure fully takes advantage of both the implicit and explicit codes. It is able to substantially speed up the calculation of damage/fracture processes in low-rate loading cases, and at the same time to give sufficiently accurate solutions. For most of fracture tests, specimens immediately fail once a crack is formed. About 99% of the loading process can be modeled with the implicit code, and hence it is unnecessary to artificially scale up mass density with the explicit code. To consistently import ABAQUS/Standard results to ABAQUS/Explicit calculation, some cares must be taken in the implementation of user defined material subroutines. These two finite element programs have different conventions and formulations. The first inconsistency relates to the order of stress and strain vectors. For example, a strain vector is stored in Standard in the order: e11, e22, e33, c12, c13, and c23 while the order of a strain vector in Explicit is e11, e22, e33, e12, e23, and e31. The second difference is that engineering shear strains are defined in Standard while tensor shear strains are used in Explicit, e.g. c12 = 2e12. 5. Calibration of material parameters for 2024-T351 Al The calibration of true stress–strain curves up to the point at fracture is not a trivial but rather time-consuming task. Localization phenomena such as shear banding and necking, which act as precursors to fracture, result in the nonuniform distribution of stresses over a cross-section. An inverse method is commonly used by adjusting the true stress–strain curve in numerical modeling to minimize the difference in the load–displacement curve between numerical prediction and experimental data. To measure damage, Lemaitre and Chaboche [8] proposed a novel type of ‘‘hourglass-like” specimens with a weakened central section, on which plastic deformation and damage would concentrate. Elastic unloading is performed at different tensile loading levels to determine the variation of Young’s moduli, and damage is calculated using Eq. (4). Meanwhile, a small strain gauge is attached on the central section to record the accumulated plastic strain. In such a way, one would be able to describe damage evolution in terms of the accumulated plastic strain and to determine damage parameters. The true stress–strain curve can then be separately calibrated in the same way as the one for uncoupled material models. It should be pointed out that using this approach the measured damage and accumulated plastic strain are actually the average values over the weakened cross-section. In the present study, the stress–strain curve and damage parameters are calibrated simultaneously from round bar tensile tests on 2024-T351 aluminum alloy. To establish the fracture locus of this alloy, Bao and Wierzbicki [28] performed a series of experiments including tensile tests on round bars, upsetting tests on short cylinders, and shear tests on butterfly specimens. These tests provide a solid foundation to develop various fracture models. In this study, the tensile test on one smooth round bar is used to determine the true stress–strain curve, the material constant S0, as well as the critical damage value Dcr. The inverse method is adopted in the calibration by forcing the numerical prediction to be as close as possible to the experimental results. The first two items are determined by matching the load–displacement curve and Dcr was obtained by means of the critical displacement to fracture. The tensile tests on the notched round bars are not used in the calibration but serve as an examiner to ensure that the material parameters are reasonably obtained. After a number of iterations, this calibration procedure yields the optimized values of S0 = 6.0 MPa and Dcr = 0.80. The true stress–strain curve of 2024-T351 aluminum alloy is shown in Fig. 1, which can be fitted well with A = 352.0 MPa, B = 853.4 MPa, and n = 0.66. For the purpose of comparison, the stress–strain curve for the uncoupled method obtained early by Bao [20] from the tensile test on the smooth round bar is superposed. It appears that the difference between them increases with the accumulated plastic strain. The uncoupled model implicitly takes damage into account in the calibration and thus the obtained stress–strain curve is of low hardening and small yield stresses at large plastic deformation. Details of the numerical modeling of the round bar tensile tests will be presented in Section 7.1.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2029

Equivalent Stress (MPa)

Coupled model (Present)

Uncoupled model (Bao [20])

Accumulated plastic strain Fig. 1. True stress–strain curves for 2024-T351 aluminum alloy.

6. Comparison with ductile fracture locus for 2024-T351 Al It has been generally agreed that ductile fracture is controlled mainly by two variables: the stress triaxiality and the accumulated plastic strain, following the pioneer work by McClintock [29] and Rice and Tracey [30] on the enlargement of a single void. A fracture locus formulated in this two-dimensional space is often employed to describe material ductile fracture properties, e.g. [31,16,28], etc. In these studies, a number of tests were performed to determine fracture loci, which are essentially semi-empirical. Lemaitre’s continuum damage model provides another way to construct ductile fracture loci, since the model contains the stress triaxiality and the plastic strain at the same time. It would be of interest to make a comparison of fracture loci for a specific material between these two procedures. This section will develop a ductile fracture locus from the continuum damage model. This type of work was first conducted by Lemaitre [9] in analogy to Rice and Tracey’s void growth model. Here, a different approach is adopted. By assuming the constant stress triaxiality during loading, the accumulated plastic strain at the point of fracture can be derived from Lemaitre’s evolution law of damage. First, Eq. (A.6) can be transformed into a nonlinear equation for rf: A2 rf þ 2AB

rnþ1 r2nþ1 ES 0 f þ B2 f ¼ ½1  ð1  Dcr Þ2  nþ1 2n þ 1 Rv

ð31Þ

Given the critical damage, the internal variable at fracture rf can be solved in terms of the stress triaxiality. Integration of Eq. (A.7) from 0 to rf gives the accumulated plastic strain pf to fracture as a function of the stress triaxiality:  ffi Z rf Z rf sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nþ1 2nþ1 1 Rv r r þ B2 1 ¼ ð1  DÞdr ¼ A2 r þ 2AB dr: ð32Þ pf ES 0 nþ1 2n þ 1 0 0 With the true stress–strain curve and the damage parameters known, the fracture locus can be generated from the above equations. A short MATLAB program was composed to fulfill this task. As mentioned earlier, Bao and Wierzbicki [28] developed an empirical fracture locus for 2024-T351 aluminum alloy based on a series of tests and corresponding numerical simulations. This fracture locus is uncoupled with the material elasto-plastic model. Fig. 2 compares this empirical fracture criterion with the newly calibrated one from Lemaitre’s damage model. It appears that Lemaitre’s law of damage evolution is able to qualitatively capture the tendency that material ductility decreases with the stress triaxiality in the high positive range but gives a completely unrealistic prediction in the negative range. Hence, Lemaitre’s law of damage evolution cannot be applied to the case dominated by compression. This deficiency had been pointed out by Lemaitre in his original paper [9]. Lemaitre’s damage rule also indicates that materials always have the

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Plastic strain at fracture

2030

CDM with D cr = 0.80

CDM with D cr = 0.18

Bao-Wierzbicki locus

Stress triaxiality Fig. 2. Comparison of ductile fracture loci between empirical Bao–Wierzbicki’s fracture model and Lemaitre’s evolution rule of damage. The circles represent the fracture points using the uncoupled method from 15 tests by Bao and Wierzbicki [28].

maximum ductility at shear. This result is questionable. McClintock [32] compared the fracture strains between shear and tension for five metals. Four of them have higher fracture strains at tension than at shear. Finally, the present calibrated damage model predicts much higher fracture strains than the Bao–Wierzbicki (BW) empirical fracture locus, though the same set of the round bar tensile tests were used in the calibration. By artificially lowering the value of the critical damage to Dcr = 0.18, the fracture locus for the coupled model becomes close to the BW one. However, this critical damage would predict premature fracture for the smooth round bar. Note, that the shape of the fracture locus is implicitly assumed in Lemaitre’s damage model as given in Eq. (32). This reduces the required number of tests, but may introduce the limitation of the applicability of the damage model. By contrast, an empirical fracture locus calibrated from a series of tests would more realistically represent ductile fracture properties. 7. Prediction of slant fracture 7.1. Round bar tensile tests The cup–cone fracture mode of a ductile round bar under tension has long been studied. Rogers [33] experimentally explored the failure micromechanism from void nucleation and coalescence to macrocrack formation and propagation. Bluhm and Morrissey [34] performed well-controlled tensile tests on round bars using a sufficiently stiff testing machine so that the fracture process could be stopped at any stage for most of metals. They found that the sharp ‘‘knee” of the load–displacement curve is closely associated with the onset of the macroscopic crack near the symmetry axis. Tvergaard and Needleman [4] first numerically studied the cup–cone fracture pattern. The paper has been widely cited, mostly for its modified Gurson’s material constitutive model but not for the prediction of the cup–cone fracture. Besson et al. [5] attacked the same type of problems using Rousselier’s and Gurson’s continuum damage model. Scheider and Brocks [35] applied a cohesive model to simulate the fracture process of a round bar under tension. Recently, Xue [36] predicted the formation of the cup–cone fracture mode by introducing the third invariant of the deviatoric stress tensor into a new evolution rule of damage. All these investigations indicate that in order to successfully capture the cup–cone fracture pattern, one has to carefully determine elements aspect ratios, mesh density, and element types, etc. By contrast, the present paper uses continuum damage mechanics to predict the fracture process and patterns of round bars under tension. Considered here are one unnotched and two notched round bars. Fig. 3 illustrates the geometry size of these round specimens. Three finite element models were built using four-node,

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2031

reduced-integration axisymmetric elements(CAX4R) with the minimum element size 0.05 mm  0.05 mm in the gauge sections. A displacement was prescribed at one end of the specimens in the ABAQUS/Standard simulation. The calculation automatically stops when the maximum value of the accumulated damage reaches the critical data, or the convergence rate becomes extremely slow. With the ‘‘importing” technique, ABAQUS/Explicit continues the calculation until complete fracture. In the ABAQUS/Explicit modeling, a low velocity of 0.01 m/s and the real mass density are defined so that dynamic effects are minimized. With the present implicit–explicit numerical scheme, only about thirty minutes were needed in a modern personal computer for each case. By contrast, it would take several days if ABAQUS/Explicit is used all the way from plastic deformation to crack growth. Note, that many research finite element codes with implicit solvers are indeed able to simulate crack propagation, e.g. see [37]. Comparison of the load–displacement curves is given in Fig. 4 showing good agreement between the numerical predictions and the test results. The critical displacements at the point of fracture are over-predicted for the cases with the notched round bars, although the same set of the damage parameters gives a rather good correlation for the case with the unnotched round bar. Fig. 5 presents the predicted fracture process of the smooth round bar. Initially, the bar deforms homogeneously in the gauge section. At a certain stage, necking deformation takes over and introduces larger triaxial tension around the symmetry axis. A macroscopic crack is formed first at the center of the necking zone, where the highest stress triaxiality is located. The crack propagates radially in the plane of the minimum cross-section. As it approaches the lateral surface, the crack deviates out of the plane at an angle of about 42°. The rest part of the round bar is finally sheared off to form the cone part of the cup–cone fracture mode. It can be observed from Fig. 5 that the load steeply drops to zero with very little further deformation from crack formation to complete separation. At the same time, necking deformation almost terminates as soon as the crack is formed. Hence, either the reduction in area or the elongation of the gauge section at the point of fracture can be used as an indicator for the onset of a macroscopic crack in round bar tensile tests. This justifies that a local effective plastic strain to fracture can be determined from numerical simulations using these data. The fracture scenarios are marked in the load–displacement curve as the circles. It appears that macrocrack formation corresponds to the sharp ‘‘knee” of the load–displacement curve. This is consistent with experimental observations by Bluhm and Morrissey [34]. As the crack grows, a series of elastic unloading stress waves propagate away from the free surfaces of the newly generated crack. The specimen completely breaks before all these stress waves arrive at the ends. Hence, the reaction force measured at the end does not immediately drop to zero at this point. Note, that the present tensile tests are assumed to be kinematically controlled. 28.3

10.0

7.20 7.50

4.47

(a) Smooth round bar. 21.67

7.50

12.0

4.0

(b) Round bar of the notch radius R = 12 mm. 26.03 4.0 7.50

4.0

(c) Round bar of the notch radius R = 4 mm.

Fig. 3. Configurations and geometry data of three round bar specimens. Unit: mm.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Force (kN)

2032

Test data 0.05 mm × 0.05 mm 0.1 mm × 0.1 mm

Relative displacement (mm) (a) Smooth round bar

Test data

Numerical results

Force (kN)

R = 4 mm

R = 12 mm

Relative displacement (mm) (b) Notched round bars

Fig. 4. Comparison of tensile load vs. displacement. The displacement was measured by a 25.4 in. extensometer. The test data were taken from Bao [20].

Comparison of the cup–cone fracture mode between the present numerical prediction and the experiment results is given in Fig. 6 showing rather good correlation. The present numerical simulation predicts the ratio of the radius of the cup portion to the radius of the whole fracture surface equal to 0.78. This value is very close to the data 0.72 measured from the post-test specimen. Many tensile tests, e.g. [34], indicate that this ratio would decrease with the increasing ductility. The out-of-plane crack growth was attributed to weak geometrical constraints on shearing as the crack propagates to the lateral surface, e.g. [33,4]. This can be clearly interpreted by means of the third invariant of the stress tensor. To distinguish various stress states, Xue [36] and Wierzbicki and Xue [38] introduced a dimensionless variable n, defined by n¼

27 J 3 3 2 r

 1 6 n 6 1;

ð33Þ

where J3 = s1s2s3 is the third invariant of the stress tensor, s1 P s2 P s3 are the three principal deviatoric stresses. It can be shown that uniaxial and plane-strain stress states are two limiting cases with jnj = 1 and n = 0, respectively. The plane-strain stress state is defined by s1 ¼ s3

s2 ¼ 0:

ð34Þ

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2033

Force (kN)

-

Relative displacement (mm)

Fig. 5. Crack formation and growth in the unnotched round bar.

Fig. 6. The cup–cone fracture mode of the unnotched round bar.

Thus in terms of the dimensionless variable n, pure shear is equivalent to the transverse plane-strain condition. The uniaxial stress state is defined by s1 ¼ s2 > 0

s3 ¼ 2s1

or

s1 ¼ 2s2

s2 ¼ s3 < 0:

ð35Þ

This condition is always satisfied at the symmetry axis of a round specimen but it loses its validity somewhere away from the axis after necking. Fig. 7 shows the variation of the normalized third invariant n through the radius of the minimum cross-section from the symmetry axis to the lateral surface. Before necking, the smooth round bar is subjected to uniaxial tension and thus n = 1 throughout the gauge section. Necking introduces triaxial loading and n decreases with the radius. The stress state of the outer annulus of the round bar moves towards the plane-strain condition. Crack formation near the symmetry axis accelerates this tendency. As the crack propagates, the region ahead of the crack tip undergoes a rapid change of the variable n. It is interesting to notice that the stress state of the most outer annulus moves away from the plane-strain condition and thus the resistance to shearing increases. To minimize energy dissipation, the crack may deviate away from the

2034

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Axisymmetric tension

Third invariant, ξ

Incipient necking Incipient Fracture Crack tip Plane-strain condition Crack tip Axisymmetric compression Normalized radial distance Fig. 7. The variation of the third stress invariant with necking deformation and crack growth for the smooth round bar.

minimum cross-section and thus the cup–cone fracture pattern is formed. For less ductile materials, a round bar breaks before deep necking is developed. The whole necking zone would still be dominated by the uniaxial stress state. Hence, a less ductile specimens always exhibits a macroscopically flat fracture surface. For more ductile materials, deep necking preceding fracture introduces the plane-strain stress state and thus the cup–cone fracture mode is easily generated. For the present notched round bars, the maximum damage is initially located at the root of the notches and gradually shifts to the symmetry axis. The macroscopic crack starts at the axis and propagates radially to the lateral surface. The cup–cone fracture mode is not so evident for the notched round bars, since the axisymmetry stress state always dominates the whole gauge section. Fig. 8 shows the numerically predicted fracture patterns of the notched round bars, which are in good agreement with that of the post-test specimens, see Fig. 9. It should be pointed out that the cup–cone fracture mode would be more easily captured, if more elements are defined in the radial direction and a larger critical damage value is assigned in the calculation. The mesh model with relatively large elements of 0.1 mm  0.1 mm in the gauge section, i.e. there are 45 elements through the radius, predicts a flat fracture surface, see Fig. 10. By comparison, the load–displacement response is not much sensitive to the mesh density, see Fig. 4a. It should be pointed out that Besson et al. [5,39] extensively studied effects of mesh discretization on the simulation of slant fracture, including element size, aspect ratios of elements, element types, and so on.

(a) Round bar of the notch radius R = 12 mm.

(b) Round bar of the notch radius R = 4mm.

Fig. 8. Numerically predicted fracture patterns of the notched round bars.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2035

Fig. 9. Fracture mode of the notched round bar of R = 12 mm.

7.2. Plane-strain specimens under tension The second example considered is slant fracture of a transverse plane-strain specimen under tension. In this type of tests shear banding often precedes crack formation and provides a path for crack growth, see, e.g. [40,41]. The conditions for the onset of shear bands were given by Hill [42], Rudnicki and Rice [43], and Rice [44]. It was found that a solid, characterized by the classical flow theory of plasticity with a smooth yield surface, is quite resistant to the bifurcation from diffuse necking into shear banding. Hutchinson and Tvergaard [45] inspected effects of material constitutive relationships on the prediction of the formation of shear bands. A few numerical studies on shear banding were also presented in the literature. Tvergaard et al. [46] performed a numerical analysis on the growth of shear bands in plane-strain tensile tests with a phenomenological corner theory of plasticity and a nonlinear elastic model. They showed that small geometrical imperfections in the form of initial thickness inhomogeneity would trigger various shear banding patterns. Using the same procedure for the simulation of the cup–cone fracture mode, Besson et al. [39] examined various aspects of numerical prediction of slant fracture of plane-strain specimens such as element size, element types, and mesh arrangement, etc. The transition of the fracture surface from a flat plane to a slant plane in a single-edge crack specimen subjected to remote Mode-I tension was studied by Mahgoub et al. [47]. Various physical mechanisms behind the formation of shear bands have been proposed, including the vertex structure of the yield surface [45] and the texture evolution for polycrystal metals [48]. By contrast, the present paper assumes that ‘‘damage softening” would induce the formation of shear bands, which is equivalent to thermal softening in dynamic cases with high strain rates. In reality, a crack would initiate and propagate along a fully developed shear band. However, in numerical modeling, the capture of shear banding does not guarantee that a crack would grow within a shear band with minimum energy dissipation. Conventional numerical procedures usually predict a flat fracture surface perpendicular to the loading axis. Here, continuum damage mechanics is used to investigate shear banding and following fracture. To evaluate effects of the strain ratio on the material ductility, Bao [20] conducted a tensile test on flat-grooved plates made of 2024-T351 aluminum alloy. This test is simulated all the way

Fig. 10. The fracture pattern of the smooth round bar predicted with the model of the element size 0.1 mm  0.1 mm.

2036

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

24.5 2.0 1.6 5.0

8.0 Fig. 11. Geometrical size of the transverse plane-strain specimens with the width 50 mm. Unit: mm.

to fracture in the present section. Fig. 11 displays the geometry size of the specimen, which is similar in shape to that designed by Clausing [49]. To ensure the plane-strain state in the center region, the ratio of the plate width (50 mm) to the length of the gauge section (8.0 mm) is sufficiently large. A two-dimensional finite element model was built using four-node, reduced-integration, plane-strain elements (CPE4R) with the minimum element size 0.05 mm  0.05 mm in the gauge section. A displacement was prescribed at one end of the specimen in the ABAQUS/Standard simulation while a constant tensile velocity of 0.04 m/s was defined within ABAQUS/Explicit. The total calculation time is about 2 h in a modern personal computer. Fig. 12 presents the predicted deformation and fracture patterns of the plane-strain specimen. When the tensile force reaches the peak value, diffuse necking occurs, see Fig. 15. At a certain stage, two narrow shear bands of highly plastic deformation takes over diffuse necking, although necking is still quite shallow. The shear bands have the angle of 40° relative to the loading axis. Further plastic deformation of the specimen concentrates on these shear bands, see Fig. 12b. A crack initiates at the center and grows along the shear band all the way to the lateral surface. This fracture mode is consistent with the experimental results, see Fig. 13. Compared with the cup–cone fracture mode, the slant fracture mode in the plane-strain specimen is relatively easily reconstructed. The fracture pattern can also be predicted using a relatively coarse mesh model with the minimum element size 0.1 mm  0.1 mm. Meanwhile, the predicted load–displacement curves are almost identical to that based on the finer mesh model. The conventional elastic–plastic material model cannot predict the formation of shear bands at a realistic level of plastic deformation, see Fig. 14. Diffuse necking continues to develop without any sign showing the

(a) Diffusion necking.

(b) Shear banding.

(c) Slant fracture.

Fig. 12. Numerically predicted deformation and fracture process of the transverse plane-strain specimen under tension.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2037

Fig. 13. Through-thickness slant fracture of the transverse plane-strain specimen.

Fig. 14. Numerically predicted diffuse necking of the plane-strain specimen using the elastic–plastic constitutive model uncoupled with damage. The stress–strain curve calibrated by Bao [20] was used in the calculation.

occurrence of shear bands. (But if the hardening rate is very low, shear bands would indeed occur. With a rigid-plastic model, Onat and Prager [50] developed a necking mode for plane-strain specimens, which is similar to that shown in Fig. 12c.) Fig. 15 shows a comparison of the load–displacement response between the test data and the present numerical simulations. It appears that the tensile load is slightly over-predicted. For transverse plane-strain specimens, there exists two transition zones near both sides, see [51]. The zones are dominated by the

Force (kN)

Uncoupled method

Test Coupled method Incipient necking Incipient shear banding

Relative displacement (mm) Fig. 15. Comparison of tensile loads vs. displacement between the numerical predictions and the test result. The relative displacement was measured using a 25.4 in. extensometer.

2038

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

plane-stress state of a lower stress. Hence, the total reaction force calculated from the two-dimensional planestrain finite element model would be over evaluated. Meanwhile, the present set of the damage parameters predicts a much larger fracture displacement than the experimental data, although the same set gives a satisfactory correlation for the case with the unnotched round bar. The incipience of necking corresponds to the peak load, which is marked as a circle in the load–displacement curve. At a ceratin point, diffuse necking shifts to through-thickness shear banding, which is marked as a triangle. A crack immediately follows the formation of the shear bands. This leads to the sharp drop of the load. By contrast, in the uncoupled procedure, the load gradually decreases with the diffuse necking deformation, since no fracture criterion was specified. 8. Discussions and conclusions A numerical algorithm for Lemaitre’s continuum damage mechanics has been formulated and implemented into ABAQUS through the user defined material subroutines. The true stress–strain curve and the damage parameters were calibrated for 2024-T351 aluminum alloy. The damage model was applied to simulate various fracture patterns encountered in common lab tests. The present coupled damage-plasticity model is able to capture the cup–cone fracture mode of the round bar and the slant fracture mode of the plane-strain specimen. The shear band preceding fracture provides the path for crack formation and propagation. It was found that these fracture patterns would be easily reconstructed with small elements and large critical damage values. The true stress–strain curve and damage parameters for 2024-T351 aluminum alloy have been calibrated from the tensile test on the smooth round bar. This set of material data gives a satisfactory correlation with the test results in terms of the load–displacement curves and the fracture patterns, but overestimates the critical displacements to fracture for other tensile tests on the notched round bars and the transverse plane-strain specimen. A comparison with the empirical Bao–Wierzbicki fracture locus indicates that the two-parameter Lemaitre’s evolution law of damage may be insufficient to describe the material ductility quantitatively. It remains to be explored whether other evolution rules of damage would improve the numerical prediction. An implicit–explicit numerical scheme was proposed within ABAQUS. The implicit code is first used to calculate the plastic deformation and the damage accumulation until fracture initiation, and the explicit code then takes over the rest of the simulation to predict crack formation and propagation. It has been demonstrated that this numerical scheme significantly reduces computational cost and at the same time provides the accurate solutions for the simulation of fracture tests under low-rate loading. The approach is also applicable to other industrial practices such as sheet-metal forming problems with fracture. A limited study on mesh size effects was performed. It was found that the predicted load–displacement response is not much sensitive to the element size. However, it would not be true if a wide range of element size was considered. Since an intrinsic length scale is lacking in the present coupled damage-plasticity model, damage softening would give rise to mesh size dependence of solutions. In the present paper, the same element size was defined in both the calibration and the prediction. The element size is thought of to represent a ceratin material length scale. Acknowledgements The author is indebted to Dr. Yingbin Bao for providing the experimental data. Thanks are also due to Professor Tomasz Wierzbicki, Dr. Liang Xue, and Mr. Yuanli Bai for many inspiring discussions. Appendix. Analytical solutions for damage development This section attempts to develop an analytical solution for the damage evolution of a solid under simple loading. With the simple closed-form solution, one will be able to quickly examine the correctness of the new algorithm and implementation. This is always desirable, especially for complex numerical procedures. Assume that the true stress–strain curve of a material is of the power-law form: ~ q ¼ A þ Brn ;

ðA:1Þ

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2039

where A is the initial yield stress; B is the hardening coefficient; and n is the material hardening exponent. By combining Eqs. (9), (11), and (12), one obtains ~2 q r_ : D_ ¼ Rv 2ES 0 1  D

ðA:2Þ

Substituting the effective von Mises stress into the above equation gives the expression for the damage accumulation rate ð1  DÞD_ ¼

Rv 2 ðA þ Brn Þ r_ ; 2ES 0

ðA:3Þ

where the expression for the stress triaxiality term Rv is given in Eq. (13). Note, that the introduction of dam~h =~ age would not change the stress triaxiality, i.e. r q ¼ rh =q. Under uniaxial tension and simple shear, we have, respectively Rv ¼ 1

with

rh 1 ¼ ; 3 q

ðA:4Þ

and 2 Rv ¼ ð1 þ mÞ 3

with

rh ¼ 0: q

ðA:5Þ

If Rv is kept constant during a loading process, both sides of Eq. (A.3) can be integrated. This yields the expression for the accumulated damage sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 2nþ1 Rv rnþ1 2 2 r þB ; ðA:6Þ D¼1 1 A r þ 2AB ES 0 nþ1 2n þ 1 where the initial condition D = r = 0 at t = 0 has been used. Substituting the above equation into Eq. (9), one is able to build the relationship between the accumulated plastic strain p and the internal variable r:  ffi Z r Z r sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nþ1 2nþ1 1 Rv r r ¼ þ B2 dr: ðA:7Þ 1 ð1  DÞ dr ¼ A2 r þ 2AB p ES 0 nþ1 2n þ 1 0 0 Similar expressions for the linear hardening and Ramberg–Osgood hardening material models were developed, respectively, by Doghri [23] and Lemaitre [9].

Theoretical solution Accumulated damage, D

ABAQUS/Standard+UMAT ABAQUS/Explicit+VUMAT

Accumulated plastic strain, p Fig. A.1. Comparison of the damage accumulation between the theoretical solution and the numerical results for one single element under simple shear.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

Effective equivalent stress, (MPa)

2040

Theoretical solution ABAQUS/Standard+UMAT ABAQUS/Explicit+VUMAT

Accumulated plastic strain, p Fig. A.2. Comparison of the accumulated plastic strain between the theoretical solution and the numerical results for one single element under simple shear.

To verify the implemented damage model, a simple finite element model with only one cubic, reduced-integration, solid element (C3D8R) was built. Four loading conditions were considered including uniaxial tension, uniaxial compression, simple shear, and combined loading. The element is assumed to be made of 2024-T351 aluminum alloy with the following material parameters: E = 72.4 GPa, m = 0.33, S0 = 6.0 MPa, A = 352 MPa, B = 853.4 MPa, and n = 0.66. A comparison between the theoretical solutions and the finite element results is given in Figs. A.1 and A.2 for the case with simple shear. It appears that the numerical predictions based on the newly implemented damage subroutines have excellent agreement with the theoretical solutions. This validates the correctness of the numerical algorithm and the programming implementation. References [1] Bonora N, Gentile D, Pirondi A, Newaz G. Ductile damage evolution under triaxial state of stress: theory and experiments. Int J Plast 2005;21:981–1007. [2] Gurson A. Continuum theory of ductile rupture by void nucleation and growth: part I—yield criteria and flow rules for porous ductile media. J Engng Mater Technol 1977;99:2–15. [3] Rousselier G. Ductile fracture models and their potential in local approach of fracture. Nucl Engng Des 1987;105:97–111. [4] Tvergaard V, Needleman A. Analysis of the cup–cone fracture in a round tensile bar. Acta Metall 1984;32:157–69. [5] Besson J, Steglich D, Brocks W. Modeling of crack growth in round bars and plane strain specimens. Int J Solids Struct 2001;38:8259–84. [6] Kachanov LM. Introduction to continuum damage mechanics. Boston, Dordrecht: Martinus Nijhoff Publisher; 1986. [7] Lemaitre J. A course on damage mechanics. New York: Springer-Verlag; 1992. [8] Lemaitre J, Chaboche J-L. Mechanics of solid materials. Cambridge: Cambridge University Press; 1990. [9] Lemaitre J. A continuous damage mechanics model for ductile fracture. J Engng Mater Technol 1985;107:83–9. [10] Chaboche JL. Continuum damage mechanics: part I—general concepts. J Appl Mech 1988;55:59–64. [11] Chaboche JL. Continuum damage mechanics: part II—damage growth, crack initiation, and crack growth. J Appl Mech 1988;55:65–72. [12] Wang TJ. Unified CDM model and local criterion for ductile fracture—I. Unified CDM model for ductile fracture. Engng Fract Mech 1992;42:177–83. [13] Bonora N. A nonlinear CDM model for ductile failure. Engng Fract Mech 1997;58:11–28. [14] Hambli R. Finite element model fracture prediction during sheet-metal blanking processes. Engng Fract Mech 2001;68:365–78. [15] Børvik T, Hopperstad OS, Berstad T, Langseth ML. A computational model of viscoplasticity and ductile damage for impact and penetration. Eur J Mech A—Solids 2001;20:685–712. [16] Johnson GR, Cook WH. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Engng Fract Mech 1985;21(1):31–48. [17] Børvik T, Hopperstad OS, Berstad T, Langseth M. Perforation of 12 mm thick steel plates by 20 mm diameter projectiles with flat, hemispherical and conical noses. Part II: numerical simulations. Int J Impact Engng 2002;27(1):37–64. [18] Andrade Pires FM, de Souza Neto EA, Owen DRJ. On the finite element prediction of damage growth and fracture initiation in finitely deforming ductile materials. Comput Methods Appl Mech Engng 2004;193:5223–56.

X. Teng / Engineering Fracture Mechanics 75 (2008) 2020–2041

2041

[19] Chaboche JL, Boudifa M, Saanouni K. A CDM approach of ductile damage with plastic compressibility. Int J Fract 2006;137:51–75. [20] Bao Y. Prediction of ductile crack formation in uncracked bodies. PhD thesis, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA, 2003. [21] Chow CL, Wang J. An anisotropic theory of continuum damage mechanics for ductile fracture. Engng Fract Mech 1987;27(5):547–58. [22] Lemaitre J, Desmorat R, Sauzay M. Anisotropic damage law of evolution. Eur J Mech A—Solids 2000;19(2):187–208. [23] Doghri I. Numerical implementation and analysis of a class of metal plasticity models coupled with ductile damage. Int J Numer Methods Engng 1995;38:3403–31. [24] de Souza Neto EA. A fast, one-equation integration algorithm for the Lemaitre ductile damage model. Commun Numer Methods Engng 2002;18:541–54. [25] Doghri I, Billardon R. Investigation of localization due to damage in elasto-plastic materials. Mech Mater 1995;19:129–49. [26] ABAQUS version 6.5 documentation. ABAQUS, Inc., Providence, Rhode Island, USA [27] Hambli R, Badie-Levet D. Damage and fracture simulation during the extrusion processes. Comput Methods Appl Mech Engng 2000;186:109–20. [28] Bao Y, Wierzbicki T. On fracture locus in the equivalent strain and stress triaxiality space. Int J Mech Sci 2004;46(1):81–98. [29] McClintock FA. A criterion for ductile fracture by the growth of holes. J Appl Mech 1968;35:363–71. [30] Rice JR, Tracey DM. On the ductile enlargement of voids in triaxial stress fields. J Mech Phys Solids 1969;17:201–17. [31] Hancock JW, Mackenzie AC. On the mechanisms of ductile failure in high strength steels subjected to multi-axial stress-states. J Mech Phys Solids 1976;24:147–69. [32] McClintock FA. Plasticity aspects of fracture. In: Liebowitz H, editor. Fracture: an advanced treatise, vol. III. New York: Academic Press; 1971. [33] Rogers HC. The tensile fracture of ductile metals. Trans Metall Soc AIME 1960;218:498–506. [34] Bluhm JI, Morrissey RJ. Fracture in a tensile specimen. In: Yokobori T, Kawasaki T, Swedlow JL, editors. Proceedings of the first conference on fracture, Sendai, Japan, September 1965, vol. 3. p. 1739–80. [35] Scheider I, Brocks W. Simulation of cup–cone fracture using the cohesive model. Engng Fract Mech 2003;70:1943–61. [36] Xue L. Damage and fracture initiation of uncracked ductile solids under triaxial loading. Int J Solids Struct 2007;44:5163–88. [37] Roychowdhury S, Roy YD, Dodds RH. Ductile tearing in thin aluminum panels: experiments and analyses using large-displacement, 3-d surface cohesive elements. Engng Fract Mech 2002;69:983–1002. [38] Wierzbicki T, Xue L. On the effect of the third invariant of the stress deviator on ductile fracture. Technical report no. 138. Impact and Crashworthiness Lab, MIT, Cambridge, MA, 2005. [39] Besson J, Steglich D, Brocks W. Modeling of plane strain ductile rupture. Int J Plast 2003;19:1514–7. [40] Anand L, Spitzig WA. Initiation of localized shear bands in plane strain. J Mech Phys Solids 1980;28:113–28. [41] Spencer K, Corbin SF, Lloyd DJ. The influence of iron content on the plane strain fracture behaviour of AA 5754 Al–Mg sheet alloys. Mater Sci Engng A 2002;325:394–404. [42] Hill R. Acceleration waves in solids. J Mech Phys Solids 1962;10:1–16. [43] Rudnicki JW, Rice JR. Conditions for the localization of deformation in pressure-sensitive dilatant materials. J Mech Phys Solids 1975;23:371–94. [44] Rice JR. The localization of plastic deformation. In: Koiter WT, editor. Theoretical and applied mechanics, Proceedings of the 14th international conference on theoretical and applied mechanics. Amsterdam: North-Holland; 1976. p. 207–20. [45] Hutchinson JW, Tvergaard V. Shear band formation in plane strain. Int J Solids Struct 1981;17:451–70. [46] Tvergaard V, Needleman A, Lo KK. Flow localization in the plane strain tensile test. J Mech Phys Solids 1981;29:115–42. [47] Mahgoub E, Deng X, Sutton M. Three-dimensional stress and deformation fields around flat and slant cracks under remote mode i loading conditions. Engng Fract Mech 2003;70(18):2527–42. [48] Inal K, Wu PD, Neale KW. Instability and localized deformation in polycrystalline solids under plane-strain tension. Int J Solids Struct 2002;39:983–1002. [49] Clausing DP. Effect of plastic strain state on ductility and toughness. Int J Fract Mech 1970;6:71–85. [50] Onat ET, Prager W. The necking of a tension specimen in plane plastic flow. J Appl Phys 1954;25:491–3. [51] McClintock FA, Zheng ZM. Ductile fracture in sheets under transverse strain gradients. Int J Fract 1993;64:321–37.