A non-local finite element based on volumetric strain gradient: Application to ductile fracture

Computational Materials Science 45 (2009) 762–767

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A non-local finite element based on volumetric strain gradient: Application to ductile fracture R. Bargellini a,*, J. Besson b, E. Lorentz a, S. Michel-Ponnelle a a b

EDF R&D, LAboratoire de Mécanique des Structures Industrielles Durables, UMR CNRS/EDF 2832, 1 Avenue du Général De Gaulle, 92141 Clamart Cedex, France Mines Paris ParisTech, Centre des Matériaux Pierre-Marie Fourt, UMR CNRS 7633, BP87, 91003 Evry Cedex, France

a r t i c l e

i n f o

Article history: Received 22 November 2007 Received in revised form 11 September 2008 Accepted 15 September 2008 Available online 1 November 2008 PACS: 02.70.Dh 62.20.mm 62.20.fq

a b s t r a c t The aim of this work is to propose a finite element formulation adapted to ductile fracture simulation using Continuum Damage Mechanics and unravelling two main encountered difficulties. First, as ductile damage represents voids nucleation and growth, constitutive behaviors constrain the plastic volumetric strain through damage evolution laws, leading to volumetric locking. A specific formulation is consequently needed. A three-field formulation is defined, in which volume change is treated as a new unknown; its relation to the displacement field is weakly enforced by mean of a Lagrange multiplier. Then, strain and damage localisation occurs due to softening, and leads to spurious mesh dependent solutions. A coupling between neighbouring material points is thus introduced through a volume change gradient term added in the three field formulation. First results show that this formulation permits to control localisation and to unravel mesh dependency. Ó 2008 Elsevier B.V. All rights reserved.

Keywords: Enhanced finite element Damage mechanics Ductile fracture Non-local model Plastic incompressibility

1. Introduction Ductile fracture prediction and control remain pivotal topics in many industrial situations such as gas and oil transportation by pipelines, airplanes maintenance. In these applications, installations cost and size render full-scale experiments difficult or even impracticable. Predictive, robust and reliable numerical analyses are thus of great importance to model crack nucleation, growth and path. In the last decade, finite element (FE) tools have been developed to perform such analyses. Three main approaches are classically used: (i) global approaches to crack resistance (e.g. [1]), (ii) Cohesive Zone Models (e.g. [2,3]) and (iii) Continuum Damage Mechanics (CDM) (e.g. [4–6]). However, global approaches are restricted to the extension of a pre-existing crack using crack growth criteria (either based on the J-integral [7], a Crack Tip Opening Displacement [8] or Crack Tip Opening Angle [9]) and do not represent the fracture process which develops ahead of growing cracks. Cohesive Zone Models constitute a tool for gradu* Corresponding author. E-mail addresses: [email protected], [email protected] (R. Bargellini), [email protected] (J. Besson). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.09.020

ally developing cracks, by lumping the fracture process zone into a zero-width band, the behavior of which is represented by a traction-separation type constitutive equation; however, despite of recent attempts (e.g. [10]), crack path prediction remains difficult because constrained by the mesh (the path coincides with the CZM element boundaries) and under-estimated structure stiffness. X-FEM techniques [11,12] can help overcoming the problem (see [13] or [14]), but criteria for crack growth and bifurcation must also be improved (overall for ductile rupture). This work thus focuses on CDM, which permit to describe the fracture process zone by a continuous constitutive model (e.g. [15] or [16], modified by [17,18]) and the crack itself as some continuous regions loosing their load-carrying capacity. Among the main difficulties encountered when using CDM for ductile fracture simulation, one may mention: (i) volumetric locking phenomena, caused by the constraint on the volumetric strain (guiding damage growth) in the same way as in von Mises isochoric plasticity; (ii) spurious mesh dependency, due to strain localisation characterizing softening constitutive laws [19]. The aim of this work is to propose some answers to these problems. In a first Section, a mixed FE formulation, adapted to quasi-incompressible behavior, is thus presented. It allows to avoid locking

R. Bargellini et al. / Computational Materials Science 45 (2009) 762–767

phenomenon. Section 2 introduces a coupling between material points in the FE formulation to control strain and damage localisation and to avoid mesh dependency. Some tests on a compact tension (CT) specimen show the capabilities of the proposed formulation. 2. A quasi-incompressible finite element formulation The first crucial task is the development of appropriate finite elements which can model the quasi-incompressible nature of the plastic flow without exhibiting volumetric locking behavior. This phenomenon, characterized by strong hydrostatic stress oscillations, is classically observed with von Mises plasticity where plastic strain should remain isochoric, and more generally when a constraint exists on the volumetric flow. In ductile damage behavior, this volumetric strain is related to voids growth, and is constrained by constitutive models (see e.g. [15,16]). Volumetric locking phenomena is consequently expected. As it is thought to be caused by kinematic constraints, two major solutions exist. The first one consists in relaxing these kinematics constraints by reducing the number of integration points (low-order elements). However, at high strain levels reached in ductile materials, spurious hydrostatic stresses can be observed. The second main solution is to enrich the kinematics by additional degrees of freedom (enhanced finite elements). This Section presents a specific enrichment using a volume change field; a three-field FE formulation is defined. Target elements are both quadrangular and triangular. As previously mentioned, ductile fracture is associated to high strain levels. In order to avoid volumetric locking, following [20,21], the deformation gradient F is enriched by a scalar new kinematical unknown G, a priori not displacement-dependant and representing the volume in the current configuration ðJ ¼ detðFÞÞ. The resulting enriched deformation gradient e F is expressed as

 13 G e F¼ F; J

so that detð e FÞ ¼ G

ð1Þ

In order for G to keep its volume change meaning, its relation to the volume J is weakly enforced by mean of a Lagrange multiplier field, namely the hydrostatic stress or pressure P. A three-field (displacement U, volume change G and pressure P) formulation of the HuWashizu type is thus obtained. Calling dL the Eulerian deformation variation, W ext the potential of external forces, s the Kirchhoff stress and sD its deviatoric part, the variational formulation of the problem reads

Z

ðsD þ PId Þ : dL dX0  Z  trs þ  P dG dX0 3 X Z 0 þ ðJ  GÞdP dX0 ¼ dW ext

8ðdU; dG; dPÞ

X0

The interpolation of the three different fields remains to be stated. Displacement field is chosen quadratic. A linear interpolation is suggested in [20] for both pressure and volume change fields; however, preliminary tests showed that this choice results in not enough constrained kinematics. Therefore, a richer interpolation is here chosen for the Lagrange multiplier: G is taken linear, whereas P is taken quadratic. Finally, the spatial integration scheme is based on three (respectively four) integration points for triangles (respectively quadrangles). One triangle and one quadrangle finite elements are thus built (see Fig. 1). On average, a mesh composed of N triangles (resp. quadrangles) contains 6.5 N (resp. 10 N) degrees of freedom. The next step is to enrich this formulation with a non-local ingredient to unravel spurious mesh dependency due to localisation.

3. Volume change gradient based non-local formulation In the first Section, an enrichment of the kinematics by a volume change unknown has been presented. A three-field FE formulation, adapted to behaviors in which volumetric strain is constrained, such as ductile damage models, has been introduced. The second main difficulty is due to the softening character (negativity of the second order work) of these behaviors, leading to strain and damage localisation and mesh dependency [19]. A coupling between material points through non-local model is needed. Different options have been proposed, either based on integral regularization tools [22,23], internal variable gradient or high-order displacement gradient [31]. A synthesis of such approaches can be found in [24] or [25]. Other kind of regularization tools have been proposed, introducing some physical time measure through viscosity (e.g. [26]) or limited damage rate [27,28]. The here proposed approach enters the non-local models framework; the idea is to used the volume change field introduced in Section 2 to build a formulation closed to high-order displacement gradient theory. The ductile damage framework is first reminded; then, the non-local formulation is presented; some simulations show its efficiency to unravel mesh dependency. 3.1. Ductile damage framework Continuum ductile damage models generally introduce a scalar internal variable [15,18] or parameter [16] d, which represents the void volume fraction (porosity) of a ductile media and describes the local loss of material integrity. d enters the definition of a local free energy Uloc also depending on the cumulative plastic strain p and the deformation gradient F:

Uloc ¼ Uloc ðF; p; dÞ ð2Þ

X0

where X0 refers to the initial configuration and Id to the second order identity tensor.

ð3Þ

Evolution of damage is linked to the volumetric strain, as it represents void growth inside the material. More precisely, a mass conservation criteria of an ideal von Mises plastic matrix containing a growing cavity leads (neglecting volumetric elastic strain) to the following relation [16]:

d_ ¼ trðDp Þ 1d

Fig. 1. Three-field quadrangle and triangle elements.

763

ð4Þ

where Dp is the eulerian plastic strain rate. Following simulations are performed in the finite strain framework introduced by Simo and Miehe [29] and modified by Lorentz and Cano [30] using the Rousselier model [15] and a hyperelastic law, without loss of generality. In this model, porosity is introduced as an internal variable; the local free energy Uloc is decomposed into an elastic Ue , a plastic Up and a damage Ud parts:

764

R. Bargellini et al. / Computational Materials Science 45 (2009) 762–767

Uloc ðF; p; dÞ ¼ Ue ðFÞ þ Up ðpÞ þ Ud ðdÞ

ð5Þ

As two processes (plasticity and damage) are here considered, two thermodynamic (driving) forces are defined: loc

A¼



loc

ð6Þ

A yield criterion f = 0, defining when damage and plasticity evolve, depends on the driving forces and on both hydrosatic sD and von Mises norm kskVM of the Kirchhoff stress: Y

sH

er1 er1 Y

1 þ er1

 RðpÞ 6 0

ð7Þ

where D and r1 are material parameters. The chosen form of Ud ðdÞ, defining the link between Y and d, and an associated normality rule with respect to the yield surface (7) permit to retrieve the mass conservation Eq. (4). In addition, a non zero initial porosity dð0Þ ¼ d0 is assumed. Plasticity is represented by a nonlinear isotropic rule based on the von Mises norm of the Kirchhoff stress; the hardening fonction RðpÞ is of power type:

"



rY RðpÞ ¼ r0 p þ r0

1n #n ð8Þ

Material parameters are given in the following table. Elasticity parameters

Hardening rule parameters

Young’s modulus

200 GPa

Poisson’s ratio

0.3

rY r0

315 MPa 300 MPa 0.1

n Damage parameters

D

ry Initial porosity d0

2 315 MPa 0.0001

As previously mentioned, Eq. (4) shows that damage evolution depends on the volumetric (plastic) strain, which is consequently constrained: volumetric locking is expected; the three-field formulation is thus used. In addition, the softening character of these models leads to spurious mesh dependency (see Fig. 4); the next Subsection defines a non-local strategy permitting to unravel this difficulty. 3.2. Non-local formulation The aim of this Section is to give the theoretical background of a non-local formulation adapted to ductile damage simulations, based on the three-field formulation presented in Section 2. 3.2.1. Variational formulation As shown by Eq. (4), damage is governed by local volumetric increase. To control damage and strain localisation, the proposed strategy consists in penalizing important volume change gradients. The starting point is the three-field variational form (2), presented in Section 2. In the spirit of second order displacement gradient methods (see e.g. [31]), it is enriched by a quadratic gradient penalty term and consequently reads:

Z

ðsD þ PId Þ : dL dX0  Z  trs þ  P dG dX0 3 X0

8ðdU; dG; dPÞ

Z

ðJ  GÞdP dX0

X0

þ

Z

X0



oU oU d ¼ RðpÞ; Y ¼  ¼ r1 ln 1d op od

f ðs; Y; RðpÞÞ ¼ kskVM þ r1 D

þ

c

oG odG : dX0 ¼ dW ext oX oX

ð9Þ

where c is a parameter (force) implicitly introducing an internal length and X denotes coordinates in the initial configuration X0 . In ductile fracture the inter-voids distance is usually taken as the internal length (see [32]). Here, the choice of the quadratic gradient of the volume change as a non-local term has been led by practical considerations rather than physical justifications. However, the resulting finite element formulation could be seen as a micromorphic theory (see [33] for a synthesis), and more precisely to microdilatation theory (see [34]), with no specific constitutive law at the micro-scale (the volumetric micro and macro strains are implicitly assumed to be equal): the coupling is introduced directly through Relation (1). Note that the added non-local term is isotropic; for ductile damage applications, this assumption seems reasonable. However, anisotropic gradient term could be used as well, but would certainly lead to complex identification of parameters. Variational form (9) gives some answer to the localisation problem. The additional quadratic term percludes the emergence of zones with high gradients of G, typically encountered when volumetric strain (and thus damage) localisation occurs. The quadratic form acts as a penalization term while minimizing the integral, leading to a control of damage and volume change. In the presented form, the volume change gradient is a Lagrangien one; the effects of an Eulerian gradient should be further investigated. 3.2.2. Discretised expressions The following spatial discretisation of the unknows is introduced, with Un , Gr and Ps the shape functions. Note that, as mentioned in Section 2, Un and Ps are taken quadratic, whereas Gr is taken linear.

8 P i;n i > > > U ðxÞ ¼ n2N U Un ðxÞ > U > > P r < G Gr ðxÞ GðxÞ ¼ r2N G > > > P > i > > P s Ps ðxÞ : P ðxÞ ¼

ð10Þ

s2N P

A derivation of the displacement field provides the discretised expression of dL:

dLij ðxÞ ¼ dU i;n

oUn ðxÞ ¼ dU i;n Ljn ðxÞ xj

ð11Þ

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 displacement, the volume change and the pressure, read (with Einstein convention on repeated indices):

8 P UðgÞ DðgÞ ðgÞ jðgÞ > > > fi;n ¼ g xg ðsij þ P dij ÞLn > > > > < PðgÞ P fs ¼ xg ðJðgÞ  GðgÞ ÞPðgÞ s g > >     > > P trsðgÞ > GðgÞ r0 oGr0 oGr ðgÞ ðgÞ > > f ¼ x þ cG  P G g : r r 3 oXg oXg g

ð12Þ

where dij is the Kronecker symbol. The stiffness matrix is of the form:

2

X0

½K ¼

X g

ðgÞ

K 6 UU ðgÞ xg 6 4 K GU ðgÞ K UP

ðgÞ

K UG ðgÞ

K GG

ðgÞ K GP

ðgÞ

K UP

3

7 ðgÞ K GP 7 5 0

ð13Þ

R. Bargellini et al. / Computational Materials Science 45 (2009) 762–767

765

Expression (13) shows that the stiffness matrix is non-symmetric (indeed K UG –K GU ) and exhibits a 0-term ðK PP Þ; a specific solver, namely MUMPS, is used to deal with these particularities. Note that the non-local term only modifies the K GG component of the stiffness matrix. 3.2.3. Energetic and balance interpretations The previous Section has introduced a non-local variational (consequently weak) formulation. The aim of this Section is to retrieve a possible strong form interpretation. Variational form (9) can be seen as a second order displacement gradient theory, in which only the volumetric part is taken into account; corresponding strong forms are thus quite similar. e which deA possibility is to define an enriched free energy U pends on the volume change gradient:

1 2

e ðF; p; d; $0 GÞ ¼ Uloc ðF; p; dÞ þ c $0 G:$0 G U

ð14Þ

where $0 refers to Lagrangian gradient. A vectorial second order stress T is consequently exhibited:

T ¼ c $0 G

ð15Þ

This second order stress modifies balance equations and boundary conditions. In order to be able to get these relations, some assumptions are made concerning external forces. It is assumed that no body nor boundary double forces (see [31]) act, this means that only classical body forces (denoted b) and surface traction forces (denoted t) are applied. Under these assumptions, integrating Relation (9) by parts and using divergence formula lead to following relations:

8 trðsÞ > > P¼  trð$0 TÞ > > > 3 < divðs  trð$0 TÞ Id Þ þ b ¼ 0 > > > ðs  trð$0 TÞ Id Þ:n ¼ 0 > > : T:n ¼ 0

on X0

ðaÞ

on X0

ðbÞ

on oX0

ðcÞ

on oX0

ðdÞ

Fig. 3. Simplified CT meshes.

ð16Þ

where n is the external normal on the boundary oX0 . System (16) highlights some consequences of the volume change gradient term. Lagrange multiplier P does no more exactly refer to the pressure, but to a so called ’augmented pressure‘ also related to the hydrostatic part of the second order stress gradient. Both classical boundary and balance equations are altered by the same ingredient. Indeed, Eqs. 16(b) and 16(c) can be seen as classical boundary and balance equations just depending on a specific stress measure C ¼ s  trð$0 TÞ Id instead of Kirchhoff stress s. Nev-

ertheless, the volume change gradient term only modifies the hydrostatic part: shear components are not concerned. A specific new boundary condition on T is also added (Eq. 16(d)). These modifications and new equation have to be ensured during simulations. In this Part, theoretical background of a non-local formulation has been presented. It is validated in the next Subsection on a simplified CT simulation. 3.3. Validation on a simplified CT specimen In this section, an academic example of a simplified CT specimen, whose characteristics are given in Fig. 2, is considered. It is

Fig. 2. Geometry of the simplified CT specimen.

766

R. Bargellini et al. / Computational Materials Science 45 (2009) 762–767

Fig. 4. Global response for the different meshes (c = 0).

composed of a vertical band of rigid material, representing the pin, and the CT specimen itself. As loading, a purely vertical displacement is imposed on one node of the rigid band; note that specific boundary conditions (Eqs. 16(c) and 16(d)) are respected. Simulations are performed using six different meshes: three of them (Q1, Q2, Q3) are composed of quadrangles (of respective size in the region of interest 0.2 mm for Q1, 0.1 mm for Q2 and 0.05 mm for Q3), and the three others (T1, T2, T3) are composed of triangles (of the same size as quadrangles). These meshes are presented in Fig. 3. The objective is to show the mesh insensitivity reached with the non-local formulation presented in this Section compared to the local results. All simulations are performed with Code_Aster [35]. Local simulations (with c = 0) are first performed with the different meshes; a simulation for an elastoplastic law is also performed. Fig. 4 presents the obtained global responses with the six different meshes, and underlines mesh dependency; the difference between the quadrangular meshes reaches and overpasses 10%; similar trend is observed for triangular meshes. The more refined the mesh, the smaller the global force. In addition, for a given element size, difference between triangular and quadrangular elements reaches 5%. The next step is to test the non-local approach presented in this Section. The influence of parameter c must be emphasized. A traction on the Q1 mesh is simulated with different c values, ranging from 0 up

Fig. 6. Global response for the different meshes (c = 1.5 N).

to 50 N. Global Force/Pin displacement responses are presented in Fig. 5. As expected, the higher c, the higher the global force necessary to reach a given pin displacement. It is related to the dissipated energy, which decreases with the number of damaged elements; for an infinitely refined mesh and a local model, it tends towards zero. On the contrary, when c tends towards infinity, the material tends to behave as an elastoplastic one. The choice of the non-local parameter c value for the following simulations has been led by practical considerations rather than physical justifications. As the other meshes are more refined than Q1, its value is chosen to have a diffuse damage on some elements for this mesh. In the following non-local tests, the value c = 1.5 N is taken. Fig. 6 presents the simulated global responses with the six different meshes. It shows that mesh dependency is substancially reduced: the difference between triangular meshes is of the percent order; it reaches 3% for quadrangular meshes, because Q1 mesh is not sufficiently refined (the difference between Q2 and Q3 is also of the percent order). In addition, the difference between quadrangular and triangular meshes is of the percent order. The mesh insensitivity is also demonstrated by the simulated damage profiles, presented in Fig. 7 for the four most refined meshes (T2, T3, Q2, Q3). In each case, damaged zone is 4.4 mm long and 0.5 mm high; the most damaged area is located on the ligament of the specimen, 0.5 mm downstream to the initial crack tip, and is 0.7 mm long. Simulated maximal damage values are about 0.6, and differ by about 3%: triangular meshes slightly underestimate damage compare to quadrangular ones. 4. Closure

Fig. 5. Global response for different c values and Q1 mesh.

A specific FE formulation, adapted to ductile damage simulation has been presented. A kinematics enrichment through a volume change field, weakly related to the volumetric strain by Lagrange multiplier, is used to define a three-field formulation of the HuWashizu type. It gives some answer to volumetric locking phenomena, by stabilizing the maximum of the stress in the media, which can be an important data when simulating brittle crack propagation in a plastic media. Then, an original non-local approach, based

R. Bargellini et al. / Computational Materials Science 45 (2009) 762–767

767

Fig. 7. Simplified CT meshes.

on the volume change gradient, permits to unravel spurious mesh dependency. The finally obtained formulation thus represents a reliable ductile damage simulation tool, with some important advantages. (i) The formulation is generic; the non-local procedure, althought altering balance equations, is independant of the constitutive behavior, as long as the volumetric strain totally controls the damage mechanism: non-local treatment may be performed at the element level, and not at the constitutive relation level. This represents an advantage over internal variable gradient approaches. Note that, as second gradient theory, variational form (9) gives the possibility of defining specific constitutive behaviors through e ðF; p; d; $0 GÞ, from which generalized stress ðs; TÞ are a potential U derived (see [36]).(ii) Its cost is limited. Even if some degrees of freedom are introduced by the three-field formulation, their number remains limited compared to full second gradient theories. In addition (iii), the non-local formulation does not introduce higher complexity nor additional variables into the computing procedure; the spatial discretisation of the volume change field is available thanks to three-field FE formulation; as a linear interpolation is used for G, its gradient is easily computed. In the future, some investigations should be made in order to increase the performance of the computations to enable full three dimensional simulations, which remain difficult to reach; remeshing tools, significantly reducing the number of degrees of freedom, seem to be an efficient way to follow. In addition, remeshing gives the possibility of simulating the smooth transition from the continuous damage zone to the discrete crack (see [37]). References [1] [2] [3] [4] [5]

L. Xia, C.F. Shih, J.W. Hutchinson, J. Mech. Phys. Solids 43 (1995) 389–413. V. Tvergaard, J. Mech. Phys. Solids 49 (2001) 2191–2207. M. Anvari, I. Scheider, C. Thaulow, Engrg. Frac. Mech. 73 (2006) 2210–2228. V. Tvergaard, A. Needelman, Acta Metall. 32 (1984) 157–169. J. Besson, D. Steglich, W. Brocks, Int. J. Plas. 19 (2003) 1517–1541.

[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

[20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

[6] M. Mashayekhi, S. Ziaei-Rad, J. Parvizian, J. Niklewicz, H. Hadavinia, Mech. Mater. 39 (2007) 623–636. J.R. Rice, G.F. Rosengren, J. Mech. Phys. Solids 16 (1968) 1–12. A.A. Wells, Br. Weld. J. 10 (1963) 563–570. H. Anderson, J. Mech. Phys. Solids 21 (1973) 337–356. I. Scheider, W. Brocks, Engrg. Frac. Mech. 70 (2003) 1943–1961. T. Belytschko, T. Black, Int. J. Numer. Meth. Eng. 45 (1999) 601–620. N. Moes, J. Dolbos, T. Belytschko, Int. J. Numer. Meth. Eng. 46 (1999) 131– 150. C. Comi, S. Mariani, U. Perego, Int. J. Num. Anal. Met. Geomech. 31 (2007) 231– 238. R. de Borst, J.C.J. Remmers, A. Needleman, Eng. Frac. Mech. 73 (2006) 160–177. G. Rousselier, J. Mech. Phys. Solids 49 (2001) 1727–1746. A.L. Gurson, J. Eng. Mater. Tech. 99 (1977) 2–15. V. Tvergaard, Adv. Appl. Mech. 27 (1990) 83–151. J. Lemaitre, J. Eng. Mater. Technol. 107 (1985) 83–89. A. Benallal, R. Billardon, G. Geymonat, Bifurcation and Localization in RateIndependant 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. R.L. Taylor, Int. J. Numer. Meth. Eng. 47 (2000) 205–227. M. Bellet, Comp. Meth. Appl. Mech. Eng. 175 (1999) 19–40. G. Pijaudier-Cabot, Z.P. Bazant, J. Eng. Mech. 113 (1987) 1243–1256. H. Baaser, V. Tvergaard, Comp. Meth. Appl. Mech. Eng. 121 (2003) 15–28. E. Lorentz, S. Andrieux, Int. J. Solids Struct. 40 (2003) 2905–2936. M. Jirasek, Int. J. Solids Struct. 35 (1998) 4133–4145. L.J. Sluys, R. De Borst, Int. J. Solids Struct. 29 (1992) 2945–2958. A. Suffis, T.A.A. Lubrecht, A. Combescure, Int. J. Solids Struct. 40 (2003) 3463– 3476. O. Allix, J.F. Deü, Eng. Trans. 45 (1997) 29–46. J.C. Simo, C. Miehe, Comp. Meth. Appl. Mech. Eng. 98 (1992) 41–104. E. Lorentz, V. Cano, Commun. Numer. Meth. Eng. 18 (2002) 851–859. P. Germain, J. Méca. 12 (1973) 235–274. A.S. Gullerud, X. Gao, R.H. Dodds Jr., R. Haj-Ali, Eng. Frac. Mech. 66 (2000) 65– 92. S. Forest, R. Sievert, Int. J. Solids Struct. 43 (2006) 7224–7245. R. Fernandes, C. Chavant, R. Chambon, Int. J. Solids Struct. (in press), doi: 10.1016/i.ijsolstr.2008.05.032. Code_Aster, Opensource finite element software available on , courtesy of EDF. R. Chambon, D. Caillerie, C. Tamagnini, C.R. Acad. Sci. Paris, Série IIb 329 (2001) 797–802. J. Mediavilla, R.H.J. Peerlings, M.G.D. Geers, Eng. Frac. Mech. 73 (2006) 895– 916.