J Besson

D Steglich

W Brocks

Modeling of crack growth in round bars and plane strain specimens

Keywords: Localization, cup-cone fracture, slant fracture Rousselier, Gurson

 and by Gurson (Tvergaard and Needleman, 1984) are used. The analysis takes into account viscoplasticity and void nucleation. Different 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 flat fracture. The crack path (flat or cupcone/slant) can be correlated to the size of the localization zone which is formed ahead of the central penny shaped crack.

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 [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF] and phenomenologically extended by [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] (so called GTN model) has been most frequently used. An approach based on continuum damage mechanics (CDM) and thermodynamics has also been proposed by [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF]. 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 modified 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. [START_REF] Xia | Ductile crack growth -I. A numerical study using computational cells with microstructurally-based length scales[END_REF]). They have also been used to model fracture of small uncracked laboratory tests samples such as smooth and notched round tensile bar [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF][START_REF] Becker | Void growth and failure in notched bars[END_REF] or plane strain specimens [START_REF] Becker | Effect of yield surface curvature on necking and failure in porous plastic solids[END_REF][START_REF] Leblond | Bifurcation effects in ductile metals with nonlocal damage[END_REF]. 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 (figure 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 first numerically analyzed in [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] 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 discontinuity surface) in elastoplastic solids have been described by [START_REF] Rice | The localisation of plastic deformation[END_REF] ; the specific case of dilatant pressure sensitive materials has been investigated in [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF]. Using a simplified model consisting of a solid body having a plane of imperfection, the susceptibility to localization of the Gurson model has been first studied by [START_REF] Yamamoto | Conditions for shear band localization in the ductile fracture of void-containing materials[END_REF] with an emphasis on the effect of heterogeneities. This model can also be used to study the evolution of the band beyond the onset of localization [START_REF] Tvergaard | Ductile failure by cavity nucleation between larger voids[END_REF]. The role of void nucleation was studied in [START_REF] Saje | Void nucleation effects on shear localization in porous plastic solids[END_REF] 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 [START_REF] Mear | Influence of yield surface curvature on flow localization in dilatant plasticity[END_REF], together with nucleation in [START_REF] Tvergaard | Effect of yield surface curvature and void nucleation on plastic flow localization[END_REF] showing a decrease of ductility with increasing kinematic hardening. The choice of the corotational stress rate was studied in [START_REF] Tvergaard | Effect of plastic spin on localization predictions for a porous ductile material[END_REF]. In particular it is shown that localization remains unaffected 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 [START_REF] Fressengeas | Inertia and thermal effects on the localization of plastic flow[END_REF][START_REF] Anand | Onset of shear localization in viscoplastic solids[END_REF] incorporating heating induced by plastic deformation, heat diffusion and inertia effects. The case of rigid-plastic dilatant solids has been treated in [START_REF] Rousselier | Application de l'analyse de stabilité d'une perturbation à la localisation de la déformation dans un matériau dilatable adoucissant[END_REF](Rousselier, , 1995b,a) ,a) 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 infinite medium in which a band-like discontinuity appears. In actual structures, this situation is indeed never met. In this case, [START_REF] Billardon | Prediction of macro-crack initiation by damage localization[END_REF]; [START_REF] Doghri | Numerical implementation and analysis of a class of metal plasticity models coupled with ductile damage[END_REF] proposed to compute Rice's condition for localization during the Finite Element 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.

In this paper, the Finite Element (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 experimental results obtained on a modern steel containing very low sulphur and phosphorus contents. In section 3, different localization indicators are presented and compared. In section 4, FE simulations of cup-cone and slant fracture are carried out. The influence of the following parameters is studied: (i) mesh refinement and element formulation, (ii) plastic and viscoplastic flow, (iii) constitutive model: Gurson and Rousselier models, effect of void nucleation. Following Billardon and Doghri, localization indicators are also evaluated 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.

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 denote matrices (M). Dots and colons are used to indicate the usual contracted products:

a. b = i a i b i , a : b = i,j a ij b ij , (A : a) ij = k,l A ijkl a kl , (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 are denoted by 1, 1 and I. The notation n.L. n represents the second order tensor A such that

A ij = k,l n k L kijl n l .
2 Material models

Constitutive equations

Models for porous materials proposed by [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF] and Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF][START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF] are used in this study. The first 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 unified framework. The plastic (or viscoplastic) flow potential φ is then written as:

φ = σ -R(p) (1)
where R is the yield stress of the undamaged material (matrix) and p an effective plastic strain representative of the matrix hardening. σ is an effective scalar stress which is a function of both the macroscopic stress tensor σ and the porosity. σ is defined by the following equations:

Gurson Φ = σ 2 eq σ 2 + 2q 1 f cosh q 2 2 σ kk σ -1 -q 2 1 f 2 def. σ = 0 (2) Rousselier Φ = σ eq (1 -f )σ + σ 1 σ f D exp σ kk 3(1 -f )σ 1 -1 def. σ = 0 (3)
where σ eq is the von Mises equivalent stress and σ kk the trace (tr) of the stress tensor. q 1 , q 2 , D and σ 1 are material coefficients which are assumed to be constant. f is a function of the porosity f which was introduced on a purely phenomenological basis to represent void coalescence [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]. In the case of the Gurson model (eq. 2) the definition of σ is implicit whereas it is explicit in the case of the Rousselier model (eq. 3).

In the case of a plastic material, one has φ = 0 whereas in the viscoplastic case φ > 0. The irreversible deformation rate εp is obtained assuming the normality rule, so that:

εp = (1 -f ) ṗ ∂φ ∂σ = (1 -f ) ṗ∂σ ∂σ = (1 -f ) ṗυ (4)
where υ is the normal to the flow potential. The total strain rate ε is expressed as εe + εp where εe is the elastic strain rate. Elastic strains and stresses are related by: σ = C : ε e . As the void volume fraction is usually small, the elasticity tensor C is assumed to be constant. υ can be computed noting that for a fixed porosity, a variation of σ induces a variation of σ such that Φ remains equal to zero (this derivation must be used in the case of the Gurson model, where the definition of σ is implicit). Therefore:

δΦ = ∂Φ ∂σ : δσ + ∂Φ ∂σ δσ = 0 ⇒ υ = ∂σ ∂σ = - ∂Φ ∂σ -1 ∂Φ ∂σ (5)
The plastic multiplier ṗ can be related to the deformation rate tensor. In the case of the Gurson model, one notes that (∂σ /∂σ) : σ = σ , so that:

εp : σ = (1 -f ) ṗσ (6)
In the case of the Rousselier model, the previous relation does not hold and one has:

ṗ = εeq ≡ 2 3 ė : ė (7)
where ė is the deviator of the strain rate tensor. In this case, ṗ corresponds to the von Mises equivalent strain rate. ṗ can be computed writing the consistency condition φ = 0 in the case of plasticity or using the viscous flow law of the dense material in the case of viscoplasticity:

ṗ = F(φ) = F(σ -R) (8)
In the following the creep function F will be chosen as a Norton law, i.e. F(φ) = φ/K n , where K and n are material parameters and where . is such than x = x if x > 0, and x = 0 otherwise.

The evolution of the porosity is given by mass conservation modified to account for strain controlled void nucleation [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF]:

ḟ = (1 -f )tr εp + A n ṗ = (1 -f ) 2 υ : 1 + A n ṗ (9)
A n is a material function used to represent nucleation. In the following, it will be assumed that it depends on p only.

Remark: In the case of the viscoplastic Rousselier model [START_REF] Barbier | Localisation et instabilité dans les matériaux élastoplastiques endommageables[END_REF], the ratio ḟ / ṗ depends on ṗ so that porosity grows faster for higher strain rate. This is due to the fact that σ 1 is taken to be rate independent so that the ratio σ kk /σ 1 will increase with strain rate. This dependence is not observed for the Gurson model as Φ can be expressed as a function of σ/σ only (in the case of the Rousselier model Φ is a function of both σ/σ and σ/σ 1 ). A solution could be to modify the Rousselier model as follows:

Φ = σ eq (1 -f )σ + 2 3 f D exp q R 2 σ kk (1 -f )σ -1 def. σ = 0 ( 10
)
where q R 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 σ 1 and D adjusted for the experimental strain rate. Parametric studies involving a strong effect of viscosity will be carried out using the Gurson model only.

Comparison: Gurson-Rousselier models

At this point, it is interesting to outline some differences between the Gurson and Rousselier models. Figure 2 compares both yield surfaces in the σ eq -σ kk plane in the case of tensile stress states (σ kk > 0). Under pure shear (σ kk = 0), damage is still generated in the case of the Rousselier model (as the normal to the yield surface does not coincide with the σ eq axis) whereas, in absence of nucleation, the Gurson model does not lead to damage growth. Under pure hydrostatic stress states (σ eq = 0), the Rousselier yield surface has a vertex which implies that at high stress triaxiality ratios (τ = 1 3 σ kk /σ eq ) the plastic deformation tensor always keeps a non-zero shear component. Note that the model proposed in [START_REF] Fleck | Yielding of metal powder bonded by isolated contacts[END_REF] for plastic metal powders has the same property.

The ratio of volumetric to shear deformation rate (strain rate triaxiality ratio) τ ε = 1 3 εkk / εeq is given for both models by:

Gurson τ ε = 1 2 q 1 q 2 f sinh q 2 2 σ kk σ σ σ eq (11)
Rousselier

τ ε = 1 3 f D exp σ kk 3(1 -f )σ 1 (12)
Taking the limits of the strain rate triaxiality ratio for σ eq → 0 such that the yield condition is met, one gets 1 :

Gurson lim

σeq→0 τ ε = q 1 f 2 3 2 q 2 σ σ eq = +∞ (13)
Rousselier lim

σeq→0 τ ε = 1 3 σ σ 1 (14)

Material parameters

The material parameters used in this study correspond to a X70 HSLA (High Strength Low Alloyed) ferritic-pearlitic steel [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF]. Due to its very low sulphur 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 µm), composed of two phases combining calcium sulfide (CaS) and aluminum (Al 2 O 3 ) or magnesium 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 f 0 .

Due to the low inclusion volume fraction, the effect 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 [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF]. In the following, it will however be assumed that the material is isotropic in 1 Using the modified Rousselier model (equation 10), the strain rate triaxiality ratio is equal to 1 3 Df q R exp(q R σ kk /2(1 -f )σ ). The limit for σ eq → 0 is equal to q R /2.

order to simplify the FE calculations. Elastic properties are: Young's modulus E = 210 GPa, Poisson's ratio ν = 0.3. Plastic hardening is described using a simple power law relationship:

R(p) = K (p + ε 0 ) n (15)
with K = 795 MPa, ε 0 = 0.002 and n = 0.13. Parameters of the Norton law are: K = 55 MPa.s 1/n , n = 5. Tests carried out at different temperatures have also shown that the behavior is almost unaffected 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 [START_REF] Mackenzie | On the influence of state of stress on ductile failure initiation in high strength steels[END_REF][START_REF] Decamp | Size and geometry effects on ductile rupture of notched bars in a C-Mn steel: Experiments and modelling[END_REF] as well as plane strain specimens [START_REF] Anand | Initiation of localized shear bands in plane strain[END_REF]. 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. Figure 1 shows examples of fracture surfaces obtained in round bars and plane strain specimens. Note the anisotropic deformation of bars. An example of cup-cone formation in another isotropic material is also given in figure 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 (σ 1 and D for the Rousselier model ; q 1 , q 2 and f 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 [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF] are D = 2 and σ 1 = 1 3 (R e + R m ) = 321 MPa, where R m is the maximum engineering stress and R e the yield limit. These values lead however to an underestimation of the ductilities. Adjusted values are equal to D = 1.4 and σ 1 = 450 MPa. In the case of the Gurson model, many studies published in the literature use q 1 = 1.5 and q 2 = 1.0 whereas the function f is simply defined as follows:

f =    f if f < f c f c + δ (f -f c ) if f > f c ( 16
)
where δ > 1 and f c (critical porosity at which coalescence starts) have to be adjusted. Using f = f leads to an overestimation of the ductilities. In section 4.3 different solutions (use of f , use of nucleation, modification of q 2 ) will be envisaged in order to adjust the experimental data.

Failure occurs for the Gurson model when f = 1/q 1 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 f R . In the following f R = 0.9 will be used.

Localization indicators

In a infinite 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 displacement jump across the band whose direction is denoted g (figure 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 different 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.

Plastic materials: bifurcation analysis

In the case of elastoplastic materials, the incremental constitutive equation can be expressed as:

σ = L t : ε (17)
where L t is the elastoplastic tangent matrix. The calculation of L t is detailed in appendix A. Writing the continuity of displacements and the stress equilibrium, it can be shown [START_REF] Rice | The localisation of plastic deformation[END_REF][START_REF] Rice | The mechanics of earthquake rupture[END_REF] 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:

∃ n, det(A t ( n)) = 0 with A t ( n) = n.L t . n (19)
g is then the eigenvector of A t ( 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 A t ( n) < 0 [START_REF] Rice | The mechanics of earthquake rupture[END_REF][START_REF] Borré | On linear versus nonlinear flow rules in strain localization analysis[END_REF]. Equation 19 can be modified 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: [START_REF] Rice | The mechanics of earthquake rupture[END_REF][START_REF] Mear | Influence of yield surface curvature on flow localization in dilatant plasticity[END_REF].

∃ n, det(A t ( n) + R) = 0 with 2R = -n⊗( n.σ)+( n.σ)⊗ n+( n.σ. n)1-σ (
The bifurcation criterion was implemented as a post-processor of the FE calculations. In practice, the condition det A t = 0 is never exactly met.

Localization will be defined as occurring when det A t < 0 for the first time.

n is then defined as the vector which minimizes det A t ( 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 A t ( n) is then computed ; it coincides with g when det A t becomes zero. In the case of the specimens and constitutive equations investigated in this work, accounting for the Jauman rate modifies 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 flow direction, expressed in its eigen-coordinate system, is supposed to be of the following form:

υ ∝        1 0 0 0 u 0 0 0 t        with t ≤ u ≤ 1 (20)
This corresponds to a general tensile situation where minor stresses can be negative. Assuming that t ≤ u ≤ 1, the vectors n and g lie in the x 1 -x 3 plane provided dilatancy is not too large [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF][START_REF] Yamamoto | Conditions for shear band localization in the ductile fracture of void-containing materials[END_REF]. This localization plane is always assumed in the following. n and the angle ψ between g and n (fig. 3) are given by:

n =   1 + νu 1 -t , 0, ± -t -νu 1 -t   cos ψ = 1 + 2νu + t 1 -t (21)
The angle ψ characterizes the type of failure which varies from pure opening fracture ψ = 0 to pure tangential fracture (ψ = π 2 ). The results are summarized in fig. 3 in the t-u plane for -2 < u < 1. The case t + νu > 0 (shaded area) corresponds to pure opening fracture ( n // g), the normal to the band being aligned with the principal strain direction. The case 1+νu < 0 corresponds, for ν = 0.3, to negative volume changes (crosshatched region) and is not relevant in the present study. Constant band orientation (constant θ) or constant fracture mode (constant ψ) are represented by straight lines in the previous diagram. Location of constant strain rate triaxiality (τ ε = 0, 1 4 , 1 2 , 1, 2) are also indicated on the diagram (grey lines). In particular, it can be seen that for τ e < 1 2 normal separation is never possible. As noted in section 2.2, τ e lies between 0 and 1 3 σ /σ 1 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 A t ( n) = 0 can then be rewritten using the previous solution for n as:

H 1 -f =      -υ 2 11 E u 2 if t + νu < 0 -υ 2 11 E u 2 + 2νut + t 2 1 -ν 2 if t + νu ≥ 0 (22) with H = H -[(1 -f ) 2 υ 11 (1 + u + t) + A n ]σ ,f
where H is the plastic hardening modulus and σ ,f = ∂σ /∂f . H incorporates both plastic hardening and softening due to porosity growth. The states most resistant to localization are those of axisymmetric extension u = t whereas localization is much easier for plane strain u = 0 [START_REF] Needleman | Limits to ductility set by plastic flow localization[END_REF].

Viscoplastic materials: linear perturbation analysis

The previous bifurcation analysis is however only valid for materials for which the elastoplastic tangent matrix is defined. 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 field which is supposed to be of the following form [START_REF] Fressengeas | Inertia and thermal effects on the localization of plastic flow[END_REF][START_REF] Anand | Onset of shear localization in viscoplastic solids[END_REF][START_REF] Rousselier | Application de l'analyse de stabilité d'une perturbation à la localisation de la déformation dans un matériau dilatable adoucissant[END_REF]Rousselier, , 1995b;;[START_REF] Barbier | Relation théorique entre la méthode de perturbation linéaire et l'analyse de bifurcation pour la prédiction de la localisation des déformations[END_REF][START_REF] Barbier | Localisation et instabilité dans les matériaux élastoplastiques endommageables[END_REF]:

g = δ u exp(iq x. n + ωt) (23)
ω 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 diffusion, dynamic loading or non-local constitutive equations. It will not be considered hereafter.

Following the generic treatment of the problem proposed in [START_REF] Barbier | Relation théorique entre la méthode de perturbation linéaire et l'analyse de bifurcation pour la prédiction de la localisation des déformations[END_REF], the material is characterized by internal variables denoted Z = (ε e , 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 differential equations:

Ż = F(Z, ε) (24)
which, applying a perturbation, leads to:

δ Ż = ∂F ∂Z .δZ + ∂F ∂ ε .δ ε (25)
The perturbation δZ can be estimated from the perturbed rate as [START_REF] Anand | Onset of shear localization in viscoplastic solids[END_REF]: δZ = δ Ż/ω. The perturbed state variables are then related to the perturbed strain by:

δZ = I - 1 ω ∂F ∂Z -1 . ∂F ∂ ε .δε = H p (ω).δε (26)
The calculation of H p is detailed in appendix B for the viscoplastic case. As shown in [START_REF] Barbier | Relation théorique entre la méthode de perturbation linéaire et l'analyse de bifurcation pour la prédiction de la localisation des déformations[END_REF], the plastic case coincides with Rice analysis for ω → +∞. The sub-matrix of H p relating the perturbed total strain to the perturbed elastic deformation is denoted H p . The perturbation on the stresses is therefore computed as:

δσ = C : H p (ω) : δε = L p (ω) : δε (27)
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 A p ( n, ω) = 0 with A p ( n, ω) = n.L p (ω). n (28)
g (or δ u) is the associated eigenvector. The criterion for localization is then applied by finding the maximum value of ω for which a vector n exists which verifies det A p ( n, ω) = 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: ω ṗ/p.

Consistent tangent matrix as localization indicator

In the framework of the FEM, material constitutive equations can be integrated using an implicit integration scheme [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF][START_REF] Chaboche | Unified cyclic viscoplastic constitutive equations : development, capabilities and thermodynamic framework[END_REF][START_REF] Aravas | On the numerical integration of a class of pressure-dependent plasticity models[END_REF][START_REF] Foerch | Polymorphic constitutive equations in finite element codes[END_REF]. The increment of the state variable ∆Z over a finite time step ∆t is given by the following set of implicit (with respect to ∆Z) nonlinear equations:

I 0 .∆Z = F (Z, ∆t, ∆ε) (29)
In the following, only fully implicit integration will be used so that Z = Z 0 + ∆Z where Z 0 is the known set of state variables at the beginning of the time increment. The matrix I 0 , defined 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 ∂F /∂Z (see appendix C). Any infinitesimal variation of the deformation increment will induce a variation of the solution state variables such that eq. 29 is still satisfied, so that:

I 0 .δZ = ∂F ∂Z .δZ + ∂F ∂ε .δε (30)
Solving the previous equation for δZ gives:

δZ = I 0 - ∂F ∂Z -1 . ∂F ∂ε .δε = H c .δε (31)
The calculation of H c is detailed in appendix C. Using the same arguments as in 3.2, the variation of the stresses is related to the variation of the total strain by: δσ = C :

H c : δε = L c : δε (32)
where H c is the sub-matrix of H c relating δε and δε e . L c corresponds to the so-called "consistent tangent matrix". If a small variation of stresses δσ develops in a planar band, the mechanical equilibrium of the band requires that: δσ. n = 0. The corresponding variation of the deformation must be of the form given by equation 18 to ensure the continuity of displacement. Using the equilibrium condition together with equation 32, a non null vector g exists if:

∃ n, det(A c ( n)) = 0 with A c ( n) = n.L c . n (33)
This condition is similar to Rice localization criterion (equation 19) with L t being replaced by L c . In the following, 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.

Comparison of the different localization indicators

In this section, the three previously described localization indicators will be compared for plane strain conditions using the Gurson model with q 2 = 1.15 and f = f . Otherwise specified, material parameters are those given in section 2.3.

Figure 4 compares, in the case of an elastoplastic material behavior, the localization indicator based on the 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 D e corresponding to a purely elastic behavior:

D e = det n.C. n = 1 4 (1 -ν)E 3 (1 -2ν)(1 + ν) 3 ∀ n (34)
It can be seen that both indicators essentially give similar results for increments of deformation ∆ε up to 0.5%. A slight difference is observed for ∆ε = 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. Figure 6 compares, in the case of elastoviscoplastic material behavior, the localization indicator based on the consistent matrix for different 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 ε → 0. As shown in table 2, the porosity at the onset of localization (defined by min n det n.L c . n = 0) increases with increasing K ; this also corresponds to an increase of θ > π 4 and to a decrease of ψ < π 2 : the normal component of the fracture mode increases. For K → 0 (or ε → 0) the viscoplastic analysis (porosity, angles θ and ψ, strain) coincides with the plastic case. As in the plastic case, the deformation increment influences the consistent tangent matrix and consequently the time step for which min n det n.L c . n starts to be negative. This effect remains limited for plasticity (figure 4) but increases with increasing viscosity. For example, for K = 50 the corresponding deformation is equal to 2.09 for ∆ε = 0.5%, 1.78 for ∆ε = 1.0% and 1.56 for ∆ε = 2.0%. This corresponds to the fact that non-converging global time steps can be divided into converging sub-steps.

In order to further evaluate the different 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 field f 0 = 1.5 10 -4 ± 1.0 10 -6 . The calculations were postprocessed in order to determine: (i) the strain ε D at which min n n.L c . n becomes negative at at least one Gauss point in the structure, (ii) the strain ε B at which a continuous elastically deforming band is formed. Using the perturbation analysis of a uniform volume element, the growth factors ωD and ωB corresponding these strains can be computed. Finally, ε L is defined as the strain at which min n n.L c . n < 0 in a uniformly deforming element. Results are gathered in table 3. For plasticity and low values of K one gets ε D ≈ ε L < ε B . This shows that the indicators underestimate the actual numerical onset of localization. Similar results are obtained for plastic behavior using equation 19. On the other hand, for high values of K, one gets ε D ≈ ε B < ε L . 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 ωD and ωB are decreasing functions of K and onset of localization does not correspond to some "critical value". As soon as localization starts the ω field becomes strongly heterogeneous with some very high values (typically over 10 4 ). 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.

Finite Element simulations -Discussion

Finite element simulations were performed using the FE softwares ABAQUS (ABAQUS, 1998) and Zébulon [START_REF] Besson | Large scale object-oriented finite element code design[END_REF][START_REF] Foerch | Polymorphic constitutive equations in finite element codes[END_REF]. The FE implementation of both Gurson and Rousselier models in ABAQUS follows the method proposed in [START_REF] Aravas | On the numerical integration of a class of pressure-dependent plasticity models[END_REF] The method used in Zébulon is detailed in appendix C. Both methods use a fully implicit integration scheme. It was shown by [START_REF] Zhang | A class of generalized midpoint algorithms for the gurson-tvergaard material model[END_REF] that more accurate integration can be obtained using a semi-implicit scheme with ζ = 0.75-0.85 (ζ = 0 fully explicitly, ζ = 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 hand, Zébulon uses the full elastic strain tensor. Although less numerically efficient, this makes it possible to easily extend the method to plastically anisotropic materials [START_REF] Grange | An anisotropic Gurson model to represent the ductile rupture of hydrided Zircaloy-4 sheets[END_REF] 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 (ABAQUS [START_REF] Hughes | Finite rotation effects in numerical integration of rate constitutive-equations arising in large deformation analysis[END_REF], Zébulon [START_REF] Ladevèze | Sur la théorie de la plasticité en grandes déformations[END_REF]). The localization analysis was implemented in Zébulon only. It is made in the rotated material frame, defined using the Jauman rate for the stresses, for each Gauss point using eq. 19, 28 or 33. Vectors n and g can then be expressed in the fixed reference frame. In the following, calculations obtained using Zébulon are only presented. Very similar results in terms of mesh dependence, occurrence of localization, formation of cup-cone were obtained with ABAQUS.

Effect of meshing

As already mentioned in [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF], mesh design plays an important role in describing localization and cup-cone formation. [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] used square elements (4 nodes) which were divided in 4 linear triangles (3 nodes). In order to numerically study the effect 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 [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]. However the loading ends were meshed ; this is sufficient 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 figure 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 [START_REF] Bonora | Identification and measurement of ductile damage parameters[END_REF]. 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) ≤ f ≤ 0.1(black). Figure 7 shows damage maps at Gauss points for meshes having a number of elements in the minimum cross section N h equal to 20, 40 or 80 and for initial element aspect ratios r h equal to 3:1, 6:1 and 12:1. Elements have quadratic shape functions with 8 nodes and reduced integration (4 Gauss points). The aspect ratio equal to 6:1 leads to approximately square elements at the onset of fracture. Too flat elements (r h = 12 : 1) always lead to flat fracture. A minimum number of elements in the cross section is required to trigger the cup-cone : it was found that N h ≥ 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 r h = 6 : 1 and N h = 40 will be used.

Figure 8 compares the damage maps obtained with different type of elements. First linear (cax4: 4 nodes, 4 Gauss points) and quadratic (cax8r: 8 nodes, 4 Gauss points) square elements are compared. The selective integration method [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and non linear media[END_REF] is used in the case of the cax4 element. It is shown that the cax4 leads to flat fracture. With this element, using 80 elements leads also to cup-cone formation. The number of degree of freedom (DOF) is then approximatively the same as in the case of the cax8r elements with N h = 40. This is also consistent with the flat fracture path obtained in the case cax8r-N h = 20 for which the number of DOF ( 4302) is about the same as for cax4-N h = 40 (4018).

The previous square elements were divided into 4 triangular elements as done in [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]: cax8r→cax6r (6 nodes, 3 Gauss points) and cax4→cax3 (3 nodes, 1 Gauss point). In both cases, cup-cone is observed. A difference 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 (figure 8). As a result, cup-cone angles obtained with quadratic and triangular elements differ. 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: Effect of symmetry Figure 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 flat crack is simulated with the first. This indicates that cup-cone formation is favored when using symmetry.

Cup-cone formation

The evolution of damage f and of the localization indicator min n n.L c . n is detailed in figure 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 figure 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/S 0 -diameter reduction ∆d/d 0 curve (F : force, S 0 initial cross section, ∆d: diameter variation, d 0 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 flat 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 deflection starts with a relatively small angle. Deflection seems to be possible when the localization wings extend over 2 and more elements. At step (5) the wings have grown sufficiently so that the cup-cone develops (6)(7). Some secondary highly damaged regions, which are no longer deforming, are left behind the main crack (figure 7). They could possibly correspond to the bands observed on the surface of the specimens in figure 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.

Figure 11 shows the localization angles θ expressed in the fixed frame as a function of the element position along the radial direction. The values correspond to the angles computed when min n det n.L c . n becomes negative for the first time. The graph is divided in different segments 1••6. Rupture remains flat 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 θ = ±60 • . From 5 to 10 elements (seg. 3), localization wings start to grow, and slightly more dispersed angles are obtained. Note that the macroscopic fracture surface remains flat 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 first moves towards the upper part of the specimen (seg. 4) corresponding to a negative localization angle θ ≈ -60 •2 . Segment 5 corresponds to the formation of two localization bands (step 5 in figure 10). After this point (seg. 6) the crack takes a new direction (going towards the bottom of the specimen, i.e. positive values of θ) with increasing values of θ. Dots on figure 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 (θ > 90 • ) or abruptly change its path with θ ≈ 10 • • • 20 • . Numerically, the crack tends to become 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, q 2 = 1.15, N h = 80 in figure 12).

Remark 1: Effect of material heterogeneities Following the method presented in [START_REF] Decamp | Size and geometry effects on ductile rupture of notched bars in a C-Mn steel: Experiments and modelling[END_REF][START_REF] Devillers-Guerville | Notch fracture toughness of a cast duplex stainless steel: modelling of experimental scatter and size effect[END_REF], the mesh was divided in square regions of size (d 0 /80)×(d 0 /80). In each square an initial value of the porosity was randomly selected according to a uniform distribution with 1.1 10 -4 ≤ f ≤ 1.9 10 -4 . Starting from this level of heterogeneity it becomes possible to generate zigzagging crack paths which are similar to what was experimentally observed (figure 1-(c)). The deformation corresponding to the first sharp drop of the load is reduced. On the other hand, the strain to failure (defined 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 effects are limited in the present case and represent about 2% of the total strain to failure.

Remark 2: Mesh size effect The mesh size effect 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 different sizes are computed. As the size of the localization 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 3: Results on notched bars It is experimentally observed that fracture in notched bars (figure 1-(g,h)) is essentially flat. 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 figure 9 showing that a flat fracture is modeled using N h = 40. Once again, the local localization angle θ is smaller than π/2 but the localization "wings" remain too small to trigger deviation. For N h = 80 a cup-cone is formed. Similar results are obtained applying rate independent plasticity.

Role of constitutive equations

Gurson model

As already mentioned, the Gurson model together with the most commonly used parameters for q 1 and q 2 leads to an overestimation of the ductility. Three different solutions were envisaged to fit the ductility: (1) use of the f function, (2) adding strain controlled nucleation, (3) adjustment of q 2 . In the first case, the adjusted parameters are: f c = 0.005, δ = 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 first case was taken. Above this value, A n is constant and equal to 0.2. In the third case, q 2 = 1.15 gives the best fit. Numerical results using the Gurson model are shown in figure 12.

Using f , flat fracture is always obtained for both cax8r or cax6r elements with 40 or 80 elements in the cross section, with plastic or viscoplastic behavior. f 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 flat 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 difference between the calculations carried out with the Rousselier model and the Gurson model using q 2 = 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 flat part on the cup-cone extends over more that 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 N h = 40 and cax8r elements a slightly zigzagging crack path is obtained. This zigzagging crack path is more pronounced with an elastoplastic behavior which is more susceptible to localization. Using a finer geometrical description (N h = 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: Effect of f R on the crack path In the case of the Rousselier model, the value of the failure porosity f R 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, f R acts as f so that the development of the localization zones is inhibited: this, again, leads to flat fracture. Similarly, using the localization criterion as a rupture criterion at each Gauss point will probably lead to the same results.

Viscosity

Figure 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.L c . n are negative correspond to zones where the normalized growth factor ω/( ṗ/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 flat fracture surfaces as shown on figure 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 σ /σ 1 thus increasing the possibility of normal fracture. The viscosity effect was therefore checked using the Gurson model with q 2 = 1.15 and cax6r elements. As necking is more pronounced in that case the mesh had to be slightly redesigned to obtain approximatively square elements at the onset of rupture. Results are also shown in figure 13 indicating that cup-cone formation is also suppressed in that case.

Plane strain specimens

Results obtained with plane strain specimens are consistent with those obtained on round bars. Some typical results are shown in figure 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 N h = 40 except when using the f 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.

For plastic behavior, the localization indicator is based on Rice's analysis of bifurcation. For time dependent 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 first 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 effects 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 effect which lead to the diminution of the LZ. Increasing viscosity will lead to flat fracture.

Strain rate: For a given viscosity, increasing the strain rate also limits localization thus promoting flat fracture. This result is in agreement with creep experiments carried out by [START_REF] Kobayashi | Nucleation and growth of creep voids in circumferentially notched specimens[END_REF] on notched aluminum bars. Under a high mean stress (30 MPa, high strain rate) specimens exhibit a flat fracture surface whereas under stress (16 MPa, low strain rate) cup-cone fracture is formed.

Use of f : The f function employed in the GTN model can cause the crack path to remain flat. As its derivative with respect to porosity is discontinuous, it may suddenly induce localization as f reaches f c 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 f c is large enough so that localization occurs for lower values of the porosity. This suggests also that it would be more appropriate to adjust q 1 and q 2 with a high value of f c used for numerical purposes only. This solution was in fact adopted by [START_REF] Gullerud | Simulation of ductile crack growth using computational cells: numerical aspects[END_REF] ; q 1 and q 2 can be adjusted from unit cell calculations as in [START_REF] Faleskog | Cell model for nonlinear fracture analysis -I[END_REF].

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 f can produce flat 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 f 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 hydrostatic 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 [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF][START_REF] Brocks | Verification of the transferability of micromechanical parameters by cell model calculations with visco-plastic materials[END_REF] 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 σ eq -σ kk /3 plane (figure 2) [START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF]. 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 f 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 [START_REF] Fleck | Yielding of metal powder bonded by isolated contacts[END_REF] (high porosity) which presents a vertex as done in [START_REF] Redanz | Analysis of shear band instabilities in sintered metals[END_REF]). An alternative solution could consist of using a combination of the Gurson model and of the straight line yield surface derived by [START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF]. In [START_REF] Redanz | Analysis of shear band instabilities in sintered metals[END_REF] and [START_REF] Fleck | Yielding of metal powder bonded by isolated contacts[END_REF] the transition porosity between both model is fixed; it could also be derived from micromechanical models of coalescence as in [START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF][START_REF] Zhang | A class of generalized midpoint algorithms for the gurson-tvergaard material model[END_REF][START_REF] Benzerga | Coalescence-controlled anisotropic ductile fracture[END_REF]. Effect of mesh refinement and element initial aspect ratio on the formation of the cup-cone (cax8r elements). 36

Effect 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 N h = 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. ε L deformation at which min det A c < 0 for the homogeneous analysis.

Table 3 Comparison of FE calculations and localization analysis (∆ε = 0.5%, ε = 2. 10 -3 ). The vector F introduced in section 3.2, is then written as: F(Z, ε) = (F e , F f , F p ) Using the following notations,

N = ∂υ ∂σ = ∂ 2 σ ∂σ 2 = 1 χ 2 ∂ 2 Φ ∂σ∂σ ⊗ ∂Φ ∂σ - 1 χ ∂ 2 Φ ∂σ 2 (B.4) υ ,f = ∂υ ∂f = ∂ 2 σ ∂σ∂f = 1 χ 2 ∂ 2 Φ ∂σ ∂f ∂Φ ∂σ - 1 χ ∂ 2 Φ ∂σ∂f (B.

   gives an example of the evolution of the growth factor ω = ω/( ṗ/p) as a function of deformation for different 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 ψ and θ is given on figure5-(b)showing a progressive change from pure tangential to pure opening fracture as deformation increases. Curves obtained for different values of K coincide as soon as the growth factor is large enough (≈ 1000).

  the localization indicator as a function of the deformation using the elastoplastic tangent matrix (dashed line) and for consistent tangent matrix (full line) with ∆ε = 0.1, 0.2, 0.5, 1.0% (plane strain, Gurson model). 33 5 (a) Localization indicator ω/( ṗ/p) as a function of the deformation for different values of the creep parameter K (b) Localization angles as a function of the deformation. ( ε = 2. 10 -3 s -1 , Gurson model). 34 6 Minimum value of the localization indicator as a function of deformation using the consistent tangent matrix for an elastoviscoplastic material (full lines) for different values of the creep parameter K (∆ε = 0.5% and ε = 2. 10 -3 , Gurson model). The elastoplastic case is indicated by dashed lines. 35 7

Fig. 1 .

 1 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.

Fig. 3 .

 3 Fig. 3. Rupture mode map.

Fig. 4 .

 4 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 ∆ε = 0.1, 0.2, 0.5, 1.0% (plane strain, Gurson model).

Fig. 5 .

 5 Fig. 5. (a) Localization indicator ω/( ṗ/p) as a function of the deformation for different values of the creep parameter K (b) Localization angles as a function of the deformation. ( ε = 2. 10 -3 s -1 , Gurson model).

Fig. 6 .

 6 Fig. 6. Minimum value of the localization indicator as a function of deformation using the consistent tangent matrix for an elastoviscoplastic material (full lines) for different values of the creep parameter K (∆ε = 0.5% and ε = 2. 10 -3 , Gurson model). The elastoplastic case is indicated by dashed lines.

Fig. 8 .

 8 Fig.8. Effect 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 N h = 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.

Fig. 9 .

 9 Fig. 9. Fracture path in notched round specimens.

Fig. 10 .

 10 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 indicator).

Fig. 11 .

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

Fig. 12 .

 12 Fig. 12. Simulated fracture of round bars using the Gurson model with different sets of parameters.

Fig. 13 .

 13 Fig. 13. Effect of viscosity on the macroscopic behavior and on the cup-cone formation.

Fig. 14 .

 14 Fig. 14. Fracture path in plane strain specimens for different models.

  needed to compute L p are given by:∂F e ∂ε e = -(1 -f )FN : C -(1 -f )F (υ : C) ⊗ υ ∂F e ∂f = F (υ -(1 -f )υ ,f ) -(1 -f )F σ ,f υ ∂F e ∂p = (1 -f )F Hυ ∂F f ∂ε e = (1 -f ) 2 FN : C : 1 + ((1 -f ) 2 υ : 1 + A n )F υ : C ∂F f ∂f = F(1 -f ) ((1 -f )υ ,f : 1 -2υ : 1) + ((1 -f ) 2 υ : 1 + A n )F σ ,f ∂F f ∂p = A n F -((1 -f ) 2 υ : 1 + A n )f ) 2 ∆pN : C : 1 ∂F f ∂f = ∆p(1 -f ) ((1 -f )υ ,f : 1 -2υ : 1) ∂F f ∂p = (1 -f ) 2 υ : 1 + A n + A n ∆p

  Angles θ and ψ (figure3), porosity and deformation (ε) at localization predicted using the consistent tangent matrix for elastoplastic and elastoviscoplastic behaviors (Gurson model, ε = 2. 10 -3 s -1 , ∆ε = 0.5%). The evolution of the corresponding localization indicator is given on figure6. D deformation at which min det A c < 0 at, at least one Gauss point in the FE calculation.ε B deformation at which a localization band is formed in the FE calculation.ωD normalized growth factor corresponding to ε D .ωB normalized growth factor corresponding to ε B .

		The arrow indicates a change in the crack direction in the
		undeformed mesh. smooth	notched	notched		plane	37
	9	Fracture path in notched round specimens. r = 4 mm r = 2 mm strain	37
			0.79 FE calculation Localization analysis 0.60 0.46 0.40
	10 Load vs. diameter reduction curve and cup-cone formation for Table 1 K ε D ε B ωD ωB ε L a round bar. Contour maps indicate damage (white: f = 0, black: f > 0.1) and the localization indicator (white: elastic Area reductions at fracture. pl. 1.28 1.66 --1.29
		unloading, gray: plastic loading, black plastic loading and negative indicator). 1 1.34 1.67 6882 61150 1.32	38
			5 1.43	1.67	2350 12800 1.44
	11 Predicted angle of localization θ as a function of the element position along the radial direction and comparison with the 10 1.53 1.69 2300 6860 1.55
		actual numerical crack orientation (dots). 25 1.67 1.71 2252 2990	1.80	39
	12 Simulated fracture of round bars using the Gurson model with 50 1.76 1.77 1850 1900 2.09
		different sets of parameters. 100 1.82 1.83	1140 1250	2.54	39
	13 Effect of viscosity on the macroscopic behavior and on the ε
		cup-cone formation.					40
	14 Fracture path in plane strain specimens for different models. K pl. 1 5 10 25 50	40
		θ ( • )	45.65 45.73 46.06 46.49 46.95 40.75
		ψ ( • ) 88.67 88.53 87.88 87.03 84.05 78.44
		f (%) 1.31	1.45	2.14	3.01	6.32	13.1
		ε	1.29	1.32	1.44	1.55	1.80	2.09
	Table 2

Preprint submitted to Int. J. Solids and Structures 2001

As slightly different parameters were used for the calculation of figure11and figure 10, crack paths do not exactly coincide. In figure10, the first crack deflection corresponds to the θ = +60 • branch. This clearly illustrates the non-uniqueness of the solution. Otherwise, results of both calculations are consistent.

Acknowledgments

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 acknowlegded for financial support and hospitality.

List of Tables 1 Area reductions at fracture. 28 2 Angles θ and ψ (figure 3), porosity and deformation (ε) at localization predicted using the consistent tangent matrix for elastoplastic and elastoviscoplastic behaviors (Gurson model, ε = 2. 10 -3 s -1 , ∆ε = 0.5%). The evolution of the corresponding localization indicator is given on figure 6. 28 3

Comparison of FE calculations and localization analysis (∆ε = 0.5%, ε = 2. 10 -3 ). 29 A Calculation of the elastoplastic tangent matrix: L t .

List of Figures

The rate equations governing the evolution of the state variables are:

The consistency condition φ = 0 yields:

H = dR/dp is the hardening modulus. The derivative σ ,f = ∂σ /∂f can be computed (in the case of the Gurson model) as:

with χ = ∂Φ/∂σ . In the following, it will be assumed that h > 0. h = 0 can be encountered for low porosities and high stress triaxiality ratios. In this case, the constitutive equations cannot be integrated for a prescribed strain rate tensor (snap back effect). Finally the elastoplastic tangent matrix is given by:

which is a special case of the more general form studied in [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF] for instance.

B Perturbation analysis: calculation of H p

In the case of viscoplasticity, the rate equations governing the evolution of the elastic strain and the porosity are still given by eq. A.1 and A.2. The evolution of p is given by the viscoplastic flow law:

C Calculation of the consistent tangent matrix: L c

The fully implicit integration of the model is obtained using a time discretization of the rate equations A.1, A.2 and A.3 or B.

3

All quantities should be considered at the end of the time increment. The equation corresponding to ∆p is written as:

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: ∆p = F(φ)∆t ≡ F p (C.4)

The matrix I 0 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 F is given by (F e , F f , F p ).

The derivative ∂F /∂ε and ∂F /∂Z can then be calculated noting that ∂∆x/∂x = 1, ∀x. One gets:

∂F ∂ε = (I, 0, 0) (C.5)

The derivatives of F are: