Numerical simulation of low stress triaxiality ductile fracture

Computers and Structures 84 (2006) 1641–1650 www.elsevier.com/locate/compstruc

Numerical simulation of low stress triaxiality ductile fracture T. Pardoen

*

De´partement des Sciences des Mate´riaux et des Proce´de´s, Universite´ Catholique de Louvain, IMAP, Place Sainte Barbe 2, B-1348 Louvain-la-Neuve, Belgium Received 17 October 2005; accepted 3 May 2006 Available online 7 July 2006

Abstract One of the main shortcoming of existing damage models when applied to the simulation of metal forming operations is their limited validity under low stress triaxiality conditions. An extended Gurson model incorporating the effect of void shape and relative void spacing on the growth and coalescence has been developed in order to encompass both low and large stress triaxiality regimes. The constitutive model has been implemented into an implicit finite element code within a finite strain set up. Identification procedures from experimental data are proposed and illustrated in an application on copper bars exhibiting different strain hardening capacity. A parametric analysis of the damage evolution in specimens deformed under uniaxial tension with necking demonstrates the importance of properly accounting for void shape and void coalescence in the case of large strains problems.  2006 Elsevier Ltd. All rights reserved. Keywords: Solid mechanics; Damage; Plasticity; Anisotropy; Fracture; Void growth

1. Introduction Metal forming operations are characterized by low stress triaxiality which allows large plastic deformation with minimum damage accumulation. On the one hand, when the factor controlling material formability is the appearance of localized plastic bands, a very sophisticated damage model is not mandatory for the prediction of failure. Indeed, the coupling between the condition for localization and the damage evolution is relatively weak because the porosity levels are usually very small in most industrial alloys. Thus, detailed considerations related for instance to the void morphology are, in that case, un-important. Furthermore, the prediction of fracture strain after localization has set in is usually not a key issue because the practical failure criterion is set by the onset of plastic localization. On the other hand, in several forming operations such as extrusion, drawing, bending, hole expansion and many others, failure can be primarily caused by the growth and coalescence of voids in the bulk or near the surface of the

deforming solid. Simulating, in those applications, the damage evolution during plastic deformation requires physically realistic damage models. Many studies have revealed a significant variation of the void aspect ratio in the low stress triaxiality regime1: whatever their initial shape, voids elongate in the direction of maximum strain [1,2]. This change of aspect ratio not only affects the void growth rate with respect to the reference rate of growth of a spherical void, but affects also the void coalescence process which determines the fracture strain [2–8]. Furthermore, non-spherical voids ranging from penny-shape cracks to highly prolate voids are common in industrial metal alloys. A primary rolling process during which globular second phases experience severe deformations leads to a preferential orientation of the second phases and to high aspect ratio. Proper modelling of a secondary forming process has to account for this initial anisotropy as well. Many forming operations also involve large amount of shear. Voids then not only elongate but rotate

1

*

Tel.: +32 10 472417; fax: +32 10 474028. E-mail address: [email protected]

0045-7949/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.05.001

The stress triaxiality is defined by the ratio of the hydrostatic stress divided by the equivalent stress. A stress triaxiality is considered to be low when it is typically smaller than 1.

1642

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

during deformation which might also affect the condition for void coalescence and fracture. The complex behaviour associated to the void orientation, shape and rotation requires to incorporate many adjustable parameters in phenomenological models in order to insure their validity under general stress state conditions, especially in the low stress triaxiality regime typical of forming operations. With the increasing number of tuning parameters, the identification procedures become more and more complex as well as the transfer from one material to another. These difficulties have prompted many research groups to work on advanced micromechanics-based constitutive models (e.g. [2,4–9]). In this work, an enhanced Gurson model [10] for dilatant plastic materials based on the contributions by Gologanu et al. [4] and Thomason [3] has been developed in order to encompass both low and large triaxiality regimes within the same model [2] (see also [7]). The advantages of such a model are (i) (in principle) a limited number of adjustable parameters; (ii) a better connection between microstructure (and thus the materials processing history) and mechanical design. The disadvantage is an increasing complexity in the numerical integration and implementation of the model. Ideally, the physical parameters of such a model have to be determined experimentally by microstructural characterization. Electron microscopes and microstructure image analysis tools are obviously not always available. In that case, the physical parameters can still be identified based on macroscopic mechanical tests using inverse procedures as for phenomenological damage models.2 The paper starts with a presentation of the constitutive model and of its implementation. The identification procedures are discussed in a second section and applied, for the sake of illustration, to the damage and fracture of copper bars. This application on copper offers an interesting validation of the model as the strain hardening capacity of this material can be modified by thermal treatment without affecting the parameters related to the damage evolution. The last section is devoted to a computational parameter study of ductile fracture during uniaxial tensile tests with necking in order to illustrate the need for advanced constitutive models in the low triaxiality regime. The current challenges involved in the development and in the transfer of micromechanically-based damage model to metal forming are outlined in the conclusions.

cence of spheroidal voids3on the behaviour of the material has been presented by Pardoen and Hutchinson [2,11]. The response of the matrix is assumed to obey J2 flow theory. Plastic anisotropy of the matrix can be introduced following the framework proposed by Benzerga et al. [7]. The version of the model presented here considers the voids present from the beginning of the deformation history which is a good approximation for many materials. Void nucleation models can be readily incorporated into the present constitutive model (see for instance [8,12,13]). The model is based on two solutions for the expansion of a void in an elastoplastic material: one solution is called ‘‘void growth’’ corresponding to diffuse plasticity in the matrix around the void, and the other is called ‘‘void coalescence’’ corresponding to localized plasticity in the ligament between the growing voids. These two solutions are cast in the form of two distinct plastic yield surfaces supplemented by the normality rule for the plastic strain increments and evolution laws for the internal variables: porosity f, void aspect ratio W (or S  ln W), relative void spacing v (void diameter divided by void spacing), and mean yield stress of the matrix material ry: Yield surfaces Growth: Ugrowth

 g C 0 r 2 g  2 kr þ grh Xk þ 2qðg þ 1Þðg þ f Þ cosh j h ry ry  ðg þ 1Þ2  q2 ðg þ f Þ2 ¼ 0

ð1Þ

Coalescence: kr0 k 3 jrh j þ  F ðW ; vÞ ¼ 0 ry 2 ry sffiffiffi# "  2 3 1  v 1 F ðW ; vÞ ¼ ð1  v2 Þ a þb 2 vW v

Ucoalescence 

ð2Þ

with

ð3Þ

Evolution law for f Growth and coalescence: f_ ¼ ð1  f Þ_epkk

ð4Þ

Evolution law for W (or S = ln W) 2. An enhanced model for void growth and coalescence Growth : An extension of the Gurson model [10] based on the contributions by Gologanu et al. [4] and Thomason [3] and accounting for the effect of the growth and coales-

2 In principle, a physical model should behave at least as well as a phenomenological model involving the same number of fitting parameters.

  3 e_ p S_ ¼ ð1 þ h1 Þ e_ p  kk d : P þ h2 e_ pkk 2 3

ð5Þ

3 A spheroid is an ellipsoid with a symmetry of revolution around one axis.

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

Coalescence: 

W ð2v2  1Þ 3 p 1 W_ ¼ e_ : P  e_ pkk 2f 2 2

 ð6Þ

Evolution laws for ry Growth: ry e_ py ð1  f Þ ¼ r : e_ p

ð7Þ

Coalescence: r_ y ¼

ory v2 p e_ oepy f e

ð8Þ

Evolution for v Coalescence:   vð3  2v2 Þ 3 p 1 p v_ ¼ e_ : P  e_ kk 6f 2 2

ð9Þ

In Eqs. (1)–(9), • The first yield surface Ugrowth is a Gurson-type yield surface derived by Gologanu–Leblond–Devaux for spheroidal voids [4] and extended to strain hardening materials [2]. The yield function Ucoalescence has been developed in the spirit of the work by Thomason [3], extended to the full coalescence response and to strain hardening materials. • r is the Cauchy stress tensor; r 0 is the deviatoric part of the Cauchy stress tensor; rgh is a generalized hydrostatic stress defined by rgh ¼ r : J. • X is a tensor associated to the void axis and defined by 2/3ez  ez  1/3ex  ex  1/3ey  ey(ez is the base vector parallel to the cavity axis); J is a tensor associated to the void axis and defined by (1  2a2)ez  ez + a2ex  ex + a2ey  ey; P is a projector tensor, defined by ez  ez; d is the Kronecker tensor. • The symbol k k represents the von Mises norm. • Gologanu et al. [4] derived analytical relationships for the dummy parameters C, g, g, j, h1, h2, a2 in terms of the state variables W and f (see [2] for full expression). • The parameter q has been calibrated as a function of f0, of W0, and of the mean strain hardening capacity of the material (exponent n in Eq. (16)) (see [2] for complete expressions); a has been calibrated as a function of the mean strain hardening index n: a = 0.1 + 0.217n + 4.83n2, see [2]; b is equal to 1.24. • The plastic strain rate tensor e_ p is defined in the reference configuration; e_ pe is the equivalent plastic strain rate (or accumulated plastic strain rate); e_ py is the equivalent plastic strain rate of the matrix material. The effective plastic strain of the matrix material epy is related to the current yield stress of the matrix material ry through the material true stress true strain curve, see later Eq. (16).

1643

The response of the material under loading is first elastic while the stress state lies within the elastic domain of the two yield surfaces. In the beginning of the deformation, the voids are small and the spacing is large resulting in a diffuse mode of plastic deformation. Except for pathological cases, the first yield surface to be reached is thus Ugrowth. With increasing deformation Ugrowth first expands due to strain hardening and then contracts due to damage softening. Void growth and ligament reduction also induce contraction of Ucoalescence. When the two yield surfaces intersect at the current loading point, transition to coalescence occurs. With increasing deformation the coalescence yield surface rapidly contracts towards the zero stress state. This model gives predictions in good agreement with finite element void cell calculations [2]. Note that the parameter v does not affect the void growth response. It is calculated during the void growth phase only to test the condition of coalescence using the following geometrical relationship: v¼

fk cg W

!13 ð10Þ

where k is equal to 2Lz/(Lx + Ly) and Lx, Ly and Lz are the mean void spacings in the three directions corresponding to the three void axis (whose evolution can directly be obtained from the strains in the corresponding directions), and cg is a geometric factor which depends on the void distribution. In the case of periodic packing p offfiffiffi voids, cg = p/6 = 0.523 for a simple cubic array, cg ¼ 3p=9 ¼ 0:605 for an hexagonal distribution and cg = 2/3 = 0.666 for a void surrounded by a cylindrical matrix. The model is completed by the normality rule for the plastic strain increment and an evolution law for the void axis e_ pij ¼ c_

oU orij

e_ z ¼ Xez

ð11Þ ð12Þ

where the spin tensor X is taken here as the material spin tensor. The voids are thus assumed to rotate with the material. This is the assumption also proposed by Gologanu [14]. A more accurate void rotation law, motivated by 3D finite element void cell calculations, as been recently developed (see [15] where a modified version of the coalescence model is also introduced). Isotropic elasticity is assumed while using the Eshelby [16]–Mori–Tanaka [17] homogenization method to approximately couple the elastic constants to the porosity (spherical voids approximation) K¼ G¼

4ð1  f ÞK 0 G0 4G0 þ 3fK 0 ð1  f ÞG0 0 þ12G0 Þ 1 þ f ð6K 9K 0 þ8G0

ð13Þ ð14Þ

1644

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

where G and K are the overall shear and bulk moduli, respectively and G0 and K0 are the shear and bulk moduli of the undamaged material, respectively. Finally, the following uniaxial response has been chosen r Ee ¼ when r < r0 ; r0 r0  n r Eep ¼ 1þ when r > r0 : r0 r0

ð15Þ ð16Þ

The constitutive model has been implemented in the finite element program ‘‘ABAQUS Standard’’ through a User defined MATerial subroutine [18]. The key aspects of the implementation are detailed hereafter by presenting first the integration of the void growth model and then the integration of the void coalescence model. 2.1. Integration of the void growth model The void growth model has been implemented using a fully implicit integration scheme based on the return mapping algorithm [19], see also [20] for an implementation of the Gurson model. The set of 10 first-order differential equations (i.e. Eqs. (1), (4), (5), (7) and the 6 equations (11); use has been made of the consistency condition _ ¼ 0 to determine the plastic multiplier appearing in U (11)) involved in the void growth model are solved using the Newton–Raphson method after proper linearization. The frame dependent internal variables appearing in the model (e.g. in (6)) are properly updated using the rotation subroutine provided by Abaqus, called ‘‘ROTSIG’’ [18]. A total formulation approach [21–23] would lead to more accurate results at very large strain, but requires extra complexity in the writing of the UMAT. Because of the complexity of the model, the theoretical tangent modulus is used as an approximation of the consistent tangent operator (see [24] for a recent work in which the consistent tangent operator has been worked out for the present model). The tangent moduli Lijkl are found by specialising the algorithm proposed by Hill [25] such that r_ ij ¼ Lijkl e_ kl

ð17Þ

with Lijkl ¼ C ijkl 

mij mkl A þ lmn mmn

ð18Þ

where the Cijkl are the elastic moduli, and the m’s, l’s and A are given by (here U  Ugrowth) lij ¼

oU ; orij

ð19Þ

mij ¼ C ijkl lkl ; ð20Þ 2      ory epy rij oU oU 3 oU 1 oU A ¼ 4 þ Þ  d ð1 þ h 1 ij P ij oepy ry ð1  f Þ orij ory 2 orij 3 orkk 3  oU oU oU oU5 þh2 þ ð1  f Þ : orkk oS orkk of

ð21Þ

Because of the use of the theoretical rather than the consistent tangent operator, quadratic convergence is not achieved. The model is highly nonlinear and only small deformation increments can be used with the Newton– Raphson scheme (typically the maximum local strain increments are about 0.001–0.003). A sub-stepping method has been implemented in order to authorize the local integration of the model while keeping the global time step large enough. 2.2. Integration of the void coalescence model The coalescence response being only controlled by the hydrostatic, rh, and the deviatoric, re  kr 0 k, components of the stress tensor (see Eq. (2)), the Aravas scheme [26] has been applied for the integration of the coalescence model in order to reduce the system of equations to be solved by the Newton–Raphson method. The yield surface has a corner at rh = 0. This corner at rh = 0 has been rounded in order to allow computation of the normal to the yield surface. This rounding is purely heuristic. For jrh/jrej < Tc, the yield surface writes Ucoalescence corner  r2h þ ðkr0 k þ ry F ðW ;vÞA1 Þ2 þ A2 ðkr0 kF Þ2 ¼ 0 ð22Þ

with A1 ¼

4T c  6 6 þ 9T

and

A2 ¼ 

13 T 2c

9 1 þ 3Tc 2 2

ð23Þ

where Tc is a numerical parameter that should be chosen small enough for not perturbing the response of the model. Numerical problems sometimes occur at the onset of coalescence due to the corner at the intersection of the yield surfaces leading to an abrupt change in strain rate direction. The very unlikely event wherein coalescence is interrupted with the resumption of diffuse void growth is not allowed in the current implementation. When the relative void spacing becomes close to one, say v = 0.99, the remaining stresses are ramped down to zero. Note that, for some materials, a second population of voids tend to significantly accelerate the coalescence process. In that case, a critical void spacing ratio, vc can be introduced, above which the element looses all its stress carrying capacity [11]. Another more advanced proposal has been recently worked out to account for the presence of a second population of voids [27]. 3. Identification and validation of the model 3.1. Identification procedures The parameters of the model which have to be identified are : • the parameters related to the initial void distribution f0, W 0, v 0;

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

  • the effective stress effective strain curve ry epy of the undamaged material; • the elastic constants (usually known); • A critical relative void spacing vc if necessary (i.e. if the material under consideration shows a premature ligament fracture during the coalescence process. Proper modeling of the energy dissipated during the coalescence process can be essential at large stress triaxiality where it might account for a significant part of the energy dissipated [2,11]. At low stress triaxiality, the energy associated to the final ligament fracture process is always negligible when compared to the total work spent to reach that point). Two identification procedures are discussed now: the first assumes that standard characterization techniques are available to determine the parameters related to the microstructure; the second only requires mechanical tests. Of course, both procedures can be combined in different manners. 3.1.1. Identification based on microstructure characterization • f0 – Voids almost always nucleate on second phase particles (inclusions, intermetallic particles, oxides, etc.). The first step is to determine which particles give rise to voids, either by interface decohesion or by fracture, as well as the fraction n of particles that nucleate voids. A simple method consists of analyzing fracture surfaces along transverse sections (after proper surface preparation). Indeed, the material just underneath the fracture surfaces is representative of the state of deformation and damage just before fracture. The volume fraction of particles fp is measured on un-deformed, metallographically prepared samples using image analysis procedures (involving Delaunay tesselation and Voronoi triangulation [28–30]). If voids nucleate by particle decohesion, then f0 = nfp. If voids nucleate by particle fracture, then the initial void shape is very flat. It has been shown that the value of the initial void shape W0 has no effect on the damage process as long as it is smaller than 0.01 and provided that f0 is scaled with W0, i.e. f0 = nfpW0 [8,12]. • W0 – If voids nucleate mechanism is by particle decohesion, then W0  Wp, where Wp is the mean void aspect ratio. Not all particle are spheroidal. One needs thus to approximate the real particle morphology by spheroids, for instance by defining an effective diameter from the void area projected in the equatorial plane. If nucleation occurs through particle fracture, the initial void shape, as explained above, is unimportant as long as is taken  smaller than typically 0.01 and f0 = nfpW0. • ry epy curve – If f0 is smaller than typically 0.1–0.2%, the effect of the porosity on the global stress–strain curve is usually very small even close to void coalescence (the porosity at void coalescence does usually not exceed

1645

1–2% for initial porosity around 0.1–0.2%). If f0 is larger than about 0.1–0.2%, then an inverse identification procedureis  recommended in order to accurately determine the ry epy curve. The idea is to simulate the test using the damage model and to find the correction to apply on the global stress–strain curve (see [31] for an example). Such procedure also delivers the flow response   after necking. Another option is to infer the ry epy curve from a torsion test in which damage can usually be neglected and in which no necking takes place. • v0 – Because particles and second phases do not distribute periodically, there is no way to estimate this parameter univocally. For most materials, the relevant v0 will be smaller than the mean relative particle spacing. Indeed, a void will always tend to coalesce, in some way, with its nearest neighbour. In a real material, the criterion will involve the distance between the voids but also the angle between the main loading direction and the direction relating the two voids. There is thus no simple way to estimate v0. This parameter will thus be considered here as the only adjustable parameter of the model that must be identified from macroscopic mechanical tests. Nevertheless, even if v0 will account for all the approximations involved in the model, it is bounded to the physics and one can, afterwards, assess whether the identified value is realistic or not, hence if the model is valid or not (see next section for the application on copper). 3.1.2. Identification from mechanical tests When characterization methods are not available or that a faster identification procedure is needed or that the identification based on microstructure characterization is not trusted, it is always possible to rely only on mechanical tests in order to adjust the parameters of the model. Literature survey and discussions with material scientists might be helpful to define a sensible first guess for f0 (e.g. do we expect a very clean material? remembering that in most industrial metals f0 is lower than 1%) and also about W0 (e.g. a primary rolling operation tends to elongate the particles). In order to supplement the classical uniaxial tests (which are difficult to interpret quantitatively, see further in Section 4), it is recommended to use specimen geometries that cover various stress triaxiality levels. Indeed, the stress triaxiality is the main factor controlling the damage evolution. The most attractive geometry is the cylindrical notched round bar specimen, as initially proposed by Hancock and McKenzie [32] and promoted, since the early 80s by A. Pineau and co-workers, e.g. [13,33,34]. By changing the notch radius, the stress triaxiality can be varied between typically 0.6 and 2. Furthermore, the stress triaxiality in such a geometry remains almost constant all along plastic deformation in the most damaged region which is located at the centre of the minimum section of the specimen (except if the notch is too sharp, then the most damaged region moves close to the tip). A robust set of param notch  eters f0, W0, v0, and ry epy can be identified from a single

1646

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

uniaxial tensile test and two tests on notched round bars with two different radii of curvature.

500 As-received Cu as-received

400 σ0

=312MPa

3.2. Application to copper 3.2.1. Identification of the parameters The first identification procedure described in the previous section has been applied on cold drawn copper bars through the nucleation of voids by interface decohesion between the copper matrix and copper oxide particles. The voids then grow by plastic deformation and finally coalesce by internal necking until final void impingement. The key microstructural parameters for modeling damage in this material are thus the parameters related to the inclusions: the volume fraction, the shape and the distribution (see [30] for details). The particle volume fraction fp is equal to about 0.3%. Microscopic observations have shown that all particles give rise to voids and thus f0  0.3%. The initial particle shape is equal to 1.6. Because voids nucleate by interface cracking, the initial void shape W0 can thus be taken as equal to 1.6. The characterization of the inclusion spacings provided a mean distance between neighbours equal to 16 lm and a mean distance between nearest neighbours equal to 6 lm.4 The mean diameter of the particles is equal to 1.3 lm. The physical bounds for v0 are thus between 1.3/6 = 0.22 and 1.3/16 = 0.08. The as-received copper bars were heavily strain hardened due to the drawing process. Part of the bars were annealed in order to recover the large strain hardening capacity of copper without affecting the shape, volume fraction and size of the particles. This system offers thus a means for assessing the ability of the model to capture the effect of the strain hardening on the damage evolution. An iterative was followed in order to identify  methods  the curve ry epy after the onset of necking, taking into account the damage evolution (see [31] for details). The uniaxial stress–strain curves are compared in Fig. 1 showing the large difference in strain hardening between the two materials for strains smaller than 0.2. Uniaxial tensile tests were also performed on notched round bars. Different radius of curvature were used in order to cover a wide range of stress triaxiality and to offer a large basis for critical assessment of the model. The specimens are noted Rx with x being equal to the radius of curvature in mm (the ligament diameter is equal to 4 mm and the specimen diameter is equal to 9 mm). Four notch radii were tested: R1, R2, R4, R8. Table 1 summarizes the fracture strains obtained for both materials and for the different geometries.

4

Each particle has several neighbours, defined using the Delaunay tessalation, and among these neighbours, the particle has one nearestneighbour. These two values are then averaged for the set of particles.

Annealed Cu

σ (MPa)

300

200

100 σ

annealed =40MPa

0

0

0

0.2

0.4

0.6

0.8

1

1.2

ε Fig. 1. True stress–true strain curves for the as-received and annealed copper bars. These curves were obtained from uniaxial tensile tests using an iterative inverse identification procedure after the onset of necking [25].

Table 1 Experimental fracture strains for the as-received and annealed copper bars and each specimen geometry

As-received Cu Annealed Cu

Uniaxial test (with necking)

R8

R4

R2

R1

0.95 1.27

0.64 1.06

0.505 0.96

0.27 0.865

0.12 0.755

3.2.2. Validation The uniaxial tensile tests on the notched and smooth copper bars were simulated using the finite element code ABAQUS Standard and the ‘‘User defined MATerial’’ option to implement the constitutive model described in the previous section. Care has been taken to insure that the meshes were sufficiently refined and that the results were independent of the degree of refinement (convergence was checked not only for the global stress–strain curve but also for the local quantities such as the stress triaxiality history in the most loaded element in which fracture initiates). Axisymmetric elements were employed while imposing, in the case of the un-notched specimens, a slight reduction of the diameter in order to trigger the necking. A first set of calculations on un-notched specimen was run for the as-received material using the parameters identified above and different values v0. The value of v0 = 0.25 leads to a fracture strain equal to the experimental fracture strain, i.e. 0.95. This value is slightly outside the physical bounds given in the previous section (where the largest value was equal to 0.22). Note that it is expected that voids will tend to coalescence with the nearest neighbours and it makes thus sense to be close to the upper bound. Note also that because the model assumes a perfect coalescence mode until void impingement and no other accelerating factors, it always tends to overestimate the fracture strain. The other tests on the notched specimens as well as the tests on the annealed material were then simulated using the values v0 = 0.25, f0 = 0.003 and W0 = 1.6 and the stress–strain curves shown in Fig. 1. Fig. 2 shows the com-

T. Pardoen / Computers and Structures 84 (2006) 1641–1650 1.5

1647

2 FE results

Δεef = 0.5

1.5 1

εef

εef

Δn = 0.15

1

Annealed 0.5

0.5

f0=0.2%, W0 =1, σ0 /E=0.001

As-received Cu 0

0 0.4

0.6

0.8

1

1.2

1.4

1.6

0

0.0 5

0.1

Fig. 2. Comparison between the experimental and predicted fracture strains as a function of the average stress triaxiality in the most damaged region of the specimens for both the as-received and annealed samples. The stress triaxiality increases with decreasing notch radius: from uniaxial to R8–R4 to R2–R1.

parison between the experimental and predicted fracture strains as a function of the mean stress triaxiality computed in the most loaded element (always in the centre of the minimum section area). The exponential variation of the fracture strain has been shown experimentally by Hancock and McKenzie [32] a long time ago and can be justified theoretically by simple void growth analysis, for instance with the Rice and Tracey model assuming a constant critical void radius [35]. The agreement between the experimental results and the simulations for the as-received materials is conspicuous. For the annealed material, the ductility is slightly overestimated at the lowest stress triaxiality and overestimated at the largest triaxiality. This result indicates that there is still room for improvement, probably at the level of the void growth flow potential, for large strain hardening exponents. Note that there are also some uncertainty on the material behaviour in the case of the annealed samples coming from variations of hardness and flow properties along the diameter of the bar. 4. Ductility in uniaxial tension – a parameter study Modelling uniaxial tension with necking involves typical features of forming processes: finite strains effects, plastic localization, non-radial loadings, low and evolving stress triaxiality. The modelling of the uniaxial loading of cylindrical specimens is thus relevant for assessing the potential of an advanced constitutive model to simulate damage evolution during, for instance, deep drawing, extrusion or rolling. Uniaxial tensile tests on cylindrical specimens is also the most common type of mechanical test for quantifying the ductility of materials. First, three simulations are reported for three ideal materials characterized by three different strain hardening exponents n = 0.05, 0.1, and 0.2, but otherwise identical properties: ry/E = 0.001, m = 0.3, f0 = 0.002, W0 = 1 assuming a uniform distribution (i.e. k0 = 1) with an hexag-

0.15

0.2

0.25

n

Mean stress triaxiality in most damaged region

Fig. 3. Variation of the ductility as a function of the strain hardening exponent.

1.2

f0=0.2%, W0 =1, σ0 /E=0.001

n = 0.1 n = 0.05

0.9 σh/σe

n = 0.2

0.6

1/3 = pure uniaxial tension

0.3

0

0

0.5

εe

1

1.5

Fig. 4. Variation of the stress triaxiality in the most damaged element (i.e. the centre of the minimum cross-section area) as a function of the average effective strain in the minimum section.

onal arrangement (giving v0 = 0.149). Fig. 3 shows the variation of the ductility, defined as the strain at cracking initiation (when full coalescence has occurred in the most loaded element), as a function of the strain hardening exponent. As expected, the ductility significantly increases with strain hardening. It is shown in Fig. 3 that this increase is larger than the change of ductility originating from the different strain at necking considering no void growth before necking (the strain at necking is almost equal to n). This result agrees also with the experiments on as-received (low strain hardening eu = 0.01 and eef = 0.96) and annealed (high strain hardening eu = 0.29 and eef = 1.3) copper described in the previous section giving Deu = 0.28 and Deef = 0.34, respectively.5 The main reason for this effect is the following. Fig. 4 shows that the stress triaxiality 5

Note that care must be taken to put to much emphasis on that comparison because the strain hardening exponent of copper changes significantly during plastic deformation – annealed and as-received Cu show closer strain hardening exponents at large strains when fracture occurs (see Fig. 1).

1648

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

increases more slowly with straining for larger n. This is a geometric finite strain effect resulting from that larger hardening capacity promotes more diffuse necks. A lower triaxiality tends to decrease the void growth rate as shown in Fig. 5 and to favor void elongation, as shown in Fig. 6. Another reason for this result is that, for identical stress triaxiality, larger n tend to decrease the void growth rate. This effect is less important than the indirect effect on the necking evolution. This discussion shows also that different fracture strains are expected when changing the specimen configuration for the tensile test: the evolution of necking and of the stress triaxiality inside the neck depends, for instance, on the width over thickness ratio in plate type specimens (e.g. [36]). These results demonstrate the complex couplings existing between the stress state (related to the applied loading configuration and to finite strain geometric effects), the flow properties (strain hardening) and the microstructure (shape, void volume fraction). Only an enhanced micromechanical model is able to properly capture these effects.

A parametric study has also been performed in order to derive trends for the ductility of metal alloys subjected again to uniaxial tension (with necking). Fig. 7 shows the variation of the ductility as a function of the initial porosity for various strain hardening exponents, keeping the initial void shape spherical. Most industrial alloys exhibit strain hardening exponent ranging between n = 0.1 and 0.3. The effect of the strain hardening is more marked at larger porosities. Indeed, the relative contribution of the straining before necking becomes more and more important when the fracture strain decreases. Fig. 8 shows the variation of the ductility as a function of the initial porosity for different initial void shapes ranging from very flat, i.e. W0 = 0.01 to very elongated, i.e. W0 = 100, with the strain hardening exponent always equal to 0.1. The effect of the initial void shape is important and explains the anisotropy in ductility observed when testing rolled plates which involves elongated second phases with a preferential orientation along the rolling direction.

3 W0 = 1, σ0 /E = 0.003

n = 0.3

0.075

f0=0.2%, W0=1, σ0/E=0.001

2.5

εef

n = 0.1 0.05

homogenous distribution

n = 0.2 n = 0.1

2

n = 0.05

1.5

f

n = 0.2

1

0.025

0.5 10 0

0

0.5

1

1.5

εe

-7

10

-6

10

-5

-4

10

-2

10

-1

n = 0.1 σy /E = 0.003

W0 =100

homogenous distribution

f0=0.2%, W0 =1, σ0 /E=0.001 n = 0.2

-3

Fig. 7. Variation of the ductility (fracture strain) as a function of the initial porosity for different strain hardening exponent.

3 2.5

10

f0

2

Fig. 5. Variation of the void volume fraction in the most damaged element (i.e. the centre of the minimum cross-section area) as a function of the average effective strain in the minimum section.

10

W0 =6

2.5

2

n = 0.1 1.5

ln(W)

ε ef

n = 0.05

W0 =1

2

W0 =1/6

1.5 W0 =0.01

1

1 0.5

0.5 0

0

0.5

1

1.5

2

εe Fig. 6. Variation of the void aspect ratio in the most damaged element (i.e. the centre of the minimum cross-section area) as a function of the average effective strain in the minimum section.

0 10

-7

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

f0 Fig. 8. Variation of the ductility (fracture strain) as a function of the initial porosity for different initial void aspect ratio.

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

Considering only materials with early void nucleation, the trends presented above should be seen as upper bounds for the fracture strain of metallic alloys. Differences with the predictions of Figs. 7 and 8 can result from various effects such as premature void coalescence due to early ligament cracking, highly heterogeneous particle distribution or contributions from other damage mechanisms. 5. Conclusions and open issues for simulating ductile fracture during forming operation Micromechanics based constitutive models for plastically deforming and damaging materials have been developed since the late 70s starting with the initial version of the Gurson model [9] (see [37] for a detailed review). Most of the researches that have followed were driven by structural integrity issues of prime interest for the nuclear, naval and aeronautic industry. The focus was mostly on predicting the resistance to crack propagation and assessing whether a structure will survive in the presence of cracks under maximum loading conditions. Most strikingly, the strong geometry dependence of crack growth resistance curves emerged without the introduction of phenomenological parameters for the characterization of crack tip constraint (e.g. [34] for a review and [38–46] for ‘‘classical references’’). Crack tip problems are characterized by large stress triaxiality. Even with significant loss of constraint, the stress triaxiality remains usually larger than 1.5–2 and void shape effects are never significant (a notable exception comes with thin plate fracture, see [47]). This is probably why extensions of Gurson based models valid in the low stress triaxiality regime accounting for void shape related effects were not much studied (except for the case of viscoplastic materials undergoing low triaxiality creep conditions, see [1]). It is only recently that, motivated by metal forming applications (e.g. [48]) and by problems related to the ductile tearing of thin sheets, a need for new progresses in that area has appeared. An example of such enhanced Gurson model has been presented in this paper with its finite element implementation, its identification, and validation. The model has been used to simulate the damage evolution during uniaxial tension with necking for a wide range of flow properties and microstructure related parameters. These simulations have demonstrated the complex couplings existing between finite strain effects, strain hardening, and the microstructure parameters, mostly the initial porosity and shape. 5.1. Perspectives and open issues Many forming operations involve significant amount of shear. Voids rotate under shear deformation. There is very limited information available in the literature about this mechanism and how it may affect the void growth rate and the coalescence condition [49–52]. Void rotation is embedded in the present model but more a more accurate formulation is discussed elsewhere [14]. A limitation of

1649

the present model is the isotropic behaviour of the matrix. First steps towards the extension of Gurson types models to plastic anisotropy have been recently undertaken, for instance, by Benzerga et al. [5,7]. Finally, although quite sophisticated, the present model does not rest on a fully realistic description of the microstructure and micromechanisms of damage in industrial metal alloys which sometimes involve (i) delayed void nucleation (e.g. [12]), (ii) accelerated void coalescence due to local cleavage or the growth of a second population of defects (e.g. [27,53]), (iii) a competition with intergranular damage (e.g. [54]), (iv) strain gradient effects when the void size is smaller than a few micrometers [55]. Acknowledgements The formulation of the present model within a rigorous finite strain framework owes much to the collaboration with Patrick Onck from the University of Groningen. His help is gratefully acknowledged. This work as been carried out in the framework of the Project IAP P5/08 supported by the ‘‘Federal Public Planning Service Science Policy’’. References [1] Budiansky B, Hutchinson JW, Slutsky S. Void growth and collapse in viscous solids. In: Hopkins HG, Sewell MJ, editors. Mechanics of solids, The Rodney Hill 60th anniversary volume. Pergamon Press; 1982. p. 13–45. [2] Pardoen T, Hutchinson JW. An extended model for void growth and coalescence. J Mech Phys Solids 2000;48:2467–512. [3] Thomason PF. Ductile fracture of metals. Oxford: Pergamon Press; 1990. [4] Gologanu M, Leblond J-B, Perrin G, Devaux J. In: Suquet P, editor. Continuum micromechanics. Berlin: Springer-Verlag; 1995. p. 61– 155. [5] Benzerga AA, Besson J, Batisse R, Pineau A. Synergistic effects of plastic anisotropy and void coalescence on fracture mode in plane strain. Mod Simul Mater Sci Eng 2002;10:73–102. [6] Benzerga AA, Besson J, Pineau A. Anisotropic ductile fracture – Part I: experiments. Acta Mater 2004;52:4623–38. [7] Benzerga AA, Besson J, Pineau A. Anisotropic ductile fracture – Part II: theory. Acta Mater 2004;52:4639–50. [8] Lassance D, Scheyvaerts F, Pardoen T. Growth and coalescence of penny-shaped voids in metallic alloys. Eng Fract Mech 2006;73: 1009–34. [9] Klo¨cker H, Tvergaard V. Growth and coalescence of non-spherical voids in metals deformed at elevated temperature. Int J Mech Sci 2003;25:1283–308. [10] Gurson AL. 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 1977;99:2–15. [11] Pardoen T, Hutchinson JW. Micromechanics-based model for trends in toughness of ductile metals. Acta Mater 2003;51:133–48. [12] Huber G, Brechet Y, Pardoen T. Void growth and void nucleation controlled ductility in quasi eutectic cast aluminium alloys. Acta Mater 2005;53:2739–49. [13] Besson J, editorLocal approach to fracture. Paris, France: Les Presses de l’Ecole des Mines; 2004. [14] Gologanu M. Etude de quelques proble`mes de rupture ductile des me´taux. PhD thesis, Universite´ Paris VI, Paris, 1997. [15] Scheyvaerts F, Onck PR, Bre´chet Y, Pardoen T. A multiscale model for the cracking resistance of 7000 Al. In: Besson J, Moinereau D,

1650

[16]

[17]

[18] [19] [20] [21]

[22]

[23] [24]

[25] [26]

[27] [28] [29]

[30]

[31]

[32]

[33] [34]

[35] [36]

T. Pardoen / Computers and Structures 84 (2006) 1641–1650

Steglich D, editors. Proceedings of the 9th European mechanics of materials conference on local approach to fracture. Euromech – Mecamat 2006. Moret-sur-Loing, France: Les Presses de l’Ecole des Mines de Paris; 9–12 May 2006. p. 447–52. Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc Roy Soc London, Ser A 1957; 241:376–96. Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21: 571–4. ABAQUS/Standard version 6.3. User’s manual. Hibbit, Karlsson & Sorensen, Inc. Providence (RI), US, 2002. Simo JC, Hughes TJR. Computational inelasticity. New York: Springer-Verlag; 1998. Kojic M, Bathe KJ. Inelastic analysis of solids and structures. Springer; 2005. Eterovid AL, Bathe KJ. A hyperelastic-based large strain elastoplastic constitutive formulation with combined isotropic-kinematic hardening using logarithmic stress and strain measures. Int J Numer Meth Eng 1990;30:1099–114. Montans FJ, Bathe KJ. Computational issues in large strain elastoplasticity: an algorithm for mixed hardening and plastic spin. Int J Numer Meth Eng 2005;63:159–96. Bathe KJ. Finite element procedures. Prentice-Hall; 1996. Kim J, Gao X. A generalized approach to formulate the consistent tangent stiffness in plasticity with application to the GLD porous material model. Int J Solids Struct 2005;42:103–22. Hill R. Eigenmodal deformations in elastic/plastic continua. J Mech Phys Solids 1967;15:371–86. Aravas N. On the numerical integration of a class of pressuredependent plasticity models. Int J Numer Meth Eng 1987;24: 1395–416. Fabre`gue D. Pardoen T., submitted for publication. Spitzig WA, Kelly JF, Richmond O. Quantitative characterization of second-phase populations. Metallography 1985;18:235–61. Joly P, Meyzaud Y, Pineau A. Micromechanisms of fracture of an aged duplex stainless steel containing a brittle and a ductile phase: development of a local criterion of fracture. Advances in local fracture/damage models for the analysis of engineering problems. AMD, vol. 137. The American Society of Mechanical Engineers; 1992. p. 151–80. Pardoen T, Doghri I, Delannay F. Experimental and numerical comparison of void growth models and coalescence criteria for the prediction of ductile fracture in copper bars. Acta Mater 1998;46: 541–52. Pardoen T, Delannay F. Assessment of void growth models from porosity measurements in cold drawn copper bars. Metall Mater Trans A 1998;29:1895–909. Hancock JW, Mackenzie AC. On the mechanisms of ductile failure in high-strength steels subjected to multi-axial stress-states J. Mech Phys Solids 1977;14:147–69. Marini B, Mudry F, Pineau A. Experimental study of cavity growth in ductile rupture. Eng Frac Mech 1985;6:989–96. Pineau A. Global and local approaches of fracture – transferability of laboratory test results to components. In: Argon AS, editor. Topics in fracture and fatigue. Springer-Verlag; 1992. p. 197– 234. Rice JR, Tracey DM. On the ductile enlargement of voids in triaxial stress fields. J Mech Phys Solids 1969;17:201–17. Zhang ZL, Hauge M, Søvik OP. Determining material true stressstrain curves from tensile specimens with rectangular cross-section. Int J Solids Struct 1999;36:3497–516.

[37] Pineau A, Pardoen T. Failure mechanisms of metals. second ed. Comprehensive structural integrity, vol. 2. Elsevier; 2006 [chapter 6]. [38] Needleman A, Tvergaard V. An analysis of ductile rupture modes at a crack tip. J Mech Phys Solids 1987;35:151–83. [39] Rousselier G, Devaux J-C, Mottet G, Devesa G. A methodology for ductile fracture analysis based on damage mechanics: an illustration of a local approach of fracture. In: Landes JD, Saxena A, Merkle JG, editors. Nonlinear fracture mechanics: vol. II – Elastic–plastic fracture, ASTM STP 995. Philadelphia: The American Society for Testing and Materials; 1989. p. 332–54. [40] Mudry F, di Rienzo F, Pineau A. Numerical comparison of global and local fracture criteria in compact tension and center-crack panel specimens. In: Landes JD, Saxena A, Merkle JG, editors. Nonlinear fracture mechanics: vol. II – Elastic–plastic fracture, ASTM STP 995. Philadelphia: The American Society for Testing and Materials; 1989. p. 24–39. [41] Bilby BA, Howard IC, Li ZH. Prediction of the first spinning cylinder test using ductile damage theory. Fatigue Fract Eng Mater Struct 1993;16:1–20. [42] Xia L, Shih CF. Ductile crack growth – I. A numerical study using computational cells with microstructurally-based length scales. J Mech Phys Solids 1995;43:233–59. [43] Xia L, Shih CF, Hutchinson JW. A computational approach to ductile crack growth under large scale yielding conditions. J Mech Phys Solids 1995;43:389–413. [44] Brocks W, Klingbeil D, Kunecke G, Sun D-Z. Application of the Gurson model to ductile tearing resistance. In: Kirk M, Bakker A, editors. Constraint effects in fracture theory and applications: vol. 2, ASTM STP 1244. Philadelphia: The American Society for Testing and Materials; 1995. p. 232–52. [45] Ruggieri C, Panontin TL, Dodds Jr RH. Numerical modeling of ductile crack growth in 3-D using computational cell elements. Int J Fract 1996;82:67–95. [46] Gao X, Faleskog J, Shih CF. Cell model for nonlinear fracture analysis – II. Fracture-process calibration and verification. Int J Fract 1998;89:374–86. [47] Pardoen T, Hachez F, Marchioni B, Blyth H, Atkins AG. Mode I fracture of sheet metal. J Mech Phys Solids 2004;52:423–52. [48] Brunet M, Morestin F, Walter H. Damage identification for anisotropic sheet-metals using a non-local damage model. Int J Damage Mech 2004;13:35. [49] Fleck NA, Hutchinson JW. Void growth in shear. Proc Roy Soc London, Ser A 1986;407:435–58. [50] Fleck NA, Hutchinson JW, Tvergaard V. Softening by void nucleation and growth in tension and shear. J Mech Phys Solids 1989;37: 515–40. [51] Bordreuil C, Boyer J-C, Salle´ E. On modelling the growth and the reorientation changes of ellipsoidal voids in a rigid plastic matrix. Mod Simul Mater Sci Eng 2003;11:365–80. [52] Xue L. Void shearing in ductile fracture of porous materials. In: Besson J, Moinereau D, Steglich D, editors. Proceedings of the 9th European mechanics of materials conference on local approach to fracture. Euromech – Mecamat 2006. Moret-sur-Loing, France: Les Presses de l’Ecole des Mines de Paris; 9–12 May 2006. p. 483–88. [53] Faleskog J, Shih CF. Micromechanics of coalescence – I. Synergistic effects of elasticity, plastic yielding and multi-size-scale voids. J Mech Phys Solids 1997;45:21–45. [54] Pardoen T, Dumont D, Deschamps A, Brechet Y. Grain boundary versus transgranular ductile failure. J Mech Phys Solids 2003;51: 637–65. [55] Hutchinson JW. Plasticity at the micron scale. Int J Solids Struct 2000;37:225–38.