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 eﬀect 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 ﬁnite element code within a ﬁnite strain set up. Identiﬁcation procedures from experimental data are proposed and illustrated in an application on copper bars exhibiting diﬀerent 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 signiﬁcant 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 aﬀects the void growth rate with respect to the reference rate of growth of a spherical void, but aﬀects 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
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 deﬁned 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.
T. Pardoen / Computers and Structures 84 (2006) 1641–1650
during deformation which might also aﬀect 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 identiﬁcation procedures become more and more complex as well as the transfer from one material to another. These diﬃculties 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  for dilatant plastic materials based on the contributions by Gologanu et al.  and Thomason  has been developed in order to encompass both low and large triaxiality regimes within the same model  (see also ). 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 identiﬁed 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 identiﬁcation 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 oﬀers an interesting validation of the model as the strain hardening capacity of this material can be modiﬁed by thermal treatment without aﬀecting 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 ﬂow theory. Plastic anisotropy of the matrix can be introduced following the framework proposed by Benzerga et al. . 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 diﬀuse 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
Coalescence: kr0 k 3 jrh j þ F ðW ; vÞ ¼ 0 ry 2 ry sﬃﬃﬃ# " 2 3 1 v 1 F ðW ; vÞ ¼ ð1 v2 Þ a þb 2 vW v
Evolution law for f Growth and coalescence: f_ ¼ ð1 f Þ_epkk
Evolution law for W (or S = ln W) 2. An enhanced model for void growth and coalescence Growth : An extension of the Gurson model  based on the contributions by Gologanu et al.  and Thomason  and accounting for the eﬀect 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 ﬁtting parameters.
3 e_ p S_ ¼ ð1 þ h1 Þ e_ p kk d : P þ h2 e_ pkk 2 3
3 A spheroid is an ellipsoid with a symmetry of revolution around one axis.
T. Pardoen / Computers and Structures 84 (2006) 1641–1650
W ð2v2 1Þ 3 p 1 W_ ¼ e_ : P e_ pkk 2f 2 2
Evolution laws for ry Growth: ry e_ py ð1 f Þ ¼ r : e_ p
Coalescence: r_ y ¼
ory v2 p e_ oepy f e
Evolution for v Coalescence: vð3 2v2 Þ 3 p 1 p v_ ¼ e_ : P e_ kk 6f 2 2
In Eqs. (1)–(9), • The ﬁrst yield surface Ugrowth is a Gurson-type yield surface derived by Gologanu–Leblond–Devaux for spheroidal voids  and extended to strain hardening materials . The yield function Ucoalescence has been developed in the spirit of the work by Thomason , 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 deﬁned by rgh ¼ r : J. • X is a tensor associated to the void axis and deﬁned 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 deﬁned by (1 2a2)ez ez + a2ex ex + a2ey ey; P is a projector tensor, deﬁned by ez ez; d is the Kronecker tensor. • The symbol k k represents the von Mises norm. • Gologanu et al.  derived analytical relationships for the dummy parameters C, g, g, j, h1, h2, a2 in terms of the state variables W and f (see  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  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 ; b is equal to 1.24. • The plastic strain rate tensor e_ p is deﬁned in the reference conﬁguration; 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 eﬀective 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).
The response of the material under loading is ﬁrst 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 diﬀuse mode of plastic deformation. Except for pathological cases, the ﬁrst yield surface to be reached is thus Ugrowth. With increasing deformation Ugrowth ﬁrst 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 ﬁnite element void cell calculations . Note that the parameter v does not aﬀect 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
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 ofﬃﬃﬃ 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_
e_ z ¼ Xez
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 . A more accurate void rotation law, motivated by 3D ﬁnite element void cell calculations, as been recently developed (see  where a modiﬁed version of the coalescence model is also introduced). Isotropic elasticity is assumed while using the Eshelby –Mori–Tanaka  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
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
The constitutive model has been implemented in the ﬁnite element program ‘‘ABAQUS Standard’’ through a User deﬁned MATerial subroutine . The key aspects of the implementation are detailed hereafter by presenting ﬁrst 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 , see also  for an implementation of the Gurson model. The set of 10 ﬁrst-order diﬀerential 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’’ . 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  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  such that r_ ij ¼ Lijkl e_ kl
with Lijkl ¼ C ijkl
mij mkl A þ lmn mmn
where the Cijkl are the elastic moduli, and the m’s, l’s and A are given by (here U Ugrowth) lij ¼
oU ; orij
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
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  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
13 T 2c
9 1 þ 3Tc 2 2
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 diﬀuse 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 signiﬁcantly 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 . Another more advanced proposal has been recently worked out to account for the presence of a second population of voids . 3. Identiﬁcation and validation of the model 3.1. Identiﬁcation procedures The parameters of the model which have to be identiﬁed are : • the parameters related to the initial void distribution f0, W 0, v 0;
T. Pardoen / Computers and Structures 84 (2006) 1641–1650
• the eﬀective stress eﬀective 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 signiﬁcant part of the energy dissipated [2,11]. At low stress triaxiality, the energy associated to the ﬁnal ligament fracture process is always negligible when compared to the total work spent to reach that point). Two identiﬁcation procedures are discussed now: the ﬁrst 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 diﬀerent manners. 3.1.1. Identiﬁcation based on microstructure characterization • f0 – Voids almost always nucleate on second phase particles (inclusions, intermetallic particles, oxides, etc.). The ﬁrst 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 ﬂat. It has been shown that the value of the initial void shape W0 has no eﬀect 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 deﬁning an eﬀective 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 eﬀect 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
1–2% for initial porosity around 0.1–0.2%). If f0 is larger than about 0.1–0.2%, then an inverse identiﬁcation procedureis recommended in order to accurately determine the ry epy curve. The idea is to simulate the test using the damage model and to ﬁnd the correction to apply on the global stress–strain curve (see  for an example). Such procedure also delivers the ﬂow 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 identiﬁed 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 identiﬁed value is realistic or not, hence if the model is valid or not (see next section for the application on copper). 3.1.2. Identiﬁcation from mechanical tests When characterization methods are not available or that a faster identiﬁcation procedure is needed or that the identiﬁcation 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 deﬁne a sensible ﬁrst 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 diﬃcult 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  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 identiﬁed from a single
T. Pardoen / Computers and Structures 84 (2006) 1641–1650
uniaxial tensile test and two tests on notched round bars with two diﬀerent radii of curvature.
500 As-received Cu as-received
3.2. Application to copper 3.2.1. Identiﬁcation of the parameters The ﬁrst identiﬁcation 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 ﬁnally coalesce by internal necking until ﬁnal 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  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 aﬀecting the shape, volume fraction and size of the particles. This system oﬀers thus a means for assessing the ability of the model to capture the eﬀect 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  for details). The uniaxial stress–strain curves are compared in Fig. 1 showing the large diﬀerence in strain hardening between the two materials for strains smaller than 0.2. Uniaxial tensile tests were also performed on notched round bars. Diﬀerent radius of curvature were used in order to cover a wide range of stress triaxiality and to oﬀer 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 diﬀerent geometries.
Each particle has several neighbours, deﬁned using the Delaunay tessalation, and among these neighbours, the particle has one nearestneighbour. These two values are then averaged for the set of particles.
ε 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 identiﬁcation procedure after the onset of necking .
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)
3.2.2. Validation The uniaxial tensile tests on the notched and smooth copper bars were simulated using the ﬁnite element code ABAQUS Standard and the ‘‘User deﬁned MATerial’’ option to implement the constitutive model described in the previous section. Care has been taken to insure that the meshes were suﬃciently reﬁned and that the results were independent of the degree of reﬁnement (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 ﬁrst set of calculations on un-notched specimen was run for the as-received material using the parameters identiﬁed above and diﬀerent 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
2 FE results
Δεef = 0.5
Δn = 0.15
f0=0.2%, W0 =1, σ0 /E=0.001
As-received Cu 0
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  a long time ago and can be justiﬁed theoretically by simple void growth analysis, for instance with the Rice and Tracey model assuming a constant critical void radius . 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 ﬂow 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 ﬂow 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: ﬁnite strains eﬀects, 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 diﬀerent 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-
Mean stress triaxiality in most damaged region
Fig. 3. Variation of the ductility as a function of the strain hardening exponent.
f0=0.2%, W0 =1, σ0 /E=0.001
n = 0.1 n = 0.05
n = 0.2
1/3 = pure uniaxial tension
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 eﬀective strain in the minimum section.
onal arrangement (giving v0 = 0.149). Fig. 3 shows the variation of the ductility, deﬁned 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 signiﬁcantly 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 eﬀect 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 signiﬁcantly during plastic deformation – annealed and as-received Cu show closer strain hardening exponents at large strains when fracture occurs (see Fig. 1).
T. Pardoen / Computers and Structures 84 (2006) 1641–1650
increases more slowly with straining for larger n. This is a geometric ﬁnite strain eﬀect resulting from that larger hardening capacity promotes more diﬀuse 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 eﬀect is less important than the indirect eﬀect on the necking evolution. This discussion shows also that diﬀerent fracture strains are expected when changing the specimen conﬁguration 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. ). These results demonstrate the complex couplings existing between the stress state (related to the applied loading conﬁguration and to ﬁnite strain geometric eﬀects), the ﬂow properties (strain hardening) and the microstructure (shape, void volume fraction). Only an enhanced micromechanical model is able to properly capture these eﬀects.
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 eﬀect 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 ﬂat, i.e. W0 = 0.01 to very elongated, i.e. W0 = 100, with the strain hardening exponent always equal to 0.1. The eﬀect 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
f0=0.2%, W0=1, σ0/E=0.001
n = 0.1 0.05
n = 0.2 n = 0.1
n = 0.05
n = 0.2
0.5 10 0
n = 0.1 σy /E = 0.003
f0=0.2%, W0 =1, σ0 /E=0.001 n = 0.2
Fig. 7. Variation of the ductility (fracture strain) as a function of the initial porosity for diﬀerent strain hardening exponent.
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 eﬀective strain in the minimum section.
n = 0.1 1.5
n = 0.05
1.5 W0 =0.01
ε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 eﬀective strain in the minimum section.
f0 Fig. 8. Variation of the ductility (fracture strain) as a function of the initial porosity for diﬀerent 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. Diﬀerences with the predictions of Figs. 7 and 8 can result from various eﬀects 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  (see  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.  for a review and [38–46] for ‘‘classical references’’). Crack tip problems are characterized by large stress triaxiality. Even with signiﬁcant loss of constraint, the stress triaxiality remains usually larger than 1.5–2 and void shape eﬀects are never signiﬁcant (a notable exception comes with thin plate fracture, see ). This is probably why extensions of Gurson based models valid in the low stress triaxiality regime accounting for void shape related eﬀects were not much studied (except for the case of viscoplastic materials undergoing low triaxiality creep conditions, see ). It is only recently that, motivated by metal forming applications (e.g. ) 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 ﬁnite element implementation, its identiﬁcation, and validation. The model has been used to simulate the damage evolution during uniaxial tension with necking for a wide range of ﬂow properties and microstructure related parameters. These simulations have demonstrated the complex couplings existing between ﬁnite strain eﬀects, strain hardening, and the microstructure parameters, mostly the initial porosity and shape. 5.1. Perspectives and open issues Many forming operations involve signiﬁcant amount of shear. Voids rotate under shear deformation. There is very limited information available in the literature about this mechanism and how it may aﬀect 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 . A limitation of
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. ), (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. ), (iv) strain gradient eﬀects when the void size is smaller than a few micrometers . Acknowledgements The formulation of the present model within a rigorous ﬁnite 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  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.  Pardoen T, Hutchinson JW. An extended model for void growth and coalescence. J Mech Phys Solids 2000;48:2467–512.  Thomason PF. Ductile fracture of metals. Oxford: Pergamon Press; 1990.  Gologanu M, Leblond J-B, Perrin G, Devaux J. In: Suquet P, editor. Continuum micromechanics. Berlin: Springer-Verlag; 1995. p. 61– 155.  Benzerga AA, Besson J, Batisse R, Pineau A. Synergistic eﬀects of plastic anisotropy and void coalescence on fracture mode in plane strain. Mod Simul Mater Sci Eng 2002;10:73–102.  Benzerga AA, Besson J, Pineau A. Anisotropic ductile fracture – Part I: experiments. Acta Mater 2004;52:4623–38.  Benzerga AA, Besson J, Pineau A. Anisotropic ductile fracture – Part II: theory. Acta Mater 2004;52:4639–50.  Lassance D, Scheyvaerts F, Pardoen T. Growth and coalescence of penny-shaped voids in metallic alloys. Eng Fract Mech 2006;73: 1009–34.  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.  Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: Part I – Yield criteria and ﬂow rules for porous ductile media. J Eng Mater Tech 1977;99:2–15.  Pardoen T, Hutchinson JW. Micromechanics-based model for trends in toughness of ductile metals. Acta Mater 2003;51:133–48.  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.  Besson J, editorLocal approach to fracture. Paris, France: Les Presses de l’Ecole des Mines; 2004.  Gologanu M. Etude de quelques proble`mes de rupture ductile des me´taux. PhD thesis, Universite´ Paris VI, Paris, 1997.  Scheyvaerts F, Onck PR, Bre´chet Y, Pardoen T. A multiscale model for the cracking resistance of 7000 Al. In: Besson J, Moinereau D,
   
  
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 ﬁeld 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 misﬁtting 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 stiﬀness 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 ﬁelds. 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.
 Pineau A, Pardoen T. Failure mechanisms of metals. second ed. Comprehensive structural integrity, vol. 2. Elsevier; 2006 [chapter 6].  Needleman A, Tvergaard V. An analysis of ductile rupture modes at a crack tip. J Mech Phys Solids 1987;35:151–83.  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.  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.  Bilby BA, Howard IC, Li ZH. Prediction of the ﬁrst spinning cylinder test using ductile damage theory. Fatigue Fract Eng Mater Struct 1993;16:1–20.  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.  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.  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 eﬀects in fracture theory and applications: vol. 2, ASTM STP 1244. Philadelphia: The American Society for Testing and Materials; 1995. p. 232–52.  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.  Gao X, Faleskog J, Shih CF. Cell model for nonlinear fracture analysis – II. Fracture-process calibration and veriﬁcation. Int J Fract 1998;89:374–86.  Pardoen T, Hachez F, Marchioni B, Blyth H, Atkins AG. Mode I fracture of sheet metal. J Mech Phys Solids 2004;52:423–52.  Brunet M, Morestin F, Walter H. Damage identiﬁcation for anisotropic sheet-metals using a non-local damage model. Int J Damage Mech 2004;13:35.  Fleck NA, Hutchinson JW. Void growth in shear. Proc Roy Soc London, Ser A 1986;407:435–58.  Fleck NA, Hutchinson JW, Tvergaard V. Softening by void nucleation and growth in tension and shear. J Mech Phys Solids 1989;37: 515–40.  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.  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.  Faleskog J, Shih CF. Micromechanics of coalescence – I. Synergistic eﬀects of elasticity, plastic yielding and multi-size-scale voids. J Mech Phys Solids 1997;45:21–45.  Pardoen T, Dumont D, Deschamps A, Brechet Y. Grain boundary versus transgranular ductile failure. J Mech Phys Solids 2003;51: 637–65.  Hutchinson JW. Plasticity at the micron scale. Int J Solids Struct 2000;37:225–38.