Damage anisotropy and its effect on the plastic anisotropy evolution under finite strains

Accepted Manuscript Damage anisotropy and its effect on the plastic anisotropy evolution under finite strains Houssem Badreddine, Khemais Saanouni, Trong Dung Nguyen PII: DOI: Reference:

S0020-7683(15)00049-9 http://dx.doi.org/10.1016/j.ijsolstr.2015.02.009 SAS 8650

To appear in:

International Journal of Solids and Structures

Received Date: Revised Date:

14 May 2014 23 January 2015

Please cite this article as: Badreddine, H., Saanouni, K., Nguyen, T.D., Damage anisotropy and its effect on the plastic anisotropy evolution under finite strains, International Journal of Solids and Structures (2015), doi: http:// dx.doi.org/10.1016/j.ijsolstr.2015.02.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

Damage anisotropy and its effect on the plastic anisotropy

2

evolution under finite strains

3 4

Houssem Badreddine*, Khemais Saanouni and Trong Dung Nguyen

5 6

ICD/LASMIS, UMR 6281, University of Technology of Troyes,

7

12 Rue Marie Curie, CS 42060, 10004 TROYES CEDEX, France

8

*Corresponding author: [email protected]

9 10 11

Abstract

12 13 14 15 16 17 18 19 20 21 22 23

In this work, non-associative finite strain anisotropic elastoplastic model fully coupled with anisotropic ductile

24

Keywords: Finite strain, isotropic hardening, kinematic hardening, ductile damage, plastic anisotropy, damage

25

induced anisotropy.

damage is considered. First, the kinematics of large deformation based on multiplicative decomposition of the total transformation gradient is briefly recalled. Using a thermodynamically-consistent framework, the effective state variables are defined to introduce the effect of the anisotropic damage through the total energy equivalence assumption. For an accurate plastic anisotropy description, non-associative plasticity framework is considered for which equivalent stresses in yield function and in plastic potential are separately defined. The evolution equations governing the overall dissipative phenomena are deduced from the generalized normality rule applied to the plastic potential, while the consistency condition is still deduced from the yield function. Applications are made to a generic material defined by selected values of material parameters and considering different typical strain controlled loading paths as uniaxial tension (UT), plane tension (PT) and equibiaxial tension (ET). For each loading path, the effect of the anisotropic damage evolution on the plastic anisotropy and vice-versa is carefully studied in the context of finite plasticity assuming small elastic strains.

26 27

1

1

1. Introduction

2 3 4 5 6 7 8 9 10

In the field of modern sheet metal forming processes, the use of high strength steels, advanced high strength

11 12 13 14 15 16 17 18 19 20 21

Concerning the mechanical constitutive models, the widely used class of models considers conventional

22 23 24 25 26

During this last decade the tendency is to use more predictive models developed using classical macro or micro-

27 28 29 30 31 32 33 34 35 36 37 38

In metal forming by large deformations, the plastic strains tend to localize inside narrow bands during the metal

steels and aluminum alloys has been steadily increased. The tendency to higher material qualities will increase due to the rising demands of customers and legislation on conflicting regulatory requirements. These are requirements on lightweight design in order to improve the energy consumption and cost efficiency while enforcing the safety requirements. This leads to enhancing the material properties to avoid the strain localization and fracture of the structural components subject to more or less severe and complex loading paths. To meet in these economic, environmental and mechanical strength requirements, the “virtualization” of the forming processes becomes a necessary and inescapable practice. This requires the use of advanced and highly predictive constitutive models together with associated numerical methods.

plasticity. They account for nonlinear isotropic hardening based on yield functions with quadratic or nonquadratic equivalent stresses in the framework of the associated flow rules. Only the initial plastic anisotropy is accounted for through the equivalent stress definition using quadratic or non-quadratic equivalent stresses (Hill, 1948; Vial et al., 1983, Barlat and Lian, 1989; Hill, 1990; Barlat et al., 1991; Karafillis and Boyce, 1993; Hill, 1993; Barlat et al., 1997; Banabic et al., 2003; Barlat et al., 2003; Cazacu and Barlat, 2004; Yoon et al., 2004; Banabic et al., 2005; Barlat et al., 2005; Barlat et al., 2011; Aretz and Barlat, 2013 ); while the yield function remains fixed in the stress space (i.e. no kinematic hardening). This type of isotropic hardening models cannot predict accurately the stress state at each point of the material deforming with large plastic strains under nonsimple local loading paths which are generally non-monotonic and non-proportional even if the motion of the loading tools is monotonic.

macro approaches. This allows taking into account initial and induced anisotropies with non-associative theory, strong coupling with the ductile damage and its effect on the other fields, the damage induced anisotropy, and the induced volume variation … A recent state of the art in this field is given in the recent book (Saanouni, 2012).

deformation often giving rise to the initiation, growth and coalescence of some microvoids and microcracks usually called ductile damage. This ductile damage occurrence is at the origin of induced softening (or negative hardening) which succeeds to the “positive” hardening exhibited by the material during the early stage of its plastic deformation. In order to accurately describe this damage occurrence as well as its effect on the material behavior, two classes of theories are widely used. The first one deals with the Gurson’s based damage theories (Rice and Tracey, 1969; Gurson, 1977; Needleman and Triantafyllidis, 1980; Beremin, 1981; Aravas, 1986; Rousselier, 1987; Onate and Kleiber, 1988; Gelin , 1990; Tvergaard, 1990; Gologanu et al.,1995; Pardoen et al., 1998; Rousselier, 2001; Benzerga et al., 2002 among many others). The second one concerns the continuum damage mechanics or CDM-based theories (Chaboche, 1978; Cordebois and Sidoroff, 1979; Cordebois and Sidoroff, 1980; Kachanov, 1980;

Murakami and Ohno, 1981; Betten, 1986; Kachanov, 1986; Chow and

Wang, 1987; Murakami, 1987; Simo and Ju, 1987; Murakami, 1988; Chow and Lu, 1989; Ju, 1989;

2

1 2 3 4 5

Chaboche, 1992; Lemaitre, 1992; Voyiadjis and Kattan, 1992; Chaboche, 1993; Ladevèze, 1993; Saanouni et

6 7 8 9 10 11 12 13

For the physically-based Gurson’s type approach, the material is assumed to be formed by a homogenous matrix

14

growth and coalescence of the microvoids, i.e. f
= f
nuc + f
gro + f
coal . In the isotropic case, the microvoids and

15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

micro-defects are of spherical or cylindrical shape (Rice and Tracey, 1969; Gurson, 1977; Needleman and

30 31 32 33 34 35 36 37 38 39

In the CDM-based theories, the ductile damage is mostly represented phenomenologically by a scalars in the

al., 1994; Chaboche et al. 1995; Voyiadjis and Kattan, 1996; Chow, 1999; Voyiadjis and Kattan, 1999; Brünig, 2002; Brünig, 2003; Saanouni and Chaboche, 2003; Brünig and Ricci, 2005; Lemaitre and Desmorat, 2005; Chaboche et al. 2006; Saanouni, 2008; Lemaitre et al., 2009; Voyiadjis et al., 2009; Murakami, 2012; Saanouni, 2012 among many others).

containing various types of inclusions (precipitates, second phases, …) sometimes together with an initial distribution of existing microvoids. Under the applied loading path, new microvoids nucleate around the existing inclusions which, together with the initially existing microvoids, grow and coalesce until the initiation of a macroscopic crack (see for example Francois et al, 1993).

Depending on the applied loading path, this

macrocrack propagates inside the deformed part by the same mechanisms of nucleation, growth and coalescence, until the final fracture of the part. From theoretical point of view, the damage is then represented by the microvoids volume fraction (f) the evolution of which is decomposed into three terms related to the nucleation

Triantafyllidis, 1980; Beremin, 1981; Aravas, 1986; Onate and Kleiber, 1988; Rousselier, 1987; Gelin, 1990; Tvergaard, 1990) and are homogenously and uniformly distributed inside the isotropic matrix. Otherwise, when the microvoids and microcracks have a privileged direction, the ductile damage is anisotropic and a microvoids shape factor is used to describe the transformation of the spherical (or cylindrical) microvoids to and elliptical ones to describe the damage induced anisotropy (Gologanu et al, 1995; Benzerga et al, 2002; Benzerga et al, 2004; Steglich et al, 2008; Gao et al, 2009; Benzerga and Leblond, 2010; Steglich et al, 2010a; Steglich et al, 2010b; Keralavarma et al, 2011). The effects of the ductile damage on plastic behavior of the material are considered by introducing the damage variable into the plastic potential. This affects the plastic flow through the normality rule giving a decrease on stresses when the void volume fraction increases under the applied loading. It is worth noting, that the ductile damage effects only the plastic flow leaving the elastic behavior unaffected by the ductile damage i.e. the elastic stiffness is insensitive to the damage occurrence. This means that inside the damaged zones, the stresses decreases but the stiffness is kept unchanged! However, the models proposed by Rousselier (Rousselier, 1987; Rousselier, 2001) and Gelin (Gelin, 1990) which are based on the Gurson’s type theory, account for the stiffness decrease when the ductile damage increases.

isotropic case (Lemaitre and Chaboche, 1985; Lemaitre, 1992; Saanouni et al., 1994; Saanouni and Chaboche, 2003; Chaboche et al. 2006; Saanouni, 2008; Badreddine, 2010; Saanouni, 2012). Also it can be represented by tensors of various ranks for the anisotropic case Cordebois and Sidoroff, 1980; Kachanov, 1980;

as in (Chaboche, 1978; Cordebois and Sidoroff, 1979;

Murakami and Ohno, 1981; Betten, 1986; Kachanov, 1986;

Chow and Wang, 1987; Murakami, 1987; Simo and Ju, 1987; Murakami, 1988; Chow and Lu, 1989; Ju, 1989; Chaboche, 1992; Voyiadjis and Kattan, 1992;

Chaboche, 1993; Ladevèze, 1993; Chaboche et al. 1995;

Voyiadjis and Kattan, 1996; Chow, 1999; Voyiadjis and Kattan, 1999; Brünig, 2002; Brünig, 2003; Brünig and Ricci, 2005; Lemaitre and Desmorat, 2005; Lemaitre et al., 2009; Voyiadjis et al., 2009; Murakami, 2012 a mong many others). For both isotropic and snisotrpic cases damage-effect on the material behavior is accounted

3

1 2 3 4 5 6 7 8

for thanks to the hypothesis of the strain equivalence, stress equivalence or energy equivalence. This kind of

9 10 11 12 13 14 15 16

Within the CDM approach, the description of the anisotropic damage assumes that the microcracks orientations

17 18 19

The more complex situation corresponds to the initially isotropic materials, where the damage directionality is

20 21 22

Various possibilities for the tensorial nature of the damage variables can be considered as presented hereafter.

23 24 25 26 27 28 29 30

•

31 32 33 34 35 36 37 38

•

CDM approaches has been used for damage prediction in metal forming assuming various types of coupling in some published works (Zhu and Cescotto, 1995; Saanouni and Chaboche, 2003; Khelifa et al., 2005; Saanouni, 2006; Saanouni, 2008; Haddag et al, 2009; Badreddine et al, 2010; Saanouni et al, 2010; Saanouni et al, 2011; Saanouni, 2012; Labergère et al, 2011; Issa et al, 2011; Issa et al, 2012). Note that a unified formulations proposed in (Chaboche et al, 2006; Hammi and Horstemeyer, 2007; Saanouni, 2012) have shown the possibilities to gather the Gurson’s type theories and the CDM-based theories to propose a unified formulation of the fully coupled damage-behavior strong coupling accounting for the damage-induced volume variation.

are highly influenced by the initial microstructure of the material and its evolution (texture) as well as by the direction of the principal stresses at each material point. Two situations are considered in the literature for the modeling of the damage anisotropy. The simplest situation corresponds to initially anisotropic material as composite materials with long fibers, where the damage directionality is entirely defined by the initial symmetries of the material. Microcracks are then oriented parallel or perpendicular to the fiber directions and the initial symmetries are preserved. In such case, it is sufficient to consider several scalar damage variables to describe various damage mechanisms.

driven by the history of the applied loading path. For a general case, the damage variable contains its own directionality as the directions of microcracks which evolve depending on the applied loading path.

Three approaches are described in the literature according to the mathematical type of the damage variable (Besson et al, 2010): Damage representation by a set of vectors depending on the main damage mechanisms (Davidson and Stevens, 1973; Krajcinovic and Fonseka, 1981; Suaris and Shah, 1984; Costin, 1985; Talreja, 1985; Costin and Stone, 1987): This choice corresponds well to the brittle damage description at microscale, taking into account the microcracks directions. In the case of materials with a well-defined initial directionality, say composite materials, the damage direction may be fixed by the initial microstructure of the material, independent of the applied loading. For initially isotropic materials, the microcracks directions are driven by the actual stress or strain states and their past history. This approach finds its limitation when the damage develops in plans of different orientations. Damage representation by second-rank tensors (Vakulenko and Kachanov, 1971; Dragon and Mroz, 1979; Cordebois and Sidoroff, 1979; Cordebois and Sidoroff, 1980; Kachanov, 1980; Murakami and Ohno, 1981; Oda, 1983; Betten,1986; Murakami, 1987; Suaris, 1987; Murakami,1988; Chaboche, 1992; Chaboche, 1993; Ladveze, 1993; Hansen and Schreyer, 1994; Voyiadjis, 1996; Desmorat, 2000; Lemaitre et al., 2000; Brünig, 2002; Brünig, 2003; Lemaitre and Desmorat, 2005; Desmorat and Cantournet, 2007; Murakami, 2012; among others.): The geometrical properties of the microcracks may be described, in a continuum sense, by a second-rank tensor. This tensor defines the reduction of the net section of material as in the initial approach of (Vakulenko and Kachanov, 1971). Based on the second-rank damage tensor, various

4

1 2 3 4 5

forms of the damage-effect fourth-rank tensors are proposed (Cordebois and Sidoroff, 1980; Chow and Wang, 1987; Murakami, 1987; Chow and Lu, 1989; Chen and Chow, 1995; Lemaitre and Desmorat, 2005; Murakami, 2012). Such theories offer an interesting possibilities and serious advantage of describing the actual state of damage by only symmetric second-rank tensors. However, the damaged elastic stiffness tensor and the damage release energy tensor are not easy to express in the general case.

6 7 8 9 10 11 12 13

•

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

In (Mazars, 1986; Chaboche, 1992; Chaboche, 1993; Ladevèze, 1993, Desmorat, 2000; Ekh et al., 2003;

29 30 31 32 33 34 35 36 37 38 39

Recently, some works as (Lehmann, 1991; Voyiadjis and Kattan, 1992a; Voyiadjis and

Damage representation by fourth-rank tensors (Chaboche, 1978; Ortiz, 1985; Simo and Ju, 1987; Simo et al., 1987; Yazdani and Schreyer, 1988; Ju, 1989; Chow, 1999): This is the most natural way to introduce the effective stress and the effective behavior of the elastic material without restrictions. It was first proposed by Chaboche (Chaboche, 1979) based on homogenization results for the cracked elastic material. In (Ortiz, 1985; Simo and Ju, 1987) this approach was extended by the use of the elasticity fourth-rank tensor as internal damage variable. This theory was applied exclusively for the simple case of initially anisotropic materials (as composite, concrete,..). However, this approach reveals several difficulties, when the material is initially isotropic and the damage evolution is driven by the applied loading path.

Lemaitre and Desmorat, 2005, Lemaitre et al., 2009) is given the modeling of the loading direction effects on the opening/closure of the microcracks and the microcavities. In fact, when the applied loading path tends to open the microcracks (tension or positive loading path), the material loses its force carrying capacity and some properties decrease due to the presence of these created microcracks. However, when the applied loading path direction is inverted (compression or negative loading paths) tending to close these microcracks, some materials properties recover their initial values due to the microcracks closure. Also the damage growth is higher in tension than in compression. Such approach is called active /passive unilateral damage theory. These theories are based on appropriate spectral decomposition of the tensorial state variables (strain-like and stress-like tensorial state variables) into positive and negative parts and the damage effect or coupling is made differently on the positive and the negative parts. This theory has been developed within the three tensorial type of damage presented above. Although the most general framework for modeling the anisotropic damage considers fourthrank damage tensors (Chaboche, 1978; Ortiz, 1985; Ju, 1989; Voyiadjis and Kattan, 1999; Chow, 1999; Lemaitre and Desmorat, 2005; Lemaitre et al 2009; Murakami, 2012), the damage representation by second-rank symmetric tensors is widely used at least for metallic materials. Kattan, 1992b,

Voyiadjis and Kattan, 1996; Steinmann and Carol, 1998; Voyiadjis and Kattan, 1999; Menzel and Steinmann, 2001, Al-Rub and Voyiadjis, 2003; Brünig, 2002; Brünig, 2003; Ekh et al., 2003; Ekh et al., 2004; Menzel et al., 2005; Voyiadjis et al., 2009, Bamman and, 2010, Brünig and Gerke, 2011a; Vignjevic et al., 2012) extrapolate the anisotropic damage model (using mainly second-rank tensor for the damage variable) to the large deformations. In (Lehmann, 1991; Steinmann and Carol, 1998; Carol et al., 2001; Menzel and Steinmann, 2001; Brünig, 2002; Brünig, 2003; Ekh et al., 2004; Menzel et al., 2005; Voyiadjis et al., 2009, Brünig and Gerke, 2011a, Vignjevic et al., 2012b) a specific kinematics for large deformations is adopted based on a nonlinear metric transformation tensor. This kinematics is characterized by a multiplicative decomposition of this metric transformation tensor into inelastic (damage and plastic) parts and elastic part. So within this approach a kinematic definition of damage is given associated to the damaged inelastic configuration. However, in this

5

1 2

works the objectivity of the damage variable and its evolution equation have not been investigated in the

3 4 5 6 7 8 9 10 11 12 13 14

The present work is an extension of our previous work (Badreddine et al, 2010) to deal with anisotropic ductile

15 16 17 18 19 20

The second part (Section 3) is dedicated to a detailed parametric study to investigate the predictive capabilities of

21

Throughout this paper the classical indicial notations for the tensorial variables are used. Accordingly

22

x, xi , xij , xijkl indicate the zero (scalar), first (vector), second and fourth rank tensors respectively. The classical

23

    inner product in the orthonormal Cartesian frame of basis e j (e1 , e2 , e3 ) is used with the Einstein

24

summation convention as for example xijkl ylk = aij .

25

2. Constitutive framework

26 27 28 29

The present modeling is performed in the framework of the thermodynamic of irreversible processes applied, at a

30

which ai indicate the strain-like variables and Ai their associated or conjugated stress-like variables. This

31 32 33

framework will be used here but limited to the isothermal processes for the sake of shortness. The extension to

framework of the finite strains; and even less the effect of the anisotropic damage on the plastic anisotropy.

damage. Furthermore the general introduction and the conclusion, the paper is organized in two main parts. In the first part (Section 2) the constitutive framework is presented starting by the outline of the kinematics of finite strains. Follow the definitions of the damage tensor and damage-effect tensor and closed with the detailed formulation of the fully coupled constitutive equations in the framework of non-associative finite strain anisotropic elastoplasticity with mixed isotropic and kinematic nonlinear hardening fully coupled with anisotropic ductile damage. Using a thermodynamically-consistent framework, the effective state variables are defined to introduce the effect of the anisotropic damage on the other variables through the equivalent total energy assumption. The equivalent stresses in yield function and in plastic potential are separately defined. The evolution equations governing the overall dissipative phenomena are deduced from the generalized normality rule applied to the plastic potential, while the consistency condition is still deduced from the yield function. Finally, the effect of microcracks closure on the damage rate under compressive loading paths is introduced.

the proposed model using a generic material with various sets of the material parameters. By selecting different values of the material parameters and considering different typical strain loading paths as uniaxial tension (UT), plane tension (PT) and equibiaxial tension (ET); the influence of the damage energy density release rate is first investigated. Then, the effect of the stress triaxiality and the initial and induced anisotropies of the plastic flow are also discussed.

macroscopic scale, to a homogenous (RVE) Representative Volume Element (Germain et al., 1983; Lemaitre and Chaboche, 1985; Saanouni, 2012). In this thermodynamically-consistent framework each physical phenomenon to be described is represented at least by a couple of state variables formally noted ( ai , Ai ) in

the anisothermal processes poses no particular difficulties (see Saanouni, 2012). Let’s begin by the recalling the kinematics of the finite strain in the framework of the rotating frame formulation.

34 35

6

1

2.1. On the kinematics of finite strains

2

The elastoplastic kinematics is based on the concept of intermediate configuration using the classical

3

multiplicative decomposition of the total gradient Fij into elastic Fij e and plastic Fij p parts, according to:

4

Fij = Fike .Fkjp = Qik Vkle . Fljp = Qik Fkj

5

where Fije = Qik .Vkje is the elastic transformation gradient in which Vije is the rotated right elastic stretch tensor

6

and Qij is the rotation tensor which represents the orientation of the rotating frame with respect to the current

7

configuration. Accordingly, the upper bar ( • ) for any tensor indicates its rotated form thanks to the orthogonal

8

rotation tensor Qij . Fijp is the plastic gradient with respect to the “isocline” configuration C p and Fij is the

9 10 11

rotated total transformation gradient as schematized in Figure 1. The satisfaction of the objectivity requirement is

12

the current configuration is transported to the configurations C p and C according to: Tij = Qki .Qlj Tkl for second-

13

rank tensor and Tijkl = Qmi .Qnj .Qrk .Qsl Tijkl for fourth-rank tensor.

14

As shown in (Badreddine et al, 2010), the use of Eq.(1) together with the small elastic strain assumption (i.e.

15

Vije = δ ij + ε ije with ε ije << 1 ) in the expression of the spatial velocity gradient ( Lij = F
ik .Fkj −1 ), gives the additive

16

decomposition of the total strain rate as well as the rotating frame rate (or spin) under the following form:

(1)

based on the rotating frame formulation (see Badreddine et al, 2010 and references given there). Using this framework, any tensorial quantity (namely Tij for second rank tensor and Tijkl for fourth-rank tensor) defined in

17

S Dij = [ Lij ]S = ε
ije + 2 ε ike . Wkj  + Dijp = ε
ijeJ + Dijp

(2)

18

WijQ = Q
ik .Q jk = Wij − Wij = Wij − K ijkl (Vmn ) Dkl

(3)

19

Where Dijp = [ Lijp ]S = [ F
ikp . Fkjp −1 ]S is the plastic the strain rate,

20

eJ strain rate and the relative total spin (with respect to Jaumann rotating frame rate) respectively; while ε
ij

21

defines the Jaumann rate of the small elastic strain tensor. The total spin is defined as Wij = [ F
ik . Fkj −1 ] A . Note

22

that the superscripts [i]S and [i ]A stand for symmetric and skew-symmetric parts of [i] respectively. The relative

23

total spin Wij is postulated a priori as a linear function of the total strain rate: Wij = K ijkl (Vmn ) Dkl in which

24

K ijkl (Vmn ) is a fourth-rank tensor depending on the rotated total right stretch tensor Vmn . Note that the Jaumann

25

rotating frame is defined by K ijkl ≡ 0 and various other cases can be defined as discussed in (Badreddine et al,

26

2010). Note that this simple kinematics is quite different from other kinematics which introduces an extended

27

multiplicative decomposition of the transformation gradient of the form Fij = Fike .Fklp .Fljd where Fijd is an

28 29

additional metric associated with damage having the sense of “damage driven deformation” as described in

Dij = [ Lij ]S and Wij are respectively the total

(Brünig, 2002; Brünig, 2003).

7

1 2

To simplify the mathematical notations, the upper bar indicating the rotated tensorial variables is deleted throughout the paper.

Ct

Vije

Fij

Qij

Fij

C0

Ct
Fijp

Vije

Fijp

Cet

p Ct

Qij Fije

3 4 5

Figure 1. Kinematic decomposition of the total deformation gradient: rotating frame concept

2.2. Damage tensor and damage-effect tensor

6 7 8 9 10 11 12 13 14

As discussed in the introduction, in continuum damage mechanics framework, the ductile damage can be

15 16 17 18 19

material parameters (Murakami, 1981; Cordebois and Sidoroff, 1980; Betten, 1986; Chaboche, 1992; Ladeveze,

20

δ S (containing some defects as microcracks and microvoids) oriented by its normal ni is replaced in a

21

homogenous way by the reduced elementary surface δ S (free from micro-defects) with a normal ni obtained

22

by applying the damage-effect operator (δ ij − d ij ) on the elementary surface δ S ni :

23 24

represented either by scalar variables or by tensorial variables of different ranks. For the textured polycrystalline metallic materials subject to complex loading paths, the microcracks develop in preferred directions controlled by the loading directions. The classical isotropic damage theory, in which damage is represented by a scalar variable (d), is not suitable. Accordingly, the use of anisotropic damage theory in which the damage variable is of tensorial nature is mandatory. This is particularely the case when the anisotropuic damage interacts with the anisotropy of plastic flow with hardening (full coupling) affecting significantly, the direction of the plastic flow and microcracks evolution. For the most metallic materials, it has been shown that a second-rank damage tensor dij is sufficient to describe the orthotropic distribution of microcracks, while leading to a smaller number of

1993; Hansen, 1994; Voyadjis and Kattan, 1992a; Voyadjis and Kattan, 1992b; Voyadjis and Kattan, 1996; Voyadjis and Kattan, 1999; Brünig, 2002; Brünig, 2003; Lemaitre and Desmorat, 2005; Voyiadjis et al, 2009; among others). In Figure 2, is given the schematization of the three orthogonal families of micro-defects represented by the second-rank damage tensor dij inside a typical RVE. In this figure the elementary surface

(δ ij − d ij ) n j δ S = ni δ S

(4)

Consequently, the second-rank damage tensor dij has the following form:

8

1

 d11  d ij  =  d12 d  13

d12 d 22 d 23

d13   n11 n21 n31   d1 0 0   n11 n12 n13       d 23  =  n12 n22 n32   0 d 2 0   n12 n22 n23  =  Pki d kl∆ Plj  d 33   n13 n23 n33   0 0 d 3   n31 n32 n33        d kl∆   

[ Pki ]

(5)

 Plj 

0 ≤ d k ≤ 1 ( k ∈ {1, 2,3} ) represent the

2

where the eigenvalues of the second-rank symmetric damage tensor

3

orthogonal families of micro-defects and

4 5 6 7

transfer matrix beteewn the observation basis and the damage principal basis. The isotropic damage d

8

condition that can be considered is the whole damage eigenvalues reach the critical value: d k = dc . To simplify,

9

in this work, we consider that the RVE is failed when the highest damage eigenvalue reaches its critical value:

( n ) defines its principal basis. The orthogonal tensor P k i

ij

defines the

corresponds to the case where micro-defects density is independent from the orientation of each surface. In this isotropic case the total failure of the RVE is achieved when the final fracture condition is reached i.e.

d = dc ≈ 0.99 . However, the failure condition in the case of anisotropic damage is not trivial. The final fracture

Max ( d k ) = d c with d c = 0.99

10

(6)

k =1..3

11 12

Figure 2. Schematic representation of the anisotropic damage in a typical RVE.

13

Since in the framework of the Continuum Damage Mechanics (CDM), the real deformed configuration Ct

14

which contains several micro-defects (microcracks and microvoids), is replaced by a fictive undamaged

15

configuration C t free from defects. The mechanical state of the material is then described in the configuration

16

Ct

17

representing the kinematic hardening; (iii)

18

representing the ductile damage. While, in the fictive configuration C t , the mechanical state is described by the

19

couples of effective state variables ( εije , σ ij ) , (α ij , X ij ) and ( r , R ) so that the total energy defined in the tow

by the following couples of state variables:

( r, R )

(ε

e ij

, σ ij ) representing the elastoplasticity;

representing the isotropic hardening, and

(α

ij

(d

, X ij )

ij

, Yij

)

9

1 2 3 4

configurations is the same (Saanouni et al, 1994; Saanouni and Chaboche, 2003; Murakami, 2012; Saanouni, 2012). If the total energy is assumed to be the sum of three parts namely: (i) The elastic enegy, (ii) the energy stored in the kinematic hardening and (iii) the energy stored in the isotropic hardening, then the energy equivalence principle is applied separately to each energy component, leading to:

εije = M ijkl ε kle

5

αij = N ijkl α kl

6

r = 1 − dij

7

γ 2

and and

−1 σ ij = M ijkl σ kl

(7)

X ij = N ijkl −1 X kl

(8)

R

r and R =

1 − dij

(9)

γ 2

8 9

Through this choice, it is worth noting that damage affects differently the elasticity, the kinematic hardening and

10

tensor M ijkl , the kinematic hardening through the fourth-rank damage-effect tensor N ijkl and the isotropic

11

hardening through the scalar valued function

12

and

13

−1 −1 −1 −1 admit inverse tensors noted M ijkl and N ijkl defined so that M ijrs M rskl = δ ij δ kl and N ijrs N rskl = δ ij δ kl . Even if

14

these operators may be quite different, they are generally taken equal (i.e. N ijkl = M ijkl ) for the sake of simplicity

15

as well as to reduce the number of material parameters (Saanouni and Chaboche, 2003; Saanouni, 2012).

16

Different forms of the damage effect operator M ijkl (or M ijkl −1 ) are considered in the literature (Murakami,

17 18 19 20 21 22

1981; Cordebois and Sidoroff, 1980; Betten, 1986; Chow and Lu, 1989; Voyadjis and Kattan, 1992a; Voyadjis

the isotropic hardening. In fact, the damage affects the elastic behavior through the fourth-rank damage-effect

1 − d ij

γ 2

where dij

2

=

1 tr (d ij2 ) 3

=

1 3

dii 2 is a damage norm

γ is a material parameter (Saanouni, 2012). These symmetric and positive-definite fourth-rank tensors

and Kattan, 1992b; Voyadjis and Kattan, 1996; Chen and Chow, 1995; Voyadjis and Kattan, 1999; Lemaitre and Desmorat, 2005; Voyiadjis et al, 2009; Murakami, 2012; among many others). Five symmetric fourth-rank damage effect tensors having the form of diagonal matrices are exhaustively discussed in (Murakami, 2012). In this work, the form proposed by (Lemaitre and Desmorat, 2005) is generalized and rewritten here, in the rotated configuration, in the following compact form:

M ijkl =

23

1 2

(h h ik

jl

+ hil h jk ) − 13 ( hij2δ kl + hkl2 δ ij ) + 91 hnn2 δ ijδ kl + 31 (1 − d rs )η δ ij δ kl 1

(10)

24

where hij = (δ ij − d ij )1 / k is the second-rank effective damage tensor defined on the rotated configuration. This

25

form M ijkl has the advantage to affect differently the deviatoric and hydrostatic parts of each tensorial state

26

D H + M ijkl variable. Indeed, M ijkl = M ijkl is decomposed into two parts:

27

•

28 29

A deviatoric part acting on the deviatoric part of the stress-like fields :

D M ijkl =

•

1 2

(h h ik

jl

+ hil h jk ) − 31 ( hij2δ kl + hkl2 δ ij ) + 91 hnn2 δ ij δ kl

(11a)

A hydrostatic part acting on the hydrostatic part of the stress-like fields:

10

1

H M ijkl = 31 (1 − d rs )η δ ij δ kl 1

(11b)

2

Note that in this paper the subscript [i]D denotes the deviatoric part of a tensor. The parameters η and k are

3

material constants governing the effect of the damage (coupling) on the hydrostatic and deviatoric stress parts

4

respectively. The form of M ijkl initially proposed by (Lemaitre and Desmorat, 2005) is recovered if k=2; while

D

5

H M ijkl

6

H M ijkl = 31 (1 − η drs )δ ijδ kl . On the other hand, the first term of the RHS of Eq.(11.a) defines the damage 1

7

is

different

from

the

one

proposed

by

(Lemaitre

and

Desmorat,

2005)

which

is

effect tensor proposed by (Chow and Wang, 1987) if k=2.

8 9 10

Here, the parameter k is fixed to the value k=4 in order de recover the classical relationships developed in our

11

d ij = d rs .δ ij and hij = (1 − d rs )1/ 4 δ ij ; while the damage-effect operator reduces to M ijkl = 1 − d δ ik δ jl with

12

η =1 so that the equations Eq. (7-9) become:

previous works for the isotropic damage case (Hammi, 2000; Saanouni and Chaboche, 2003; Badreddine et al, 2010; Saanouni, 2012). Recall that, in the fully isotropic damage case, the damage tensor reduces to 1

1

13

 e σ ij e εij = 1 − d ε ij and σij = 1− d   X ij αij = 1 − d αij and X ij = 1− d   R γ  r = 1 − d r and R = 1− d γ 

14

Using the Voigt notations, the damage-effect operator M ijkl has the following matrix form in the damage

15

eigenvectors basis:

16

 49 h12 + 19 ( h2 2 + h3 2 ) 19 h3 2 − 92 ( h12 + h2 2 ) 91 h2 2 − 92 ( h1 2 + h3 2 )  2 2 2 2 2 2 4 h + 19 ( h1 + h3 ) 19 h1 − 29 ( h2 + h3 ) 9 2  2 2 2 4  h + 91 ( h1 + h2 ) 9 3 M ijkl =    sym    1 1 0 0 0 1   1 1 0 0 0    1 0 0 0 η 1 + 3 (1 − drs 1 )   0 0 0   sym 0 0   0  

(12)

1 2

0 0

0 0

0

0 0

h1h2 1 2

h2 h3

0   0  0   0  0   1 hh  2 1 3 

(13)

11

1

where hi are the eigenvalues of the effective damage second-rank tensor hij . This form is symmetric but not

2

diagonal shaped matrix compared to forms discussed in (Murakami, 2012).

3

2.3. Fully coupled constitutive equations

4

2.3.1. Stress-like state variables

5

The Helmholtz free energy

6

space is taken as a state potential (Germain et al., 1983). To introduce the effect of damage, this potential is built

7

in the fictitious rotated configuration using the effective state variables ( ψ (εije , α ij , r ) ). Note that in this fictitious

8

space, the kinematics defined through Eq. (1) to Eq. (3) remain unchanged. In this framework of large plastic

9

strains, the rotated damage tensor d ij = Qki .Qlj d kl is used which has the same eigenvalues as dij but is

10

characterized by rotated eigenvectors with the rotation tensor of the rotating frame Qij . For the sake of the

11

notation simplicity, and as mentioned above, the tensor dij represents the damage tensor with respect to the

12

rotating frame.

13 14

Assuming that the plastic strain and the hardening do not affect the elastic properties of the material (no

ψ , positive and convex function of all the strain-like state variables in the strain

texture), the state potential can be additively decomposed as following:

1 2

1 3

1 2

ρψ (εije , α ij , r) = ρΨ e (εije ) + ρψ p (αij , r) = Λ ijkl εijeεkle + Cαijαij + Qr 2 γ 1 1 1 e M ijmn Λ ijkl M klrsε mn ε rse + CM ijmn M ijrsα mnα rs + Q(1 − duv 2 )r 2 2 3 2 1  1 1 e  2 ε rse + C mnrsα mnα rs + Qr = Λ mnrsε mn 2 3 2

=

15

(14)

16

where Λ ijkl = 2µe δ ik δ jl + λe δ ij δ kl is the symmetric fourth-rank tensor defining the material elastic properties of

17

the undamaged RVE for the isotropic elasticity where µe and λe are the classical Lame’s constants; the

18

parameters C and Q are the kinematic and the isotropic hardening modules respectively. Their damage-

19

γ D D 2    ) respectively. affected values are given by: Λ mnrs = M ijmn Λ ijkl M klrs , C mnrs = 3 CM ijmn M ijrs and Q = Q (1 − d ij

20

  Clearly, the operators Λ mnrs and Cmnrs lose their initial isotropy due to the damage-induced anisotropy through

21

the operator M ijkl . It is worth noting that, since the back stress tensor X ij has a deviatoric nature so its

22

associated effective tensor X ij is also deviatoric. Eq. (14) becomes as flows with considering the isotropic

23

elasticity and using Eq.(11):

24

2

1 2

1 3

1 2

2 2 ρψ = µε rkeD hkm ε mleD hlr2 + ke (1 − 13 d nn )η ε rre 2 + Cα rk hkm α ml hlr2 + (1 − ( 13 d nn2 )γ / 2 )Qr 2

(15)

25

where ke = λe + 2 µe / 3 is the elastic bulk modulus. The stress-like variables are easily obtained using Eq.(13) in

26

their compact form:

12

1

σ ij = ρ

∂ψ e  = Λ ijkl ε kle = M ijmn Λ mnrs M rskl ε kleD + ke (1 − d uv 1 )η (ε tte )δ ij ∂ε ije

(16)

∂ψ p  2 D D α kl = Cijklα kl = CM ijmn M mnkl ∂α ij 3

(17)

X ij = ρ

2

R(r) = ρ

3

Yij = − ρ 4

5

(18)

∂ψ = Yije + Yijα + Yijr ∂dij

D 1 ∂M mnkl 1 ∂M mnkl e D M klrs ε rse − C α mnα rs + 16 γ duv =− Λ klpq M pqrsε mn 2 ∂dij 3 ∂dij

(19)

γ −2

2

2

Qr dij

By concedering the relation Eq.(15), the relations Eq.(16) to Eq.(19) become:

∂ψ e eD 2 D = 2 µ ( him2 ε ml hlj ) + k e (1 − 13 d nn )η ε rre δ ij ∂ε ije

(20)

D ∂ψ p 2 = C ( him2 α ml hlj2 ) ∂α ij 3

(21)

∂h 2 1 ∂ψ e 2 ε mleD lr + keη (1 − 13 d nn )η −1 ε rre 2δ ij = −2 µeε rkeD hkm ∂d ij ∂d ij 6

(22)

σ ij = ρ

6

X ij = ρ

7

8

γ ∂ψ p  = Q r = Q(1 − dij 2 ) r ∂r

Yije = − ρ

9

Yijα = − ρ

∂h 2 ∂ψ pα 2 2 = − Cα rk hkm α ml lr ∂dij ∂d ij 3

(23)

10

Yijr = − ρ

∂ψ pr 1 = 6 γ d kl ∂d ij

(24)

γ −2

2

Qr 2 ⋅ d ij

11

In Eq. (22) and Eq. (23) appears the derivative of the effective damage function ∂hlr2 / ∂dij which is developed in

12

detail in the Appendix 1.

13

The Eq.(19) and Eq.(22) to Eq.(24) represent the damage energy density release rate, symmetric and second-rank

14

tensor which is divided into three contributions coming from elasticity Yije , kinematic hardening Yijα and

15

isotropic hardening Yijr (Saanouni et al, 1994; Saanouni, 2012). The close examination of the elastic contribution

16

Yije given by Eq. (22) shows that this latter contains separately the effect of the deviatoric strain (the first term of

17 18

the RHS of Eq. (22)) which is expected to govern the shear damage growth; while the second term of the RHS of Eq. (22), which is spherical in nature, controls the isotropic growth of voids under hydrostatic stress state.

13

1

The kinematic and isotropic hardening contributions given by Eq. (23) and Eq. (24) reflect the effect of the

2

loading history. Particularly, the isotropic hardening contribution Yijr (Eq. (24)) is dependent on the damage

3

directions due to the choice of the damage norm d ij

4

isotropic character of the hardening.

5

Finally, let’s note that if the damage is isotropic d ij = d δ ij and η = 1 leading to M ijkl = 1 − d δ ik δ jl , the damage

6 7

energy density release rate reduces to its classical form developed in our previous works (Saanouni et al. 1994;

8

2.3.2. Strain-like flux variables

2

giving an anisotropic damage evolution despite the

Saanouni and Chaboche, 2003; Badreddine et al, 2010; Saanouni, 2012).

9 10

The thermodynamical admissibility of the model is ensured by fulfilling the following Clausius - Duhem

11

Φ = σ ij Dijp − X ijα
ij − Rr
+ Yij d
ij ≥ 0

12 13 14 15

It is worth noting, that the plastic spin does not contribute in the plastic dissipation defined by Eq. (25) since the

16

The flux variables ( Dijp , α
ij , r
,

17 18 19

be defined in such a manner that the inequality Eq.(25) is a posteriori identically satisfied. Assuming the rate

20 21

F , both scalar valued, positive and convex functions of their main arguments in the stress space. It is worth noting that two yield functions and two dissipation potentials can be introduced relative to the plasticity

22

( f p ,F

23 24

2009; Saanouni, 2012, Vignjevic et al., 2012). In this study, a single surface model is used to describe the

25

f σ ij , X ij , R = σ ijD − X ijD − R − σ y ≤ 0

inequality which is the result of combining the first and the second thermodynamics conservation laws:

small elastic strain assumption is adopted. Two other cases for which this happens concern either if the elasticity is isotropic (which is not the case here due to the coupling with the anisotropic damage) or if the Jaumann derivatives is considered i.e. K ijkl ≡ 0 (see Section 3 Eqs. (17) to (28) in Badreddine et al. 2010).

d
ij ) associated to the stress-like variables defined by Eq.(16) to Eq.(19), should

independent plasticity dissipation, the generalized standard material framework (Germain et al., 1983) is used assuming the non-associative plasticity theory. This leads to introduce a yield function f and a plastic potential

p

) and to the damage ( f d , F

d

) (Zhu and Cescotto, 1995; Saanouni and Chaboche, 2003; Voyiadjis,

damaged elastoplastic behavior, with the same yield function and plastic potential chosen as:

(

26

F

27

F

28

(25)

(σ p

ij

(σ

F

d

)

)

 ,d = F , X ij ,R;Y ij ij ij

p

)

, X ij ,R = σ ij − X ij

(Y ,d ) = ij

ij

(26)

c

(σ p

)

(Y ,d )

(27)

 b R − 1 ) + 43Ca X ij X ij + R( 2Q

(28)

ij

, X ij ,R + F

Yij − Y0

S

( s + 1) (1 −

d

3

d rs

1

)

β

S

ij

ij

s +1

(29)

14

1

where σ y is the initial yield stress, a and b characterize the nonlinearity of the kinematic and the isotropic

2

hardening, respectively. Similarly, the material parameters ( β , S , s , Y0 ) characterize the nonlinear ductile

3

damage evolution. The notation

4

x ≤ 0 . The norm Yij

5 6

stresses characterizing the yield criterion and the plastic potential respectively. In this work the well-known

3

x indicates the positive value of x i.e. x = x if x > 0 and

is defined by: Yij 3 = 3Yij Yij while the stress norms i c and i

p

x = 0 if

are the equivalent

anisotropic Hill48 quadratic equivalent stress is adopted for both:

σ ijD − X ijD =

7

i

(σ

D ij

)

(

)

(σ

I − X ijD H ijkl σ klD − X klD =

D ij

I − X ijD ) H ijkl (σ klD − X klD ) with I = {c , p} (30)

c p H ijkl and H ijkl are the initial anisotropic operators of the undamaged material. Each one of these two

8

where

9

operators is characterized by 6 martial constants ( F i , G i , H i , Li , M i , N i ) to describe the initial orthotropic

10

plasticity. The following remarks are highlighted:

11

•

The classical associative theory is recovered when taking

12

•

c p The classical isotropic von Mises equivalent stress is recovered when taking H ijkl = H ijkl = δ ik δ jl − 13 δ ijδ kl .

13

c D −1 c D −1 p D −1 p D −1 The operators H ijkl = M ijmn H mnrs M rskl and H ijkl = M ijmn H mnrs M rskl are the actual anisotropic operators of the

14

damaged RVE, indicating that the plastic anisotropies evolve with the anisotropic damage growth. Note that,

15

using

16 17

character of the plastic strain tensor by ensuring a purely deviatoric tensor for the normal to the plastic potential

18

In this non-associative formulation, if the generalized normality rule is applied to the plastic potential F instead

19

of the yield function f , the following evolution equations governing the overall dissipative phenomena are

20

obtained:

21

22

23

24

D −1 M ijkl

instead of

M ijkl −1

c p H ijkl = H ijkl ,

c p in the definition of H ijkl and H ijkl above, aims to guaranty the deviatoric

(see Eq. (31)).

∂F D −1 p Dijp = λ
nkl with nijp = = λ
nijp = λ
M ijkl ∂σ ij

α
ij = −λ


rs

rs

∂F = Dijp − λ
aα ij ∂X ij

∂F 1 r
= −λ
= λ
( ∂R 1 − dij ∂F λ
d
ij = λ
= ∂Yij 1 − drs

(

(

p H ijkl σ klD − X klD σ D − X D

γ

1

)

(31)

p

(32)

− br )

(33)

2

Ymn − Y0 β

)

S

3

3Yij

S

Ykl

(34) 3

15

1

Note that nijp = ∂F ∂σ ij is the outward normal to the plastic potential F

2

nijp = ∂F ∂σ ij is the outward normal to F

3 4

The expressions Eq. (31) and Eq. (32) clearly show the effect of damage induced anisotropy on the evolution of

5

directions) is governed by the damage energy density release rate tensor Yij . This allows accounting for the

6

contribution of the stress state and loading history on the damage evolution through the various terms of damage

7

energy (elastic and hardening contributions). The threshold energy Y0 is introduced to account for the fact

8 9

observed for wide range of materials assuming that the evolution of damage takes place only if the energy stored

10

Remark: The necessary has been made in the definition of the effective (damage affected) properties (see the

11

p expressions of Cijkl and H ijkl ) in order to have a purely deviatoric normal nijp and consequently a purely

12

deviatoric plastic strain rate tensor defined by Eq. (31).

13

For this case of time-independent plasticity, the plastic multiplier λ
is deduced from the classical consistency

14

condition applied to the yield function i.e. f
= 0 if f = 0.

15

expression of λ
can be obtained as (see Appendix 2):

in the effective stress space ( σ ij ) and λ
is the plastic multiplier.

plastic flow and kinematic hardening. Eq. (34) shows that the damage evolution (in terms of norm and

during the loading reaches a given value.

After some algebraic calculations, the final

 1 ∂f ∂σ ik  (δ δ − ϒklmn ) Dmn if
λ =  ℵ ∂σ ij ∂ε ejl km ln  0 otherwise 

16

17

in the stress space ( σ ij ), while

where

f
= 0, f = 0

(35)

ℵ is the scalar valued tangent damaged elastoplastic hardening modulus given by: ℵ=

18

∂f ∂σ ij ∂F ∂f ∂X ij ∂F ∂f ∂R ∂F + + e ∂σ ij ∂ε kl ∂σ kl ∂X ij ∂α kl ∂X kl ∂R ∂r ∂R

∂F − ∂Yij

 ∂f ∂σ kl ∂f ∂X kl ∂f ∂R ∂f + + +   ∂σ kl ∂d ij ∂X kl ∂dij ∂R ∂dij ∂dij

  

(36)

19

The mathematical expressions of the different derivatives in Eq.(36) are given in the Appendix 2. The fourth-

20

rank tensor ϒijkl appearing in Eq. (35), is due to large strain assumption and represents the rotating frame choice

21

via the fourth-rank operator Kijkl (see Eq. (2) and Eq. (3)) . By using the indicial notations, this operator can be

22

expressed function of the elastic strain tensor by:

23

e e K mjkl − ε mj K imkl ϒijkl = ε im

24 25

Using the state relations Eq. (16) to Eq. (19) and the evolution relations Eq. (31) to Eq. (34), the relations Eq.

(37)

(35) and Eq. (36) write as:

16

λ
=

1

1 c nij Λ ijkl ( δ kmδ ln − ϒ klmn ) Dmn ℵ

 + C ) n p − an c X + Q − ℵ = nijc (Λ ijkl ijkl kl ij ij

bR 1 − d rs

2

Ymn − Y0

3

−

(

3

2 1 − d rs 3 4 5 6

with

)

S

γ 2

  Yij  c ∂Λ ∂C klpq γ Qr klpq e c n n ε − α pq + kl pq kl  ∂d ij ∂d ij Yuv  3 1 − d ot 3 

D −1 c nijc = ∂f / ∂σ ij = M ijkl nkl is the outward normal to the yield surface

(

)

c nijc = H ijkl σijD − X ij / σ klD − X kl

c

f

γ 2

d xz

γ −2 2

 (39)  d ij   

in the stress space and

is the outward normal in the effective stress space.

The accumulated plastic strain rate is then given using plastic work definition and by combining the relations Eq. (30) and Eq. (31) :

(σ ijD − X ij ) Dijp = σ klD − X kl

7 8

1

β

s

(38)

p

. p
⇒ p
= λ
and

t


p = ∫ pdt 0

(40)

2.3.3. Effect of the microcracks closure on the damage growth

9 10 11 12 13 14 15 16 17 18 19 20 21 22

The phenomenon of microcracks closure is generally observed for the loading paths alternating tension and

23

xij + =

compression phases. Indeed, if compressive load succeeds to tensile load during which some microcracks have been created, these latter are progressively closed as the compressive load increases. The microcracks can be fully closed when the amount of the compressive load is enough high. Accordingly, the damage becomes inactive in compression while it was active in tension. This damage deactivation induces a partial or complete recovery of the material mechanical properties (elastic modulus, strain hardening modules, etc ...) and reduces or even completely cancels the damage evolution under compression. This activation (under tension) and deactivation (under compression) of damage has been largely studied in the literature for various types of materials (Marigo, 1985; Ortiz, 1985; Mazars, 1986; Simo et al., 1987; Ju, 1989, Mazars et al., 1990; Chaboche, 1992; Lemaitre, 1992; Chaboche, 1993; Ladeveze, 1993. Chaboche, 1995, Desmorat, 2000; Lemaitre and Desmorat, 2005; Desmorat, 2007; Lemaitre, 2009; Voyiadjis, 2009 ; among many others). The basic idea for the modeling of this phenomenon consists in introducing the damage effect (coupling) differently into the positive and negative parts of each tensorial state variable. Recall that the decomposition into positive and negative parts of any symmetric second-rank tensor x is performed through a spectral decomposition of the form r≤ n

∑x

k

Eijk and xij − = xij − xij + with xk are the r distinct eigenvalues of the tensor xij and Eijk = eik ekj

k =1

24

 (without summation on k) are the eigenprojection tensors in which e k are the eigenvectors associated to the

25

eigenvalues xk . Such kind of modeling poses no particular difficulties in the case of an isotropic damage with

26 27

isotropic elasticity (see for example Ladeveze, 1993; Lemaitre and Desmorat, 2005; Lemaitre et al, 2009; Issa et al, 2010; Saanouni 2012). However, with anisotropic damage, serious difficulties arise because of the continuity

17

1 2

and convexity loss of the state and dissipation potentials due to this spectral decomposition of the state variables

3

To avoid these difficulties, in this work, we will follow the same method but applied only to the damage energy

4

density release rate Yij . Indeed, we wish to consider the microcracks closure effect only on the damage evolution

5 6 7

without trying to recover the undamaged material properties under negative loading paths. This preserves the

8

e α changing only the components Yij and Yij of the damage energy density release rate Yij (Eq. (19) , Eq. (22)

9

e r α and Eq. (23)). In other words, we define a new energy density release rate noted Yij * = Yij * +Yij * +Yij with

10

respect to which the damage rate is derived. We keep the same damage potential defined by Eq. (29) in which

11

the norm

12

“microcracks closure” has the following form in place of Eq. (34):

(Ju, 1989; Chaboche, 1992, 1993; Ladeveze, 1993; Desmorat, 1999).

continuity and convexity of the state and dissipation potentials while giving a different damage rate under positive and negative loading paths. To this end, we keep unchanged the state relations Eq. (13) to Eq. (15) by

Yij

3

is simply replaced by Yij * . Accordingly, the damage evolution equation accounting for the 3

d
ij =

13

Ymn * − Y0

λ


(

1 − d rs

S

3

1

)

β

S

3Yij *

Yuv *

(41)

3

14

The new form of Yij * is derived from an appropriate form of the state potential accounting for the spectral

15

decomposition of the effective small elastic strain tensor εije and the effective kinematic strain tensor αij

16

defined by:

17

D eD εije = M ijkl ε kl + +

1 1 e D eD (1 − d rs )η ε nn .δ ij + M ijkl (1 − d pq )η −ε tte .δ ij − ε kl − − 1 1 3 3

(42)

18

D D αij = M ijkl α kl + + M ijkl −α kl −

19

with ε ijeD+ and ε ijeD− are the positive and negative parts of the deviatoric elastic strain tensor classically defined

20

e δ ij = ε ijeD+ + ε ijeD− . The damage-effect operator by ε ijeD = ε ije − 13 ε kk

21

deviatoric elastic and kinematic hardening tensors and it has exactly the same form as Eq. (11a) with the

22

effective damage operator hij − = δ ij − d ij

23 24

related to the microcracks closure effect and allows reducing the effect of damage on the negative parts of the

25 26

leading to a lower damage growth under compressive loading, with two extreme cases: (i) = 1 makes no

27 28 29

effect; and (ii) = 0 for which the damage growth take place only under the positive loads while no crack

30

the state potential from which derives the new damage energy density release rate Yij * is defined by:

(

1/ 4

)

(43)

D M ijkl is applied to the negative part of the

(

instead of hij = δ ij − d ij

1/ 4

)

. The parameter ( 0 ≤ ≤ 1 ) is

state variables. For metallic materials, the value of this parameter is usually taken = 0.2 (Lemaitre, 1992)

difference between the damage growth under positive and negative loading paths i.e. no microcracks closure

growth under the compressive phase of the loading paths. In this latter case, the microcracks close as soon as the local loading path is compressive termed here as a full microcrack closure effect. Accordingly, the new form of

18

1 2

2

2 2 eD 2 ρψ * = µε rkeD+ hkm ε mleD+ hlr2 + µε rkeD− hkm ke (1 − 13 d nn )η ε rre − ε ml − hlr − +

1 1 1 2 α + ml hlr2 + Cα − rk h−2kmα − ml h−2lr + (1 − ( 13 d nn2 )γ / 2 )Qr 2 Cα + rk hkm 3 3 2

1 2 3

1 + ke (1 − 13 d nn )η −ε rre 2

2

(44)

from which the new form of the elastic and kinematic contributions to the damage energy density release rate is easily obtained:

Yije * = − ρ

4

∂h 2 1 ∂ψ e * 2 ε mleD+ lr + keη (1 − 13 d nn )η −1 ε rre = Yije + * +Yije− * = −2 µeε rkeD+ hkm ∂dij ∂d ij 6

2 eD − 2 µeε rkeD− hkm −ε ml −

Yijα * = − ρ

5

∂hlr2 − 1 + keη (1 − 13 d nn )η −1 −ε rre ∂dij 6

2

2

δ ij (45)

δ ij

∂ψ pα * 2 ∂h 2 2 ∂hlr2 − 2 2 = Yijα + * +Yijα − * = − Cα rk + hkm α ml + lr − Cα rk − hkm α − ml − ∂d ij ∂d ij 3 ∂dij 3

(46)

6

where Yije + * and Yije − * (respectively Yijα + * and Yijα − * ) are the elastic (respectively kinematic hardening)

7 8

contributions depending respectively on positive and negative parts of the elastic strain tensor (respectively

9

Depending on the selected choice, the expressions of Eq. (45) and Eq. (46) give rise to various particular cases:

kinematic hardening strain tensor).

10 11

•

12 13 14

•

recovered leading to the damage growing identically in tension and in compression,

kinematic hardening are: 2  e 1 e+ eD 2 eD ∂hlr Y * Y * 2 h = = µ ε ε + keη (1 − 13 d nn )η −1 ε rre ij e rk + km ml +  ij ∂dij 6   2 Y α * = Y α + * = 2 Cα h 2 α ∂hlr ij rk + km ml +  ij 3 ∂dij 

•

2

δ ij

(a ) (b) (47)

Isotropic damage without microcracks closure ( = 1 ): the damage evolves identically in tension and in compression driven by:

 e 1 1 e eD eD η −1 e 2  1 Yij * = Yij = 3  µeε kl ε kl + 6 keη (1 − 3 d nn ) ε rr  δ ij     Y α * = Y α = 1  1 Cα α  δ ij kl kl  ij  ij 3  3 

18

19 20

Anisotropic damage with full microcracks closure effect ( = 0 ): in this case the damage grows only under positive load and the associated damage energy density release rate related to elasticity and

15 16 17

Anisotropic damage without microcracks closure ( = 1 ): the damage evolution equations Eq. (34) is

•

(a) (48)

(b )

Isotropic damage with full microcracks closure effect ( = 0 ): the damage evolves only in tension driven by the positive parts of the damage energy density release rate:

19

1

 e 1 1 e+ eD eD η −1 e 1 Yij * = Yij * = 3  µeε kl +ε kl + + 6 keη (1 − 3 d nn ) ε rr    Y α * = Y α + * = 1  1 Cα α  δ ij kl + kl +  ij  ij 3  3 

2

 δ ij 

(a) (49)

(b )

2

3. Parametric investigation of the proposed fully coupled model

3 4 5 6 7

This section is devoted to a relatively detailed parametric study aiming to show the capability of the proposed

8 9

( = 0 ) and without ( = 1 )microcracks closure effect, as well as the interactions between the plastic and

10 11 12

First, attention is focused on the effect of the various elastic and hardening components of the damage energy

13

Mises equivalent stress: ϖ = σ kk / 3 3σ ijDσ ijD / 2  , on the damage and plastic strain evolution is investigated.  

14

Particularly, the effects of the microcrack closure parameter and the hydrostatic coupling parameter η are

15

investigated.

16 17

Finally, the effects of the plastic anisotropy on the damage evolution and vice-versa are studied. Only the Hill

18

Table 1: Elastic and plastic parameters of the generic material.

model to describe damage growth in anisotropic materials under finite plastic strains with different loading paths in a typical RVE. The effects of various parameters of the model on the anisotropic damage growth and on the plastic flow with hardening are investigated. The elastic and hardening parameters, of the generic or virtual plastically isotropic material given in Table 1, are kept unchanged. The damage induced anisotropy with full

damage anisotropies in that cases, are investigated using various damage parameters shown in Table 2.

density release rate, as well as the effect of the coupling parameters governing them for the generic plastically isotropic material. Also the effect of the stress triaxiality, defined as the ratio of the hydrostatic stress to the von

anisotropy parameters are introduced assuming the non-associative normality rule for the sake of shortness.

Parameters

σy

E

ν

Q

b

C

a

Unit

MPa

GPa

-

MPa

-

MPa

-

300.0

210.0

0.3

2600.0

1.50

7000.0

170.0

19 20

  The plane loading paths imposed to the RVE, based on plane kinematics defined in plane (e1 , e2 ) , are:(i)

21 22

Uniaxial tension (UT), (ii) Plane tension (PT) and (iii) Equibiaxial tension (EB). In this plane case, the whole

23

mechanical fields represented by symmetric second-rank tensors have four non-zero components, including the    damage tensor which has the following form, with respect to the orthogonal observation triade (e1 , e2 , e3 ) :

20

 d11  d ij  =  d12   0

1

d12 d 22 0

 n11 n21 0   d I 0 0 0   n11 n12 0      0  0   n12 n22 0  =  n12 n22 0   0 d II d 33  ( e ,e ,e )  0 0 1   0 0 d III   0 0 1  1 2 3         dkl∆   

[ Pki ] 1/ 2

(50)

 Plj 

2

where d I = (d11 + d 22 ) / 2 + ( ( d11 − d 22 ) 2 − 4d122 )

3

  plane (e1 , e2 ) ,

4

with

5

orthogonal rotation tensor P defining the base change is represented by the following matrix:

/ 2 ,is the maximum eigenvalue of the damage tensor in the 1/ 2

d II = (d11 + d 22 ) / 2 − ( (d11 − d 22 )2 − 4 d122 )

/ 2 is the minimum eigenvalue and d 33 = d3 = d III ,

d1 , d2 and d3 are the eigenvalues of the damage tensor defined in Eq. (5). Here, the second-rank

 cos φ Pij =  − sin φ  0

6

sin φ cos φ 0

0 d12 0  with φ = arctg ( − ) d 22 − d I 1 

(51)

7

The angle φ , governed by the state and history of the applied loading, is the gap between the loading direction

8

    e1 and the maximum damage direction n1 (in the plane (e1 , e2 ) ).

9

3.1. Effect of the damage anisotropy on the initially isotropic plasticity

10

c p Let’s consider an RVE made from an initially isotropic plastic material (i.e. H ijkl = H ijkl = δ ik δ jl − 13 δ ijδ kl with

11

F i = Gi = H i = 0.5 and Li = M i = N i = 1.5 ) with mixed isotropic and kinematic hardening defined by the

12 13

parameters given in Table 1. The ductile damage is supposed either isotropic with fixed material parameters or

Stress J2(σ) [100 MPa]

anisotropic with four different sets of the material parameters as defined in Table 2. 16 14 12 10

isotropic damage γ=10 Anisotropic damage γ=2 γ=3 γ=6 γ=10

8 6 4 2 0 0,0

0,2

0,4

0,6

0,8

1,0

Equivalent plastic strain p

14 15 16

Figure 3. Uniaxial tension results in terms of stress versus equivalent plastic strain with isotropic plastically

17

Table 2: Damage parameters of the generic material given for the isotropic and anisotropic damage cases.

generic material.

Parameters

S

s

Y0

β

γ

η



Dc

21

Isotropic damage

Anisotropic damage

6.8

1.3

0.0

1.0

10.0

1.0

0.0 / 1.0

0.999

60.0

1.3

0.0

2.0

2.0

1.0

0.0 / 1.0

0.999

40.0

0.7

0.0

2.0

3.0

1.0

0.0 / 1.0

0.999

10.5

1.3

0.0

1.0

6.0

1.0

0.0 / 1.0

0.999

10.5

1.3

0.0

1.0

10.0

1.0

0.0 / 1.0

0.999

1 2

For the isotropic damage case, only γ=10.0 is considered; while for the anisotropic damage case the four

3 4 5 6

different values of γ ={2.0, 3.0, 6.0 , 10.0} are considered. These parameters are chosen in such a manner that all the cases give the same equivalent plastic strain at the final failure (or material ductility) for the UT loading path as shown in Figure 3. In Table 3 are presented the imposed total strain rate tensor as well as the resulting

8

damage tensor with and without microcracks closure effects for the three loading paths under concern. For the  uniaxial tension the maximum damage component is driven by the loading direction e1 , while two identical   eigenvalues are defined in the orthogonal plane (e2 , e3 ) normal to that loading direction. For =1, this double

9

eigenvalue is equal to d II = 14 d I and for =0 is equal to zero. Note that for this last case ( =0) the damage

7

10 11

  don’t evolve in the negative strain directions (i.e. plane (e2 , e3 ) ).

Table 3: The loading paths considered and their corresponding damage tensor with plastically isotropic material

Uniaxial Tension (UT)

Plane Tension (PT) Equi-Biaxial Tension (EB)

d ij

Dij

Loading

2 0 0  0 −1 0  2 0 0 −1

ε


0 d I 0 0 d  0 II    0 0 d II  with d I > d II = 14 d I

1 0

0 0   0 0 −1

dI 0   0

ε
0 0

1 0 0  ε
 0 1 0   0 0 −2

dI 0   0

0 dI 0

0 0 0 0  0 d I  0  0  d III 

with d III > d I = 14 d III

Scheme

=0

=1

dI 0   0

0 0 0 0  0 0 

dI 0   0

0 0 0 0  0 0 

dI 0   0

0 dI 0

0 0  0 

12 13

  For the PT loading path with =1 a maximum damage is observed in the plane (e1 , e2 ) , while for =0 exactly

14

the same damage state as the UT path is obtained with the maximum damage along the positive strain direction  e1 . Finally, for the EB loading path a major difference is observed between the two cases =1 (without

15

22

1 3

microcracks closure) and =0 (full microcracks closure). For =1 the maximum damage component is driven    by the normal direction e3 to the plane of the kinematics (e1 , e2 ) . For the case =0 the maximum damage   component is in the plane of the kinematics (e1 , e2 ) .

4

3.1.1. Effects of the damage energy density release rate

5 6 7 8 9

In Figure 4 are summarized, for the uniaxial tension (UT) path, the evolutions of the maximum components of

2

10

the damage energies density release rate: elastic (Figure 4a), kinematic hardening (Figure 4b), isotropic hardening (Figure 4c) and total (Figure 4D); performed using the different sets of material parameters defined in Table 2. In this figure, with different scales for the ordinates axis due to the different scales of these components, are compared the whole maximum components of the damage energy density release rate in the maximum  damage direction e1 (i.e. YIe , YIα and YIr ) and their summation (i.e. YI = YIe + YIα + YIr ) for both isotropic

11 12

damage and anisotropic damage. As expected, we observe that for the whole components the energies predicted

13

On the other hand, for the anisotropic damage, we observe that YIe (Figure 4a) and YIα (Figure 4b) are much

14

less sensitive to the values of the parameter γ than YIr (Figure 4c) since this parameter affects only the damage

15

effect on the isotropic hardening according to Eq. (9). In fact YIr decreases when γ increases as shown in Figure

16

4c influencing by the way the variation of the total energy YI as clearly shown in Figure 4d. In conclusion, lower

17

values of the parameter γ allow an important contribution of the isotropic hardening component of the damage

18

energy YIr since the early stage of the plastic loading and higher values of γ delays this higher contribution to

19

the final fracture stage. For example for γ =10.0, the damage evolution is governed from the beginning of the

20

plastic yield by only the elastic YIe and kinematic hardening YIα damage energies while the isotropic hardening

21

damage energy YIr contributes significantly when the final fracture is imminent. 8

0,5

Isotropic damage case with γ=10

7 6

Kinematic hardening α Energy density release Y I

Elastic Energy density release YI

e

by the anisotropic damage are higher than those predicted by the isotropic damage.

Anisotropic damage case

5 4

γ=2 γ=3 γ=6 γ=10

3 2 1 0 0,0

0,2

0,4

0,6

Equivalent plastic strain p

(a) Elastic energy density release

0,8

0,4

Anisotropic damage case α YI

0,3

γ=2 γ=3 γ=6 γ=10

0,2

Isotropic damage case with γ=10 α Y

0,1 0,0 0,0

0,2

0,4

0,6

0,8

Equivalent plastic strain p

(b) Kinematic hardening energy density release

23

Total Energy density release YI

Isotropic hardening r Energy density release YI

2,0 Isotropic damage case with γ=10 r Y

1,8 1,6

Anisotropic damage case r YI

1,4 1,2

γ=2 γ=3 γ=6 γ=10

1,0 0,8 0,6 0,4 0,2 0,0 0,0

0,2

0,4

0,6

0,8

8 Anisotropic damage case YI

7 6

γ=2 γ=3 γ=6 γ=10

5 4

Isotropic damage case with γ=10

3 2 1 0 0,0

0,2

Equivalent plastic strain p

(c) Isotropic hardening energy density release

1

0,4

0,6

0,8

Equivalent plastic strain p

(d) Total energy density release

Figure 4. Uniaxial tension results in terms of damage energies versus the equivalent plastic strain.

0,5

0,5

0,4

0,4

0,4

0,3 0,2 0,1

Damage dIII

0,5

Damage dII

Damage dI

2

0,3 0,2

1,0

1,2

1,4 1,6 1,8

Equivalente plastic strain p

(a) d I .

0,0 0,0

0,2 0,1

0,1

0,0 0,0 0,2 0,4 0,6 0,8

0,3

0,2 0,4

0,6

0,8

1,0 1,2

1,4

1,6

1,8

Equivalente plastic strain p

(b) d II .

0,0 0,0 0,2 0,4 0,6 0,8 1,0 1,2

1,4 1,6 1,8

Equivalente plastic strain p

(c) d III .

3

Figure 5. Evolution of the damage tensor eigenvalues for the different loading paths.

4 5

In Figure 5 are shown the evolution of the damage eigenvalues versus the accumulated plastic strain for the three

6

and without the microcracks closure effect. As expected from Table 3, for UT only d I is non-zero while

7

d II = d III =0 when the full microcrack closure ( =0) is accounted for; however when the microcrack closure is

8

disabled ( =1) we have d II = d III d I / 4 . For the PT the main difference with the UT is that d II =0 and

9

d III = d I when =1 with a lower damage rate (Figure 5a and c). Finally, for the EB-T a different situation is

10

applied loading paths defined in Table 3 for both the isotropic damage as well as the anisotropic damage with

observed with d II = d I = d III / 4 when =1 while d II = d I and d III = 0 for =0.

24

1,5 1,0

 e2 '

π − plane ∆l

ε II

D

0,5

π θ=

6

0,0

 e1 '

max damage dI=0.1

-0,5 -1,0 -1,5 -1,5

 e3 '

isotropic damage Anisotropic damage γ=2 γ=3 γ=6 γ=7 γ=8 γ=9 γ=10

max damage d I=0.9

-1,0

-0,5

0,0

εI

1

0,5

1,0

1,5

D

2

Figure 6. Isotropic and anisotropic isodamage curves plotted in the deviatoric strain space (π-plane) for

3

d I = 0.1 , d I = 0.9 and ϖ = 0

4

In Figure 6 are displayed the isodamage curves (or surfaces) in the diviatoric (π-plane) of the strain space for two

5

constant maximum damage values d I = 0.1 (moderate damage value) and d I = 0.9 (high damage value near the

6 7 8 9 10

final fracture of the RVE) while the stress triaxiality is zero ( ϖ = 0 ). As expected, the isotropic damage gives

11

damage level d I = 0.1 , the isodamage surfaces have the same polygonal shape independently from the values of

12

the parameter γ. However, for the higher damage level d I = 0.9 (i.e. near the final fracture of the RVE), the

13

shape of the isodamage surfaces degenerates from a circular shape for lower values of γ (γ=2.0 ) to polygonal

14

shape with rounded edges for higher values of γ (γ=10.0) . This change in the isodamage surfaces shape is due to

15

the contribution of the isotropic component YIr of the damage energy density release rate as discussed above

16

from Figure 4. In fact, for higher level of the maximum damage, the lower values of γ generate, as was

17

mentioned above, higher values of YIr more important than YIe , inducing by this ways the circular isodamage

18

surfaces. For higher values of γ the elastic damage energy YIe is much lower than YIr inducing the hexagonal

19

shape of the isodamage surfaces. Recall that the elastic damage energy YIe represents the loading path

20

directionality effect on both the damage state and evolution.

perfectly circular isodamage curves for the tow maximum damage values, due to the fact that isotropic damage is independent from the loading directions.

We observe also that the isodamage surface size increase with

increasing the damage level. For the anisotropic damage, a significant change on the isodamage surfaces is observed depending on the maximum damage level and on the value of the parameter γ. For the lower maximum

25

2,0 1,5

 e2 '

h=0 h=0.2 h=0.4 h=0.6 h=0.8 h=1

1,0

ε II

D

0,5 0,0

θ = π3

2π 3

-0,5

 e1 '

-1,0 -1,5 -2,0 -2,0

 e3 ' -1,5

π − plane -1,0

-0,5

0,0

0,5

1,0

1,5

2,0

D

εI

(a) Isotropic damage 2,5 2,0

h=0 h=0.2 h=0.4 h=0.6 h=0.8 h=1

 e2 '

1,5 1,0

θ = π3

εII

D

0,5

2π 3

0,0 -0,5

 e1 '

-1,0 -1,5 -2,0

 e3 '

π − plane

-2,5 -2,5 -2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 2,5 D

εI

(b) Anisotropic damage

1 2

Figure 7. Effect of the microcracks closure parameter on the isodamage surfaces for isotropic and anisotropic

3 4

Otherwise, we remark that the circular shape of the isotropic damage envelop is inscribed on the anisotropic

5

angle θ = 13 cos −1 ( 272 det(σ ijD ) J 23 (σ ijD ) ) . From Figure 6, it is worth noting that the maximum strain values are

6 7

situated on the hexagonal shape corners, which are defined by a Lode angle θ = π / 6 . This Lode angle

8

uij = ωeik ekj − ωeil elj (without summation with i ≠ j = {1, 2, 3} and ω = ±1 ). For this loading path, the difference

9

between the isotropic and the anisotropic damage cases as can be seen in Figure 6 in terms of the radius ratio

damage plotted in the diviatoric strain space (π-plane)

damage hexagonal envelop. This last shows a significant damage evolution change with respect to the Lode

corresponds to the PT loading paths (i.e. pure shear loadings) with fixed direction on the strain space

10

' ∆l / ε iso = 2 − 1 is induced by the expression of the energy norm Yij

11

potential (see Eq. (29)). This quantity is fixed for a given material.

12 13 14

When considering the microcracks closure effect, the isodamage surfaces are summarized in Figure 7 for six

3

= 3 Yij Yij governing the damage

different values of the parameter varying between 0.0 (full microcracks closure effect) and 1.0 (no microcracks closure effect). For the isotropic damage, the isodamage surfaces degenerate from the perfectly

26

1 2

initial circular shape in the absence of the microcracks closure effect (i.e., = 1 ) to a three corners star shaped

3 4 5 6 7

For the anisotropic damage the isodamage surfaces degenerate from the hexagonal shape for = 1 (without

8

direction uij = eik ekj + eil elj − 2eim emj (without summation and with i ≠ j ≠ k = {1, 2, 3} ).

9

3.1.2. Effect of the stress triaxiality

surfaces when the full microcrack closure effect is accounted for (i.e., = 0 ) as is clearly seen in Figure 7.a.

microcracks closure effect) to a three corners star shaped surfaces for = 0 (with full microcracks closure effect) as shown in Figure 7.b. Note that, the corners of the star shaped surface are sharper for anisotropic damage (Figure 7.b) than for isotropic damage (Figure 7.a). Also, for both isotropic and anisotropic damage, the direction of the star corners are given for a Lode angle θ = π / 3 corresponding to (EB) loading path defined by the

10

In this section the effects of the microcrack closure parameter and the parameter η affecting the hydrostatic

11

H part of the damage-effect tensor M ijkl (see Eq. (11b)), on the damage growth with various values of the stress

12

triaxiality ( ϖ )using the same plastically isotropic generic material defined in Table 1 and Table 2.

13

Let’s begin with examining the effect of the microcracks closure parameter on the evolution of the equivalent

14

plastic strain at fracture of the RVE (or the material ductility p fr ) for the two extreme cases: =1 and =0.

15

Figure 8 gives the evolution of the material ductility p fr with the stress triaxiality ϖ based on the plane stress

16

assumption for isotropic and anisotropic damage using the two extreme cases =1 (without microcracks closure

17

effect) and =0 (with full microcracks closure effect) for the maximum damage value d I = 0.9 . From this

18

Figure we can see the significant difference between the results issued from isotropic and anisotropic damage.

19

For the isotropic damage, p fr increases monotonically when ϖ decreases from ϖ =2/3 (EB loading path) to

20

ϖ = 0 (PC loading path) and lower when the microcracks closure is fully accounted for (i.e., =0). However,

21

when the microcracks closure is not accounted for (i.e., =1), the maximum of plastic strain at fracture

22

p fr bifurcates from the one given by =0 at ϖ = 1 / 3 (UT path) reaches a maximum value of p fr = 0.82 for

23

ϖ = 0 (PC path) and decreases for the negative values of the stress triaxiality to reach p fr = 0.78 for ϖ = −0.3 .

24

Note that, for this isotropic damage case, the p fr versus ϖ curves are the same for ϖ ≥ 1 / 3 .

25

For the anisotropic damage the equivalent plastic strain at fracture of the RVE evolves differently compared to

26

the above discussed isotropic damage. In this anisotropic damage case, the p fr versus ϖ curves are the same

27

for =1 and =0 for ϖ lying between 0.1 and 0.55. While for =0 (full microcracks closure effect)

28

p fr increase for ϖ ≥ 5.5 and reach p fr = 1.25 for ϖ = 2 / 3 ; and also increase for ϖ < 0.1 to reach a maximum

29

of p fr = 1.62 for ϖ = −0.3 . Finally, it is noted that whole curves intersect for the positive value of the stress

30

triaxiality around ϖ = +0.3 .

27

Equivalent Plastic strain at rupture ρ fr

1,8

1 3

− 13

2 3

3 3

Anisotropic damage h=1 h=0 Isotropic damage h=1 h=0

0

1,6 1,4 1,2 1,0 0,8 0,6 0,4 -0,2

0,0

0,2

Triaxiality

1

0,4

0,6

0,8

ϖ

2 3

Figure 8. Effect of the microcracks closure parameter on the plane stress damage envelops defined for a

4

H Now the effect of the parameter η affecting the hydrostatic part of the damage-effect tensor M ijkl (see Eq.

5

(11b)) on the material ductility p fr for the single case when the full microcracks closure effect ( =0) is

6

considered only for the anisotropic damage with d I = 0.9 .

maximum damage level d I = 0.9 .

Equivalent Plastic strain at rupture ρ fr

1,8

− 13

1,6

1 3

2 3

3 3

η=0.2 η=0.6 η=0.8 η=1.0 η=1.2 η=1.5 η=2.0 η=3.0

0

1,4 1,2 1,0 0,8 0,6 0,4 -0,2

7

0,0

0,2

Triaxiality

0,4

0,6

0,8

ϖ

8

Figure 9. Effect of the hydrostatic damage coupling parameter η on the material ductility versus stress triaxiality

9

for anisotropic damage with =0 for d I = 0.9

10

Figure 9 summarizes the p fr versus ϖ curves for eight different values of η varying between η =0.2 and

11

η =3.0 for the maximum damage value d I = 0.9 , obtained with anisotropic damage and = 0 (full microcracks

12

closure effect). These results show that p fr is still insensitive to the parameter η for ϖ ≤ 0.1 while they are

28

1

highly sensitive to the values of η for ϖ > 0.1 . In fact, p fr starts from p fr = 1.7 for ϖ = −0.3 decreases

2

monotonically to reach p fr = 0.57 for ϖ = 2 / 3 with η = 3.0 . However, for the whole other values of η , p fr

3

increases for ϖ > 1 / 3 to reach different levels depending on the value of η and this maximum is higher for

4

lower values of η as shown in Figure 9.

(a) =1.

(b) =0.

5

Figure 10. Effect of the microcracks closure parameter on the 3D damage surfaces defined for a maximum

6

damage level d I = 0.9 (near the final rupture).

7 8

A 3D representation of the isodamage surfaces are plotted in the 3D triade defined by the deviatoric strain plane

9 10 11

( =1) and with full microcracks closure ( =0) as shown in Figure 10. This 3D representation shows the

(π-plane) and the stress triaxiality as a third normal axis in the two cases: without microcracks closure effect

variation, with respect to the triaxiality, of the isodamage curves already displayed in Figure 7 for ϖ = 0 . This figure shows that the increasing of the triaxiality decreases the size of the isodamage surfaces in the deviatoric

29

1 2 3 4

strain π-plane and changes the shape of the isodamage surfaces from polygonal shape to perfectly circular shape

5 6 7 8 9

Note that the varying shape of the iso-damage surfaces seen in figures 6, 7 and 10 show a clear sensitivity to the

for high triaxiality level. This change is due to the fact that the increasing of the triaxiality favors the hydrostatic damage energy (see Eq. (19) and Eq. (22) to (24)). Recall that this part of the elastic damage energy leads to the isotropic damage growth which is at the origin of the circular isodamage surfaces for high triaxiality.

Lode angle. This sensitivity is naturally accounted for, in this model, due to the strong coupling between the elastoplastic behavior and the anisotropic ductile damage, without explicitly introducing the Lode angle in the expression of the damage potential (Eq. 29) as made in the fracture surface proposed in (Wierzbicki et al., 2005; Teng, and Wierzbicki, 2006).

10

3.2. Effect of the initial plastic anisotropy on the anisotropic damage

11

In this section we consider plastically anisotropic material characterized by the following anisotropic material

12

parameters considering associative plasticity for the sake of simplicity: F p = F c = 0.2 , G p = Gc = 0.4 ,

13 14 15

H p = H c = 0.6 , N p = N c = 1.0 . These anisotropic parameters allow a strong plastic anisotropy as can be seen in the Lankford ratio evolution curves given in Figure 11, Figure 12 and Figure 13. The remaining material

16

plane of the kinematics, are defined by the orientation ϕ 0 (see Table 4). Recall that the Lankford ratios have the

17

following forms for the quadratic Hill48 equivalent stresses:

18

parameters are given in Table 1 and 2 with γ = 10.0 . The initial directions of the material anisotropy, in the

rϕ

( 4H =

p

+ G p + F p − 2 N p ) c os 4 ϕ0 + ( 2 N p − 4H p − F p − G p ) c os 2 ϕ0 + H p

(G

p

− F p ) c os2 ϕ0 + F p

(43)

19

It should be mentioned that Eq. (43) is given for the initial undamaged configuration. So, when damage starts

20 21

and evolves, this relation will be function of the damaged parameters ( F p , G p , H p and N p ). This highlights

22 23

For the UT loading path, as shown in Table 4, the damage tensor has the same loading eigenvectors as for the

the additional plastic anisotropy induced by the anisotropic character of the ductile damage.

24

isotropic plasticity except when the initial orientation of the anisotropic direction ϕ 0 is not equal to 0° and 90°.  For these two orientations the maximum value of damage d I is driven by the loading direction e1 (i.e. φ = 0° ).

25

For the case ϕ 0 = 45° and in general when ϕ 0 ≠ (0°,90°) the damage tensor has not the same directions as the

26 27

loading but the angle φ remains small ( φ = 0.67° ) so the maximum damage d I is mostly driven by the loading  direction e1 . For = 1 the transverse damage components d II = d 22 and d III = d 33 differ strongly and their ratio

28

depends on plastic anisotropy and the orientation ϕ 0 (see Table 5). However, for = 0 the transverse damage

29

components remain zero ( d II = d III = 0 ) as for the plastically isotropic material.

30 31 32

30

1

Table 4: The loading paths considered and their corresponding damage tensor with anisotropic plasticity

ϕ0 = 0° and ϕ0 = 90°

Uniaxial Tension (UT)

d ij

Dij

Loading

2 0 0  0 −1 0  2 0 0 −1

ε


ϕ0 = 45°

0 d11 0 0 d 0  22   0 0 d 33  with d11 > d33 > d22

dI 0   0

 d11 d12 0  d d 22 0   12   0 0 d33  with φ = 0.674°

 d11 d  12  0

1 0 0  ε
0 0 0   0 0 −1

1 0 0  ε
 0 1 0   0 0 −2

Equi-Biaxial Tension (EB)

d11 0   0

0 d22 0

0 0 0 0  0 0  d12 d 22

0

0 0  0 

with φ = 0.673°

d I > d33 ≈ d II

Plane Tension (PT)

Scheme

=0

=1

0 0  d 33 

 d11 0   0

0 d 22

0

0 0  0 

with d33 > d11 > d22

with d11 >> d 22

d11 0   0

 d11 0   0

0 d22 0

0 0  d 33 

with d33 > d11 > d22

0 d 22

0

0 0 0

with d11 > d 22

2

5 6

depends on the initial plastic anisotropy, except for = 0 for which the third normal damage component remains identically zero as for UT loading path (see Table 3).

Damage dI

7 8

physically wrong behavior. However, for = 0 the maximum damage remains in the plane of the kinematics  driven by the direction e1 . The transverse damage components are strongly different for all cases and their ratio

0,5

0,25

0,5

0,4

0,20

0,4

0,3 0,2 0,1 0,0 0,0

0,15 0,10 0,05

0,2

0,4

0,6

0,8

1,0

Equivalente plastic strain p

(a) d I .

9

Damage dIII

4

For the PT and EB loading paths, we observe a fundamental difference between the two cases = 0 and = 1 .  For = 1 , the maximum damage is driven by the direction e3 normal to the plane of kinematics which is

Damage dII

3

1,2

0,00 0,0

0,3 0,2 0,1

0,2

0,4

0,6

0,8

1,0

Equivalente plastic strain p

(b) d II .

1,2

0,0 0,0

0,2

0,4

0,6

0,8

1,0

1,2

Equivalente plastic strain p

(c) d III .

Figure 11. Evolutions of the eigenvalues of the damage tensor for the different loading paths.

31

1 2 3 4

In Figure 11 are shown the evolutions of the three principal damage components versus the accumulated plastic

5 6

In Figure 12 is shown the effect of the initial plastic anisotropy on the isodamage surfaces plotted in the

7 8 9 10

reduced size compared to the isodamage surface for isotropic plasticity given in Figure 7.b, and evolves toward

11

characterized by the same initial yield stress σ y since

12

  surfaces coincide for this direction as shown by Figure 12. For the two other directions ( e2 and e3 ) the elastic

13

stress limits are respectively σ y / G c + F c

14 15 16 17 18

radiuses of the isodamage surfaces shrink along these loading directions as seen in Figure 12. For the anisotropic

strain for the anisotropic damage with and without the microcracks closure under the three applied loading paths. Compared to the isotropic plasticity given in Figure 5, the ratios of the non-zero damage components are dependent on the plastic anisotropy as shown in Table 4.

deviatoric strain π-plane for d I = 0.9 . For the isotropic damage, the isodamage surface becomes elliptical with

the three corners star shape. This shape variation is due to the elastic limit change governed by the anisotropic parameter. Recall that, an increasing of the elastic damage energy is proportional to the increasing of plastic  yield stress. Accordingly, for the loading direction e1 both isotropic and anisotropic plastic cases are

G c + H c = 1 leading to the fact that all the isodamage

and σ y / H c + F c

which are greater than

σy

. For this reason the

damage case, the same effects are observed i.e. the change of the isodamage surfaces and the decrease of their size. Furthermore, we can observe a slight rotation of the isodamage surfaces accompanying the shape variation. This rotation can be recognized through the change of the Lode angle associated with the three corners star shaped isodamage surfaces.

2,5 2,0

 e2 '

1,5 1,0

εII

D

0,5

2π 3

0,0 -0,5

 e1 '

Isotropic plasticity Anisotropic damage h=0 Anisotropic damage h=1 Isotropic damage h=1 Anisotropic plasticity Anisotropic damage h=0 Anisotropic damage h=1 Isotropic damage h=1

-1,0 -1,5 -2,0

 e3 '

π − plane

-2,5 -2,5 -2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 2,5 D

19 20 21 22 23 24 25

εI

Figure 12. Effect of plastic anisotropy on the isodamage surface using isotropic and anisotropic damage plotted in the diviatoric strain π-plane for d I = 0.9 When an RVE is pre-strained with the one of the plane loading paths given in Table 4, the evolving damage tensor allows the change of the plastic anisotropy. Hence, after the pre-straining and the creation of some microcracks in the RVE, the anisotropic plasticity is changed depending on the previous pre-straining directions and amplitude. Concretely, this leads to an important variation of the effective plastic anisotropy parameters

32

1 2 3

( F p , G p , H p and N p ) . In the following, we will examine the effect of an initial damage induced by a prestraining state on the evolution of the plastic anisotropy in terms of plastic flow and yielding limits for tow maximum damage levels d I = 0.1 and d I = 0.5 . 6

5

initial case of h=0 dI=0.1

4

dI=0.5

3

case of h=1 dI=0.1

Lankford ratio rψ

Lankford ratio rψ

6

Increasing damage

dI=0.5

2 1

10

20

30

40

50

60

70

80

dI=0.5

4

case of h=1 dI=0.1

3

dI=0.5

2 1 Increasing damage 0 0 10 20

0 0

initial case of h=0 dI=0.1

5

90

ψ0[°]

2 anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

Normalized Stress σ90/σy

Normalized Stress σ90/σy

8

dI=0.5

4

dI=0.5

-8 -8

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

50

60

70

80

90

(b) Lankford ratio for ϕ0 = 90°

initial surface isotropic damage case dI=0.1

6

40

ψ0[°]

(a) Lankford ratio for ϕ0 = 0° 8

30

6

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case with h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case with h=1 dI=0.1

-4 -6

dI=0.5

-8

8

-8

(c) Yield surface for ϕ0 = 0°

-6

-4 -2 0 2 4 Normalized Stress σ 0/σy

6

8

(c) Yield surface for ϕ0 = 90°

4 5

Figure 13. Damage induced change on Lankford ratio curves and yield surface due to UT path with two material

6 7

In Figure 13 are plotted the evolutions of the Lankford ratio of the pre-strained RVE with UT loading path. We

8

the maximum damage evolves along the direction of tension and consequently the Lankford ratio rϕ = D22p / D33p

9

increases strongly along the transverse direction which is ϕ0 = 90° (see Figure 13a). Indeed, after pre-damaging

10

in the ϕ0 = 0° direction, the resulting softening increases significantly the plastic flow in this direction. So, when

11

an imposed loading occurs in the transverse direction (i.e. ϕ0 = 90° ), the plastic strain rate component D22p shall

12

increases strongly compared to the component D33p . The increase is even larger as the level of the maximum

13

damage is important. The inverse situation is observed if the pre-straining direction is ϕ0 = 90° as shown in

14

Figure 13b. It is observed that the Lankford ratio increases significantly for the orientation ϕ0 = 0° which is not

orientations ϕ 0 .

observe in this figure that the ϕ 0 orientation induces an important change in the shape of the curves. For ϕ0 = 0° ,

33

1

the case for the orientation ϕ0 = 90° . Note that the accounting of the microcracks closure ( = 0 ), favors the

2

increasing of the plastic anisotropy for all cases.

3

The effect of the initial damage state due to UT loading (with material orientation ϕ0 = 0° and ϕ0 = 90° ) on the

4

yield surfaces plotted in the normalized stress space is showed in Figure 13c for ϕ0 = 0° and in Figure 13d for

5

ϕ0 = 90° . In these figures are compared the results of the isotropic and anisotropic damage through two

6

maximum damage levels d I = 0.1 and d I = 0.5 . It is observed that the change induced on the yield surface

7 8 9 10 11 12

increases with the pre-damaging created by the previous UT loading path. This change is related to expansion

13 14 15

microcracks closure effect) induces significant difference compared to = 1 (without microcracks closure

16

for = 1 ( d II = 0.08 and d III = 0.140 see Table 5). Also from Figure 13c and Figure 13d we observe that the

17

anisotropic yield surfaces have grater size than the isotropic damage ones.

18 19

For the PT loading path are given, in Figure 14, the Lankford ratio evolutions (Figure 14a) and the yield surface

20

loading we observe an important difference between the two cases = 0 and = 1 for anisotropic damage. For

21

= 0 , the curves of the Lankford ratio increase with the maximum damage level in the side of the transverse

22

direction (for ϕ0 = 90° ) similar to the UT loading path. However, for = 1 , the Lankford ratio decreases with

23

the maximum damage level in the side of the loading direction (ψ 0 = 0° ). This is evidently due to the different

24

damage state obtained for the two cases (see Table 5). Indeed, for = 0 the damage state is similar to the one

25 26

obtained with the UT path. However, for the case = 1 the maximum damage component is carried by the

27

plane of kinematics. Consequently, the plastic flow component in the normal direction D33p increases more than

28 29

the transverse component which explains the decreasing of the Lankford ratio. It should be mentioned that the

30

damage component d I in the direction of the pre-straining and the normal component d III approaches 1 (final

31

fracture) as we can see in Table 5.

32 33 34 35

Furthermore, we observe also for this TP loading path, an important effect of the damage anisotropy on the yield

and rotation of the elliptic yield surface compared to the isotropic damage case. Note that the size of the yield surface is mostly reduced for the isotropic damage case. This is due to the fact that the entire stress components are identically affected when considering isotropic damage. However, with the anisotropic damage the transverse damage components (see Table 5) are significantly smaller than the maximum component and consequently the softening in these directions is smaller than with the isotropic damage model. For this same reason, = 0 (full

effect). This difference appears particularly in the increasing of the yield stress in biaxial expansion due to the fact that the transverse damage components are identically zero for = 0 ( d II = d III = 0. ) so lower than those

(Figure 14b) for d I = 0.1 and d I = 0.5 , obtained with the different versions of the model. For this pre-straining

normal to the plane of kinematics. So the softening in the normal direction is more enhanced than those in the

values of the Lankford ratio don’t vary when ϕ 0 approaches 90°. This is due to the fact that the ratio of the

surface plotted in the stress space. As for the UT path, the effects induced by damage anisotropy consist on rotation and expansion of the yield surface compared to the isotropic one. For = 0 allows the increasing of the yield stress in the side of the biaxial tension.

34

1 2

Table 5: Damage eigenvalues obtained with the different loading paths and for a 0.5 maximum damage level Loading

ϕ 0 = 0° ϕ 0 = 45° ϕ0 = 90°

UT

PT EB

dI

= 1.0 d II

d III

dI

= 0.0 d II

d III

0.500

0.038

0.140

0.500

0.

0.

0.500

0.065

0.095

0.500

0.0002

0.

0.500 0.203 0.116

0.022 0.019 0,050

0,188 0.500 0.500

0.500 0.500 0.500

0. 0.083 0.217

0. 0. 0.

3 4 5

Finally, the results obtained with the EB loading paths are shown in Figure 15 in terms of the Lankford ratio

6

d I = 0.5 . From Figure 15a we find that the microcracks closure parameter influences the evolution of the

7 8 9 10 11 12 13

Lankford ratio. Indeed, for = 1 , the maximum damage component is carried by the direction normal to the

variation (Figure 15a) and the yield surfaces (Figure 15b) for the same maximum damage values d I = 0.1 and

plane of kinematics (see in Table 4 and Table 5). So the plastic flow in this direction increases highly than the components in the plane of kinematics. In addition, the maximum damage component is at least twice larger than the components in the plane of kinematics. So for these raisons the various curves of the Lankford ratio evolution decrease for all directions with the increasing of the maximum damage level. The inverse situation is observed for the case = 0 for which the Lankford ratio evolution curves increase with damage due to the fact that the softening in the plane of kinematics is more important than in the transverse direction.

5

initial case of h=0 dI=0.1

4

dI=0.5

3

case of h=1 dI=0.1

8

Increasing damage h=0

dI=0.5

2 1

Increasing damage h=1

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

dI=0.5

-8

0 0

10

20

30

40

50

60

70

80

90

ψ0[°]

(a) Damage induced change on the Lankford ratio evolution.

14 15

Normalized Stress σ90 /σy

Lankford ratio rψ

6

-8

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

6

8

(b) Damage induced change on yield surface displayed in the normalized stress space.

Figure 14. Effect of the induced damage anisotropy on the yield surface and the Lankford evolution after (PT) loading.

35

initial case of h=0 dI=0.1

5 4

8 Normalized Stress σ90/σy

Lankford ratio rψ

6

dI=0.5

Increasing damage h=0

case of h=1 dI=0.1

3

dI=0.5

2 1

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

Increasing damage h=1

0 0

10

20

dI=0.5

-8 -8

30

40

50

60

70

80

90

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

6

8

ψ0[°] (a) Damage induced change on the Lankford ratio evolution

(b) Damage induced change on the yield surface displayed in the normalized stress space

1 2

Figure 15. Effect of the induced damage anisotropy on the yield surface and the Lankford evolution after (EB)

3 4 5 6

Regarding the yield surface, we find also for this pre-straining with EB loading path, an important change

7

σ 90 > 0 ). With = 0 however, we obtain an increase of the yield surface in the side of biaxial tension domain

8

( σ 0 > 0 and σ 90 > 0 or σ 0 < 0 and σ 90 < 0 ).

loading.

induced by the damage anisotropy compared to the isotropic damage model. We observe an important increase in the size of the yield surface depending on the value of the microcracks closure parameter . For = 1 , we obtain an increase of the yield surface in the side of pure shear domain ( σ 0 > 0 and σ 90 < 0 or σ 0 < 0 and

9 10 11 12 13

In conclusion we can see that the anisotropy of damage induces an important change in plastic anisotropy in

14

4. Conclusions

15 16 17 18 19

A thermodynamically-consistent model accounting for large anisotropic plastic strains, mixed non-linear

20 21 22 23

The results obtained from a parametric study on the RVE with generic anisotropic non-associative plasticity,

terms of plastic flow and yield stress. The change on these surfaces is function of the history of the loading, i.e. the pre-straining loading paths. Moreover, the consideration of the microcracks closure effect (through the parameter ) allows significant changes on the evolution of the plastic anisotropy (induced anisotropy). When the value of this parameter approaches 0, the induced plastic anisotropy evolution is stronger.

hardening fully coupled with anisotropic ductile damage in which the microcracks closure effect on the damage growth is accounted for, is developed in framework of non-associative plasticity theory. The ductile damage anisotropy is described by a symmetric second-rank tensor and its effect on the mechanical fields is described by symmetric fourth-rank tensor.

enabled as to show the predictive capabilities of the proposed model in the description of the evolution of the plastic anisotropy induced by the anisotropic damage. Indeed, the comparison of the result of the anisotropic damage model with various values of the microcracks closure parameter revealed an important change in the

36

1 2 3 4 5

plastic anisotropy depending on the damage state created by pre-straining loading paths. The changes affect the

6 7 8 9 10 11

All these aspects are expected to have an important influence on the numerical simulation of metal forming

12

Acknowledgement.

13 14

Thanks are due to the Region Champagne – Ardenne (France) for the financial support.

plastic flow anisotropy, the isodamage surfaces as well as the yield surface. It is found that the iso-damage surfaces are sensitive to the microcracks closure (Figure 7, Figure 10) and the stress triaxiality (Figure 8, Figure 9 and Figure 10) where the shape of the isodamage surfaces change significantly with respect to Lode angle due to the strong coupling with the anisotropic damage.

processes where the initial and induced anisotropies play an important role. However, before any application to metal forming simulations, the material parameters entering the constitutive equation have to be identified by comparison with experimental results conducted until the final fracture of the specimens subject to various loading paths. This task is under progress based on available experimental works completed by specific experimental results as we will show in a forthcoming paper.

15

37

1

Appendix 1

2

Calculation of the damage effect tensor derivatives

3

In this Appendix is given in details of the calculation of the derivatives of the effective damage tensor

4

∂hij2 / ∂d kl . The damage tensor can be rewrite from the equation Eq.(5) in the following form: r ≤n

r≤n

k =1

k =1

dij = ∑ d k Eijk and hij = ∑ hk Eijk

5

d k are the r distinct eigenvalues of the damage tensor and Eijk are their associated eigenprojection

6

with

7

tensors.

8

associated eigenprojection tensors. The derivative of the effective damage tensor is then given by:

hk = (1 − d k )1/ 4 are the eigenvalues of the effective damage tensor h which has E k are also their ij ij

2 ∂hpq

9

(1.1)

∂d rs

r≤n r≤n

= ∑∑ k =1 l =1

r≤n r≤n r ≤n r ≤n ∂hk2 k l ∂ (1 − d k )1/ 2 k l 1 E pr Eqs = ∑∑ E pr Eqs = − ∑∑ δ Ek El 1/ 2 kl pr qs ∂d l ∂d l k =1 l =1 k =1 l =1 2(1 − d k )

(1.2)

10 11

Appendix 2

12

Calculation of the plastic multiplier

13

The constancy condition is given by:

f
(σ ij , X ij , R; dij ) = 0 ⇒

14 15

∂f ∂f
∂f
∂f
σ
ij + X ij + R+ dij = 0 ∂σ ij ∂X ij ∂R ∂dij

(2.1)

Using the state relation given by the equation Eq.(16-18) and the evolution equations Eq.(32-34), we obtain :

σ
ij =

16

∂σ ij

∂ε

e kl

ε
kle +

∂σ ij
∂σ ij ∂σ ij ∂F d kl = e ε
kle + λ
∂d kl ∂ε kl ∂d kl ∂Ykl

(2.2)

17

∂X ij ∂X ∂X ∂F ∂X ij ∂F X
ij = α
kl + ij d
kl = −λ
ij + λ
∂α kl ∂d kl ∂α kl ∂X kl ∂d kl ∂Ykl

(2.3)

18

∂R ∂R
∂R ∂F ∂R ∂F R
= r
+ dij = −λ
+ λ
∂r ∂dij ∂r ∂R ∂dij ∂Yij

(2.4)

19

Note that the elastic strain rate is defined by the additive decomposition (see Eq.(2))

20

ε
ije

21

Dij so we shall rewrite : 2 ε ik . Wkj 

S

S

= Dij − Dij − 2 ε ik . Wkj  in which the quantity 2 ε ik . Wkj  is a linear function of the total strain rate tensor p

e

e

e

S

= ϒijpq Dpq with ϒijkl = ε ime K mjkl − ε mje Kimkl is a symmetric fourth-

22

e p rank tensor. Hence, we thus arrive at ε
ij = Dij − Dij − ϒ ijpq D pq and

23

σ
ij =

24

(2.2) becomes:

∂σ ij ∂ε

e kl

Dkl −

∂σ ij ∂ε

e kl

ϒ klpq Dpq −

∂σ ij ∂ε

e kl

Dklp +

∂σ ij
d kl .Using the evolution equation Eq.( 31), the equation ∂d kl

38

∂σ ij

(δ

δ ql − ϒ klpq ) Dpq − λ


1

σ
ij =

2 3

Hence, the consistency condition Eq.(1.2) becomes :

4

∂ε

e kl

pk

 ∂f ∂σ ij ∂F ∂f ∂X ij ∂F ∂f ∂R ∂F + +  e ∂f ∂σ ij  ∂σ ij ∂ε kl ∂σ kl ∂X ij ∂α kl ∂X kl ∂R ∂r ∂R δ δ − ϒ klpq ) Dpq − λ
 f
= e ( pk ql ∂σ ij ∂ε kl  − ∂F  ∂f ∂σ kl + ∂f ∂X kl + ∂f ∂R + ∂f  ∂Yij  ∂σ kl ∂d ij ∂X kl ∂d ij ∂R ∂d ij ∂d ij   where f
= a pq D pq − λ
ℵ = 0

6

with a pq =

7

 ∂f ∂σ ij ∂F ∂f ∂X ij ∂F ∂f ∂R ∂F + +  e  ∂σ ij ∂ε kl ∂σ kl ∂X ij ∂α kl ∂X kl ∂R ∂r ∂R ℵ=   − ∂F  ∂f ∂σ kl + ∂f ∂X kl + ∂f ∂R + ∂f  ∂Yij  ∂σ kl ∂d ij ∂X kl ∂d ij ∂R ∂dij ∂d ij  

8

The expressions of the different derivatives :

9

∂ σ klD − X klD ∂f = ∂σ ij ∂σ ij

11

12

13

14

∂σ ij ∂ε kle

=

(2.5)

    = 0 (2.6)    

1 ⇒ λ
= a pq D pq ℵ

5

10

∂σ ij ∂F ∂σ ij ∂F + λ
e ∂ε kl ∂σ kl ∂d kl ∂Ykl

∂f ∂σ ij (δ pkδ ql − ϒklpq ) and ∂σ ij ∂ε kle

 εe ) ∂ (Λ ijkl kl

∂ε kle

∂ σ klD − X klD ∂F = ∂σ ij ∂σ ij

c

=

∂

(σ

D pq

       

D − X pq ) H cpqrs (σ rsD − X rsD )

∂σ ij

=

c H ijrs (σ rsD − X rsD ) = nc ij σ D − X D kl

kl c

 =Λ ijkl p

=

∂

(σ

D pq

D p − X pq ) H pqrs (σ rsD − X rsD )

∂σ ij

=

p H ijrs (σ rsD − X rsD ) = n p ij σ D − X D kl

kl

p

∂ σ − X ∂f ∂f = =− = −nijc ∂X ij ∂X ij ∂σ ij ∂X ij ∂ (Cijklα kl )  = = Cijkl ∂α kl ∂α kl ∂ σ klD − X klD ∂( 3a X X ) ∂F p D −1  X mn = + 4C kl kl = −nijp + 23Ca M ijmn ∂X ij ∂X ij ∂X ij D kl

D kl c

D −1 2 D p = −nijp + 23Ca M ijmn 3 CM mnpqα pq = − nij + aα ij

15

 ∂f ∂R ∂  R =− =−  ∂R ∂R ∂R 1 − d  ij 

γ 2

 1  =− 1 − d ij  

γ 2

39

1

Q(1 − d γ ) r  ∂ ij 2 ∂R   = Q(1 − d γ ) =  ij 2 ∂r ∂r  b R − 1 ) ∂  R( ∂F 1 2Q  = ( b R − 1 ) ∂R = =  Q ∂R ∂R ∂R 1 − dij

2 =

3 4

5

∂σ ij ∂d kl

∂X ij ∂dkl

=

=

∂d kl ∂C

ijpq

∂d kl

1 − dij

2

γ

1 − dij  ∂Λ ijpq

γ

 1  ( b 1 − dij r − 1 ) = −  2  1 − d

1 γ 2

γ 2

γ

( Qb Qr − 1 )

2

  − br  

ε epq

α pq

γ γ ∂ Q (1 − dkl 2 ) r  ∂  dkl 2  ∂R     =  = −Qr  = −γ Qr d kl ∂dij ∂dij ∂dij

=−

1

( Qb R − 1 ) =

γ Qr

d kl

3

γ −2 2

γ −1 2

∂  d rs 2    = −γ Qr d kl ∂dij

γ −1

∂ 

2

drr 2   ∂dij 1 3

dij

6  ∂F ∂  S =  ∂Yij ∂Yij ( s + 1) 1 − d kl 

(

=

7

Yrs 3 − Y0 1

)

β

Yuv − Y0

1

(1 − d )

s

∂ Yrs ∂Yij

S

kl 1

=

Yuv 3 − Y0

1

(1 − d )

β

S

8 9

3

=

1

)

β

Yuv − Y0

1

(1 − d )

Yrs − Y0

3

β

S

kl 1

s

3

∂ Yrs

S

∂Yij

3

s

∂ 3YrsYrs ∂Yij

s

3Yij Yrs

kl 1

 1 =   1 − d kl

(

S

3

β

s +1

3

Finally, using the expressions of these last derivatives in the equation Eq.(6.2), we obtain

(

 a pq = nijc Λ ijkl δ pk δ ql − ϒ klpq

10

 + C ) n p − an c X + Q − ℵ = nijc (Λ ijkl ijkl kl ij ij

bR 1 − d rs

11 −

Ymn − Y0

3

(

2 1 − d rs

3

1

)

β

S

s

)

γ 2

  Yij  c ∂Λ ∂C klpq γ Qr klpq e c n n ε − α pq + kl pq kl  ∂d ij ∂d ij Yuv  3 1 − d ot 3 

γ 2

d xz

γ −2 2

  d ij   

12 13

40

1

References

2 3

Al-Rub R.K.A., and Voyiadjis G.Z. , 2003, On the coupling of anisotropic damage and plasticity models for

4 5

Aravas, N., 1986, The analysis of void growth that leads to central burst during extrusion. J. Mech. Phys. Solids

6 7

Aretz, H., Barlat, F., 2013. New convex yield functions for orthotropic metal plasticity. Int. J. Non-Linear

8 9

Badreddine, H., Saanouni K. and Dogui A., 2010, On non-associative anisotropic finite plasticity fully coupled

10 11

Bamman D.J., Solanki K.N., 2010, On kinematics, Thermodynamics and kinetic coupling of a damage theory for

12 13

Banabic, D., Kuwabara, T., Balan, T., Comsa, D.S., Julean, D., 2003, Non-quadratic yield criterion for

14 15

Banabic, D., Aretz, H., Comsa, D.S., Paraianu, L., 2005, An improved analytical description of orthotropy in

16 17

Barlat, F., Lian, J., 1989, Plastic behavior and stretchability of sheet metals. Part I: a yield function for

18 19

Barlat, F., Lege, D.J., Brem, J.C., 1991, A six-component yield function for anisotropic materials. Int. J.

20 21 22

Barlat, F., Becker, R.C., Brem, J.C., Lege, D.J., Murtha, S.J., Hayashida, Y., Maeda, Y., Yanagawa, M., Chung,

23 24

Barlat, F., Ferreira Duarte, J., Gracio, J.J., Lopes, A.B., Rauch, E.F., 2003, Plastic flow for non-monotonic

25 26

Barlat, F., Aretz, H., Yoon, J.W., Karabin, M.E., Brem, J.C., Dick, R.E., 2005. Linear transformation-based

27 28

Barlat, F., Gracio, J.J., Lee, M-G., Rauch, E. F., Vincze, G., 2011, An alternative to kinematic hardening in

29 30

Benzerga AA, Besson J, Batisse R, Pineau A., 2002, Synergistic effects of plastic anisotropy and void

31 32

Benzerga, A. A., Besson, J., Pineau, A. 2004a. Anisotropic ductile fracture: Part I: Experiments. Acta Materialia,

ductile materials, Int. J. Solids Structures, 40, 11, 2611–2643.

34, 55–79.

Mechanics 51, 97-111. doi:10.1016/j.ijnonlinmec.2012.12.007.

with isotropic ductile damage for metal forming, International Journal of Plasticity, 26, pp 1541–1575.

polycrystalline materials, International Journal of Plasticity, 26, pp. 775-793.

orthotropic sheet metals under plane-stress conditions. Int. J. Mech. Sci. 45, 797–811.

metallic sheets. Int. J. Plasticity 21, 493–512.

orthotropic sheets under plane stress conditions. Int. J. Plasticity 5, 51–66.

Plasticity 7, 693–712.

K., Matsui, K., Hattori, S., 1997, Yielding description for solution strengthened aluminum alloys. Int. J. Plasticity 13 (4), 385–401.

loading conditions of an aluminum alloy sheet sample, Int. J. Plasticity 19, 1215–1244.

anisotropic yield functions. Int. J. Plast. 21, 1009–1039.

classical plasticity. Int. J. Plasticity 27, 1309–1327.

coalescence on fracture mode in plane strain. Modell Simul Mater Sci Eng;10:73–102

52, pp. 4623–4638.

33

41

1 2

Benzerga A. A., Besson J., Pineau A., 2004b, Anisotropic ductile fracture Part II: theory, Acta Materialia 52,

3 4

Benzerga, A. A., Leblond, J.B., 2010, Ductile fracture by void growth to coalescence. Advances in Applied

5 6

Beremin FM., 1981, Three-dimensional constitutive relations of damage and fracture, In: Nemat-Nasser S,

7 8

Besson J., Cailletaud G., Chaboche J.L., Forest S., Bletry M., 2010, Non-linear mechanics of Materials,

9 10

Betten, J., 1986, Applications of tensor functions to the formulation of constitutive equations involving damage

11 12

Bonora, N., Gentile, D., Pirondi, A., Newaz, G., 2005. Ductile damage evolution under triaxial state of stress:

13 14

Boudifa M., Saanouni K;, Chaboche J.L., 2009, A micromechanical model for inelastic ductile damage

15 16

Brünig, M., 2002. Numerical analysis and elastic–plastic deformation behavior of anisotropically damaged

17 18

Brünig, M., 2003. An anisotropic ductile damage model based on irreversible thermodynamics. Int. J. Plast. 19,

19 20

Brünig, M., Ricci S., 2005, Nonlocal continuum theory of anisotropically damaged metals International Journal

21 22

Brünig M., Gerke S., 2011a, Simulation of damage evolution in ductile metals undergoing dynamic loading

23 24

Brünig, M., Albrecht, D., Gerke, S., 2011b, Modeling of ductile damage and fracture behavior based on

25 26 27

Carol, I., Rizzi, E., Willam, K., 2001. On the formulation of anisotropic elastic degradation. Part I: Theory based

28 29

Cazacu, O., Barlat, F., 2004, A criterion for description of anisotropy and yield differential effects in pressure

30 31

Chaboche, J.L., 1978, “Description thermodynamique et phénoménologique de la viscoplasticité cyclique avec

4639–4650.

Mechanics, 44, pp. 169–305.

editor. New York: Pergamon Press.

Springer, Dordrecht, Germany.

and initial anisotropy. Engineering Fracture Mechanics, vol. 25, no. 5-5, pages 573–584.

Theory and experiments. International Journal of Plasticity, 21, 981–1007.

prediction in polycrystalline metals, Int. Journal of Mechanical Sciences, 51, pp :453-464.

solids. Int. J. Plast. 18, 1237–1270.

1679–1713.

of Plasticity 21, p. 1346–1382.

conditions, , International Journal of Plasticity, Volume 27, Issue 10, Pages 1598–1617.

different micromechanisms. International Journal of Damage Mechanics 20 (4), 558–577.

on a pseudo-logarithmic damage tensor rate, Part II: Generalized pseudo-Rankine model for tensile damage. Int. J. Solids Structures 38, 491–518, 519–546.

insensitive metals. Int. J. Plasticity 20, 2027–2045.

endommagement”, Thèse université Paris VI et publication ONERA n° 1978-3.

42

1 2

Chaboche, J.L., 1992, “Damage Induced Anisotropy: On the Difficulties Associated with the Active/Passive

3 4

Chaboche, J.L., 1993, “Development of Continuum Damage Mechanics for Elastics Solids Sustaining

5 6

Chaboche , J.L., Lesne, P.M. and Maire, J.F., 1995, “Continuum Damage Mechanics, anisotropy and damage

7 8

Chaboche J.L., Boudifa M., Saanouni M., 2006, “A CDM approach of ductile damage with plastic

Unilateral Condition” Int. J. of Damage Mechanics, 1(2), pp. 148-171.

anisotropic andunilateral damage” Int. J. Damage Mech., vol. 2, pp. 311-329.

deactivation for brittle materials like concrete and ceramic composites”, Int. J. Damage Mechanics, 4, 5-22.

compressibility", International Journal of fracture, N°137, pp: 51-75.

9 10

Chen ,XF., and Chow, CL. , 1995, On damage strain energy release rate Y. International Journal of Damage,

11 12

Chow, C.L. and Wang, June, 1987, An anisotropic theory of continuum damage mechanics for ductile fracture.

13 14

Chow, C.L and Lu, T.J. , 1989, On evolution laws of anisotropic damage. Engineering Fracture Mechanics, vol.

15 16

Chow, C.L. , 1999, Constitutive Modeling of Material Damage for Fatigue Failure Prediction. International

17 18

Cordebois, J.P., Sidoroff, F., 1979, “Anisotropie élastique induite par endommagement” Colloque Euromech

19 20

Cordebois, J.P., Sidoroff, F., 1980. Anisotropic damage in elasticity and plasticity. Journal of Theory and

21 22

Desmorat, R., 2000, Quasi-unilateral conditions in anisotropic elasticity.Comptes Rendus de l’Académie des

23 24

Desmorat, R., Cantournet, S., 2007,Modeling Micro-defects Closure Effect with Isotropic/Anisotropic Damage.

25 26

Dragon, A., and Mroz, Z., 1979, A continuum model for plastic-brittle behavior of rock and concrete, Int. J.

27 28

Ekh, M., Menzel, A., Runesson, K., Steinmann, P., 2003, Anisotropic damage with the MCR-effect coupled to

29 30

Ekh, M., Lillbacka, R., Runesson, K., 2004, A model framework for anisotropic damage coupled to crystal

31 32

François D., Pineau A., Zaoui A., 1993, Comportement mécanuique des matériaux. Vol2 : Endommagement,

Mechanics 4, pp. 251–263.

Engineering Fracture Mechanics, vol. 27, no. 5, pages 547 – 558.

34, no. 3, pages 679 – 701.

journal of DAMAGE MECHANICS, vol. 8, pages 355–375.

115, Grenoble, Editions du CNRS, n° 395.

Applied Mechanics, Numero Special, pp. 45–60.

Sciences - Series IIB - Mechanics, vol. 328, no. 6, pages 445 – 450.

International Journal of damage mechanics, vol. 00, pages 1–32.

Engrg. Sci., 17, 121-137.

plasticity. International Journal Engineering Sciences 41, 1535–1551.

(visco)plasticity. International Journal of Plasticity, 20 (12), 2143–2159.

mécanique de la rupture , mécanique du contact, Hermes, Paris .

43

1 2

Gao X., Zhang T., Hayden M., Roe C., 2009. Effects of the stress state on plasticity and ductile failure of an

3 4

Gelin J.C., 1990, Finite Element Analysis of Ductile Fracture and defects Formations in Cold and Hot Forging,

5

Germain, P., Nguyen, Q.S., Suquet, P., 1983, Continuum thermodynamics. J. Appl. Mech. 50, 1010–1020.

6 7

Gologanu, M., Leblond, J.-B., Perrin, G., Devaux, J., 1995, Recent extensions of Gurson’s model for porous

8 9

Gurson, A.L., 1977, Continuum theory of ductile rupture by void nucleation and growth: Part I – Yield criteria

aluminum 5083 alloy. International Journal of Plasticity 25, 2366–2382.

Annals of the CIRP, 39, pp.215-218.

ductile metals, in: Suquet, P. (Ed.), Continuum Micromechanics, Springer, pp. 61–130.

and flow rules for porous ductile media, J. Eng. Mat. Tech. 99, 2–15.

10 11 12

Haddag, B. , ABED-Meraim, F., Balan, T., 2009, Strain localization analysis using a large deformation

13 14 15

Hammi, Y., Horstemeyer, M.F., 2007, A physically motivated anisotropic tensorial representation of damage

16 17 18

Hansen, N.R., 1994, and Schreyer, H.L., A thermodynamically consistent framework for theories of

19 20

Hill, R., 1948. A theory of the yielding and plastic flow of anisotropic metals. Proc. Roy. Soc. London A193,

21 22

HILL, R., 1990, Constitutive modelling of orthotropic plasticity in sheet metals,.J. Mech. Phys. Solids, 38,

23 24

HILL, R., 1993, A user-friendly theory of orthotropic plasticity in sheet metals, Int. J. of Mech. Sci., 35, Issue 1,

25 26 27

Issa, M., Saanouni, K., Labergère, C., Rassineux, A., 2011. Prediction of serrated chip formation in orthogonal

28 29

Issa M., Labergère, C., Saanouni, K., Rassineux, A., 2012. Numerical prediction of thermomechanical fields

30 31

Ju., J.W., 1989, On energy-based coupled elastoplastic damage theories : Constitutive modeling and

32

Kachanov, M., 1980, Continuum model of medium with cracks, J. Engrg. Mech., ASCE 106(5), 1039-1051.

anisotropic elastic–plastic model coupled with damage, International Journal of Plasticity, 25, (10), pp. 19701996.

with separate functions for void nucleation, growth, and coalescence, Int. J. of Plas., Volume 23, Issues 10–11, Pages 1641-1678.

elastoplasticity coupled with damage. International Journal of Solids and Structures, vol. 31, no. 3, pages 359 – 389.

281–297.

pp.405-417.

19-25.

metal cutting by advanced adaptive 2D numerical methodology. Int. J. Machining and Machinability of Materials, 9, Nos. 3/4.

localization in orthogonal cutting, CIRP Journal of Manufacturing Science and Technology, 5, pp.175–195.

computational aspects. International Journal of Solids and Structures, vol. 25, no. 7, pages 803 – 833.

44

1 2

Kachanov, L. M., 1986, Introduction to Continuum Damage Mechanics , Netherlands, Martinus Nijhoff

3 4

Karafillis, A.P., Boyce, M.C., 1993 , A general anisotropic yield criterion using bounds and a transformation

5 6

Keralavarma S.M., Hoelscher S., Benzerga A.A., 2011, Void growth and coalescence in anisotropic plastic

7 8 9

Khelifa, M., Badreddine, H., Gahbiche, M.-A., Saanouni, K., Cherouat, A., Dogui, A., 2005, Effect of

Publishers.

weighting tensor. J. Mech. Phys. Solids 41, 1859–1886.

solids, International Journal of Solids and Structures 48 1696–1710.

anisotropic plastic flow on the ductile damage evolution in sheet metal forming. Application to the circular bulging test, International Journal of Forming Processes. Vol. 8, N°:2, pp. 271-289.

10 11

Labergère, C., Rassineux, A., Saanouni K., 2011., 2D adaptive mesh methodology for the simulation of metal

12 13

Ladevèze P. , 1993, On an Anisotropic Damage Theory, Failure Criteria of Structured Media, Ed. Boehler,

14 15

Lee, H., Peng, K., and Wang, J., 1985, An anisotropic Damage Criterion for Deformation Instability and its

16 17

Lehmann, Th. , 1991, Thermodynamical foundations of large inelastic deformations of solid bodies including

18 19

Lemaitre, J. and Desmorat, R., 2005, Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle

20

Lemaitre, J., 1992. A Course on Damage Mechanics. Springer, Berlin.

21

Lemaitre, J., Chaboche, J.L., 1985, Mécanique des matériaux solides. Dunod, Paris, 1ème éditions.

22 23

Lemaitre, J., Chaboche, J.L., BENALLAL, A., et Desmorat., R., 2009, Mécanique des matériaux solides.

24 25

Lemaitre, J., Desmorat, R., et Sauzay, M., 2000, Anisotropic damage law of evolution. Eur. J. Mech. A/Solids,

26 27

Marigo, J.J., 1985, Modelling of brittle and fatigue damage for elastic material by growth of microvoids,

forming processes with damage, Int. J. Mater. Forming, 4, pp.317–328.

Balkema,Rotterdam, 355-363.

Application to Forming Limit Analysis of Metal Plates, Engng. Fract. Mechs., 21(5):1031-1054.

damage. Int. J. Plasticity 7, 79–98.

Failures, Springer.

Dunod, Paris, 3ème éditions.

vol. 19, pages 187–208.

Engineering Fracture Mechanics, Volume 21, Issue 4, Pages 861-874.

45

1 2 3

Mazars, J., 1986, "A Model of a Unilateral Elastic Damageable Material and its Application to concrete",

4 5

Mazars, J., Berthaud, Y., Ramtani, S., 1990, Engineering Fracture Mechanics, Volume 35, Issues 4–5, Pages

6 7

Menzel A., and Steinmann P., 2001, A theoretical and computational framework for anisotropic continuum

8 9

Menzel A., Ekh M., Runesson K., Steinmann P., 2005, A framework for multiplicative elastoplasticity with

RILEM Int. Conf. Fracture Mechanics of Concrete, Lausanne, 1985, Fracture Toughness Fracture Energy of Concrete, New York : Elsevier.

629-635.

damage mechanics at large strains, Int. J. Solids Structures, 38, 9505–9523.

kinematic hardening coupled to anisotropic damage, International Journal of Plasticity 21, p. 397–434.

10 11

Murakami , S., and Ohno, N., 1981, A continuum theory of creep and creep damage. In A.R.S.Ponter, editeur,

12 13 14

Murakami, S., 1987, Anisotropic damage theory and its application to creep crack growth analysis, Constitutive

15

Murakami, S., 1988. Mechanical modeling of material damage. J. Appl. Mech. 55 (2), 280–286.

16 17

Murakami, S., 2012, Continuum Damage Mechanics – A Continuum Mechanics Approach to the Analysis of

18 19

Needleman A., and Triantafyllidis N., 1980, Void growth and local necking in biaxially stretched sheets, J. Eng.

20 21

Oda, M., 1983, A method for evaluating the effect of crack geometry on the mechanical behavior of cracked

22 23

Onate, E. and Kleiber M., 1988, Plastic and Viscoplastic Flow of Void Containing Metal - Applications to

24 25

Ortiz, M., 1985, A constitutive theory for the inelastic behavior of concrete. Mech. Mater, vol. 4, no. 1, pages

26 27

Pardoen, T., Doghri, I., Delannay, F., 1998, Experimental and numerical comparison of void growth models and

28 29

Rice, J.R., Tracey, D.M., 1969, On the enlargement of voids in triaxial stress fields, J. Mech. Phys. Sol. 17, 201–

30 31

Rousselier, G.,1987, Ductile fracture models and their potential in local approach of fracture, Nuclear

32 33

Rousselier, G., 2001, Dissipation in porous metal plasticity and ductile fracture, Journal of the Mechanics and

Creep of Structure IUTAM symp., pages 422– 444, Berlin, Federal Republic of Germany, Springer-Verlag.

laws for engineering materials: Theory and applications. C S. Desai et al., eds., Elsevier, New York, N.Y., 187194.

Damage and Fracture, Springer: Dordrecht, Heidelberg, London, New York.

Mat. Technol., 100, p:164-172.

rock masses, Mech. of Mat., 2, 163-171.

Axisymmetric Sheet Forming Problem , Int. J. Num. Meth. In Engng. 25:237-251.

67–93.

void coalescence criteria for the prediction of ductile fracture in copper bars. Acta Mater. 46, 541–552.

217.

Engineering and Design, Volume 105, Issue 1, Pages 97-111.

Physics of Solids, Volume 49, Issue 8, Pages 1727-1746.

46

1 2

Saanouni, K., Forster, Ch., and Ben Hatira, F., 1994. «On the anelastic Flow with Damage», J.Damage

3 4

Saanouni, K., and Hammi, Y., 2000, Numerical simulation of damage in metal forming processes, Continuous

5 6

Saanouni, K., Cherouat, A., and Hammi, Y., 2001, Numerical aspects of finite elastoplasticity with isotropic

7 8 9

Saanouni, K., and Chaboche, J.L., 2003, Computational damage mechanics. Application to metal forming.

10 11

Saanouni, K., 2006. Virtual metal forming including the ductile damage occurrence: Actual state of the art and

12 13

Saanouni, K., 2008. On the numerical prediction of the ductile fracture in metal forming. Engineering Fracture

14 15 16

Saanouni, K., Lestriez, P., Labergère, C., 2011. 2D adaptive simulations in finite thermo-elasto-viscoplasticity

17 18

Saanouni, K., 2012. Damage Mechanics in metal forming: Advanced modeling and numerical simulation,

19 20

Sidoroff, F., The geometrical concept of intermediate configuration and elastic-plastic finite strain. Arch. Mech.,

21 22

Sidoroff, F., 1982, Incremental constitutive equation for large strain elastoplasticity, Int. J. Engng Sc., vol.20, No

23 24

Sidoroff, F., and Dogui, A., 2001, Some issues about anisotropic elastic-plastic models at finite strain », Int. J.

25 26 27

Simo, J.C., and Ortiz, M., 1985, A Unified Approach to Finite Element Deformation Elastoplastic Analysis

28 29

Simo, J.C., and Ju, J.W., 1987, Stress and strain based continuum damage models. International Journal of Solids

30

Simo, J.C., Hughes, T.J.R., 1998, Computational Inelasticity, Springer – Verlag, New York, Inc.

Mechanics, vol. 3, pp. 140-169.

Damage and Fracture, Edition A. Benallal, Elsevier, ISBN. 2-84299-247-4, pp.353-363.

ductile damage for metal forming, Revue européenne des éléments finis, 2-3-4, pp. 327-351.

Chapter 7 of Vol. 3 “Numerical and Computational methods”, Ed. I. Milne, R.O. Ritchie and B. Karihaloo, ISBN 0-08-043749-4, Elsevier Oxford (UK), p. 321-376.

main perspectives. Journal of Materials Processing Technology, 177, pp. 19–25.

Mechanics, 75, pp. 3545-3559.

with ductile damage: Application to orthogonal metal cutting by chip formation and breaking. Int. J. of Damage Mechanics, 20, pp. 23–61.

ISTE/Wiley, London.

vol. 25, no. 2, pages 299–308, 1973.

1, pp.19-26.

Sol. Str. 38, 9569-9578.

Based on the Use of Hyperelastic Constitutive Equations", Comp. Meth. Appl. Meth. Engng., vol. 49, pp. 221245.

and Structures, vol. 23, no. 7, pages 841 – 869.

47

1 2

Steglich, D., Brocks, W., 1997. Micromechanical modeling of the behavior of ductile materials including

3 4

Steglich D, Brocks W, Heerens J, Pardoen T., 2008, Anisotropic ductile fracture of Al 2024 alloys. Engng Fract

5 6

Steglich D., Wafai H., Besson J., 2010a, Interaction between anisotropic plastic deformation and damage

7 8

Steglich D, Wafai H, Brocks W., 2010b, Anisotropic deformation and damage in aluminium 2198 T8 Sheets. Int

9 10

Steinmann, P., and Carol, I. , 1998, "A framework for geometrically nonlinear continuum damage mechanics",.

11 12

Suaris, W., 1987, A damage theory for concrete incorporating crack growth characteristics, Constitutive laws for

13 14

Talreja, R. , 1985, A Continuum Mechanics Characterization of Damage composite Material, Proc. Royal

15 16

Teng, X., and Wierzbicki , T., 2006, Evaluation of six fracture models in high velocity perforation,

17

Tvergaard, V., 1990, Material failure by void growth to coalescence. Adv. Appl. Mech., 27, 83-151.

18 19

Vakulenko, A. A., and Kachanov, M. L., 1971, Continuum theory of cracked media (in Russian), Mekh

20 21

Vignjevic , R., Djordjevic, N., Campbell, J., Panov, V., 2012a, Modelling of dynamic damage and failure in

22 23

Vignjevic , R., Djordjevic, N., Panov, V., 2012b, Modelling of dynamic behaviour of orthotropic metals

24 25

Vial, C., Hosford, W.F., Caddell, R.M., 1983, Yield loci of anisotropic sheet metals, Int. J. of Mech. Sci., 25,

26 27

Voyiadjis, G.Z., and Kattan, P.I., 1992a. A plasticity-damage theory for large deformations of solids. Part I:

28 29

Voyiadjis, G.Z., and Kattan, P.I., 1992b. Finite strain plasticity and damage in constitutive modeling of metals

particles. Computational Material Science, 9, p. 7–17.

Mech;75:3692–706.

evolution in Al 2198 sheet metal, Engineering Fracture Mechanics, Volume 77, Issue 17, Pages 3501–3518.

J. Damage Mech;19(2):131–52.

Int. J. Eng. Sci. 36, 1793–1814.

engineering materials: theory and applications, C S. Desai et al., eds., Elsevier, New York, N.Y.

Society, London, pp. 195-216.

Engineering Fracture Mechanics, Volume 73, Issue 12, Pages 1653-1678.

Tverdogo Tela. U.S.S.R ., 6,159-166.

aluminium alloys, International Journal of Impact Engineering, Volume 49, Pages 61-76.

including damage and failure, International Journal of Plasticity, Volume 38, Pages 47-85.

Issue 12, 899-915.

theoretical formulation. Int. J. Eng. Sci. 30, 1089–1108.

with spin tensors. App. Mech. Rev. 45, S95–S109.

48

1 2

Voyiadjis, G.Z. and Kattan, P.I., 1996, On the Symmetrization of the Effective Stress Tensor in Continuum

3 4

Voyiadjis, G.Z., Kattan, P.I., 1999. Advances in Damage Mechanics: Metals and Metal Matrix Composites.

5 6 7

Voyiadjis, G.Z., Taqieddin, Z.N., et Kattan, P.I., Theoretical Formulation of a coupled elastic-plastic Anisotropic

8 9

Wierzbicki, T., Bao, Y., Lee , Y-W., Bai, Y., 2005, Calibration and evaluation of seven fracture models,

Damage Mechanics. Journal of the Mechanical Behavior of Materials, vol. 7, no. 2, pages 139–165.

Elsevier Oxford, UK.

Damage Model for concrete using the strain Energy Equivalence Concept. International journal of Damage mechanics, pages 603–637, 2009.

International Journal of Mechanical Sciences, Volume 47, Issues 4–5, Pages 719-743 .

10 11

Yazdani, S., and Schreyer, H. L., 1988, An anisotropic damage model with dilatation for concrete, Mech. of

12 13

Yoon, J.W., Barlat, F., Dick, R.E., Chung, K., Kang, T.J., 2004, Plane stress yield function for aluminum alloy

14 15

Zhu, Y. Y., Cescotto, S., 1995. A fully coupled elasto-visco-plastic damage theory for anisotropic materials.

Mat., 7, 231-244.

sheets – Part II: FE formulation and its implementation. Int. J. Plasticity 20, 495–522.

Int.J.Molids and Structures, 32, pp. 1607–1641.

16

49

1

Ct

Vije

Fij

Qij

Fij

C0

Ct
Fijp

Vije

Fijp

Cet

p Ct

Qij

2 3

Fije

Figure 1. Kinematic decomposition of the total deformation gradient: rotating frame concept

4

50

1 2

3 4

Figure 2. Schematic representation of the anisotropic damage in a typical RVE.

5

51

1 2

Stress J2(σ) [100 MPa]

3 16 14 12 10

isotropic damage γ=10 Anisotropic damage γ=2 γ=3 γ=6 γ=10

8 6 4 2 0 0,0

4 5 6

0,2

0,4

0,6

0,8

1,0

Equivalent plastic strain p Figure 3. Uniaxial tension results in terms of stress versus equivalent plastic strain with isotropic plastically generic material.

7

52

1

8

0,5

Isotropic damage case with γ=10

7 6

Kinematic hardening α Energy density release Y I

Elastic Energy density release YI

e

2

Anisotropic damage case γ=2 γ=3 γ=6 γ=10

5 4 3 2 1 0 0,0

0,2

0,4

0,6

0,8

0,4

Anisotropic damage case α YI γ=2 γ=3 γ=6 γ=10

0,3 0,2 0,1 0,0 0,0

0,2

Isotropic hardening r Energy density release YI

2,0

1,4 1,2 1,0 0,8 0,6

Isotropic damage case with γ=10 r Y Anisotropic damage case r YI γ=2 γ=3 γ=6 γ=10

0,4 0,2 0,0 0,0

0,2

0,4

0,6

0,8

Equivalent plastic strain p

(c) Isotropic hardening energy density release

3

0,6

0,8

(b) Kinematic hardening energy density release Total Energy density release YI

(a) Elastic energy density release

1,6

0,4

Equivalent plastic strain p

Equivalent plastic strain p

1,8

Isotropic damage case with γ=10 α Y

8 7

Anisotropic damage case YI

6

γ=2 γ=3 γ=6 γ=10

5 4 3

Isotropic damage case with γ=10

2 1 0 0,0

0,2

0,4

0,6

0,8

Equivalent plastic strain p

(d) Total energy density release

Figure 4. Uniaxial tension results in terms of damage energies versus the equivalent plastic strain.

4

53

0,5

0,5

0,4

0,4

0,4

0,3 0,2 0,1

0,3 0,2 0,1

0,0 0,0 0,2 0,4 0,6 0,8

1,0

1,2

1,4 1,6 1,8

Equivalente plastic strain p

(a) d I .

2

Damage dIII

0,5

Damage dII

Damage dI

1

0,0 0,0 0,2

0,3 0,2 0,1

0,4

0,6 0,8

1,0

1,2

1,4 1,6

Equivalente plastic strain p

(b) d II .

1,8

0,0 0,0

0,2

0,4

0,6

0,8

1,0

1,2

1,4

1,6

1,8

Equivalente plastic strain p

(c) d III .

Figure 5. Evolution of the damage tensor eigenvalues for the different loading paths.

3

54

1 1,5 1,0

 e2 '

π − plane ∆l

εII

D

0,5

π θ=

6

0,0

 e1 '

max damage dI=0.1

-0,5 -1,0 -1,5 -1,5

2

 e3 ' -1,0

isotropic damage Anisotropic damage γ=2 γ=3 γ=6 γ=7 γ=8 γ=9 γ=10

max damage dI=0.9

-0,5

0,0

εI

0,5

1,0

1,5

D

3

Figure 6. Isotropic and anisotropic isodamage curves plotted in the deviatoric strain space (π-plane) for

4

d I = 0.1 , d I = 0.9 and ϖ = 0

5

55

1 2,0 1,5

 e2 '

h=0 h=0.2 h=0.4 h=0.6 h=0.8 h=1

1,0

ε II

D

0,5 0,0

θ = π3

2π 3

-0,5

 e1 '

-1,0 -1,5 -2,0 -2,0

 e3 ' -1,5

π − plane -1,0 -0,5

0,0

0,5

1,0

1,5

2,0

D

εI

(c) Isotropic damage 2,5 2,0

h=0 h=0.2 h=0.4 h=0.6 h=0.8 h=1

 e2 '

1,5 1,0

θ = π3

εII

D

0,5

2π 3

0,0 -0,5

 e1 '

-1,0 -1,5 -2,0

 e3 '

π − plane

-2,5 -2,5 -2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 2,5 D

εI

(d) Anisotropic damage

2 3

Figure 7. Effect of the microcracks closure parameter on the isodamage surfaces for isotropic and anisotropic damage plotted in the diviatoric strain space (π-plane)

4

56

1 2 3

Equivalent Plastic strain at rupture ρ fr

1,8

− 13

1,6

1 3

2 3

3 3

0

1,4

Anisotropic damage h=1 h=0 Isotropic damage h=1 h=0

1,2 1,0 0,8 0,6 0,4 -0,2

4

0,0

0,2

Triaxiality

0,4

0,6

0,8

ϖ

5

Figure 8. Effect of the microcracks closure parameter on the plane stress damage envelop defined for a

6

maximum damage level d I = 0.9 .

7

57

1 2 3 4

Equivalent Plastic strain at rupture ρ fr

1,8

− 13

1,6

1 3

2 3

3 3

η =0.2 η =0.6 η =0.8 η =1.0 η =1.2 η =1.5 η =2.0 η =3.0

0

1,4 1,2 1,0 0,8 0,6 0,4 -0,2

5

0,0

0,2

0,4

0,6

0,8

Triaxiality ϖ

6

Figure 9. Effect of the hydrostatic damage coupling parameter η on the material ductility versus stress triaxiality

7

for anisotropic damage with =0 for d I = 0.9

8

58

1 2 3

(a) =1.

(b) =0.

4

Figure 10. Effect of the microcracks closure parameter on the 3D damage surfaces defined for a maximum

5

damage level d I = 0.9 (near the final rupture).

6

59

0,25

0,5

0,4

0,20

0,4

0,3 0,2

0,2

0,4

0,6

0,8

1,0

Equivalente plastic strain p

(a) d I .

2

0,15 0,10 0,05

0,1 0,0 0,0

Damage dIII

0,5

Damage dII

Damage dI

1

1,2

0,00 0,0

0,3 0,2 0,1

0,2

0,4

0,6

0,8

1,0

Equivalente plastic strain p

(b) d II .

1,2

0,0 0,0

0,2

0,4

0,6

0,8

1,0

1,2

Equivalente plastic strain p

(c) d III .

Figure 11. Evolutions of the eigenvalues of the damage tensor for the different loading paths.

3

60

1 2 3 2,5 2,0

 e2 '

1,5 1,0

εII

D

0,5

2π 3

0,0 -0,5

 e1 '

Isotropic plasticity Anisotropic damage h=0 Anisotropic damage h=1 Isotropic damage h=1 Anisotropic plasticity Anisotropic damage h=0 Anisotropic damage h=1 Isotropic damage h=1

-1,0 -1,5 -2,0

 e3 '

π − plane

-2,5 -2,5 -2,0 -1,5 -1,0 -0,5 0,0 0,5 1,0 1,5 2,0 2,5 D

εI

4 5

Figure 12. Effect of plastic anisotropy on the isodamage surface using isotropic and anisotropic damage plotted

6

in the diviatoric strain π-plane for d I = 0.9

7

61

1 2 3 6

5

initial case of h=0 dI=0.1

4

dI=0.5

3

case of h=1 dI=0.1

Lankford ratio rψ

Lankford ratio rψ

6

Increasing damage

dI=0.5

2 1

10

20

30

40

50

60

70

80

dI=0.5

4

case of h=1 dI=0.1

3

dI=0.5

2 1 Increasing damage 0 0 10 20

0 0

initial case of h=0 dI=0.1

5

90

ψ0[°]

2 anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

Normalized Stress σ90/σy

Normalized Stress σ90/σy

8

dI=0.5

4

dI=0.5

-8 -8

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

60

70

80

90

6

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case with h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case with h=1 dI=0.1

-4 -6

dI=0.5

-8

8

-8

(c) Yield surface for ϕ0 = 0°

4 5

50

(b) Lankford ratio for ϕ0 = 90°

initial surface isotropic damage case dI=0.1

6

40

ψ0[°]

(a) Lankford ratio for ϕ0 = 0° 8

30

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

6

8

(c) Yield surface for ϕ0 = 90°

Figure 13. Damage induced change on Lankford ratio curves and yield surface due to UT path with two material orientations ϕ0 .

6

62

1 2

5

initial case of h=0 dI=0.1

4

dI=0.5

3

case of h=1 dI=0.1

8

Increasing damage h=0

dI=0.5

2 1

Increasing damage h=1

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

dI=0.5

-8

0 0

10

20

30

40

50

60

70

80

90

ψ0[°]

(a) Damage induced change on the Lankford ratio evolution.

3 4

Normalized Stress σ90/σy

Lankford ratio rψ

6

-8

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

6

8

(b) Damage induced change on yield surface displayed in the normalized stress space.

Figure 14. Effect of the induced damage anisotropy on the yield surface and the Lankford evolution after (PT) loading.

5

63

1

initial case of h=0 dI=0.1

5 4

8

Normalized Stress σ90/σy

Lankford ratio rψ

6

dI=0.5

Increasing damage h=0

case of h=1 dI=0.1

3

dI=0.5

2 1

initial surface isotropic damage case dI=0.1

6

dI=0.5

4 2

anisotropic damage case h=0 dI=0.1

0 -2

dI=0.5 anisotropic damage case h=1 dI=0.1

-4 -6

Increasing damage h=1

0 0

10

20

dI=0.5

-8 -8

30

40

50

60

70

80

90

-6

-4 -2 0 2 4 Normalized Stress σ0/σy

6

8

ψ0[°] (a) Damage induced change on the Lankford ratio evolution

2 3

(b Damage induced change on the yield surface displayed in the normalized stress space

Figure 15. Effect of the induced damage anisotropy on the yield surface and the Lankford evolution after (EB) loading.

4 5

64

1 2 3 4 5 6 7 8 9

Table 1: Elastic and plastic parameters of the generic material. Parameters

σy

E

ν

Q

b

C

a

Unit

MPa

GPa

-

MPa

-

MPa

-

300.0

210.0

0.3

2600.0

1.50

7000.0

170.0

10 11

65

1 2 3

Table 2: Damage parameters of the generic material given for the isotropic and anisotropic damage cases. Parameters Isotropic damage

Anisotropic damage

S

s

Y0

β

γ

η



Dc

6.8

1.3

0.0

1.0

10.0

1.0

0.0 / 1.0

0.999

60.0

1.3

0.0

2.0

2.0

1.0

0.0 / 1.0

0.999

40.0

0.7

0.0

2.0

3.0

1.0

0.0 / 1.0

0.999

10.5

1.3

0.0

1.0

6.0

1.0

0.0 / 1.0

0.999

10.5

1.3

0.0

1.0

10.0

1.0

0.0 / 1.0

0.999

4 5

66

1 2 3 4

Table 3: The loading paths considered and their corresponding damage tensor with plastically isotropic material

Uniaxial Tension (UT)

dij

Dij

Loading

2 0 0  0 −1 0  2 0 0 −1

ε


dI 0   0

0 0  d II 

0 d II 0

with d I > d II = 14 d I

Plane Tension (PT) Equi-Biaxial Tension (EB)

1 0 0  ε
0 0 0   0 0 −1

1 0 0  ε
 0 1 0   0 0 −2

dI 0   0 dI 0   0

0 dI 0

0 0  0 d I  0

0

0  0  d III 

with d III > d I = 14 d III

Scheme

=0

=1

dI 0   0

0 0 0 0  0 0 

dI 0   0

0 0 0 0  0 0 

dI 0   0

dI

0 0

0 0  0 

5 6

67

1 2

Table 4: The loading paths considered and their corresponding damage tensor with anisotropic plasticity

ϕ0 = 0° and ϕ0 = 90°

Uniaxial Tension (UT)

d ij

Dij

Loading

2 0 0  0 −1 0  2 0 0 −1

ε


ϕ 0 = 45°

0  d11 0 0 d  0 22    0 0 d33  with d11 > d33 > d 22

 d11 d  12  0

d12 d 22 0

0 0  d33 

with φ = 0.674° d I > d 33 ≈ d II

Plane Tension (PT)

1 0 0  ε
0 0 0   0 0 −1

Equi-Biaxial Tension (EB)

1 0 0  ε
 0 1 0   0 0 −2

Scheme

=0

=1

dI 0   0

0 0 0 0  0 0 

 d11 d12 0  d   12 d 22 0   0 0 0  with φ = 0.673°

0  d11 0 0 d 0  22   0 0 d33  with d33 > d11 > d 22

 d11 0 0  0 d 0  22   0 0 0  with d11 >> d 22

 d11 0   0

 d11 0   0

0 d 22 0

0 0  d33 

with d33 > d11 > d 22

0 d22 0

0 0  0 

with d11 > d 22

3 4

68

1 2 3

Table 5: Damage eigenvalues obtained with the different loading paths and for a 0.5 maximum damage level Loading

UT

ϕ 0 = 0° ϕ 0 = 45° ϕ 0 = 90° PT EB

dI

= 1.0 d II

d III

dI

= 0.0 d II

d III

0.500

0.038

0.140

0.500

0.

0.

0.500

0.065

0.095

0.500

0.0002

0.

0.500 0.203 0.116

0.022 0.019 0,050

0,188 0.500 0.500

0.500 0.500 0.500

0. 0.083 0.217

0. 0. 0.

4 5 6

69

1 2

Highlights

3 4

•

The anisotropic ductile damage is discrebeb by second-rank tensor.

5 6

•

The effect of anisotropic damage on the mechanical properties is introduced through the equivalent total

7 8

•

9 10

•

energy assumption. For an accurate plastic anisotropy description, non-associative plasticity framework is considered for which equivalent stresses in yield function and in plastic potential are separately defined. The damage anisotropy inducse an important change on the shape of the iso-damage envelop which appears with a dependence on lode angle.

11

70