International Journal of Solids and Structures 38 (2001) 8259±8284

[Link]/locate/ijsolstr

Modeling of crack growth in round bars and plane

strain specimens
J. Besson a,b,*, D. Steglich b, W. Brocks b
a
Ecole des Mines de Paris, Centre des Materiaux Pierre-Marie, UMR CNRS 7633, BP 87, 91003 Evry Cedex, France
b
Institute of Materials Research, GKSS Research Center, Geesthacht 21502, Germany
Received 26 September 2000

Abstract
The formation of cup-cone fracture in round bars and of slant fracture in plane strain specimens is studied using the
®nite element (FE) method. Constitutive models proposed by Rousselier [Nucl. Engng. Des. 105 (1987) 97] and by
Gurson [Acta Metall. 32 (1984) 157] are used. The analysis takes into account viscoplasticity and void nucleation.
Di€erent indicators of localization are computed during FE calculations. The analysis shows that cup-cone is more
easily formed using the Rousselier model than the Gurson model. Cup-cone fracture is inhibited in highly viscous
materials. The use of the f function in the Gurson model favors ¯at fracture. The crack path (¯at or cup-cone/slant)
can be correlated to the size of the localization zone which is formed ahead of the central penny shaped crack. Ó 2001
Published by Elsevier Science Ltd.

Keywords: Localization; Cup-cone fracture; Slant fracture Rousselier; Gurson

1. Introduction

Models able to represent the strength and toughness of ductile materials have found increasing interest
and application. The micromechanically based model proposed by Gurson (1977) and phenomenologically
extended by Tvergaard and Needleman (1984) (so called GTN model) has been most frequently used. An
approach based on continuum damage mechanics (CDM) and thermodynamics has also been proposed by
Rousselier (1987). Both models modify the von Mises yield potential by introducing a single scalar damage
quantity, namely the void volume fraction of cavities, f. Both models can be modi®ed to describe the
nucleation of cavities at inclusions or to account for viscoplastic behavior.
These models have been successfully applied to model crack propagation in precracked structures (e.g.
Xia et al., 1995). They have also been used to model fracture of small uncracked laboratory tests samples
such as smooth and notched round tensile bar (Tvergaard and Needleman, 1984; Becker et al., 1988) or

*
Corresponding author. Address: Ecole des Mines de Paris, Centre des Materiaux Pierre-Marie, UMR CNRS 7633, BP 87, 91003
Evry Cedex, France. Tel.: +33-160-76-30-37; fax: +33-160-76-31-50.
E-mail address: besson@[Link] (J. Besson).

0020-7683/01/$ - see front matter Ó 2001 Published by Elsevier Science Ltd.

PII: S 0 0 2 0 - 7 6 8 3 ( 0 1 ) 0 0 1 6 7 - 6
8260 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Fig. 1. Examples of experimentally observed cup-cone and slant fractures: (a,b) Isotropic material. (c,d) Anisotropic material; the
arrow indicates a secondary deformation band. (e,f) Plane strain specimen; secondary bands can also be observed on the surface. (g,h)
Notched bar.

plane strain specimens (Becker and Needleman, 1986; Leblond et al., 1994). Rupture of such specimens
involves both initiation and crack propagation. Cracks propagate in regions where deformation and
damage are localized leading to either ``cup-cone'' fracture in round bars or to ``shear lips'' in plane strain
samples (Fig. 1). This last fracture mode will be referred to as ``slant fracture'' in the following as the stress
state does not correspond to pure shear. Cup-cone formation has been ®rst numerically analyzed in
Tvergaard and Needleman (1984) but has received little attention since then.
As previously noted, cup-cone or slant fracture formation results from the localization of damage in
narrow bands. Conditions for localization (expressed as the possibility of forming a strain rate disconti-
nuity surface) in elastoplastic solids have been described by Rice (1976); the speci®c case of dilatant
pressure sensitive materials has been investigated in Rudnicki and Rice (1975). Using a simpli®ed model
consisting of a solid body having a plane of imperfection, the susceptibility to localization of the Gurson
model has been ®rst studied by Yamamoto (1978) with an emphasis on the e€ect of heterogeneities. This
model can also be used to study the evolution of the band beyond the onset of localization (Tvergaard,
1982). The role of void nucleation was studied in Saje et al. (1982) showing, in particular, that localization
can be growth or nucleation controlled depending on the stress state and material parameters. The role of
kinematic hardening was studied in Mear and Hutchinson (1985), together with nucleation in Tvergaard
(1987) showing a decrease of ductility with increasing kinematic hardening. The choice of the corotational
stress rate was studied in Tvergaard and van der Giessen (1991). In particular it is shown that localization
remains una€ected by the choice of the objective stress rate for purely isotropic hardening. Instability
in solids has also been studied using the linear perturbation analysis in the case of viscous non-voided
materials (Fressengeas and Molinari, 1985; Anand et al., 1987) incorporating heating induced by plastic
deformation, heat di€usion and inertia e€ects. The case of rigid±plastic dilatant solids has been treated in
Rousselier (1991, 1995a,b) showing that ``shear lips'' fracture dominates for small porosities whereas
``normal'' fracture occurs at high porosities.
The above analyses of instability are however limited to the ideal situation of an in®nite medium in
which a band-like discontinuity appears. In actual structures, this situation is indeed never met. In this case,
Billardon and Doghri (1989); Doghri and Billardon (1995) proposed to compute Rice's condition for lo-
calization during the ®nite element (FE) calculation; macro-crack initiation is assumed to occur when the
condition for localization is met; the calculation is then stopped as stability and uniqueness of the solution
are no longer insured.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8261

In this paper, the FE method is used to numerically investigate the conditions for the formation of cup-
cone and shear fracture. In Section 2, a common framework describing Gurson and Rousselier models is
presented together with an extension to viscoplastic materials. Model parameters are adjusted to experi-
mental results obtained on a modern steel containing very low sulfur and phosphorus contents. In Section
3, di€erent localization indicators are presented and compared. In Section 4, FE simulations of cup-cone
and slant fracture are carried out. The in¯uence of the following parameters is studied: (i) mesh re®nement
and element formulation, (ii) plastic and viscoplastic ¯ow, (iii) constitutive model: Gurson and Rousselier
models, e€ect of void nucleation. Following Billardon and Doghri, localization indicators are also evalu-
ated during the FE calculations. They are however not used as failure criterion; calculations are carried out
beyond localization. On the other hand, simulated crack paths are compared with predictions obtained
from the indicators.

1.1. Notations

Tensorial notation is used for convenience. First order tensors are denoted as ~ a, second order tensors as
a and fourth order tensors as A. Boldface symbols P denote matrices
P (M). Dots andPcolons are used to in-
dicate the usual contracted products: ~
a
~
b ˆ a b , a : b ˆ a b , A : a† ˆ
P i i i i;j ij ij ij k;l Aijkl akl , M
N†ij ˆ
k M ik N kj , etc. The symbol denotes the tensorial product. Voigt or standard notations for the tensors
will be used depending on the context. Identity tensors and matrices P are denoted by 1, 1 and I. The notation
~
n
L
~
n represents the second order tensor A such that Aij ˆ k;l nk Lkijl nl .

2. Material models

2.1. Constitutive equations

Models for porous materials proposed by Rousselier (1987) and Gurson (Gurson, 1977; Tvergaard,
1990) are used in this study. The ®rst one was developed based on the thermodynamical considerations
whereas the second one was derived from a micromechanical description of the porous material. In both
cases, damage is represented by a single scalar variable: the porosity f. In the following, both models will be
described using a uni®ed framework. The plastic (or viscoplastic) ¯ow potential / is then written as
/ ˆ rH R p† 1†
where R is the yield stress of the undamaged material (matrix) and p an e€ective plastic strain representative
of the matrix hardening. rH is an e€ective scalar stress which is a function of both the macroscopic stress
tensor r and the porosity. rH is de®ned by the following equations:
 
r2eq q2 rkk def:r
Gurson U ˆ 2 ‡ 2q1 fH cosh 1 q21 fH2 ˆ H 0 2†
rH 2 rH
 
req r1 rkk def:r
Rousselier Uˆ ‡ fD exp 1 ˆH 0 3†
1 f †rH rH 3 1 f †r1
where req is the von Mises equivalent stress and rkk the trace (tr) of the stress tensor. q1 , q2 , D and r1 are
material coecients which are assumed to be constant. fH is a function of the porosity f which was in-
troduced on a purely phenomenological basis to represent void coalescence (Tvergaard and Needleman,
1984). In the case of the Gurson model (Eq. (2)) the de®nition of rH is implicit whereas it is explicit in the
case of the Rousselier model (Eq. (3)).
8262 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

In the case of a plastic material, one has / ˆ 0 whereas in the viscoplastic case / > 0. The irreversible
deformation rate e_ p is obtained assuming the normality rule, so that:
o/ orH
e_ p ˆ 1 f †p_ ˆ 1 f †p_ ˆ 1 _
f †pv 4†
or or
where v is the normal to the ¯ow potential. The total strain rate e_ is expressed as e_ e ‡ e_ p , where e_ e is the
elastic strain rate. Elastic strains and stresses are related by: r ˆ C : ee . As the void volume fraction is
usually small, the elasticity tensor C is assumed to be constant. v can be computed noting that for a ®xed
porosity, a variation of r induces a variation of rH such that U remains equal to zero (this derivation must
be used in the case of the Gurson model, where the de®nition of rH is implicit). Therefore:
  1
oU oU orH oU oU
dU ˆ : dr ‡ drH ˆ 0 ) v ˆ ˆ 5†
or orH or orH or
The plastic multiplier p_ can be related to the deformation rate tensor. In the case of the Gurson model, one
notes that orH =or† : r ˆ rH , so that:
e_ p : r ˆ 1 _ H
f †pr 6†
In the case of the Rousselier model, the previous relation does not hold and one has:
q
p_ ˆ e_eq  23e_ : e_ 7†

where e_ is the deviator of the strain rate tensor. In this case, p_ corresponds to the von Mises equivalent
strain rate. p_ can be computed writing the consistency condition /_ ˆ 0 in the case of plasticity or using the
viscous ¯ow law of the dense material in the case of viscoplasticity:
p_ ˆ F /† ˆ F rH R† 8†
n
In the following the creep function F will be chosen as a Norton law, i.e. F /† ˆ h/=Ki , where K and n are
material parameters and where h
i is such than hxi ˆ x if x > 0, and hxi ˆ 0 otherwise.
The evolution of the porosity is given by mass conservation modi®ed to account for strain controlled
void nucleation (Chu and Needleman, 1980):
 
f_ ˆ 1 f †tr e_ p ‡ An p_ ˆ 1 f † v : 1 ‡ An p_
2
9†

An is a material function used to represent nucleation. In the following, it will be assumed that it depends on
p only.

Remark 1. In the case of the viscoplastic Rousselier model (Barbier, 1999), the ratio f_ =p_ depends on p_ so
that porosity grows faster for higher strain rate. This is due to the fact that r1 is taken to be rate inde-
pendent so that the ratio rkk =r1 will increase with strain rate. This dependence is not observed for the
Gurson model as U can be expressed as a function of r=rH only (in the case of the Rousselier model U is a
function of both r=rH and r=r1 ). A solution could be to modify the Rousselier model as follows:
 
req 2 qR rkk def:r
Uˆ ‡ fD exp 1 ˆH 0 10†
1 f †rH 3 2 1 f †rH
where qR is a new adjustable material parameter. In that case the strain rate dependence is suppressed. The
study of this new model is however out of the scope of the present work. In the following, the original
model will be used with the model parameters r1 and D adjusted for the experimental strain rate. Para-
metric studies involving a strong e€ect of viscosity will be carried out using the Gurson model only.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8263

Fig. 2. Comparison of the Gurson and Rousselier yield surfaces.

2.2. Comparison: Gurson±Rousselier models

At this point, it is interesting to outline some di€erences between the Gurson and Rousselier models. Fig.
2 compares both yield surfaces in the req ±rkk plane in the case of tensile stress states (rkk > 0). Under pure
shear (rkk ˆ 0), damage is still generated in the case of the Rousselier model (as the normal to the yield
surface does not coincide with the req axis) whereas, in absence of nucleation, the Gurson model does not
lead to damage growth. Under pure hydrostatic stress states (req ˆ 0), the Rousselier yield surface has a
vertex which implies that at high stress triaxiality ratios s ˆ 1=3†rkk =req † the plastic deformation tensor
always keeps a non-zero shear component. Note that the model proposed in Fleck et al. (1992) for plastic
metal powders has the same property.
The ratio of volumetric to shear deformation rate (strain rate triaxiality ratio) se ˆ 1=3†e_kk =e_eq is given
for both models by
 
1 q2 rkk rH
Gurson se ˆ q1 q2 fH sinh 11†
2 2 rH req
 
1 rkk
Rousselier se ˆ fD exp 12†
3 3 1 f †r1
Taking the limits of the strain rate triaxiality ratio for req ! 0 such that the yield condition is met, one
gets: 1
 3=2
q1 fH rH
Gurson lim se ˆ q2 ˆ ‡1 13†
req !0 2 req

1
Using the modi®ed Rousselier model (Eq. (10)), the strain rate triaxiality ratio is equal to 1=3†DfqR exp qR rkk =2 1 f †rH †. The
limit for req ! 0 is equal to qR =2.
8264 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

1 rH
Rousselier lim se ˆ 14†
req !0 3 r1

2.3. Material parameters

The material parameters used in this study correspond to a X70 HSLA (high strength low alloyed)
ferritic±pearlitic steel (Rivalin, 1998). Due to its very low sulfur and phosphorus content, the inclusion
volume fraction is equal to 1:5  10 4 . Inclusions consist of small, globular particles of low mean diameter
(about 1 lm), composed of two phases combining calcium sul®de (CaS) and aluminum (Al2 O3 ) or mag-
nesium oxides (MgO). It is assumed that immediate debonding between the matrix and the inclusions takes
place so that the inclusion volume fraction corresponds to the initial porosity f0 .
Due to the low inclusion volume fraction, the e€ect of damage on the overall behavior remains very
limited up to elevated plastic strains. This implies that parameters relative to plastic hardening and viscosity
can be directly determined from tensile tests. The material was provided as hot rolled sheets; its plastic
behavior is therefore anisotropic as shown in (Rivalin, 1998). In the following, it will however be assumed
that the material is isotropic in order to simplify the FE calculations. Elastic properties are: Young's
modulus E ˆ 210 GPa, Poisson's ratio m ˆ 0:3. Plastic hardening is described using a simple power law
relationship:
n0
R p† ˆ K 0 p ‡ e0 † 15†
0 0 1=n
with K ˆ 795 MPa, e0 ˆ 0:002 and n ˆ 0:13. Parameters of the Norton law are: K ˆ 55 MPa s , n ˆ 5.
Tests carried out at di€erent temperatures have also shown that the behavior is almost una€ected by
temperature up to 300 °C so that heating due to plastic deformation can be neglected.
Ductile rupture was characterized using smooth and notched round bars (Mackenzie et al., 1977;
Decamp et al., 1998) as well as plane strain specimens (Anand and Spitzig, 1980). Round bars have an
initial minimum diameter equal to 10 mm. The notch radius is equal to 4 and 2 mm. Plane strain specimens
have a thickness of 5 mm. The area reductions at fracture are given on Table 1 for all specimens. Fig. 1
shows examples of fracture surfaces obtained in round bars and plane strain specimens. Note the aniso-
tropic deformation of bars. An example of cup-cone formation in another isotropic material is also given in
Fig. 1. This material has a similar composition as the material of this study but was subjected to a thermal
treatment leading to a ferritic±bainitic microstructure and to isotropic plastic properties.
Damage parameters (r1 and D for the Rousselier model; q1 , q2 and fH for the Gurson model) were
adjusted to represent the experimental area reductions at fracture. In the case of the Rousselier model, the
values recommended in Rousselier (1987) are D ˆ 2 and r1 ˆ 1=3† Re ‡ Rm † ˆ 321 MPa, where Rm is the
maximum engineering stress and Re the yield limit. These values lead however to an underestimation of
the ductilities. Adjusted values are equal to D ˆ 1:4 and r1 ˆ 450 MPa. In the case of the Gurson model,
many studies published in the literature use q1 ˆ 1:5 and q2 ˆ 1:0 whereas the function fH is simply de®ned
as follows:

f if f < fc
fH ˆ 16†
fc ‡ d f fc † if f > fc

Table 1
Area reductions at fracture
Smooth Notched Notched Plane
r ˆ 4 mm r ˆ 2 mm strain
0.79 0.60 0.46 0.40
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8265

where d > 1 and fc (critical porosity at which coalescence starts) have to be adjusted. Using fH ˆ f leads to
an overestimation of the ductilities. In Section 4.3 di€erent solutions (use of fH , use of nucleation, modi-
®cation of q2 ) will be envisaged in order to adjust the experimental data.
Failure occurs for the Gurson model when fH ˆ 1=q1 ' 66%. In the case of the Rousselier model, the
material loses its stress carrying capacity when f reaches 1. In order to obtain a more realistic failure
porosity, the material is usually considered as broken when f is larger than a critical value fR . In the fol-
lowing fR ˆ 0:9 will be used.

3. Localization indicators

In a in®nite homogeneous medium, localization is assumed to occur when it becomes possible to form a
strain rate discontinuity in a planar band. This band is characterized by its unit normal ~ n and the dis-
placement jump across the band whose direction is denoted ~ g (Fig. 3). Note that the magnitude of the jump
remains unknown. In the case of voided materials, ~ n and ~g are not necessarily orthogonal (Rousselier,
1995a).
In the following, Rice's condition for bifurcation is presented together with analytical results concerning
localization angles and critical hardening modulus allowing to further compare Rousselier and Gurson
models (Section 3.1). The perturbation analysis, which can be applied to viscoplastic materials is then
presented (Section 3.2). In Section 3.3, the consistent tangent matrix is computed. It is proposed to use this
matrix instead of the elastoplastic tangent matrix in Rice's condition for localization. In Section 3.4 the
di€erent localization indicators are compared; in particular it is checked that the indicator using the
consistent tangent matrix gives predictions in agreement with both other indicators.

Fig. 3. Rupture mode map.

8266 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

3.1. Plastic materials: bifurcation analysis

In the case of elastoplastic materials, the incremental constitutive equation can be expressed as
r_ ˆ Lt : e_ 17†

where Lt is the elastoplastic tangent matrix. The calculation of Lt is detailed in Appendix A. Writing the
continuity of displacements and the stress equilibrium, it can be shown (Rice, 1976; Rice and Rudnicki,
1980) that the jump of the deformation tensor is expressed as
1
2
~
g ~
n ‡~
n ~
g† 18†
and that the condition for bifurcation is written as
9~
n; det At ~
n†† ˆ 0 with At ~
n† ˆ ~
n
Lt
~
n 19†

~
g is then the eigenvector of At ~n† corresponding to the eigenvalue equal to zero. This condition corresponds
to continuous bifurcation (plastic yielding on each side of the band). Discontinuous bifurcation (plastic
yielding on one side and elastic unloading on the other side) corresponds to det At ~ n† < 0 (Rice and
Rudnicki, 1980; Borre and Maier, 1989). Eq. (19) can be modi®ed in the case of large deformations.
The stress rate in Eq. (17) being considered as the Jauman rate, the bifurcation criterion is now written
as: 9~n; det At ~n† ‡ R† ˆ 0 with 2R ˆ ~ n ~ n
r† ‡ ~n
r† ~n‡ ~n
r
~
n†1 r (Rice and Rudnicki, 1980;
Mear and Hutchinson, 1985).
The bifurcation criterion was implemented as a post-processor of the FE calculations. In practice, the
condition det At ˆ 0 is never exactly met. Localization will be de®ned as occurring when det At < 0 for the
®rst time. ~n is then de®ned as the vector which minimizes det At ~n†. The corresponding minimum value will
be used as a localization indicator. Minimization is done using the simplex algorithm. The eigenvector
corresponding to the minimum eigenvalue of At ~ n† is then computed; it coincides with ~
g when det At be-
comes zero. In the case of the specimens and constitutive equations investigated in this work, accounting
for the Jauman rate modi®es the results only slightly. In the following, results using Eq. (19) will be shown
only.
The bifurcation criterion can be used to further compare the localization behavior of Gurson and
Rousselier models. The plastic ¯ow direction, expressed in its eigen-coordinate system, is supposed to be of
the following form:
0 1
1 0 0
v / @ 0 u 0 A with t 6 u 6 1 20†
0 0 t
This corresponds to a general tensile situation where minor stresses can be negative. Assuming that
t 6 u 6 1, the vectors ~ n and ~
g lie in the x1 ±x3 plane provided dilatancy is not too large (Rudnicki and Rice,
1975; Yamamoto, 1978). This localization plane is always assumed in the following. ~ n and the angle w
between ~ g and ~ n (Fig. 3) are given by
r r!
1 ‡ mu t mu 1 ‡ 2mu ‡ t
~
nˆ ; 0;  cos w ˆ 21†
1 t 1 t 1 t

The angle w characterizes the type of failure which varies from pure opening fracture w ˆ 0 to pure tan-
gential fracture w ˆ p=2††. The results are summarized in Fig. 3 in the t±u plane for 2 < u < 1. The case
t ‡ mu > 0 (shaded area) corresponds to pure opening fracture (~ nk~
g), the normal to the band being aligned
with the principal strain direction. The case 1 ‡ mu < 0 corresponds, for m ˆ 0:3, to negative volume changes
(crosshatched region) and is not relevant in the present study. Constant band orientation (constant h)
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8267

or constant fracture mode (constant w) are represented by straight lines in the previous diagram. Location
of constant strain rate triaxiality (se ˆ 0; 1=4; 1=2; 1; 2) are also indicated on the diagram (grey lines). In
particular, it can be seen that for se < 1=2 normal separation is never possible. As noted in Section 2.2,
se lies between 0 and 13 rH =r1 for the Rousselier model. So that, in some cases, normal separation will be
impossible. One can therefore generally expect a higher tendency to pure opening fracture for the Gurson
model than for the Rousselier model.
The condition for localization det At ~ n† ˆ 0 can then be rewritten using the previous solution for ~ n as
0
 2 2
H v11 Eu if t ‡ mu < 0 2
ˆ with H 0 ˆ H ‰ 1 f † v11 1 ‡ u ‡ t† ‡ An ŠrH;f
1 f v211 Eu2 ‡ 2mut ‡ t2 =1 m2 if t ‡ mu P 0
22†
0
where H is the plastic hardening modulus and rH;f ˆ orH =of . H incorporates both plastic hardening and
softening due to porosity growth. The states most resistant to localization are those of axisymmetric ex-
tension u ˆ t whereas localization is much easier for plane strain u ˆ 0 (Needleman and Rice, 1978).

3.2. Viscoplastic materials: linear perturbation analysis

The previous bifurcation analysis is however only valid for materials for which the elastoplastic tangent
matrix is de®ned. For viscoplastic materials, the linear perturbation analysis can be applied. The method
consists in analyzing, inside an homogeneous volume element, the stability of a perturbation ~ g of the
displacement ®eld which is supposed to be of the following form (Fressengeas and Molinari, 1985; Anand
et al., 1987, Rousselier, 1991, 1995b; Barbier et al., 1998; Barbier, 1999):
~
g ˆ d~ x
~
u exp iq~ n ‡ xt† 23†
x is the growth rate of the perturbation and 1=q its characteristic length. Note that q plays only a role in
processes involving length scales such as thermal di€usion, dynamic loading or non-local constitutive
equations. It will not be considered hereafter.
Following the generic treatment of the problem proposed in Barbier et al. (1998), the material is
characterized by internal variables denoted Z ˆ ee ; z†. In the present case, z represents the plastic strain
and the porosity. The evolution laws of these variables are written as a set of di€erential equations:
Z_ ˆ F Z; e_ † 24†
which, applying a perturbation, leads to:
oF oF
dZ_ ˆ
dZ ‡
d_e 25†
oZ o_e
_
The perturbation dZ can be estimated from the perturbed rate as (Anand et al., 1987): dZ ˆ dZ=x. The
perturbed state variables are then related to the perturbed strain by
  1
1 oF oF
dZ ˆ I :
de ˆ Hp x†
de 26†
x oZ o_e

The calculation of Hp is detailed in Appendix B for the viscoplastic case. As shown in Barbier et al. (1998),
the plastic case coincides with Rice analysis for x ! ‡1. The sub-matrix of Hp relating the perturbed total
strain to the perturbed elastic deformation is denoted Hp . The perturbation on the stresses is therefore
computed as
dr ˆ C : Hp x† : de ˆ Lp x† : de 27†
8268 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Writing the stress equilibrium, as in the case of the bifurcation analysis, leads to the following condition for
the appearance of a localization band:
n; det Ap ~
9~ n; x† ˆ 0 with Ap ~
n; x† ˆ ~
n
Lp x†
~
n 28†

~
g (or d~
u) is the associated eigenvector. The criterion for localization is then applied by ®nding the maximum
value of x for which a vector ~ n exists which veri®es det Ap ~
n; x† ˆ 0. Localization occurs when the rate of
variation of the perturbation is much larger than the rate of variation of the unperturbed solution. This
condition can be expressed as: x  p=p. _

3.3. Consistent tangent matrix as localization indicator

In the framework of the FEM, material constitutive equations can be integrated using an implicit in-
tegration scheme (Simo and Taylor, 1985; Chaboche and Cailletaud, 1996; Aravas, 1987; Foerch et al.,
1997). The increment of the state variable DZ over a ®nite time step Dt is given by the following set of
implicit (with respect to DZ) non-linear equations:
I0
DZ ˆ FH Z; Dt; De† 29†
In the following, only fully implicit integration will be used so that Z ˆ Z0 ‡ DZ where Z0 is the known set
of state variables at the beginning of the time increment. The matrix I0 , de®ned in Appendix C, is used to
describe plasticity and viscoplasticity in the same framework. Eq. (29) is solved using a Newton±Raphson
method which requires the calculation of oFH =oZ (see Appendix C). Any in®nitesimal variation of the
deformation increment will induce a variation of the solution state variables such that Eq. (29) is still
satis®ed, so that:
oFH oFH
I0
dZ ˆ
dZ ‡
de 30†
oZ oe
Solving the previous equation for dZ gives:
  1
oFH oFH
dZ ˆ I0

de ˆ Hc
de 31†
oZ oe
The calculation of Hc is detailed in Appendix C. Using the same arguments as in Section 3.2, the variation
of the stresses is related to the variation of the total strain by
dr ˆ C : Hc : de ˆ Lc : de 32†

where Hc is the sub-matrix of Hc relating de and dee . Lc corresponds to the so-called ``consistent tangent
matrix''. If a small variation of stresses dr develops in a planar band, the mechanical equilibrium of the
n ˆ~
band requires that: dr
~ 0. The corresponding variation of the deformation must be of the form given by
Eq. (18) to ensure the continuity of displacement. Using the equilibrium condition together with Eq. (32), a
non null vector ~g exists if:
n; det Ac ~
9~ n†† ˆ 0 with Ac ~
n† ˆ ~
n
Lc
~
n 33†

This condition is similar to Rice localization criterion (Eq. (19)) with Lt being replaced by Lc . In the fol-
lowing, it is proposed to use it to post-process the FE calculation in a similar way as in Section 3.1. A
potential advantage of using the consistent tangent matrix is that the same criterion can be applied to both
elastoplastic and elastoviscoplastic materials.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8269

Fig. 4. Minimum value of the localization indicator as a function of the deformation using the elastoplastic tangent matrix (dashed
line) and for consistent tangent matrix (full line) with De ˆ 0:1%, 0:2%, 0:5%, 1:0% (plane strain, Gurson model).

3.4. Comparison of the di€erent localization indicators

In this section, the three previously described localization indicators will be compared for plane strain
conditions using the Gurson model with q2 ˆ 1:15 and fH ˆ f . Otherwise speci®ed, material parameters are
those given in Section 2.3.
Fig. 4 compares, in the case of an elastoplastic material behavior, the localization indicator based on the
tangent matrix to the one based on the consistent matrix using constant deformation increments in the
loading direction. Values of the indicator are normalized with respect to the value De corresponding to a
purely elastic behavior:
1 1 m†E3
n
C
~
De ˆ det~ nˆ 8~
n 34†
4 1 2m† 1 ‡ m†3
It can be seen that both indicators essentially give similar results for increments of deformation De up to
0.5%. A slight di€erence is observed for De ˆ 1:0%. For larger increments, convergence is not always
obtained. In that case sub-stepping must be used.
Indicators based on the perturbation method or the consistent tangent matrix give qualitatively similar
trends: stability is increased with increasing K and strain rate and with decreasing n. Fig. 5(a) gives an
example of the evolution of the growth factor x e ˆ x= p=p†
_ as a function of deformation for di€erent values
of K. An interesting feature of the perturbation analysis, is that it provides a continuous evaluation of the
localization angles. An example of the evolution of w and h is given on Fig. 5(b) showing a progressive
change from pure tangential to pure opening fracture as deformation increases. Curves obtained for dif-
ferent values of K coincide as soon as the growth factor is large enough (1000).
Fig. 6 compares, in the case of elastoviscoplastic material behavior, the localization indicator based on
the consistent matrix for di€erent values of the creep parameter K. It can also be seen that the results
converge to the results obtained for an elastoplastic material (dashed line) as K ! 0. Similar results are
obtained for e_ ! 0. As shown in Table 2, the porosity at the onset of localization (de®ned by
min~n det~n
Lc
~
n ˆ 0) increases with increasing K; this also corresponds to an increase of h > p=4 and to a
decrease of w < p=2: the normal component of the fracture mode increases. For K ! 0 (or e_ ! 0) the
viscoplastic analysis (porosity, angles h and w, strain) coincides with the plastic case. As in the plastic case,
the deformation increment in¯uences the consistent tangent matrix and consequently the time step for
which min~n det~ n
Lc
~
n starts to be negative. This e€ect remains limited for plasticity (Fig. 4) but increases
with increasing viscosity. For example, for K ˆ 50 the corresponding deformation is equal to 2:09 for
8270 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

_
Fig. 5. (a) Localization indicator x= p=p† as a function of the deformation for di€erent values of the creep parameter K (b) Local-
ization angles as a function of the deformation. (e_ ˆ 2  10 3 s 1 , Gurson model).

Fig. 6. Minimum value of the localization indicator as a function of deformation using the consistent tangent matrix for an elasto-
viscoplastic material (full lines) for di€erent values of the creep parameter K (De ˆ 0:5% and e_ ˆ 2  10 3 , Gurson model). The
elastoplastic case is indicated by dashed lines.

Table 2
Angles h and w (Fig. 3), porosity and deformation (e) at localization predicted using the consistent tangent matrix for elastoplastic and
elastoviscoplastic behaviors (Gurson model, e_ ˆ 2  10 3 s 1 , De ˆ 0:5%). The evolution of the corresponding localization indicator is
given on Fig. 6
K h (°) w (°) f (%) e
pl. 45.65 88.67 1.31 1.29
1 45.73 88.53 1.45 1.32
5 46.06 87.88 2.14 1.44
10 46.49 87.03 3.01 1.55
25 46.95 84.05 6.32 1.80
50 40.75 78.44 13.1 2.09

De ˆ 0:5%, 1.78 for De ˆ 1:0% and 1.56 for De ˆ 2:0%. This corresponds to the fact that non-converging
global time steps can be divided into converging sub-steps.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8271

Table 3
Comparison of FE calculations and localization analysis (De ˆ 0:5%, e_ ˆ 2  10 3 )
K FE calculation Localization analysis
eD eB eD
x eB
x eL
pl. 1.28 1.66 ± ± 1.29
1 1.34 1.67 6882 61150 1.32
5 1.43 1.67 2350 12800 1.44
10 1.53 1.69 2300 6860 1.55
25 1.67 1.71 2252 2990 1.80
50 1.76 1.77 1850 1900 2.09
100 1.82 1.83 1140 1250 2.54
eD , deformation at which min det Ac < 0 at, at least one Gauss point in the FE calculation; eB , deformation at which a localization
e D , normalized growth factor corresponding to eD ; x
band is formed in the FE calculation; x e B , normalized growth factor corresponding
to eB ; eL , deformation at which min det Ac < 0 for the homogeneous analysis.

In order to further evaluate the di€erent localization indicators, comparisons were made with simple FE
calculations. The mesh consists in a 10  10 plane strain small deformation elements square computed
using periodic boundary conditions. The imposed deformation step is 0:5% and sub-stepping was allowed.
In practice, sub-stepping is needed when numerical localization starts. An initial imperfection is introduced
as a uniformly distributed random porosity ®eld f0 ˆ 1:5  10 4  1:0  10 6 . The calculations were post-
processed in order to determine: (i) the strain eD at which min~n ~
n
Lc
~
n becomes negative at least one Gauss
point in the structure, (ii) the strain eB at which a continuous elastically deforming band is formed. Using
the perturbation analysis of a uniform volume element, the growth factors x e D and x e B corresponding these
strains can be computed. Finally, eL is de®ned as the strain at which min~n ~ n
Lc
~
n < 0 in a uniformly
deforming element. Results are gathered in Table 3. For plasticity and low values of K one gets
eD  eL < eB . This shows that the indicators underestimate the actual numerical onset of localization.
Similar results are obtained for plastic behavior using Eq. (19). On the other hand, for high values of K, one
gets eD  eB < eL . In that case, the indicator computed for a uniform material overestimates the actual
localization. On the other hand, the indicator computed during the FE calculation at each Gauss point
provides in that case pertinent informations about the onset of localization. Values of x e D and xe B are
decreasing functions of K and onset of localization does not correspond to some ``critical value''. As soon
e ®eld becomes strongly heterogeneous with some very high values (typically over
as localization starts the x
104 ). Similar conclusions were drawn from another set of FE calculations carried out on a 10  10 square
with an initial geometrical imperfection, homogeneous initial porosity and common boundary conditions.
From this part, it can be concluded that localization indicators can be used in FE calculations (Section 4)
to determine where localization is currently occurring. An overestimation of the zones where localization
occurs is to be expected for plastic and slightly viscous materials.

4. Finite element simulationsÐdiscussion

Finite element simulations were performed using the FE softwares A B A Q U S (A B A Q U S , 1998) and
Zebulon (Besson and Foerch, 1997; Foerch et al., 1997). The FE implementation of both Gurson and
Rousselier models in A B A Q U S follows the method proposed in Aravas (1987). The method used in Zebulon
is detailed in Appendix C. Both methods use a fully implicit integration scheme. It was shown by Zhang
(1995) that more accurate integration can be obtained using a semi-implicit scheme with f ˆ 0:75±0.85
(f ˆ 0 fully explicit, f ˆ 1 fully implicit). The Aravas method uses a reduced set of integration variables
(scalar deviatoric and volumetric components of the plastic strain) to describe deformation. On the other
8272 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

hand, Zebulon uses the full elastic strain tensor. Although less numerically ecient, this makes it possible to
easily extend the method to plastically anisotropic materials (Grange et al., 2000) and to directly compute
the consistent tangent matrix. In the present study the fully implicit algorithm is used. Finite strains were
treated using corotational reference frames (A B A Q U S (Hughes and Winget, 1980), Zebulon (Ladeveze,
1980)). The localization analysis was implemented in Zebulon only. It is made in the rotated material frame,
de®ned using the Jauman rate for the stresses, for each Gauss point using Eqs. (19) and (28) or Eq. (33).
Vectors ~n and ~
g can then be expressed in the ®xed reference frame. In the following, calculations obtained
using Zebulon are only presented. Very similar results in terms of mesh dependence, occurrence of local-
ization, formation of cup-cone were obtained with A B A Q U S .

4.1. E€ect of meshing

As already mentioned in Tvergaard and Needleman (1984), mesh design plays an important role in
describing localization and cup-cone formation. Tvergaard and Needleman (1984) used square elements
(four nodes) which were divided in four linear triangles (three nodes). In order to numerically study the
e€ect of meshing, calculations were performed for round bars using the Rousselier model as this model was
found to more easily lead to cup-cone formation as the Gurson model.
Specimens were meshed without using an initial geometrical imperfection as in Tvergaard and Nee-
dleman (1984). However the loading ends were meshed; this is sucient to generate stress and strain
heterogeneities so that necking and subsequent failure always occur in the middle of the specimen. An
example of mesh is shown on Fig. 7. Symmetry was not enforced so that the entire specimen and not only
one half is meshed. This choice was motivated by the experimental observation that one single crack (and
not two symmetric cracks) is generated during cup-cone formation. Ductility was characterized using the
diameter reduction as the tensile elongation is highly sensitive to strain path and calculation parameters
(Bonora, 1999). Simulations were carried out with a macroscopic strain rate equal to 2  10 3 s 1 .
In the following porosity maps will be presented using the following scale 0 white† 6 f 6 0:1 black†.
Fig. 7 shows damage maps at Gauss points for meshes having a number of elements in the minimum cross

Fig. 7. E€ect of mesh re®nement and element initial aspect ratio on the formation of the cup-cone (cax8r elements).
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8273

Fig. 8. E€ect of the element type on the formation of the cup-cone (cax8r, cax6r, cax4 and cax3 elements). A cup-cone is formed
using cax4 and Nh ˆ 80. A calculation with initial heterogeneous distribution of the porosity is also shown: a cup-cone with a ``zigzag''
is formed. The crack path in the undeformed mesh is also shown for cax8r and cax6r elements. The arrow indicates a change in the
crack direction in the undeformed mesh.

section Nh equal to 20, 40 or 80 and for initial element aspect ratios rh equal to 3:1, 6:1 and 12:1. Elements
have quadratic shape functions with eight nodes and reduced integration (4G points). The aspect ratio
equal to 6:1 leads to approximately square elements at the onset of fracture. Too ¯at elements (rh ˆ 12:1)
always lead to ¯at fracture. A minimum number of elements in the cross section is required to trigger the
cup-cone: it was found that Nh P 30 is needed here. Provided this condition is met, elements with an initial
aspect ratio of 3:1 (thus leading to elongated elements at fracture) can also generate cup-cone. In the
following and otherwise stated, elements with rh ˆ 6:1 and Nh ˆ 40 will be used.
Fig. 8 compares the damage maps obtained with di€erent type of elements. First linear (cax4: four
nodes, 4 G points) and quadratic (cax8r: eight nodes, 4 G points) square elements are compared. The
selective integration method (Hughes, 1980) is used in the case of the cax4 element. It is shown that
the cax4 leads to ¯at fracture. With this element, using 80 elements leads also to cup-cone formation. The
number of degree of freedom (DOF) is then approximately the same as in the case of the cax8r elements
with Nh ˆ 40. This is also consistent with the ¯at fracture path obtained in the case cax8r±±Nh ˆ 20 for
which the number of DOF (4302) is about the same as for cax4±±Nh ˆ 40 (4018).
The previous square elements were divided into four triangular elements as done in Tvergaard and
Needleman (1984): cax8r ! cax6r (six nodes, 3 G points) and cax4 ! cax3 (three nodes, 1 G point).
In both cases, cup-cone is observed. A di€erence can however be noted between the cax8r and the cax6r/
cax3 elements. In the case of the triangular elements the highly damaged band will tend to stay along a
preferred mesh direction. The crack appears to have more freedom to follow a direction not related to the
mesh in the case of cax4 and cax8r elements (Fig. 8). As a result, cup-cone angles obtained with qua-
dratic and triangular elements di€er. Note that the zigzagging crack path obtained with triangles is a
numerical artifact generated when the crack comes close to the coarse elements zone of the mesh.

Remark 2. (E€ect of symmetry). Fig. 9 compares calculations for a notched bar carried out using an entire
mesh or a half mesh with symmetry conditions. A zigzagging crack is obtained with the later whereas a ¯at
crack is simulated with the ®rst. This indicates that cup-cone formation is favored when using symmetry.
8274 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Fig. 9. Fracture path in notched round specimens.

Fig. 10. Load vs. diameter reduction curve and cup-cone formation for a round bar. Contour maps indicate damage (white: f ˆ 0,
black: f > 0:1) and the localization indicator (white: elastic unloading, gray: plastic loading, black plastic loading and negative in-
dicator).

4.2. Cup-cone formation

The evolution of damage f and of the localization indicator min~n ~ n
Lc
~

n is detailed in Fig. 10 during
cup-cone formation. After the maximum load has been reached, progressive necking is observed and the
specimen can be divided into two regions: elastic unloading and plastic loading ((1) in Fig. 10). Damage is
concentrated at the center of the necked region so that the localization indicator becomes negative (2). This
also corresponds to the sharp slope change on the normalized load F =S0 -diameter reduction Dd=d0 curve
(F: force, S0 initial cross section, Dd: diameter variation, d0 initial diameter). The diameter reduction is
monitored in the initial symmetry plane. The highly damaged zone grows and leaves, behind its tips, an
elastically unloaded zone which correspond to the ¯at portion of the cup-cone (3) and forms a penny-
shaped crack. Ahead of this zone, two ``wings'' where the localization criterion is met, develop and become
larger as the central crack grows (3)±(5). Crack de¯ection starts with a relatively small angle. De¯ection
seems to be possible when the localization wings extend over 2 and more elements. At step (5) the wings
have grown suciently so that the cup-cone develops (6) and (7). Some secondary highly damaged regions,
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8275

Fig. 11. Predicted angle of localization h as a function of the element position along the radial direction and comparison with the actual
numerical crack orientation (dots).

which are no longer deforming, are left behind the main crack (Fig. 7). They could possibly correspond to
the bands observed on the surface of the specimens in Fig. 1. At some point, the diameter does not vary any
more. This occurs when the point of measure leaves the active plastic zone which is ahead of the highly
damaged band.
Fig. 11 shows the localization angles h expressed in the ®xed frame as a function of the element position
along the radial direction. The values correspond to the angles computed when min~n det~ n
Lc
~
n becomes
negative for the ®rst time. The graph is divided in di€erent segments 1±6. Rupture remains ¯at over 10
elements (seg. 2). This part can however be sub-divided. Below 5 elements localization occurs only on one
single row of Gauss points (seg. 1). In this region the localization angles are h ˆ 60°. From 5 to 10 ele-
ments (seg. 3), localization wings start to grow, and slightly more dispersed angles are obtained. Note that
the macroscopic fracture surface remains ¯at whereas locally an inclined crack path is predicted: this could
explain the rough fracture surface observed in the center of the specimens. Above 10 elements, an inclined
fracture path is formed. The crack ®rst moves towards the upper part of the specimen (seg. 4) corre-
sponding to a negative localization angle h  60°. 2 Segment 5 corresponds to the formation of two lo-
calization bands (step 5 in Fig. 10). After this point (seg. 6) the crack takes a new direction (going towards
the bottom of the specimen, i.e. positive values of h) with increasing values of h. Dots on Fig. 11 show the
actual numerical orientation of the crack showing a relative good agreement up to 30 elements. Above this
value, the localization analysis indicates that the crack should continuously turn toward the top of the
specimen (h > 90°) or abruptly change its path with h  10±20°. Numerically, the crack tends to become

2
As slightly di€erent parameters were used for the calculation of Fig. 11 and Fig. 10, crack paths do not exactly coincide. In Fig. 10,
the ®rst crack de¯ection corresponds to the h ˆ ‡60° branch. This clearly illustrates the non-uniqueness of the solution. Otherwise,
results of both calculations are consistent.
8276 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Fig. 12. Simulated fracture of round bars using the Gurson model with di€erent sets of parameters.

horizontal. This discrepancy could be due to the fact that the mesh is already highly deformed in the outer
region or that post-localization evolution could modify the crack orientation. Note that an abrupt crack
path change can sometimes be observed (e.g. with the Gurson model, q2 ˆ 1:15, Nh ˆ 80 in Fig. 12).

Remark 3. (E€ect of material heterogeneities). Following the method presented in Decamp et al. (1998) and
Devillers-Guerville et al. (1997) the mesh was divided in square regions of size d0 =80†  d0 =80†. In each
square an initial value of the porosity was randomly selected according to a uniform distribution with
1:1  10 4 6 f 6 1:9  10 4 . Starting from this level of heterogeneity it becomes possible to generate zig-
zagging crack paths which are similar to what was experimentally observed (Fig. 1(c)). The deformation
corresponding to the ®rst sharp drop of the load is reduced. On the other hand, the strain to failure (de®ned
by the vertical drop of the load) is increased as more plastic strain in the central region of the specimen is
generated due to the zigzag. These e€ects are limited in the present case and represent about 2% of the total
strain to failure.

Remark 4. (Mesh size e€ect). The mesh size e€ect observed in Section 4.1 is also related to the size of the
localization zone ahead of the highly damaged region: when it extends over 1 or 2 elements only, cup-cone
fracture will not occur. Many authors consider that mesh size is a material parameter which should be
adjusted and be kept constant when specimens of di€erent sizes are computed. As the size of the local-
ization zone scales with the size of the specimen, this means that cup-cone fracture should be more easily
observed in large specimens than in small ones.

Remark 5. (Results on notched bars). It is experimentally observed that fracture in notched bars (Fig. 1(g,h))
is essentially ¯at. A narrow shear lip is formed at the very end of rupture. FE calculations for a notch radius
equal to 4 mm are shown in Fig. 9 showing that a ¯at fracture is modeled using Nh ˆ 40. Once again, the
local localization angle h is smaller than p=2 but the localization ``wings'' remain too small to trigger de-
viation. For Nh ˆ 80 a cup±cone is formed. Similar results are obtained applying rate independent plas-
ticity.

4.3. Role of constitutive equations

4.3.1. Gurson model

As already mentioned, the Gurson model together with the most commonly used parameters for q1 and
q2 leads to an overestimation of the ductility. Three di€erent solutions were envisaged to ®t the ductility: (1)
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8277

use of the fH function, (2) adding strain controlled nucleation, (3) adjustment of q2 . In the ®rst case, the
adjusted parameters are: fc ˆ 0:005, d ˆ 3. For the second case, one considers that after a critical porosity
is reached a second population of inclusions (e.g. NbC) acts as a new source for void nucleation. The same
value for the critical porosity as in the ®rst case was taken. Above this value, An is constant and equal to 0:2.
In the third case, q2 ˆ 1:15 gives the best ®t. Numerical results using the Gurson model are shown in Fig.
12.
Using fH , ¯at fracture is always obtained for both cax8r or cax6r elements with 40 or 80 elements in
the cross section, with plastic or viscoplastic behavior. fH is a function whose derivative is discontinuous, so
that any localization indicator is also discontinuous. This implies that localization can occur at one Gauss
point whereas the neighborhood remains stable. It follows that the previously described ``localization
wings'' completely vanish leading to ¯at fracture although locally an inclined crack path is predicted. In the
case of strain controlled nucleation and although a discontinuous localization function is used, a zigzagging
crack path is obtained. This is due to the fact that damage is then enhanced in the inclined highly deformed
regions formed ahead of a penny-shaped crack.
The di€erence between the calculations carried out with the Rousselier model and the Gurson model
using q2 ˆ 1:15 are outlined in the following. (i) The localization angle predicted at the center of the bar is
higher and equal to  80°. (ii) The ¯at part on the cup-cone extends over more than 35% of the initial
radius (25% for the Rousselier model). (iii) The size of the localization ``wings'' are, for a given material
behavior, smaller than for the Rousselier model. Consequently, with Nh ˆ 40 and cax8r elements a
slightly zigzagging crack path is obtained. This zigzagging crack path is more pronounced with an elas-
toplastic behavior which is more susceptible to localization. Using a ®ner geometrical description
(Nh ˆ 80±±cax8r or cax6r) leads to a clear cup-cone formation. These numerical results are in agreement
with the simple analysis of Section 3.1 which indicates that the Rousselier model can lead more easily to
inclined fracture than the Gurson model.

Remark 6. E€ect of fR on the crack path. In the case of the Rousselier model, the value of the failure porosity
fR was chosen equal to 0.9 which corresponds to a very high level of damage (Section 2.3). Similar results
are obtained with smaller values down to 0.4±0.5. Below this value, fR acts as fH so that the development of
the localization zones is inhibited: this, again, leads to ¯at fracture. Similarly, using the localization cri-
terion as a rupture criterion at each Gauss point will probably lead to the same results.

4.3.2. Viscosity
Fig. 13 compares the macroscopic response and the development of localization in a round bar for which
the material is either plastic or viscoplastic.
In the case of plasticity, it can be seen that indicators based on tangent and consistent matrix give
very similar results in terms of zones where localization occurs. At the onset of the cup-cone formation,
the localization ``wings'' extend over about 10 elements. However based on the results of Section 3.4, it is
very likely that the actual localization zone is overestimated. In the case of viscoplasticity, zones where
min~n ~
n
Lc
~ _
n are negative correspond to zones where the normalized growth factor x= p=p† is larger than
1000 showing also a good correspondence between both criteria. These zones only extend over 5 elements.
Calculations were also carried out with a much higher viscosity (e.g. K ˆ 1000) or a higher strain rate
(e.g. 1000 s 1 ). The localization zones shrink and eventually vanish; this leads to ¯at fracture surfaces as
shown on Fig. 13. As noted in Section 2.1 changing the viscosity using the Rousselier model leads to a
higher porosity growth rate and a higher ration rH =r1 thus increasing the possibility of normal fracture.
The viscosity e€ect was therefore checked using the Gurson model with q2 ˆ 1:15 and cax6r elements. As
necking is more pronounced in that case the mesh had to be slightly redesigned to obtain approximately
square elements at the onset of rupture. Results are also shown in Fig. 13 indicating that cup-cone for-
mation is also suppressed in that case.
8278 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Fig. 13. E€ect of viscosity on the macroscopic behavior and on the cup-cone formation.

Fig. 14. Fracture path in plane strain specimens for di€erent models.

4.4. Plane strain specimens

Results obtained with plane strain specimens are consistent with those obtained on round bars. Some
typical results are shown in Fig. 14. Necking occurs, so that failure is initiated at the center of the specimen.
As it is easier to form localization bands in this case, slant fracture is obtained with the Gurson model for
Nh ˆ 40 except when using the fH function.

5. SummaryÐconcluding remarks

In this study, the formation of cup-cone and slant fracture has been analyzed using the FE method.
Constitutive equations are based on the Rousselier and Gurson models including viscosity and damage
nucleation. Indicators are also computed to detect zones where localization of deformation and damage can
occur.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8279

For plastic behavior, the localization indicator is based on Rice's analysis of bifurcation. For time de-
pendent plasticity, the indicator uses the perturbation analysis. It is also proposed to use the consistent
tangent matrix to derive a localization criterion similar to Rice's condition for bifurcation but which can be
used for both plastic and viscoplastic materials. Comparison of the three indicators shows that they give
consistent descriptions of the cup-cone development in bars.
Cup-cone formation can be analyzed considering the size of the zone (hereafter referred to as LZ) where
localization occurs and where the material is not yet highly damaged (f > 20%, 30%) or broken. After
necking, localization occurs ®rst at the center of the specimen. However the LZ remains small due to the
axisymmetric deformation state prevailing at the center of the specimen. As the central crack extends, the
LZ can grow so that the cup-cone can be formed. Based on the size of the LZ, most of the computed e€ects
can be interpreted:
Mesh size: When the mesh size is too coarse to capture the LZ, cup-cone cannot be formed. In addition,
assuming that the mesh size is a material characteristic parameter should imply that cup-cone is more likely
to appear in large specimens.
Viscosity: Viscosity has a regularizing e€ect which lead to the diminution of the LZ. Increasing viscosity
will lead to ¯at fracture.
Strain rate: For a given viscosity, increasing the strain rate also limits localization thus promoting
¯at fracture. This result is in agreement with creep experiments carried out by Kobayashi et al. (1998)
on notched aluminum bars. Under a high mean stress (30 MPa, high strain rate) specimens exhibit a ¯at
fracture surface whereas under low stress (16 MPa, low strain rate) cup-cone fracture is formed.
Use of fH : The fH function employed in the GTN model can cause the crack path to remain ¯at. As its
derivative with respect to porosity is discontinuous, it may suddenly induce localization as f reaches fc at one
single Gauss point whereas the surrounding Gauss points remain in a state far from instability. This process
inhibits the formation of a large LZ and therefore of crack deviation unless fc is large enough so that lo-
calization occurs for lower values of the porosity. This suggests also that it would be more appropriate to
adjust q1 and q2 with a high value of fc used for numerical purposes only. This solution was in fact adopted
by Gullerud et al. (2000); q1 and q2 can be adjusted from unit cell calculations as in Faleskog et al. (1998).
Nucleation: Strain controlled nucleation favors cup-cone fracture, as this mechanism induces damage in
the plastic ``wings'' formed ahead of the central penny shaped crack.
Similar conclusions can be drawn from the study of round notched bars. In plane strain specimens slant
fracture can be more easily obtained; once again using fH can produce ¯at fracture.
The comparison of constitutive models shows that cup-cone is formed more easily when employing the
Rousselier model that the Gurson model even when the fH function is not used. This is attributed to the
shape of the yield function; in the case of the Rousselier model, the presence of a vertex for a tensile hy-
drostatic stress state implies that the deformation rate tensor always keeps a shear component. This inhibits
pure opening mode fracture. It is interesting to note that calculations (Koplik and Needleman, 1988;
Brocks et al., 1995) on axisymmetric unit cells indicate that during void coalescence the macroscopic strain
rate tensor has its radial and hoop component equal to zero. Assuming isotropy and using the normality
rule implies that the yield surface is then a straight line of slope 3=2 in the req ±rkk =3 plane (Fig. 2)
(Thomason, 1985). Such a behavior can be accounted by the Rousselier model although it was not explicitly
designed for this purpose. On the other hand, the Gurson model does not represent this behavior; moreover
the use of fH tends to increase the volumetric part of the deformation. An other solution could be to use a
combination of the Gurson model (low porosity) and of the model proposed by Fleck et al. (1992) (high
porosity) which presents a vertex as done in Redanz and Tvergaard (1999). An alternative solution could
consist of using a combination of the Gurson model and of the straight line yield surface derived by
Thomason (1985). In Redanz and Tvergaard (1999) and Fleck et al. (1992) the transition porosity between
both model is ®xed; it could also be derived from micromechanical models of coalescence as in Thomason
(1985), Zhang and Niemi (1995) and Benzerga et al. (1999).
8280 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Acknowledgements

The authors would like to thank Prof. A. Pineau for providing the material of this study and Dr. S.
Forest and Dr. G. Barbier for commenting the manuscript. This work was performed during the sabbatical
leave of JB at GKSS which is acknowledged for ®nancial support and hospitality.

Appendix A. Calculation of the elastoplastic tangent matrix: Lt

The rate equations governing the evolution of the state variables are:
e_ e ˆ e_ 1 _
f †pv A:1†

f_ ˆ 1
2
_ : 1 ‡ An p_
f † pv A:2†

p_ ˆ consistency A:3†
The consistency condition /_ ˆ 0 yields:

1
p_ ˆ v : C : e_ with h ˆ 1 f †v : C : v ‡ H ‰1 f †2 v : 1 ‡ An ŠrH;f A:4†
h

H ˆ dR=dp is the hardening modulus. The derivative rH;f ˆ orH =of can be computed (in the case of the
Gurson model) as
orH 1 oU
ˆ A:5†
of v of
with v ˆ oU=orH . In the following, it will be assumed that h > 0. h ˆ 0 can be encountered for low po-
rosities and high stress triaxiality ratios. In this case, the constitutive equations cannot be integrated for a
prescribed strain rate tensor (snap back e€ect). Finally the elastoplastic tangent matrix is given by

1 f
Lp ˆ C C : v† v : C† A:6†
h

which is a special case of the more general form studied in Rudnicki and Rice (1975) for instance.

Appendix B. Perturbation analysis: calculation of Hp

In the case of viscoplasticity, the rate equations governing the evolution of the elastic strain and the
porosity are still given by Eqs. (A.1) and (A.2). The evolution of p is given by the viscoplastic ¯ow law:
e_ e ˆ e_ 1 f †Fv  Fe B:1†

f_ ˆ
2
1 f † v : 1 ‡ An †F  Ff B:2†

p_ ˆ F /†  Fp B:3†
The vector F introduced in Section 3.2, is then written as: F Z; e_ † ˆ Fe ; Ff ; Fp † Using the following nota-
tions,

ov o2 rH 1 o2 U oU 1 o2 U
Nˆ ˆ ˆ B:4†
or or2 v2 ororH or v or2
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8281

ov o2 rH 1 o2 U oU 1 o2 U
v;f ˆ ˆ ˆ 2 B:5†
of orof v orH of or v orof

dAn
A0n ˆ B:6†
dp

oF
F0 ˆ B:7†
o/
the partial derivative needed to compute Lp are given by:

oFe
ˆ 1 f †FN : C 1 f †F0 v : C† v
oee
oFe

ˆF v 1 f †v;f 1 f †F0 rH;f v
of
oFe
ˆ 1 f †F0 H v
op
oFf
ˆ 1 f †2 FN : C : 1 ‡ 1 f †2 v : 1 ‡ An †F0 v : C
oee
oFf
2
ˆ F 1 f † 1 f †v;f : 1 2v : 1 ‡ 1 f † v : 1 ‡ An †F0 rH;f
of
oFf 2
ˆ A0n F 1 f † v : 1 ‡ An †F0 H
op
oFp
ˆ F0 v : C
oee
oFp
ˆ F0 rH;f
of
oFp
ˆ F0 H
op

One also gets:

oF
ˆ I; 0; 0† B:8†
oe

Appendix C. Calculation of the consistent tangent matrix: Lc

The fully implicit integration of the model is obtained using a time discretization of the rate equations
(A.1)±(A.3) or B.3
Dee ˆ De 1 f †Dpv  FH
e C:1†
2
Df ˆ ‰ 1 f † v : 1 ‡ An ŠDp  FH
f C:2†
All quantities should be considered at the end of the time increment. The equation corresponding to Dp is
written as

0 ˆ rH R  FH
p C:3†
8282 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

in the case of plasticity. This condition expresses the fact that the material lies on the yield surface at the end
of the time increment. In the case of a viscoplastic material, the previous equation is replaced by
Dp ˆ F /†Dt  FH
p C:4†
The matrix I0 appearing in Eq. (29), is therefore the unity matrix in the case of viscoplasticity. In the case of
plasticity, its diagonal term corresponding to p is set to zero. The function FH is given by FH H H
e ; Ff ; Fp †.
H H
The derivative oF =oe and oF =oZ can then be calculated noting that oDx=ox ˆ 1, 8 x. One gets:
oFH
ˆ I; 0; 0† C:5†
oe
The derivatives of FH are:
oFH
e
ˆ 1 f †DpN : C
oee
oFHe


ˆ Dp v 1 f †v;f
of
oFHe
ˆ 1 f †v
op
oFHf 2
ˆ 1 f † DpN : C : 1
oee
oFHf


ˆ Dp 1 f† 1 f †v;f : 1 2v : 1
of
oFH 2
f
ˆ 1 f † v : 1 ‡ An ‡ A0n Dp
op
(
oFHp F0 Dtv : C vp
ˆ
oee v:C p
(
oFHp F0 DtrH;f vp
ˆ
of rH;f p
(
oFHp F0 DtH vp
ˆ
op H p

References

A B A Q U S , 1998. Version 5.8. Tech. Rep., H.K.S. Inc., Pawtucket, USA.

Anand, L., Kim, K., Shawki, T., 1987. Onset of shear localization in viscoplastic solids. J. Mech. Phys. Solids 35 (4), 407±429.
Anand, L., Spitzig, W., 1980. Initiation of localized shear bands in plane strain. J. Mech. Phys. Solids 28, 113±128.
Aravas, N., 1987. On the numerical integration of a class of pressure-dependent plasticity models. Int. J. Num. Math. Engng. 24, 1395±
1416.
Barbier, G., 1999. Localisation et instabilite dans les materiaux elastoplastiques endommageables. Ph.D. Thesis, Universite Paris VI.
Barbier, G., Benallal, A., Cano, V., 1998. Relation theorique entre la methode de perturbation lineaire et l'analyse de bifurcation pour
la prediction de la localisation des deformations. Comptes rendus de l'Academie des Sciences 326 (3), 153±158.
Becker, R., Needleman, A., 1986. E€ect of yield surface curvature on necking and failure in porous plastic solids. J. Appl. Mech. 53,
491±499.
 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284 8283

Becker, R., Needleman, A., Richmond, O., Tvergaard, V., 1988. Void growth and failure in notched bars. J. Mech. Phys. Solids 36,
317±351.
Benzerga, A., Besson, J., Pineau, A., 1999. Coalescence-controlled anisotropic ductile fracture. J. Engng. Mater. Engng. 121, 121±229.
Besson, J., Foerch, R., 1997. Large scale object-oriented ®nite element code design. Comput. Meth. Appl. Mech. Engng. 142, 165±187.
Billardon, R., Doghri, I., 1989. Prediction of macro-crack initiation by damage localization. C.R. Acad. Sci. Paris 308 (Serie II),
347±352.
Bonora, N., 1999. Identi®cation and measurement of ductile damage parameters. J. Strain Anal. Engng. Des. 34 (6), 463±478.
Borre, G., Maier, G., 1989. On linear versus nonlinear ¯ow rules in strain localization analysis. Meccanica 24, 36±41.
Brocks, W., Sun, D.Z., H onig, A., 1995. Veri®cation of the transferability of micromechanical parameters by cell model calculations
with visco-plastic materials. Int. J. Plast. 11, 971±998.
Chaboche, J., Cailletaud, G., 1996. Integration methods for complex plastic constitutive equations. Comput. Meth. Appl. Mech.
Engng. 133 (1±2), 125±155.
Chu, C., Needleman, A., 1980. Void nucleation e€ects in biaxially stretched sheets. J. Engng. Mater. Technol. 102, 249±256.
Decamp, K., Bauvineau, L., Besson, J., Pineau, A., 1998. Size and geometry e€ects on ductile rupture of notched bars in a C±Mn steel:
experiments and modelling. Int. J. Fract. 88 (1), 1±18.
Devillers-Guerville, L., Besson, J., Pineau, A., 1997. Notch fracture toughness of a cast duplex stainless steel: modelling of
experimental scatter and size e€ect. Nucl. Engng. Des. 168, 211±225.
Doghri, I., Billardon, R., 1995. Investigation of localization due to damage in elasto-plastic materials. Mech. Mater. 19, 129±149.
Faleskog, J., Gao, X., Shih, C., 1998. Cell model for nonlinear fracture analysis. I. Micromechanics calibration. Int. J. Fract. 89,
355±373.
Fleck, N., Kuhn, L., McMeecking, R., 1992. Yielding of metal powder bonded by isolated contacts. J. Mech. Phys. Solids 40 (5),
1139±1162.
Foerch, R., Besson, J., Cailletaud, G., Pilvin, P., 1997. Polymorphic constitutive equations in ®nite element codes. Comput. Meth.
Appl. Mech. Engng. 141, 355±372.
Fressengeas, C., Molinari, A., 1985. Inertia and thermal e€ects on the localization of plastic ¯ow. Acta Metall. 33 (3), 387±396.
Grange, M., Besson, J., Andrieu, E., 2000. An anisotropic Gurson model to represent the ductile rupture of hydrided Zircaloy-4 sheets.
Int. J. Fract. 105 (3), 273±293.
Gullerud, A., Gao, X., Dodds Jr., R., Haj-Ali, R., 2000. Simulation of ductile crack growth using computational cells: numerical
aspects. Engng. Fract. Mech. 66, 65±92.
Gurson, A., 1977. Continuum theory of ductile rupture by void nucleation and growth. Part IÐYield criteria and ¯ow rules for porous
ductile media. J. Engng. Mater. Technol. 99, 2±15.
Hughes, T., 1980. Generalization of selective integration procedures to anisotropic and non linear media. Int. J. Numer. Meth. Engng.
15, 1413±1418.
Hughes, T., Winget, J., 1980. Finite rotation e€ects in numerical integration of rate constitutive-equations arising in large deformation
analysis. Int. J. Numer. Meth. Engng. 15 (12), 1862±1867.
Kobayashi, K., Imada, H., Majima, T., 1998. Nucleation and growth of creep voids in circumferentially notched specimens. JSME Int.
J. Series A: Solids Mech. Mater. Engng. 41 (2), 218±224.
Koplik, J., Needleman, A., 1988. Void growth and coalescence in porous plastic solids. Int. J. Solids Struct. 24, 835±853.
Ladeveze, P., 1980. Sur la theorie de la plasticite en grandes deformations. Tech. Rep., Rapport interne No. 9, LMT, ENS Cachan.
Leblond, J., Perrin, G., Devaux, J., 1994. Bifurcation e€ects in ductile metals with nonlocal damage. J. Appl. Mech. 61, 236±242.
Mackenzie, A., Hancock, J., Brown, D., 1977. On the in¯uence of state of stress on ductile failure initiation in high strength steels.
Engng. Fract. Mech. 9, 167±188.
Mear, M., Hutchinson, J., 1985. In¯uence of yield surface curvature on ¯ow localization in dilatant plasticity. Mech. Mater. 4,
395±407.
Needleman, A., Rice, J., 1978. Limits to ductility set by plastic ¯ow localization. In: Koistinen, D.P. (Ed. ), Mechanics of Sheet Metal
Forming. Plenum Publishing Corporation, New York.
Redanz, P., Tvergaard, V., 1999. Analysis of shear band instabilities in sintered metals. Int. J. Solids Struct. 36, 3661±3676.
Rice, J., 1976. The localisation of plastic deformation. In: Koiter, W. (Ed.), Proceedings of the 14th International Conference on
Theoretical and Applied Mechnics, Delft. North-Holland, Amsterdam.
Rice, J., Rudnicki, J., 1980. A note on some features of the theory of localization of deformation. Int. J. Solids Struct. 16, 597±605.
Rivalin, F., 1998. Developpement d'aciers pour gazoducs a haute limite d'elasticite et tenacite elevee: Mecanique et mecanismes de la
rupture ductile  
a grande vitesse. Ph.D. Thesis, Ecole des Mines de Paris.
Rousselier, G., 1987. Ductile fracture models and their potential in local approach of fracture. Nucl. Engng. Des. 105, 97±111.
Rousselier, G., 1991. Application de l'analyse de stabilite d'une perturbation a la localisation de la deformation dans un materiau
dilatable adoucissant. C.R. Acad. Sci. Paris 313 (Serie II), 1367±1373.
Rousselier, G., 1995a. Localisation de la deformation dans des materiaux plastiques et viscoplastiques, dilatables, soumis a des
sollicitations thermomecaniques et dynamiques tridimensionnelles. C.R. Acad. Sci. Paris 320 (Serie IIb), 265±270.
8284 J. Besson et al. / International Journal of Solids and Structures 38 (2001) 8259±8284

Rousselier, G., 1995b. Stabilite locale et modes de rupture ductile. C.R. Acad. Sci. Paris 320 (Serie IIb), 69±75.
Rudnicki, J., Rice, J., 1975. Conditions for the localization of deformation in pressure-sensitive dilatant materials. J. Mech. Phys. Sol.
23, 371±394.
Saje, M., Pan, J., Needleman, A., 1982. Void nucleation e€ects on shear localization in porous plastic solids. Int. J. Fract. 19, 163±182.
Simo, J., Taylor, R., 1985. Consistent tangent operators for rate-independent elastoplasticity. Comput. Meth. Appl. Mech. Engng. 48,
101±118.
Thomason, P.F., 1985. Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous
solids. Acta Metall. 33 (6), 1079±1085.
Tvergaard, V., 1982. In¯uence of void nucleation on ductile shear fracture at a free surface. J. Mech. Phys. Solids 30 (6), 399±425.
Tvergaard, V., 1987. E€ect of yield surface curvature and void nucleation on plastic ¯ow localization. J. Mech. Phys. Solids 35 (1),
43±60.
Tvergaard, V., 1990. Material failure by void growth to coalescence. Adv. Appl. Mech. 27, 83±151.
Tvergaard, V., Needleman, A., 1984. Analysis of cup-cone fracture in a round tensile bar. Acta Metall. 32, 157±169.
Tvergaard, V., van der Giessen, E., 1991. E€ect of plastic spin on localization predictions for a porous ductile material. J. Mech. Phys.
Solids 39 (6), 763±781.
Xia, L., Shih, C., Hutchinson, J., 1995. A computational approach to ductile crack growth under large scale yielding conditions.
J. Mech. Phys. Solids 43 (3), 389±413.
Yamamoto, H., 1978. Conditions for shear band localization in the ductile fracture of void-containing materials. Int. J. Fract. 14 (4),
347±365.
Zhang, Z., 1995. On the accuracies of numerical integration algorithms for Gurson-based pressure-dependent elastoplastic constitutive
models. Comput. Meth. Appl. Mech. Engng. 121, 15±28.
Zhang, Z., Niemi, E., 1995. A new failure criterion for the Gurson±Tvergaard dilational constitutive model. Int. J. Fract. 70, 321±334.