Numerical simulation of ductile fracture with the Rousselier constitutive law

Available online at www.sciencedirect.com

Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982 www.elsevier.com/locate/cma

Numerical simulation of ductile fracture with the Rousselier constitutive law E. Lorentz a,c,*, J. Besson b, V. Cano c a

Laboratoire de Me´canique des Structures Industrielles Durables, UMR CNRS/EDF 2832, 1 Av. Ge´ne´ral de Gaulle, 92141 Clamart cedex, France b Ecole des Mines de Paris, Centre des Mate´riaux Pierre et Marie Curie, UMR CNRS 7633, BP 87, 91003 Evry cedex, France c EDF R&D, De´partement Analyses Me´caniques et Acoustique, 1 Av. Ge´ne´ral de Gaulle, 92141 Clamart cedex, France Received 14 January 2007; received in revised form 22 November 2007; accepted 21 December 2007 Available online 3 January 2008

Abstract This contribution aims at developing robust and reliable tools dedicated to predictive numerical simulation of ductile fracture. Description of damage mechanisms relies here on the Rousselier constitutive law. To achieve robustness and objectivity, several points need then to be addressed: reduction of the set of internal variables so as to recast the model into the framework of generalised standard materials; proper expression of the normal plastic flow rule at the vertex of the yield surface; finite strain formulation based on a multiplicative split of the deformation gradient; fully implicit integration of the constitutive equations leading to a single scalar equation that admits a unique root; use of mixed finite elements to avoid volumetric locking; strain gradient-based nonlocal model to control strain localisation and remedy spurious mesh dependency. The capabilities of the resulting numerical formulation are demonstrated through simulating the cup–cone fracture of a notched tensile specimen. Ó 2008 Elsevier B.V. All rights reserved. Keywords: Constitutive law; Plasticity; Ductile fracture; Rousselier; Mixed finite element; Nonlocal model

1. Introduction Predictive numerical simulations of ductile fracture may be of great interest in industrial situations for which fullscale experimental approaches are either too costly or even impracticable. Consider for instance ductile tearing of gas pipelines over several hundred meters or interaction between ductile and brittle fracture to assess the lifespan of aging nuclear vessels in order to operate in full safety. For such applications, simulations should predict crack paths, stability, stress states, etc. On one hand, crack resistance approaches based on the J integral are not applicable because they are strongly geometry-dependent and restricted to limited extensions of pre-existing cracks [1]. * Corresponding author. Address: Laboratoire de Me´canique des Structures Industrielles Durables, UMR CNRS/EDF 2832, 1 Av. Ge´ne´ral de Gaulle, 92141 Clamart cedex, France. E-mail address: [email protected] (E. Lorentz).

0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2007.12.015

On the other hand, cohesive zone models do not seem mature to predict crack paths because they exhibit strong mesh dependency and over-estimate cracked areas [2]. Therefore, continuum damage mechanics (CDM) is generally acknowledged as the best approach to perform such analyses. CDM relies on specific constitutive models to describe damage processes; unfortunately, all of them exhibit a common property, strain-softening, that leads to ill-posed problems, characterised by pathological strain localisation [3,4]. The practical remedy consisting in setting the mesh size does not solve ill-posedness. Namely, spurious dependencies to mesh orientation and finite element types still remain [5,6] and preclude predictive simulations. On a physical ground, ill-posedness is related to the fact that the assumption of scale separation used to derive the constitutive relations is no more valid when strain localisation occurs: a constitutive coupling between neighbouring material points should emerge, resulting in a nonlocal constitutive law. Several

1966

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

phenomenological approaches are available, see [6] for a review; they all require fine meshes with finite elements several times smaller than the width of the localisation bands, which result in large numbers of degrees of freedom. Therefore, key objectives when designing ductile models should be robustness and objectivity to ensure at best convergence of the computations and objectivity to provide reasonably predictive results in various situations ranging from specimen simulations to adjust material parameters up to fullscale industrial structures. Achieving robustness and objectivity is the goal of the article; its fulfilment will be demonstrated through the simulation of the cup–cone fracture of a notched tensile specimen, so that pioneering work results [7] are retrieved. The starting point is the choice of the constitutive law. To achieve robustness, a simple model is selected: that namely implies damage isotropy (nucleation, growth and coalescence of a single population of spherical cavities) and plastic isotropy (hence a reduced dependence on the choice of finite strain formulation). Though popular, the Gurson one [7,8] is disregarded to the benefit of the Rousselier one [9–11]. In spite of acknowledged drawbacks, namely its incorrect behaviour in compression or the lack of straightforward extension to viscoplasticity [5], the latter enables transition from flat to slant fracture [5,12] without requiring specific nucleation nor coalescence laws, thanks to the vertex of its yield surface which preserves shear flow even for high triaxiality. Moreover, it is shown in Section 2 how it may be recast into the framework of generalised standard materials (GSM) [13,14] and benefit from a finite strain extension with multiplicative split of the strain into plastic and elastic parts as expressed in [15]. Note that the GSM framework also provides the exact flow rule at the vertex point. The numerical integration of the constitutive equations is described in Section 3. It is shown that it can reduce to the solution of a single scalar equation which always admits a single root. The accuracy of the integration scheme is also analysed, regarding some approximations that have been made. Finally, Section 4 deals with structural computations. The question of appropriate finite elements to avoid volumetric locking phenomena is given some answer. Rather than selecting low order finite elements [16,17], the authors turned to mixed high order triangles and tetrahedrons, following the formulation of [18,19]. For convenience, the choice of finite elements also suggests the choice of a nonlocal model based on the second gradient of the deformation, which has already been put into practice for brittle as well as ductile materials [20–24]. As announced, the applicability of the whole formulation is demonstrated by modelling a notched tensile specimen. 2. The Rousselier model revisited 2.1. Finite strain notations First, some notations are introduced. X and x are the positions of a material point of the body, respectively, in

the initial configuration X0 and in the current one X. F is the deformation gradient and D the Eulerian strain rate: F¼

ox ¼ rx; oX

1 D ¼ symðLÞ ¼ ðL þ LT Þ; 2 ð1Þ

_ 1 ; L ¼ FF

where the superscript T denotes transpose of a second order tensor and a dot denotes differentiation with respect to time. dU is the displacement variation and dL the Eulerian deformation variation: dL ¼

odU ¼ dFF1 : ox

ð2Þ

A (virtual) intermediate relaxed configuration XP is introduced by unloading every material point [15]: Fp and Fe are, respectively, the (plastic) deformation gradient from X0 to XP and the (elastic) deformation gradient from XP to X, with the corresponding multiplicative split: F ¼ Fe Fp :

ð3Þ

The elastic, plastic and total volumetric strains are, respectively, Je, Jp and J: J e ¼ det Fe ;

J p ¼ det Fp ;

J ¼ det F;

J ¼ J eJ p:

ð4Þ

The Eulerian plastic strain rate is defined as 1

Dp ¼ symðF_ p Fp Þ

ð5Þ

while symmetrical elastic and plastic strain measures are introduced as: Be ¼ Fe FeT ;

1 e ¼ ðId  Be Þ; 2

Be ¼ FGp FT :

1

Gp ¼ ðFpT Fp Þ ; ð6Þ

Note that e is the elastic Almansi strain tensor. Finally, r is the Cauchy stress and s the Kirchhoff stress, related by s ¼ J r:

ð7Þ

2.2. The Rousselier model Rousselier introduced in [9,10] an elasto-plastic damage constitutive law with isotropic hardening to describe ductile fracture of metals. Evolution of damage is set so as to ensure the following relation between volumetric plastic growth and porosity: tr Dp ¼

f_ : 1f

ð8Þ

This relation provides a physical ground to the model since it reflects the growth of a cavity inside an incompressible plastic matrix [8]. More recently, Rousselier gave a new presentation of the constitutive law, equivalent to the former one and in which porosity itself is the damage variable [11]. The material state is described by the strain and a set of internal variables: the plastic strain, the hardening variable p (cumulated plastic strain) and the porosity f. The free Helmholtz energy is composed of three terms:

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

the elastic strain energy, the stored energy due to hardening and the stored energy due to damage. It reads URS ðF; Fp ; p; f Þ ¼ Uel ðeÞ þ Ust ðpÞ þ Udam ðf Þ:

ð9Þ

The state equations are obtained by derivation of the free Helmholtz energy with respect to the state variables. Thus, the driving forces A and Y, respectively, associated to p and f read Z p oURS ¼ RðpÞ; ð10Þ Ust ðpÞ ¼ RðsÞds ! A ¼  op 0 Udam ðf Þ ¼ r1 ½f ln f þ ð1  f Þ lnð1  f Þ ! Y   oURS f ¼ r1 ln ¼ ; 1f of

ð11Þ

where R is the hardening function and r1 a material characteristic. In the context of finite strain, the stress–strain relation is given a rate expression through a hypoelastic law: r

D ¼ Dp þ E1 : r;

ð12Þ

where E is Hooke’s tensor and the superscript triangle denotes some objective stress rate. Such an additive decomposition is known to be admissible as long the elastic strain remains small compared to the total strain. The evolution of the internal variables is governed through a yield threshold G which defines the elasticity domain: Gðs; A; Y Þ ¼ seq þ r1 D

expðY =r1 Þ expðsH =r1 Þ þ A  ry ; 1 þ expðY =r1 Þ ð13Þ

where ry denotes the yield stress while D is an additional material parameter. It is interesting to note that the stress measure involved in the yield threshold is the Kirchhoff one. It has been split in its deviatoric part sD (von Mises norm seq) and its hydrostatic part sH, responsible for the volumetric plastic strain rate: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 3 D D D H H s :s : s ¼ s þ s Id; s ¼ tr s ¼ s : Id; seq ¼ 3 3 2 ð14Þ The flow rules are obtained by associated normality with respect to the yield surface: Dp ¼ k

oG ; os

p_ ¼ k

oG ; oA

oG f_ ¼ k : oY

ð15Þ

The plastic multiplier k is given by the consistency equation: Gðs; A; Y Þ 6 0;

k P 0;

kGðs; A; Y Þ ¼ 0:

ð16Þ

It can be observed from (15) that a straightforward calculation yields again the relation (8). Eqs. (10)–(16) describe the cavity growth regime. Nucleation is simply modelled by assuming that the initial porosity results from second phase particles or nonmetallic inclusions of volume fraction f0:

f ð0Þ ¼ f0 :

1967

ð17Þ

In particular, no continuous nucleation is assumed in the study. Coalescence of cavities is roughly modelled through strain localisation: no specific critical porosity is introduced on the contrary of what is proposed for the Gurson law in [7]. In spite of these assumptions, the Rousselier model is able to describe complex structural damage [2,5]. 2.3. Reformulation of the Rousselier model as an elastoplastic constitutive law Thanks to associated normal flow rules, the Rousselier model should benefit from the attractive properties of generalised standard materials, ranging from Thermodynamics to numerical integration [14,25,26]. Unfortunately, a straightforward calculation shows that the yield function G is not convex with respect to (sH, Y). Hence, it is neither convex with respect to its whole set of arguments (s, A, Y), a point seemingly missed in [11]. This hinders the abovementioned properties. Therefore, we aim at recasting the Rousselier model into the framework of GSM. The main benefit is the existence and uniqueness of the solution to the time-integration of the constitutive law, an essential property regarding robustness. The corresponding constitutive equations are summarised in Box 1, so that readers uninterested by the different steps can skip what follows and go directly to Section 3. Box 1: Constitutive equations of the revisited Rousselier model. Elastic and plastic strain definition T F ¼ Fe Fp ; Be ¼ Fe Fe ; e ¼ 12 ðId  Be Þ; T Gp ¼ ðFp Fp Þ1 ; Be ¼ FGp FT : Stress–strain relation s = (E:e)Be. State equation A = R(p). Yield criterion F(s, A;f) = seq + r1D f exp(sH/r1) + A  ry. Flow rule p_ ¼ k:   _ p FT ¼ k  2 Df expðsH =r1 ÞId  3N Be ; FG 3 ( D N ¼ sseq if seq > 0; N 2 fd; tr d ¼ 0 and d eq 6 1g if seq ¼ 0: Consistency condition F(s, A;f) 6 0; k P 0; kF(s, A;f) = 0. Porosity equation   0 f ¼ max f0 ; 1  1f : J

1968

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

First, we propose a modification of the porosity growth law (8). In the latter relation, the elastic strain has been neglected. This is an approximation and (8) may as well be rewritten as f_ : tr D ¼ 1f

ð18Þ

The advantage of expression (18) is the existence of an integrated expression: tr D ¼

1  f0 f_ () J ¼ : 1f 1f

ð19Þ

In this equation, no lower bound is taken into account for the current porosity f. However, it cannot go lower than the initial porosity f0 since the cavities are actually filled with particles. Therefore, Eq. (19) has to be complemented by a closure condition, leading to   1  f0 f ¼ max f0 ; 1  : ð20Þ J In that way, the porosity has become an explicit function of the volumetric strain and the initial porosity. The set of internal variables reduces to the plastic strain and the hardening variable, as already pointed out in [5]. The free Helmholtz energy is reduced to the elastic strain energy plus the stored energy due to hardening: UðF; Fp ; pÞ ¼ Uel ðeÞ þ Ust ðpÞ:

ð21Þ

Moreover, by substituting (11) in (13), the usual Rousselier criterion parameterised by the porosity (hence by the volumetric strain) is retrieved, see Fig. 1: F ðs; A; f Þ ¼ seq þ r1 Df expðsH =r1 Þ þ A  ry

ð22Þ

The evolution of the internal variables is still governed by the flow rule (15) and the consistency condition (16), which now read

oF oF ; p_ ¼ k ; os oA F ðs; A; f Þ 6 0; k P 0;

Dp ¼ k

ð23Þ kF ðs; A; f Þ ¼ 0:

As F is a convex function with respect to (s, A), the Rousselier model is recast into the framework of GSM. Compared to a classic von Mises elasto-plastic law with isotropic hardening, its only specificities are the dependence of the yield threshold on the hydrostatic stress and on the volumetric strain. Actually, it should be noted that the dependence of the yield threshold on the strain (through the porosity) leads to a slightly extended definition of generalised standard materials, see [27]. However, the important fact is that the properties related to the evolution of the internal variables for a given strain field are preserved (normality defined by means of convex analysis, energetic interpretation of the set of time-discretised equations, existence and uniqueness of the incremental plastic evolution). 2.4. Interpretation of the flow rule at the vertex point The dependence of the yield function F on the volumetric strain does not raise any special difficulty. But its dependence on the hydrostatic stress does. Indeed, F is not differentiable when seq = 0 (hydrostatic axis). On the contrary of von Mises and Gurson criteria, there exists a vertex (sv, Av) on the yield surface where the flow rule cannot be expressed as the differential expression (23), see again Fig. 1:  y  r þ RðpÞ ð25Þ sv ¼ r1 ln Id; Av ¼ RðpÞ: Df r1 The evolution equations (23) and (24) have to be extended. Thanks to the GSM framework, the flow rule actually derives from the maximum dissipation principle which does not rely on differentiability but on convexity only (sub-differentiability). More precisely, let us denote Kf the (convex) elasticity domain for a given porosity f: b F ð^s; A; b f Þ 6 0g: K f ¼ fð^s; AÞ;

Y

ð24Þ

ð26Þ

For any point (s, A) inside Kf, the maximum dissipation principle states that the corresponding flow direction, say (d, k), maximises the dissipation, that is it belongs to the following convex cone:

Z X

b 2 K f s : d þ Ak P ^s : d þ Akg: b Cðs; A; f Þ ¼ fðd; kÞ; 8ð^s; AÞ ð27Þ Extension of the evolution equations (23) and (24) then simply reads _ 2 Cðs; A; f Þ: ðDp ; pÞ

Fig. 1. Rousselier yield surface in the principal stress space (p = 1; J = 1.1).

ð28Þ

After some calculations detailed in Appendix A, C(s, A; f) can be expressed as the set of flow directions (d, k) that satisfy:

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

3 d ¼ kN with 2

(

D

tr d ¼ kDf exp

D

N ¼ sseq

if seq > 0;

N 2 fd; tr d ¼ 0 and d eq 6 1g if seq ¼ 0; ð29Þ

 H s ; r1

F ðs; A; f Þ 6 0;

k P 0;

kF ðs; A; f Þ ¼ 0:

ð31Þ

2.5. Finite strain formulation The elastic part of the constitutive relation (12) is based on a hypoelastic law involving some objective stress rate. This has some drawbacks: spurious elastic dissipation, arbitrary stress rate, increased complexity to achieve incremental objectivity. In the case of yield thresholds that depend on the Kirchhoff stress, Simo and Miehe [28] have brought an elegant answer to these issues. They introduced a hyperelastic law instead of (12) and deduced the flow rule from the maximal dissipation principle. This finite strain formulation seems well-suited to the Rousselier model as revisited in Sections 2.3 and 2.4; it is applied hereafter. In agreement with (21), the free Helmholtz energy is expressed as a function of the elastic strain e defined in (6) and of the hardening variable p: Z p 1 UðF; Fp ; pÞ ¼ Uel ðeÞ þ Ust ðpÞ ¼ e : E : e þ RðsÞds: 2 0 ð32Þ A quadratic expression of the elastic strain energy is reasonable as long as the elastic strain remains small. Besides, it is assumed that Hooke’s tensor is isotropic, with K and l the compressibility and shear moduli: E : e ¼ Kðtr eÞId þ 2leD :

ð33Þ

Under the assumption of no thermo-mechanical coupling, the mechanical dissipation reads

þ

  oUel e sþ B :D oe

st 1 oUel _ p FT Þ  oU p: _ : ðFG 2 oe op

ð34Þ

As the dissipation should be zero for elastic evolution, the stress–strain relation can be deduced s¼

ð36Þ

d

ð30Þ

The first alternative in (29) corresponds to regular points of the yield surface, for which the differential flow rule (23) is retrieved. The second alternative in (29) corresponds to the vertex of the yield surface, for which a whole set of flow directions are allowed.

D ¼ s : D  U_ ¼

And the dissipation reduces to   1 _ p T e1 _ þ Ap: D ¼ s :  FG F B 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

1969

oUel e B ¼ ðE : eÞBe : oe

ð35Þ

_ p is naturally It should be noticed that the time derivative G p objective since G is a material tensor, i.e. expressed in the initial configuration. Thanks to the assumption of maximal dissipation with respect to the convex yield criterion F, the results of Section 2.4 can be applied. From (28), we deduce that the evolution equations read   1 _ p T e1  FG F B ; p_ 2 Cðs; A; f Þ: ð37Þ 2 Taking into account the integrated expression of the porosity (Section 2.3), the generalised expression of the flow rule (Section 2.4) and the present finite strain formulation leads to a new expression of the Rousselier model. The whole set of constitutive equations is gathered in Box 1. 3. Numerical integration 3.1. Stress approximation We focus our attention on solving numerically the set of constitutive equations in Box 1. A straightforward discretisation with respect to time using an implicit Euler scheme is applied. As pinpointed by Simo and Miehe [28], it leads automatically to an incrementally objective formulation _ p is expressed on a material tensince the only tensor rate G sor. However, the resulting system of equations is highly nonlinear so that some approximations are suggested in [28] to enable a practical efficient solution. Here, another approximation is done, following [29]. The discrete equations are still incrementally objective while their structure is identical to that which would be obtained under the assumption of infinitesimal strain. In particular, the solution can be interpreted as the minimum of a convex energy, thanks to the normal and associated character of the evolution equations. Therefore, the system of discrete equations admits a unique solution, a highly desirable property when applying some numerical algorithm to solve the system. It should be noticed that these properties are obtained without the introduction of logarithmic strains or exponential-map as in [30,31] for instance. This is thought to be an advantage regarding the numerical solution since no diagonalisation of the strain tensor is required, a cumbersome operation namely to derive the tangent operator. However, the price to pay is a correction a posteriori of the volumetric part of the flow rule in order to avoid instabilities which may appear with relatively large strain increments, as will be shown later. This is in contrast with finite strain integration schemes based on an exponential-map for which the volume conservation is implicitly taken care of. The above-mentioned stress approximation introduced in [29] is based on the replacement of the Kirchhoff stress

1970

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

s by another stress measure s inside the yield function and the evolution equations: s ¼ sBe : ð38Þ e As long as the elastic strain remains small, B  Id so that both stress measures are close to each other. The dissipation (36) then reads   1 _p T _ F þ Ap: D ¼ s :  FG ð39Þ 2 Therefore, the multiplication by Be in the flow rule disappears (the source of the high nonlinearity) when applying the principle of maximal dissipation with respect to s instead of s. That results in a change (simplification) in the stress–strain relation, the flow rule and the consistency condition only: s ¼ E : e; ð40Þ

 H  _ p FT ¼ k  2 Df exp s Id  3N ; FG 3 r1 ( D if seq > 0; N ¼ sseq ð41Þ

N 2 d; tr d ¼ 0 and d eq 6 1 if seq ¼ 0; F ðs; A; f Þ 6 0;

k P 0;

kF ðs; A; f Þ ¼ 0:

ð46Þ ð47Þ

The elastic trial is introduced as the state for which no evolution of the internal variables (plastic strain and hardening variable) takes place during the considered time-step. This state is referenced by the subscript E: pE ¼ pn1 ; GpE ¼ Gpn1 ; 1 eE ¼ ðId  BeE Þ: 2

BeE ¼ Fn Gpn1 FTn ¼ fBen1 f T ; ð48Þ

Consequently, the evolution of the plastic strain can be expressed as Fn ðGpn  Gpn1 ÞFTn ¼ Ben  BeE ¼ 2ðen  eE Þ:

ð49Þ

For the sake of simplicity, the subscript n is henceforth omitted. The set of discrete equations is gathered in Box 2. Thanks to the introduction of the relative deformation gradient f, the plastic strain disappears from the discrete equations: only the hardening variable and the elastic strain need be stored from one time-step to the next.

ð42Þ Box 2: Discrete equations of the Rousselier constitutive law.

3.2. Discrete equations In the context of pressure-dependent yield surface, Zhang [32] showed that a h-method may be more accurate than the backward Euler scheme in the case of large timesteps with 0.75 6 h 6 0.85. However, the important fact in his demonstration is that the consistency equation is always evaluated at the end of the time-step (h = 1). In order to benefit from the above-mentioned properties of existence and uniqueness of the solution, the other equations have to be evaluated with the same value of h as the consistency equation, hence h = 1. Therefore, for the sake of robustness, the Euler backward scheme is adopted: robustness is favoured over accuracy. The set of continuous constitutive equations including the stress approximation (Section 3.1) are now discretised. Consider a time-step ranging from tn1 to tn; qn1 (resp. qn) denotes the value of a quantity q at the beginning (resp. end) of the time-step and Dq = qn  qn1. In addition, the (relative) deformation gradient f from configuration Xn1 to configuration Xn is introduced Fn ¼ fFn1 : ð43Þ The set of discrete equation reads sn ¼ E : en ; An ¼ Rðpn Þ; pn  pn1 ¼ k; Fn ðGpn  Gpn1 ÞFTn

 H  2 sn ¼ k  Dfn exp Id  3Nn 3 r1 ( sD n Nn ¼ sneq with Nn 2 fd; tr d ¼ 0 and d eq 6 1g

F ðsn ; An ; fn Þ 6 0; k P 0; kF ðsn ; An ; fn Þ ¼ 0;   1  f0 fn ¼ max f0 ; 1  : Jn

Porosity equation   0 f ¼ max f0 ; 1  1f ; J ¼ detðFn1 Þ detðfÞ: J Elastic trial pE ¼ pn1 ; eE ¼ fen1 f T þ 12 ðId  ff T Þ: Yield criterion F(s, A; f) = seq + D f exp(sH/r1) + A  ry. State equations A = R(p); s = E:e = (K tr(e) Id + 2leD). Flow rule p  pE = k,  e  eE ¼ k 13 Df expðsH =r1 ÞId þ 32 N ; ( D N ¼ sseq if seq > 0; N 2 fd; tr d ¼ 0 and d eq 6 1g if seq ¼ 0: Consistency condition F(s, A; f) 6 0; k P 0; kF(s, A; f) = 0.

ð44Þ

Stress definitions r ¼ J1 s; s ¼ sðId  2eÞ; s ¼ ðE : eÞ:

ð45Þ

Volumetric flow rule correction e J p ¼ J pn1 exp½ðp  pn1 ÞDf expðsH =r1 Þ; e J e ¼ Jp ; eJ D ~e ¼ e þ ~eH Id with ~eH solution of ð2~eH  1Þ3  13 e2eq ð2~eH  1Þþ ðdet eD  e J 2e Þ ¼ 0:

if sneq > 0; if sneq ¼ 0;

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

3.3. Solution strategy Integration of the constitutive equations consists in computing the evolution of the internal variables and the stress at the end of the time-step by considering that the deformation gradient is known at the beginning and the end of the time-step and the internal variables are known at the beginning of the time-step:  inputs: Fn1; f; pn1; en1;  outputs: rn; pn; en. Therefore, it appears that the discrete equations in Box 2 can be split in three sets. The first one consists of straightforward formulas to compute the porosity and the elastic trial. The second one is specific to the constitutive behaviour and requires a nonlinear solution which is the purpose of this section. The third one is merely a post-treatment to compute the stress tensors, functions of the elastic strain. Thus, the objective of the nonlinear solution is the following computation:  inputs: f; pE; eE;  outputs: p; e. The starting point is the consistency condition. As usual, two alternatives have to be considered: either k = 0 and FE 6 0 (elastic solution), or F = 0 and k P 0 (plastic solution). Checking that the solution is elastic does not raise any difficulty. If it is not elastic, it will be showed that the governing equations reduce to a single scalar equation the unknown of which is x ¼ tr e  tr eE :

ð50Þ

Indeed, the volumetric part of the flow rule reads  H s x ¼ trðe  eE Þ ¼ kDf exp r1     K tr eE Kx ¼ kDf exp  exp  r1 r1 

m

1 K x exp x kðxÞ ¼ GE r1



 Vertex: This point is characterised by: eeq ¼ 0 ) eD ¼ 0

and

N eq ðxÞ ¼

2eEeq : 3kðxÞ

ð54Þ

On the other hand, the solution is licit if and only if x fulfils: 2 N eq ðxÞ 6 1 () kðxÞ P eEeq : ð55Þ 3  Regular point: In that case, a classical computation shows that 3 eD 3 eD  eD () eeq  kðxÞ ¼ eEeq E ¼ kðxÞ 2 2 eeq eD ¼ eeq

and

eD E : eEeq

ð56Þ

Both alternatives can be gathered as   3 eeq ðxÞ ¼ eEeq  kðxÞ ; 2

ð57Þ

where the MacCauley bracket denotes the positive part of a real number. Finally, whatever the alternative, the consistency equation reduces to

   K F ðxÞ ¼ 2leEeq  3lkðxÞ þ r1 GE exp  x r1  RðpE þ kðxÞÞ  ry ¼ 0:

F ð0Þ ¼ F E > 0;

ð58Þ

lim F ðxÞ ¼ ½ lim RðpÞ þ ry  < 0:

x!þ1

p!þ1

ð59Þ

 K tr eE with GE ¼ Df exp  : r1 ð51Þ

ð52Þ

In order to express the yield function as a function of x only, consider the deviatoric part of the flow rule 3 eD  eD E ¼ kðxÞN; (2 D N ¼ eeeq

At that point, a distinction whether the solution is reached on the vertex point of the yield surface or a regular one is necessary.

F(x) is a strictly decreasing function with bounds



And the yield function can be expressed as

 K F ¼ 2leeq þ r1 GE exp  x  RðpE þ kðxÞÞ  ry : r1

1971

if eeq > 0;

N 2 fd; tr d ¼ 0 and d eq 6 1g if eeq ¼ 0: ð53Þ

Therefore, there exists a unique positive solution that can be found by means of Newton’s method [33] or any other solution algorithm dedicated to scalar functions. Then, the other quantities of interest can be computed from x. However, it should be noticed that introducing the unknown x (instead of the hardening variable p for instance) requires that volumetric plastic strain occurs, which precludes a smooth transition from the Rousselier model to von Mises plasticity (f0 = 0). The solution strategy is summarised in Box 3. In addition to the strict integration of the constitutive law, the consistent tangent operator is required to get a quadratic convergence rate for the solution of the equilibrium equations [34]. Its expression is obtained through the derivation of analytical expressions and applications of the chain rule, see Appendix B. Though cumbersome, it does not raise any difficulty.

1972

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

Box 3: Solution strategy for the numerical integration of the Rousselier constitutive law. Objective f ; p E ; eE

! p; e.

compute

 Elastic trial (FE 6 0) h i e ¼ eE ; Ktr eE y F E ¼ 2leEeq þ Df exp  r1  RðpE Þ  r 6 0 ! elastic solution : p ¼ pE : Plastic solution (FE > 0) Definitions:   E GE ¼ Df exp  Ktre ; r  1 kðxÞ ¼ G1E x exp rK1 x ;     F ðxÞ ¼ 2leEeq  3lkðxÞ þ r1 GE exp  rK1 x  RðpE þ kðxÞÞ  ry : Solution of the scalar function: Find x* 2 ]0 + 1[ such as: F(x*) = 0. Straightforward computations of the quantities of interest: D   e eeq ¼ eEeq  32 kðx Þ ; tr e ¼ tr eE þ x ; e ¼ 13 trðeÞId þ eeq eEE ; eq * p = pE + k(x ).

3.4. Volumetric flow rule correction An approximation of the yield function has been introduced in order to simplify the set of discrete constitutive equations. The finite strain geometrical aspects are removed from the core of the constitutive law. Moreover, existence and uniqueness of the incremental solution is achieved. However, the flow rule is normal to the yield surface so that it is also distorted by the approximation. In the case of von Mises plasticity for instance, it was observed that the incompressibility of the plastic strain was not preserved [28,29]. The deviation is of the order of the elastic strain and cumulated from one time-step to the other. Therefore, it was proposed in [28] to correct a posteriori the volumetric part of the flow rule to fulfil exactly the plastic incompressibility. In ductile plasticity, plastic volumetric expansion is a characteristic feature of the model that should also been dealt with care since it measures the degree of damage. Therefore, a correction of the same kind as for von Mises plasticity should be introduced. First, the plastic volume expansion has to be derived from the continuous constitutive equations (Box 1). After some manipulations based on the definitions (3), (5) and (6), the left-hand side d of the flow rule (29)–(31) can also be expressed as the plastic strain rate convected in the current configuration: 1 _ p T e1 1 d ¼  FG F B ¼ Fe Dp Fe : ð60Þ 2 In particular, its volumetric part reads J_ p ð61Þ tr d ¼ tr Dp ¼ p : J

Therefore, without any approximation, the flow rule (30) leads to d _ ðln J p Þ ¼ pDf expðsH =r1 Þ: dt

ð62Þ

The plastic volumetric strain is an increasing function of time, a property which is not necessarily preserved with the approximated discrete equations. Then, discretisation of (62) with the implicit Euler scheme and replacing s by s (since it is the latter that appears in the discrete consistency condition (46)) leads to

 H  e Jp s ¼ exp ðp  pn1 ÞDf exp ¼ expðtr e  tr eE Þ: p J n1 r1 ð63Þ This is a new equation which provides corrected values for the volumetric plastic strain (the tilde stays for ‘‘corrected”), compared to (49). However, in expression (63), the value of the elastic strain in the right-hand side term is computed with the approximated equations. Therefore, a posteriori corrected values are obtained for the plastic volumetric strain e J p and the elastic volumetric strain e p ðe J ¼ J =e J Þ. In the same way as in [28], we propose to correct a posteriori the trace of the elastic strain (the deviatoric component remains unchanged) so that it is consistent with the corrected volumetric strain e Je 2 ~e ¼ eD þ ~eH Id ) e J e ¼  det½X Id þ 2eD  with

X ¼ 2~eH  1:

ð64Þ

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

with the integrated form (20) are compared to those with the differential one (8). Regarding the latter case, the treatment of the porosity variable is explicit: during integration of the constitutive law, the yield function F is evaluated with the value of the porosity at the beginning of the time-step fn1; following (61) and (62), the porosity is then updated at the end of the time-step through the relation:

Table 1 Material parameters Elastic domain E = 198 000 MPa; m = 0.3; ry = 495 MPa Hardening function

  n y 1=n þ p  ry RðpÞ ¼ r0 rr0 r0 ¼ 1015 MPa; n ¼ 0:15 Damage parameters r1 = 490 MPa; D = 2

1  fn ¼ ð1  fn1 Þ expðx Þ:

Initial porosity Single point : f0 ¼ 102 NT4 : f0 ¼ 5  104 Characteristic length c = 5.2 N

A development of the determinant in (64) leads to the following cubic equation for ~eH : 4 J 2e Þ ¼ 0: X 3  e2eq X þ ð8 det eD þ e 3

1973

ð65Þ

Its solutions can be found with Cardano’s formula. If several real solutions exist, it is suggested to choose the closest to tr e/3. 3.5. Numerical validation A numerical experiment is performed which consists in the pure tension of a specimen under the assumption that the mechanical state remains homogeneous during the whole loading stage: it corresponds to the behaviour of a single material point. Since strain localisation is artificially precluded, the test is not representative of a real structure. However, it enables a focus on the integration of the constitutive law and its accuracy, namely quantitative comparisons between variants of the model and convergence observations with respect to the external load discretisation. The material parameters are those used in [35] with a powerlaw interpolation of the hardening curve, see Table 1. They are representative of A508 steel except for the initial porosity f0 which is greatly over-estimated in order to trigger the porosity growth in spite of the lack of stress concentration. The initial and current length of the specimen are denoted, respectively, l0 and l. The load consists of the stretching (l  l0)/l0; it is an increasing prescribed function of time. The response of the material point is characterised by the component of the Cauchy stress in the loading direction (the other components are zero) and the Jacobian determinant J of the deformation gradient which measures the volumetric strain. The original Rousselier model recalled in Section 2.2 is compared to its revisited formulation, the numerical integration of which relies on the stress approximation introduced in Section 3.1. To demonstrate the interest of the a posteriori correction of the volumetric flow rule (Section 3.4), some computations are led with the correction, the other without it. Moreover, in order to measure the influence of the porosity evolution law, results obtained

ð66Þ

For each of the five variants of the model, three discretisations of the load are applied: stretching increments of 10%, 1% and 0.1%. For the latter discretisation, the results can be considered converged since a finer discretisation (0.01%) leads to very close results. It should be noted that integration of the original Rousselier model with the algorithm described in [36] was not possible for the coarsest discretisation; that highlights the robustness of the algorithm dedicated to the new formulation. The results are gathered in Figs. 2 and 3. Several observations can be made:  in the revisited formulation, the flow rule correction proves necessary because it avoids unphysical results when making use of a coarse discretisation. Even though the version with the integrated porosity law appears more sensitive, this is also true for the version with the differential (explicit) porosity law since the Jacobian is under-estimated (lesser than one) in the first stage of the loading;  for the three remaining formulations (original one and those with stress approximation and volumetric flow rule correction), an accurate solution is already obtained with the intermediate discretisation (stretching increment equal to 1%). Besides, the new formulation enables use of coarser discretisations (up to 10%) which lead to reasonable results (less than 10% error relative to the converged solutions for both stress and volumetric strain responses). Both models with the differential porosity law (original one and new one) converge towards the same response which appears to be more ductile than that obtained with the integrated porosity law. This can be explained by the fact that the integrated porosity law over-estimates the porosity evolution (contribution of all the elastic volumetric strain) while the differential one under-estimates it (no contribution at all of the elastic volumetric strain). Therefore, the question of which evolution law should be used is a matter of choice, even though the integrated form leads to more conservative results. 4. Structural computations As soon as structural computations are concerned, two major issues arise. First, the volumetric strain is highly constrained by the constitutive behaviour, in quite the same manner as in von Mises plasticity where the strain should remain almost isochoric. Therefore, the locking phenome-

1974

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982 differential porosity + stress approximation + volumetric correction

1000

800

Cauchy stress (MPa)

800

Cauchy Stress (MPa)

integrated porosity + stress approximation + volumetric correction

1000

600

400

600

400 time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

200

200

0

0 0

1

2

0

3

1

3

integrated porosity + stress approximation

1000

1000

800

800

Cauchy stress (MPa)

Cauchy stress (MPa)

differential porosity + stress approximation

600

400

600

400 time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

time-step length : 0.1 % time-step length : 1 % time-step length : 10 % 200

200

0

0

0

1

2

3

0

1

Original Rousselier model

800

800

Cauchy stress (MPa)

1000

600

400

600

400

time-step length : 0.1 % time-step length : 1 %

200

integ. poro. + stress approx. diffe. poro. + stress approx. diffe. poro. + stress approx. + hydr. corr. original Rousselier model integ. poro. + stress approx. + hydr. corr.

200

0 1

3

summary for the finest discretisations

1000

0

2

Stretching

Stretching

Cauchy stress (MPa)

2

Stretching

Stretching

2

3

0 0

1

Stretching

Fig. 2. Single material point (SMP): stress–strain response.

2

Stretching

3

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982 differential porosity + stress approximation + volumetric correction

integrated porosity + stress approximation + volumetric correction 1.6

1.6

1.5

1.5 time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

1.3

1.2

1.3

1.2

1.1

1.1

1

1

0

1

time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

1.4

Volume change J

Volume change J

1.4

0.9

2

0.9

3

0

1

Stretching differential porosity + stress approximation

integrated porosity + stress approximation

1.5 time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

time-step length : 0.1 % time-step length : 1 % time-step length : 10 %

1.4

Volume change J

1.4

Volume change J

3

1.6

1.5

1.3

1.2

1.3

1.2

1.1

1.1

1

1

0.9

0.9

0

1

2

3

0

1

Stretching

2

3

Stretching

Original Rousselier model

1.6

summary for the finest discretisations

1.6

1.5

integ. poro. + stress approx. + hydr. corr. diffe. poro. + stress approx. + hydr. corr. original Rousselier model diffe. poro. + stress approx. integ. poro. + stress approx.

1.5 time-step length : 0.1 % time-step length : 1 %

1.4

Volume change J

1.4

Volume change J

2

Stretching

1.6

1.3

1.2

1.3

1.2

1.1

1.1

1

1

0.9

1975

0

1

2

Stretching

3

0.9

0

1

2

Stretching

Fig. 3. Single material point (SMP): volume change response (deformation Jacobian).

3

1976

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

non observed with von Mises plasticity is also expected with the Rousselier law. Special finite elements have to be introduced. Besides, because of softening, strain localisation occurs. As long as a the constitutive law remains purely local, that is without any coupling between neighbouring material points, it results in an ill-posed problem exhibiting spurious mesh dependency, hence nonobjective results. The solution which consists in adjusting the mesh size as a new material parameter is not sufficient to get predictive results because the dependence to the mesh is not only related to the mesh size but also to its orientation, to the kind of shape functions, to the number of Gauss points, etc. Therefore, a nonlocal formulation is required to control strain localisation in a satisfactory way.

For sake of completeness, we give the expressions of the corresponding internal nodal forces. More details regarding the derivation of the expressions and of the tangent operator are given in [37]. The following spatial discretisation of the unknowns is introduced, with Un , Gr and Ps shape functions to be specified: X X U i ðxÞ ¼ U i;n Un ðxÞ; GðxÞ ¼ Gr Gr ðxÞ; P ðxÞ ¼

r2N G s

P Ps ðxÞ:

ð68Þ

s2N P

A straightforward derivation provides the discretisation expression for dL: dLij ðxÞ ¼ dU i;n

oUn ðxÞ ¼ dU i;n Ljn ðxÞ: oxj

H = 50 mm

Mean stretch :

E=

2U a H

Ua Prescribed displacement Computed domain and boundary conditions

Fig. 4. Notched tensile specimen NT4: geometry, computed domain and boundary conditions.

The announced locking phenomenon is known to be caused by the constraints that restrict the kinematics. The solution thus consists in enriching the kinematics. Following [18,19], the volumetric strain is treated as a new unknown. Its relation to the displacement field is weakly enforced by means of a Lagrange multiplier, the pressure P. Actually, a slight modification compared to [18,19] is introduced: the new kinematical unknown G is the logarithmic volumetric strain (G  lnJ) instead of J itself, which leads to simpler expressions. Namely, the variational formulation of the problem reads (with Wext the potential of external forces):   Z

tr s D P 8dU; dG; dP dL : ðs þ pIdÞ þ dG 3 X0  þ dP ðln J  GÞ dX0 ¼ dW ext : ð67Þ

n2N U

Ua

Specime geometry

4.1. Mixed finite elements to avoid volumetric locking phenomena

X

Prescribed displacement

ð69Þ

A set of integration points is also introduced in the reference configuration, with coordinates xg and weight xg. The value of any function q evaluated at xg is denoted q(xg) = q(g). Then, the internal forces, respectively corresponding to

the displacements, the volumetric strain and the pressure, simply read (with summation over repeated indices): X DðgÞ U ¼ xg ðsij þ P ðgÞ dij ÞLjðgÞ ð70Þ fi;n n ; g

frG ¼

X g

fsP ¼

X

 ðgÞ  tr s  P ðgÞ GðgÞ r ; 3

ð71Þ

xg ðln J ðgÞ  GðgÞ ÞPðgÞ s :

ð72Þ

xg

g

The interpolation of the fields U (displacements), G (logarithmic volumetric strain) and P (pressure) remains to be stated. The target finite element is a triangle (2D)/ tetrahedron (3D) because of its convenience when meshing complex structures or making use of adaptive remeshing. The displacement is chosen quadratic (P2 interpolation). A linear (P1) interpolation is suggested in [18] for G and P. However, it has been shown in [38] that this choice results in spurious plastic localisation in 3D and axisymmetrical 2D. These oscillations are thought to be caused by the fact that the kinematics is not enough constrained. Therefore, a richer interpolation is chosen for the Lagrange multiplier: G is P1-interpolated, while P is P2-interpolated. Finally, the spatial integration scheme is classically based on three (respectively, four) integration points for triangles (resp. tetrahedrons). In order to illustrate the above-mentioned effects of the spatial interpolation, a notched tensile specimen (NT4) is considered, see Fig. 4. The material obeys von Mises plasticity with characteristics representative again of A508 steel, see Table 1. The simulation is performed with Code_Aster [39] under the assumption of axial symmetry (2D computation). The results are plotted in terms of tension stress and cumulative plastic strain, without any graphical interpolation (true integration point visualisa-

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982 Standard P2 triangle: axial stress (MPa)

Standard P2 triangle: cumulated plastic strain

3.0

3.0

2.0

2.0

1.0

0.0 0.0

1977

1.0

1.0

2.0

-1000

0

3.0 1000

4.0 2000

5.0

0.0

0.0

1.0

3000

2.0 0

P2-P1-P1 triangle: axial stress (MPa)

3.0 1

4.0 2

5.0 3

P2-P1-P1 triangle: cumulated plastic strain

3.0

3.0

2.0

2.0

1.0

1.0

0.0

0.0 0.0

1.0

2.0

-1000

0

3.0 1000

4.0 2000

5.0

0.0

3.0

2.0

2.0

1.0

1.0

-1000

2.0 0

3.0 1000

4.0 2000

3.0 1

4.0 2

5.0 3

P2-P1-P2 triangle: cumulated plastic strain

3.0

1.0

2.0 0

P2-P1-P2 triangle: axial stress (MPa)

0.0 0.0

1.0

3000

5.0 3000

0.0 0.0

1.0

2.0 0

3.0 1

4.0 2

5.0 3

Fig. 5. Notched tensile specimen NT4: von Mises elasto-plastic response (mean stretch: 8%).

tion), see Fig. 5. Use of standard quadratic triangles exhibits stress checkerboard patterns, characteristic of the volumetric locking phenomenon. The enriched (U, G, P) triangle with, respectively, P2, P1 and P1 interpolations exhibits numerous plastic localisation bands parallel to each other. Finally, only a P2–P1–P2 interpolation provides satisfactory results as well for the stress distribution as for the plastic strain. As will be seen below, the latter choice also seems robust with the Rousselier law.

4.2. Nonlocal formulation To control strain localisation, a coupling between material points has to be introduced. This can be achieved in several ways, leading to various nonlocal models. Here, a strain gradient model is adopted where only some of the strain gradient components are involved in the constitutive law. Actually, the variable that drives the damage mechanism is the volumetric strain (at least as long as void nucleation

1978

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

is not taken into account). A control of the volumetric strain gradient should be sufficient to avoid spurious localisation. Therefore, the variational formulation (67) is altered in the following way:   Z

tr s D  p þ dP ðln J  GÞ dL : s þ pIdÞ þ dG 3 X0  oG odG þc  ð73Þ dX0 ¼ dW ext ; oX oX where c is a new material parameter (a force). Thus, an internal length is implicitly introduced through c. The new quadratic gradient term in G actually corresponds to a penalty term which precludes the emergence of area with high gradients of G, typical of uncontrolled localisation band. More refined interpretations of such a variational formulation have been provided in the references mentioned in introduction. Namely, if the distinction between G and lnJ is temporarily forgotten, a higher order stress tensor T can be exhibited, which implies in particular that the equilibrium equation is also altered: Z Z . odL oG odG oG  dX0 ¼ dX0 with T ¼ cId : c T.. oX oX oX oX X0 X0 ð74Þ As can be seen, the choice of the nonlocal formulation has been led by practical considerations rather than physical justifications in order to obtain a coupling between neighbouring material points. Indeed, the nonlocal model does not introduce higher complexity nor additional variables into the computing procedure since a spatial discretisation of the field G is already available, thanks to the finite elements introduced in Section 4.1: the quadratic gradient term in (73) is easily computed. The only difference regarding the expressions of the internal nodal forces is the introduction of a new term into the nodal force relative to the volumetric strain:   X tr sðgÞ oGr 0 oGr0 frG ¼  P ðgÞ GrðgÞ þ cGr ðxg Þ ðxg Þ : xg 3 oX oX g

The different features that have been presented are triggered in the simulation. The material parameters are recalled in Table 1, including the constant c that weights the nonlocal term, which is set so that the localisation band width is about 0.3 mm in the reference configuration (which corresponds to 1 mm in the final deformed configuration). To reach a proper accuracy, the mesh size goes down to 0.03 mm per element in the refined zone where damage is expected. Because of the refinement, only a 2D simulation is conducted, benefiting from the axial symmetry. But the plane symmetry (with regard to the plane z = 0) is not assumed, in order to enable symmetry loss during the process. As for the former von Mises simulation, P2–P1–P2 mixed finite elements are used. It results in 30,000 triangles, corresponding to 200,000 degrees of freedom. The convergence analysis summed up in Figs. 2 and 3 showed that a reasonable accuracy can be reached with local increments of about 5% cumulated plastic strain (Dp 0.05) in the case of the stress approximation along with the volumetric flow rule correction and the integrated expression of the porosity. Therefore, the displacement increments DUa prescribed for each time-step on the upper and lower walls of the specimen are adjusted so that the plastic increment Dp does not exceed 5% anywhere in the structure. This requirement leads to very small increments as soon as plastic localisation occurs: in fact, 269 increments were necessary to perform the simulation up to a mean stretch of 5%, which corresponds to almost complete fracture of the specimen. The global response of the specimen is plotted in Fig. 6 in terms of normalised load (F/S0) vs. relative diameter reduction (D///0). During the first stage, the response coincides with what would be observed with von Mises law. Beyond point (1), localisation occurs and the effect of increasing porosity becomes apparent at the structural scale. The relative diameter reduction tends to reduce because of the expansion due to increasing porosity. Afterwards, the diameter reduction becomes constant which means that the 600

ð75Þ

(1)

Normalised load F/S0 (MPa)

This is thought to be an advantage over integral nonlocal formulations [30,40]. Moreover, the nonlocal procedure is independent of the constitutive behaviour, as long as the volumetric strain totally controls the damage mechanism (namely, plastic shear bands are not regularised since they do not trigger any volume change). However, variants of the model should be investigated in the future: for instance, the effect of an Eulerian gradient of G instead of the Lagrangian one, which is expected to have an impact on the localisation band width.

(2)

400 (3)

200

(4)

4.3. Numerical simulation of a notched tensile specimen A numerical simulation of a notched specimen NT4 (Fig. 4) is now performed with Code_Aster [39] to show the ability of the Rousselier law to predict ductile fracture.

0

0

0.05

0.1 0.15 Relative diameter reduction

Fig. 6. Global response of the NT4 specimen.

0.2

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

1979

Fig. 7. NT4 specimen: evolution of porosity.

localisation band has deviated outside the plane (z = 0) where the diameter reduction is measured. This analysis is confirmed by the plots of porosity in Fig. 7. The model can fully reproduce the cup–cone fracture of the specimen [5]. The localisation phenomenon is indeed controlled beyond its inception (1), so that the porosity band spreads over ten finite elements. Therefore, the mechanical variables are accurately described, which enables the successive detections of the branching of the band (2), the bifurcation stage (3) when one of the bands is selected (loss of symmetry) and finally the ultimate fracture (4) when the band reaches the boundary of the specimen.

5. Closure The objective consisting in the numerical prediction of ductile damage, from its inception up to final fracture, seems to be reached. To ensure the necessary robustness and objectivity, several topics have been investigated:  Constitutive law To keep it as simple as possible, the Rousselier model has been chosen. It has been recast into the framework of associated plasticity with isotropic hardening. Special attention has been paid to the vertex

1980









E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

point of its yield surface, where the plastic flow direction is not unique but only constrained to belong to a convex cone. Integration of the constitutive law A backward Euler scheme has been applied to discretised the constitutive equations, which ensures unconditional stability. This is not only true for the plasticity variables but also for the porosity. As the model belongs to the class of generalised standard materials, there exists a unique solution to the corresponding set of discrete equations. In addition, the set of equations can be reduced to a single scalar equation. Finite strain A variant of Simo and Miehe’ formalism for finite strains has been adopted in order to preserve the nice structure of the constitutive equations (existence and uniqueness of the local solution), to provide a consistent thermodynamic framework (namely hyperelasticity for the stress– strain relation) and to ensure frame invariance. Simultaneously fulfilling all these properties needs some assumption on the stress measure, which requires a posteriori correction of the flow rule. An acceptable accuracy has been reached, even for quite large strain increments. Finite elements In order to avoid the well-known locking phenomenon related to plastic incompressibility (or rather to the constraint on the volumetric strain), mixed finite elements have been introduced, the degrees of freedom of which are the displacement, the (logarithmic) volumetric strain and the pressure. Relying on recent observations, the pressure field has been discretised by quadratic shape functions, so as to avoid spurious plastic instabilities. Nonlocal formulation A coupling between neighbouring material points is necessary to retrieve a well-posed problem and thus avoid mesh dependency caused by strain-softening. This has been achieved by a strain gradient formulation which limits high gradients of the volumetric strain. It should be noted that this choice combined with the previous mixed finite elements does not result in additional computational cost.

In the future, two different leads should be investigated. On one hand, the physical validity of the chain of components should be confirmed by quantitative comparisons between experiments and computations. On the other hand, steps forward should be made regarding the performance of the computations to enable full 3D simulations. In particular, some attention should be paid to mesh adaptivity in order to significantly reduce the number of finite elements. Appendix A. Generalised flow rule for the rousselier criterion The convex cone that defines the flow directions reads: b 2 K f s : d þ Ak P ^s : d þ Akg: b Cðs; A; f Þ ¼ fðd; kÞ; 8ð^s; AÞ ð27Þ

It means that for a given flow direction (d, k), the driving forces (s, A) have to be solution of the following optimisation problem: b ðs; AÞ ¼ arg max½^s : d þ Ak: ð^s;b A Þ2K f

ð76Þ

Let us proceed step by step with the solution of this problem. First, the following parameterisation of the stress tensor is introduced: s ¼ sH Id þ seq n with sH 2 R;  tr n ¼ 0; n 2 S 1 i:e: neq ¼ 1:

seq 2 Rþ ; ð77Þ

Thus, the maximisation problem can be rewrite as ðsH ; seq ; n; AÞ solution of

max b

b ½^sH trd þ ^seq ^n : dD þ Ak

ð^sH ;^seq ; A Þ2K f seq P0 ^ n2S 1

¼

max b

ð^sH ;^seq ; A Þ2K f

b þ ^seq max ^n : dD : ½^sH trd þ Ak ^ n2S 1

ð78Þ

seq P0

The second maximum is computed by means of a standard tensor property: ( if deq 6¼ 0; dD =deq 2 max ^n : dD ¼ deq ; n ¼ ^ n2S 1 3 undetermined if deq ¼ 0: ð79Þ The maximisation problem then becomes: ðsH ; seq ; AÞ solution of "  # 2 b þ deq^seq max max ^sH tr d þ Ak 3 Qð^sH ;b A ÞP0 06^seq 6Qð^sH ;b AÞ with QðsH ; AÞ ¼ ry  r1 Df expðsH =r1 Þ  A: The second maximum in (80) is readily obtained:   2 2 b with seq deq^seq ¼ deq Qð^sH ; AÞ max 3 3 06^seq 6Qð^sH ;b AÞ ( b Qð^sH ; AÞ if deq 6¼ 0; ¼ b undetermined in ½0 Qð^sH ; AÞ if deq ¼ 0: Finally, the maximisation problem is:

  2 H b ðs ; AÞ solution of max A k  deq þ ^sH tr d 3 Qð^sH ;b A ÞP0  H  ^s 2 2  deq r1 Df exp þ ry deq : 3 3 r1

ð80Þ

ð81Þ

ð82Þ

This is a regular optimisation problem with a single inequality constraint. Its solution involves a Lagrangian L with Lagrange multiplier K P 0:

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

   H 2 2 s LðsH ; A; KÞ ¼ A k  deq þ sH tr d  deq r1 Df exp 3 3 r1 2 y þ r deq þ KQðsH ; AÞ: ð83Þ 3 The solution is characterised by the following set of equations:   oL H 2 d ðs ; A; KÞ ¼ tr d  þ K Df expðsH =r1 Þ ¼ 0; eq osH 3 ð84Þ oL H 2 ðs ; A; KÞ ¼ k  deq  K ¼ 0; oA 3 K P 0; QðsH ; AÞ P 0; KQðsH ; AÞ ¼ 0:

ð85Þ ð86Þ

Three alternatives should be considered to solve this nonlinear system:  K>0



 s tr d ¼ kDf exp ; r1 H

2 deq < k; 3

H

Qðs ; AÞ ¼ 0;

seq  QðsH ; AÞ 6 0;

n undetermined: ð88Þ

 K = 0 and deq > 0  H s tr d ¼ kDf exp ; r1 n¼

2 deq ¼ k; 3

ð89Þ

The three alternatives can be expressed the other way to define the convex cone. C(s, A; f) is the set of flow directions (d, k) that satisfy (

ð94Þ

where the superscript T denotes the inverse of the transpose. The variation of the porosity is deduced from (94) (with H the Heavyside function):  0 if x < 0; dJ df ¼ H ðJ  1Þð1  f0 Þ 2 with H ðxÞ ¼ 1 if x P 0: J ð95Þ And the variation of the elastic strain trial is

D

N ¼ sseq

deEeq ¼

3 eD E : deE : 2 eEeq

Regarding the variation of the Cauchy stress, the following differentiations successively apply: dJ 1 r þ ds; J J ds ¼ E : de:

dr ¼ 

ds ¼ dsðId  2eÞ  2sde; ð97Þ

de ¼ deE :

N 2 fd; trd ¼ 0 and d eq 6 1g if seq ¼ 0;

ð98Þ

In the case of a plastic solution, a distinction should be made between the solution at the vertex point and a regular solution. To that end, we introduce the following indicator:

if seq > 0; v ¼ H ð2leEeq  3lkðx ÞÞ ¼



0 1

for a singular solution; for a regular solution:

ð90Þ H

tr d ¼ kDf expðs =r1 Þ; F ðs; A; f Þ 6 0;

k P 0;

ð91Þ kF ðs; A; f Þ ¼ 0:

ð92Þ

Appendix B. Consistent tangent operator The consistent tangent operator L is a fourth order tensor that relates an infinitesimal variation of the relative deformation gradient input df to the resulting variation of the Cauchy stress dr dr ¼ L : df:

ð96Þ

It remains now to determine the variation of the elastic strain de; it relies on the discrete equations specifically relative to the constitutive law, see Box 3. If the solution is elastic, one has of course:

seq  QðsH ; AÞ ¼ 0;

dD : deq

3 dD ¼ kN with 2

dJ ¼ J f T : df;

seq ¼ 0:

 K = 0 and deq = 0 k ¼ 0;

strain specificities and a last one to compute the Cauchy stress as a post-treatment of the elastic strain. Therefore, the computation of the tangent operator is split in the same way. First, the variation of the deformation Jacobian reads

deE ¼ sym½dfð2en1  IdÞf T ;

ð87Þ

tr d ¼ 0;

1981

ð93Þ

It has been observed that the discrete equations in Box 2 are split in three sets: one relative to the elastic trial, another specific to the constitutive law but without finite

ð99Þ Besides, whatever the solution type, one has

   K K tr eE dGE ¼ Ddf  Df trðdeE Þ exp  : r1 r1

ð100Þ

It is now possible to express the variation of the root x*   2lvdeEeq þ r1 exp  rK1 x dGE   dF ¼ 0 ) dx ¼ : ð101Þ dk 3lv dk þ KGE exp  rK1 x þ dR dx dp dx Then, the variation of deeq is   3 dk  dx : deeq ¼ v deEeq  2 dx

ð102Þ

1982

E. Lorentz et al. / Comput. Methods Appl. Mech. Engrg. 197 (2008) 1965–1982

And finally the variation of e reads

!

deE eD 1 eD deD E de ¼ ðtrðdeE Þ þ dx ÞId þ deeq E þ eeq  2eq E : 3 eEeq eEeq eEeq ð103Þ Through applications of the chain rule, the full expression of the consistent tangent operator may finally be obtained but this does not prove necessary for a numerical implementation since it is based on the evaluation of each contributing quantities (94)–(103). Besides, it can be noticed that the consistent tangent operator is not modified by the volumetric correction of the flow rule (Section 3.4) since the latter is applied a posteriori, namely after computation of the stress tensor. References [1] L. Xia, C.F. Shih, J.W. Hutchinson, A computational approach to ductile crack growth under large scale yielding conditions, J. Mech. Phys. Solids 43 (3) (1995) 389–413. [2] I. Scheider, W. Brocks, Simulation of cup–cone fracture using the cohesive model, Engrg. Fract. Mech. 70 (2003) 1943–1961. [3] J. Rudnicki, J. Rice, Conditions for the localization of deformation in pressure-sensitive dilatant materials, J. Mech. Phys. Solids 23 (1975) 371–394. [4] A. Benallal, R. Billardon, G. Geymonat, Bifurcation and localization in rate-independent materials. Some general considerations, in: Q.S. Nguyen (Ed.), Bifurcation and Stability of Dissipative Systems CISM Courses and Lectures, vol. 327, Springer-Verlag, 1993, pp. 1– 44. [5] J. Besson, D. Steglich, W. Brocks, Modeling of crack growth in round bars and plane strain specimens, Int. J. Solids Struct. 38 (2001) 8259– 8284. [6] S. Forest, E. Lorentz, Localization and regularization, in: J. Besson (Ed.), Local Approach to Fracture, Presses de l’Ecole des Mines de Paris, 2004. [7] V. Tvergaard, A. Needleman, Analysis of cup–cone fracture in a round tensile bar, Acta Metall. 32 (1984) 157–169. [8] A.L. Gurson, Continuum theory of ductile rupture by void growth: part I – yield criterion and flow rules for porous ductile media, J. Engrg. Mater. Technol. 99 (1977) 2–15. [9] G. Rousselier, Finite deformation constitutive relations including ductile fracture damage, in: Nemat-Nasser (Ed.), Three Dimensional Constitutive Relations and Ductile Fracture, North Holland Publishing Company, 1987, pp. 331–355. [10] G. Rousselier, Ductile fracture models and their potential in local approach of fracture, Nucl. Engrg. Des. 105 (1987) 97–111. [11] G. Rousselier, Dissipation in porous metal plasticity and ductile fracture, J. Mech. Phys. Solids 49 (2001) 1727–1746. [12] J. Besson, D. Steglich, W. Brocks, Modeling of plane strain ductile rupture, Int. J. Plasticity 19 (2003) 1517–1541. [13] B. Halphen, Q.S. Nguyen, Sur les mate´riaux standard ge´ne´ralise´s, J. de Me´canique 14 (1975) 39–63. [14] P. Germain, Q.S. Nguyen, P. Suquet, Continuum thermodynamics, J. Appl. Mech. 50 (12) (1983) 1010–1020. [15] E. Lee, Elastic–plastic deformation at finite strains, J. Appl. Mech. 36 (1969) 1–6. [16] J.C. Simo, F. Armero, Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 33 (1992) 1413–1449. [17] E.A. de Souza Neto, F.M. Andrade Pires, D.R.J. Owen, F-bar-based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking, Int. J. Numer. Methods Engrg. 62 (2005) 353–383.

[18] R.L. Taylor, A mixed-enhanced formulation for tetrahedral finite elements, Int. J. Numer. Methods Engrg. 47 (2000) 205–227. [19] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, Solid Mechanics, fifth ed., vol. 2, Butterworth-Heinemann, Oxford, 2000, pp. 328–332. [20] P. Germain, The method of virtual power in continuum mechanics. Part 2: microstructure, SIAM J. Appl. Math. 25 (1973) 556–575. [21] S. BardenHagen, N. Triantafyllidis, Derivation of higher order gradient continuum theories in 2, 3-D nonlinear elasticity from periodic lattice models, J. Mech. Phys. Solids 42 (1) (1994) 111– 139. [22] R. Chambon, D. Caillerie, T. Matsuchima, Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies, Int. J. Solids Struct. 38 (2001) 8503– 8527. [23] N.A. Fleck, J.W. Hutchinson, A phenomenological theory for strain gradient effects in plasticity, J. Mech. Phys. Solids 41 (12) (1993) 1825–1857. [24] M. Gologanu, J.-B. Leblond, G. Perrin, J. Devaux, Recent extensions of Gurson’s model for porous ductile metals, in: INRIA, Ecole EDFCEA-INRIA, Fracture et fatigue, Rocquencourt, France, 1998. [25] E. Lorentz, S. Andrieux, Analysis of non-local models through energetic formulations, Int. J. Solids Struct. 40 (2003) 2905–2936. [26] G. Francfort, A. Mielke, Existence results for a class of rateindependent material models with nonconvex elastic energies, J. Reine Angew. Math. 595 (2004) 55–91. [27] Q.S. Nguyen, Bifurcation and stability in dissipative media (plasticity, friction, fracture), Appl. Mech. Rev. 47 (1) (1994) 1–31. [28] J.C. Simo, C. Miehe, Associative coupled thermoplasticity at finite strains formulation, numerical analysis and implementation, Comput. Methods Appl. Mech. Engrg. 98 (1992) 41–104. [29] E. Lorentz, V. Cano, A minimisation principle for finite strain plasticity: incremental objectivity and immediate implementation, Commun. Numer. Methods Engrg. 18 (2002) 851–859. [30] H. Baaser, V. Tvergaard, A new algorithmic approach treating nonlocal effects at finite rate-independent deformation using the Rousselier damage model, Comput. Methods Appl. Mech. Engrg. 192 (2003) 107–124. [31] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, SpringerVerlag, New York, 1998. [32] Z.L. Zhang, On the accuracies of numerical integration algorithms for Gurson-based pressure-dependent elastoplastic constitutive models, Comput. Methods Appl. Mech. Engrg. 121 (1995) 15–28. [33] V. Cano, E. Lorentz Mode`le de Rousselier en grandes de´formations. Documentation de Code_Aster R5.03.06. , 2002. [34] J.C. Simo, R.L. Taylor, Consistent tangent operators for rateindependent elastoplasticity, Comput. Methods Appl. Mech. Engrg. 48 (1985) 101–118. [35] A. Parrot, Instrumental Charpy test ‘‘Group of the SF2M”: summary of the numerical Round Robin on the simulation of simplified compact tension specimen in the transition. EDF internal report HT26/02/022/A, 2003. [36] G. Rousselier, R. Masson, G. Barbier, Mode`le de Rousselier pour la rupture ductile. Documentation de Code_Aster R5.03.07. , 2004. [37] S. Michel-Ponnelle, E. Lorentz, Ele´ments finis adapte´s aux comportements quasi-incompressibles. Documentation de Code_Aster R3.06.08, forthcoming on , 2006. [38] J. Besson, E. Lorentz, Calcul de dommages: simuler et pre´dire le vieillissement me´canique des installations de production et de transport d’e´nergie. CdM – Ecole des Mines de Paris, rapport annuel ACI/ARMINES 042240, 2006. [39] Code_Aster. Opensource finite element software , courtesy of EDF. [40] T. Drabek, H.J. Bo¨hm, Damage models for studying ductile matrix fracture in composites, Comput. Mater. Sci. 32 (2005) 329–336.