Effect of compression–tension loading reversal on the strain to fracture of dual phase steel sheets

International Journal of Plasticity 72 (2015) 21e43

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Effect of compressionetension loading reversal on the strain to fracture of dual phase steel sheets
phane J. Marcadet a, Dirk Mohr a, b, * Ste a

Impact and Crashworthiness Laboratory, Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge MA, USA b
Solid Mechanics Laboratory (CNRS-UMR 7649), Department of Mechanics, Ecole Polytechnique, Palaiseau, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 October 2014 Received in revised form 3 May 2015 Available online 16 May 2015

The effect of loading direction reversal on the onset of ductile fracture of DP780 steel sheets is investigated through compressionetension experiments on flat notched specimens. A finite strain constitutive model is proposed combining a Swift-Voce isotropic hardening law with two FrederickeArmstrong kinematic hardening rules and a YoshidaUemori type of hardening stagnation approach. The plasticity model parameters are identified from uniaxial tension-compression stressestrain curve measurements and finite element simulations of compressionetension experiments on notched specimens. The model predictions are validated through comparison with experimentally-measured load edisplacement curves up to the onset of fracture, local surface strain measurements and longitudinal thickness profiles. In addition, the model is used to estimate the local strain and stress fields in monotonic fracture experiments covering plane stress states ranging from pure shear to plane strain tension. The extracted loading paths to fracture show a significant increase in ductility as a function of the compressive pre-strain. A HosfordCoulomb damage indicator model is presented to provide a phenomenological description of the experimental results for monotonic and reverse loading. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Ductile fracture In-plane compression Kinematic hardening Reverse loading Hosford-Coulomb

1. Introduction Predicting the onset of ductile fracture has been an active field of research for more than 50 years. In particular, the fracture initiation after monotonic proportional loading paths has been investigated intensively (e.g. Brunig et al. (2008), Bai and Wierzbicki (2008, 2010), Sun et al. (2009), Li et al. (2011), Gruben et al. (2011), Chung et al. (2011), Lecarme et al. (2011), Khan and Liu (2012), Luo et al. (2012), Huespe et al. (2012), Malcher et al. (2012), Lou et al. (2014)). In industrial practice, in particular during sheet metal forming, ductile fracture often initiates after complex non-proportional loading histories. Among these, reverse loading is an important non-proportional loading condition which prevails for instance when a sheet is bent and unbent as it is drawn over a die radius. Simulating the mechanical response of ductile materials up to the point of fracture initiation requires the accurate modeling and identification of the hardening behavior of the material at large strains. Many plasticity models for reverse loading have been developed for life-cycle analysis. As a result, most experimental procedures are designed for characterizing


*Corresponding author. Solid Mechanics Laboratory (CNRS-UMR 7649), Department of Mechanics, Ecole Polytechnique, Palaiseau, France. E-mail address: [email protected] (D. Mohr). http://dx.doi.org/10.1016/j.ijplas.2015.05.002 0749-6419/© 2015 Elsevier Ltd. All rights reserved.

22

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

the small strain response only. One of few exceptions are the reverse shear experiments of Barlat et al. (2003) on 3 mm thick 1050-O aluminum sheets. Using wide shear specimens with a narrow gage section of reduced thickness, they achieved shear strains of up to 0.22 prior to loading direction reversal. Yoshida et al. (2002) presented an experimental study on the kinematic hardening response of sheet materials involving a finite strain compression phase. They bonded several flat specimens together and inserted the stack of specimens in an anti-buckling device during testing. Other examples of the use of antibuckling devices for testing sheet materials under in-plane compression can be found in Dietrich and Turski (1978), Kuwabara (1995), Yoshida et al. (2002), Boger et al. (2005), Cao et al. (2009) and Beese and Mohr (2011). The large strain compressionetension experiments by Yoshida et al. (2002) show that DP steels feature a Bauschinger effect, transient behavior, permanent softening and work hardening stagnation. Recall that the Bauschinger effect corresponds to an early yield after load reversal (Fig. 1b and d), transient behavior corresponds to a high hardening rate in the elastoeplastic transition regime resulting from load reversal (Fig. 1d); permanent softening prevails when the stress level after loading reversal remains below that for monotonic loading for the same equivalent plastic strain (Fig. 1b); work hardening stagnation causes a significantly reduced hardening rate after the transient hardening regime (Fig. 1c). Detailed reviews of kinematic hardening models are found in Chaboche (2008) and Eggertsen and Mattiasson (2009, 2010, 2011). Prager (1956) type of kinematic hardening, also referred to as linear kinematic hardening, describes both the Bauschinger effect and permanent softening. The main shortcoming of this model is the intrinsic coupling of both effects, i.e. a material exhibiting a Bauschinger effect without any permanent softening cannot be described with Prager's model. Furthermore, it describes neither transient behavior nor work hardening stagnation. Also, this type of hardening is unbounded and results in a persistent and often unrealistic rate of hardening at large strains. The Amstrong and Frederick (1966) kinematic hardening rule, also referred to as non-linear kinematic hardening model, describes the Bauschinger effect and transient behavior. The governing differential equation includes a recall term which activates the so-called dynamic-recovery.

Fig. 1. Comparison of the von Mises stress versus equivalent plastic strain curve for monotonic loading with that after loading reversal at a strain of 0.1 to illustrate different hardening model assumptions: (a) isotropic hardening, (b) permanent softening, (c) work hardening stagnation, (d) transient behavior.

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

23

The recall term is co-linear to the back stress tensor and is proportional to the increment in equivalent plastic strain. As a result, the evolution of the back stress is no longer linear unbounded and converges towards a saturation value under monotonic loading. Two parameters are used, one to control the Bauschinger effect and one for the transient behavior. However, the ArmstrongeFrederick model describes neither permanent softening nor work hardening stagnation. For improved approximations, several non-linear kinematics hardening models can be added with different recall constants characterizing the back stress evolution (Chaboche et al., 1979; Chaboche and Rousselier, 1983). These models give good predictions in the case of cyclic loading in the range of small strains, as they are able to describe the Bauschinger effect with great accuracy. The special case of coupling linear kinematic hardening with non-linear kinematic hardening provides good predictions in the case of moderate and large strains as it describes the permanent softening behavior during reverse deformation, especially with advanced high strength steels free from work hardening stagnation (Yoshida et al., 2002). However, the permanent softening effect is only represented through the linear kinematic term. As a result, an increase of the amount of permanent softening in the model always results in an increase in the strain hardening at large strains. Mroz (1967) proposed a multi-surface model framework to describe strain hardening. This idea was developed further by Dafalias and Popov (1976) by making use of a bounding surface in addition to the yield surface, with the distance between these two evolving surfaces defining the rate of strain hardening. Chaboche (2008) argues that a model featuring a combination of linear and non-linear kinematic hardening terms in addition to isotropic hardening can replicate the performance of a Dafalias-Popov type of model. However, one advantage of the Dafalias-Popov type of formulations is that the material response to monotonic loading can be identified independently from its response to reverse loading. The Dafalias-Popov model has been developed further by Geng and Wagoner (2000) to account for permanent softening. Yoshida and Uemori (2002) enriched the model even further and incorporated work-hardening stagnation. To the best of the authors' knowledge, experimental results on the effect of loading direction reversal on the strain to fracture are scarce and only found for bulk materials in the open literature. Bao and Treitler (2004) performed reverse loading experiments on notched axisymmetric bar aluminum 2024-T351 specimens with compression followed by tension all the way to fracture. They observed a substantial increase in ductility due to pre-compression. Papasidero et al. (2015) made use of a biaxial testing machine to subject tubular fracture specimens to non-proportional loading. Their results also demonstrated that pre-compressing aluminum 2024-T351 increases the strain to fracture for subsequent loading at higher stress triaxialities. The present paper investigates the effect of loading direction reversal on the strain to fracture in dual phase steel sheets. Section 2 presents an experimental procedure for the large strain pre-compression of 1 mm thick notched specimens prior to fracture testing. A finite strain plasticity model is presented in Section 3 combining elements of a non-associated plasticity model for advanced high strength steels (Mohr et al., 2010) with the non-linear kinematic hardening models of Chaboche (2008) and a Yoshida and Uemori (2002) type of work hardening stagnation approach. The plasticity model is then used in Section 4 to simulate all fracture experiments, before characterizing the effect of loading reversal on the strain at fracture initiation in Section 5. In the latter, a Hosford-Coulomb fracture initiation model with a non-linear damage accumulation rule is used to model the observed loading path effect.

2. Experiments 2.1. Material and specimens All experiments are performed on specimens extracted from 1.06 mm thick DP780 dual-phase steel sheets provided by US Steel. Under uniaxial tension, the material exhibits an initial yield stress of about 450 MPa and an ultimate strength of about 800 MPa. The axial stressestrain response for uniaxial tension is approximately isotropic (since the red, black and blue curves lie almost on top of each other in Fig. 2a), while the Lankford ratios are mildly loading direction-dependent: r0 ¼ 0.75, r45 ¼ 0.77, and r90 ¼ 0.95. Stocky dog-bone shaped specimens (UTC-specimens, Fig. 2b) are designed for uniaxial tension and compression testing. A gage section width of 10 mm is chosen to allow for the use of an anti-bucking device with a 5 mm wide central window for DIC strain measurement (Fig. 3c). As far as the choice of the gage section length (14 mm) is concerned, a compromise is sought between two competing effects: firstly, a large gage length-to-width ratio is desired to guarantee the validity of the assumption of uniaxial stress fields; secondly, a short gage section length is preferred to delay out-of-plane buckling under compressive loading. To investigate the plastic material response at very large strains, flat notched specimens with a 20 mm notch radius (NCTspecimens, Fig. 3d) are employed to subject the material to a compressionetension loading sequence. Notches guarantee that the localized neck will form at the specimen center perpendicular to the loading axis. This facilitates the numerical simulation of the experiments (as compared to conventional uniaxial tension specimens where the position of the through-thickness neck is not known a priori). All specimens feature 50 mm wide and 10 mm long shoulder areas. The specimens are extracted from sheet metal using water-jet cutting with the specimen axis aligned with the sheet rolling direction. Note that the edge quality of a waterjet cut is sufficient for the present experiments since fracture always initiates at the specimen center, i.e. away from the cut edges.

24

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

Fig. 2. (a) Uniaxial stressestrain response of DP780 steel for different directions under uniaxial tension; (b) Specimen geometry for uniaxial compressionetension experiments.

2.2. Experimental procedure A hydraulic testing machine (Instron, Model 8802) is used to perform all experiments. Custom-made high pressure clamps are employed to attach the specimen to the testing frame (Fig. 3a and b). Unlike conventional wedge grips, the clamps work equally well under compression and tension. Accurate alignment of the upper and lower specimen grips is critically important to delay buckling (we recommend a tolerance of less than 0.1 mm for the parallelism of the top and bottom clamping surfaces). Furthermore, it is important to tighten the specimen clamps under active force control to avoid any plastic pre-compression due to the Poisson effect during clamping of the specimen. Fig. 3a shows a photograph of the assembly of the floating anti-buckling device. The design is similar to that proposed by Beese and Mohr (2011) except for a change in type and number of springs and bolts to apply an increased average lateral pressure of about 3 MPa. Thin Teflon sheets are placed between the specimen and the device to minimize friction. Note that the applied lateral stress is very small (about 3 MPa) as compared to the axial stress in the specimen (about 800 MPa); its effect on the material response is thus neglected when processing the experimental results. The extensometer function of the digital image correlation software VIC2D (Correlated Solutions) is used to determine the relative vertical displacement of two points positioned on the longitudinal specimen axis at an initial distance of 12.5 mm and 10 mm for the uniaxial and notched tension specimens, respectively (as highlighted by blue dots in Figs. 2b and 3d). For the notched specimens, the average axial surface strain is also determined in the area of localized necking using a DIC extensometer of 2 mm length (red dots in Fig. 3d).

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

25

Fig. 3. Experimental set-up for compressionetension testing: (a) front view and (b) side view of the specimen with anti-buckling device and high pressure clamps, (c) displacement measurement using Digital Image Correlation (DIC), and (d) notched specimen (NCT) for compressionetension fracture testing.

2.3. Experimental results The experimental program includes:  Tension followed by compression on UTC-specimens for tensile pre-strains of 0.05 and 0.10 (axial engineering strains), as well as  compression followed by tension on NCT-specimens for compressive pre-strains of 0.026, 0.064, 0.091 and 0.132 (axial engineering strain at specimen center as measured with the 2 mm long virtual surface extensometer). All experiments are performed at a cross-head velocity of about 1 mm/min. Monotonic tension experiments with and without anti-buckling device yielded identical results. This partially confirms that the effect of the lateral friction due to the anti-buckling device is negligible. True compressive strains, as high as 0.2, could be achieved without noticeable buckling. In Fig. 2a, the true stressestrain curve for monotonic compression along the rolling direction (green curve) is shown next to that for tension in Fig. 2a. The curve is only shown for strains of up to 0.10 as the effect of barreling makes the assumption of uniaxial stress fields invalid at larger compressive strains. A summary of all measured forceedisplacement curves for experiments performed with loading direction reversal are shown in Fig. 4. For both stocky dogbone specimens (Fig. 4a) and notched specimens (Fig. 4b), all curves lie on top of each other during the initial phase of loading (tension phase for the UTC experiments, compression phase for the NCT experiments) which demonstrates the repeatability of the experimental procedure. 3. Combined Chaboche-Yoshida (CCY) plasticity model A new plasticity model is presented combining elements of the non-associated plasticity model for advanced high strength steels of Mohr et al. (2010), the non-linear kinematic hardening models of Chaboche (2008) and a Yoshida and

26

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

Fig. 4. Experimental results: (a) engineering stressestrain curves as obtained from uniaxial tension-compression experiments, (b) forceedisplacement curves as obtained from notched compressionetension experiments.

Uemori (2002) type of work-hardening stagnation approach. The model is embedded into the standard finite strain framework of the commercial finite element software Abaqus/explicit (Abaqus, 2012). 3.1. Yield surface The center of the yield surface in the deviatoric Cauchy stress space is described through the back stress tensor X. The tensor x is introduced to describe the relative position of a point in stress space with respect to the yield surface center

x ¼ dev½s  X:

(1)

Note from Fig. 2a that the stressestrain curves for uniaxial tension along the rolling, transverse and diagonal directions lie approximately on top of each other. We therefore use the isotropic von Mises equivalent stress measure to define the yield surface

f ðs; kÞ ¼ xVM  k ¼ 0; with

xVM

rffiffiffiffiffiffiffiffiffiffiffiffi 3 x : x: :¼ 2

(2)

(3)

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

27

Note that the above yield surface corresponds to the isotropic von Mises yield surface if X ¼ 0, whereas it is anisotropic otherwise. 3.2. Non-associated flow rule Despite the isotropic stressestrain response under uniaxial tension, the measurements of the Lankford coefficients indicate some in-plane anisotropy in the material. An anisotropic non-associated flow rule (see Stoughton (2002), Cvitanic et al. (2008), Mohr et al. (2010)) is therefore chosen to describe the evolution of the plastic strain tensor. It is assumed that the increment in plastic strains dεp is aligned with the derivative of the Hill’48 potential function in the x-space,

dεp ¼ dl

vxHill vx

(4)

with

xHill ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x211 þ G22 x222 þ ð1 þ 2G12 þ G22 Þx233 þ 2G12 x11 x22  2ð1 þ G12 Þx11 x33  2ðG22 þ G12 Þx22 x33 þ G33 x212 þ 3x213 þ 3x223 (5)

In (4), dl  0 is the plastic multiplier, while the coefficients G12, G22 and G33 are directly linked to the Lankford ratios,

G12 ¼ 

r0 r 1 þ r90 1 þ 2r45 r0 þ r90 ; G22 ¼ 0 ; and G33 ¼ : 1 þ r0 r90 1 þ r0 r90 1 þ r0

(6)

For ease of notation, we rewrite Hill's equivalent stress (5) as quadratic form:

xHill ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x:G:x

(7)

with the fourth-order tensor G. 3.3. Definition of the equivalent plastic strain The equivalent plastic strain increment is defined as work-conjugate to the von Mises equivalent stress in x-space,

dεp ¼

1 x : dεp : xVM

(8)

Combining Eqs. (4) and (8), the relationship between the equivalent plastic strain and the plastic multiplier is obtained,

dεp ¼

xHill dl: xVM

(9)

Note that in the absence of a back stress, the above definition of the equivalent plastic strain is the same as that proposed in Stoughton (2002), Mohr et al. (2010) and Mohr and Marcadet (2015). 3.4. Isotropic hardening The isotropic hardening law describes the evolution of the deformation resistance k during plastic loading. Formally, we write

  dk ¼ bH εp dεp ;

(10)

where dεp is the increment in the equivalent plastic strain. 0  b  1 is a multiplier associated with work hardening stagnation (to be detailed below), and H ¼ H½εp  is the isotropic hardening modulus. In Chaboche (2008), a Voce law is used for the isotropic hardening function. A limitation of such an evolution is its saturation at large strains. Here, we use instead a linear combination of the Voce and Swift hardening laws (Sung et al., 2010) to parameterize the isotropic hardening function,

   n1 Hε εp ¼ ð1  wÞQ tetεp þ wAn ε0 þ εp with the Voce parameters {Q,t}, the Swift parameters {A,ε0,n}, and the weighting factor 0  w  1.

(11)

28

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

3.5. Non-linear kinematic hardening The evolution law for the back stress is described through the sum of two non-linear kinematic hardening rules,

dX ¼ da1 þ da2 :

(12)

with the initial conditions ai(t ¼ 0) ¼ 0, X(t ¼ 0) ¼ X0. The corresponding evolution equations read

da1 ¼ g1



 2 G:x C1  a1 dl 3 xHill

(13)

da2 ¼ bg2

  2 G:x C2  a2 dl 3 xHill

(14)

and

with the work hardening stagnation multiplier b, and the material model parameters Ci and gi. Note that in the non-linear kinematic hardening law, the back stress evolution direction is aligned with the non-associated plastic flow. The two non-linear kinematic hardening rules serve different purposes. a1 is introduced to model the Bauschinger effect. According to the FrederickeArmstrong differential equation, the evolution of the back stress is no longer unbound under monotonic loading; it converges towards a saturation stress instead. The second term, a2, is introduced to model apparent softening. A linear kinematic hardening term is usually introduced which results in permanent softening. This has proven as a useful assumption when modeling the response of metals at moderate strains (Eggertsen and Mattiasson, 2011; Zang et al., 2011; Geng and Wagoner, 2000; Chun et al., 2002). However, according to our experimental observations, the softening effect appears to fade away for very large strains. In other words, instead of permanent softening, we introduce transient softening through a non-linear kinematic hardening rule. The key difference between transient softening and the Bauschinger effect is the strain scale, i.e. the Bauschinger effect fades away rapidly after a few percent of strain after loading reversal, whereas the strain scale associated with the transient softening effect is at least one order of magnitude higher. In the model, this is reflected in the choice of parameters (e.g. g1 z 10g2). By introducing b in the evolution of the second backstress term only, the fast recovery from the Bauschinger effect remains fully active during the work hardening stagnation phase. 3.6. Work hardening stagnation This part of our model is inspired by the work of Yoshida and Uemori (2002). Originally developed as a two surface model (Dafalias-Popov framework), we borrow some ideas from Yoshida and Uemori (2002) to define a constitutive equation for the work hardening stagnation multiplier b. The activation of work hardening stagnation depends on the loading history. To separate the effect of work hardening stagnation from the effect of the non-associated plastic flow, we characterize the loading path in terms of the strain-like path tensor

Zεp rffiffiffi 3 ndεp 2

p¼

(15)

0

With n denoting the normal to the yield surface. p is equal to the plastic strain tensor in the case of associated plastic flow. A sphere of radius r is defined around a central point q in a way that p lies always inside (Fig. 5d). Denoting the distance between q and p as d,

d¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðp  qÞ : ðp  qÞ;

(16)

we have

d  r  0:

(17)

The work hardening stagnation multiplier is then defined as

b :¼

d r

with 0  b  1. The position and size of the sphere are permitted to change according to the evolution equations

(18)

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

29

Fig. 5. Illustration of the work hardening stagnation model in the plastic strain space. The sequence shows (a) the initial values, (b) the evolution during monotonic loading, (c) the point of loading reversal, (d) the transient stagnation after loading reversal, and (e)e(f) the evolution after stagnation.

dq ¼ ð1  hÞdp dr ¼

h ðp  qÞ : dp d

(19) (20)

if b ¼ 1 (i.e. p is located on the sphere boundary, see Fig. 5b, e and f) and (p  q):dp > 0 (i.e. the loading direction dp points outwards). Note that (18) was chosen such that the consistency condition db ¼ 0 is readily fulfilled. The case b ¼ 1 with db ¼ 0 corresponds to loading with no work hardening stagnation (Fig. 5b), i.e. isotropic and kinematic hardening are both fully

30

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

active. b ¼ 0 represents the opposite limiting case of maximum work hardening stagnation. Different from the model proposed by Yoshida and Uemori (2002), we allow for partial work hardening stagnation, i.e. 0 < b < 1. 3.7. Summary of model parameters The proposed plasticity model includes sixteen material model parameters:  Three anisotropy coefficients {G12,G22,G33} defining the non-associated flow rule  Seven isotropic hardening parameters including the initial deformation resistance k0, the Swift parameters {ε0,A,n}, the Voce parameters {Q,t}, and the weighting factor w  the initial back stress tensor X0  the Bauschinger parameters C1 and g1  the softening parameters C2 and g2  for work hardening stagnation parameter h It is worth noting that the proposed Combined Chaboche-Yoshida (CCY) model reduces to a conventional Chaboche model when deactivating the work hardening stagnation option (h ¼ 0) and using G12 ¼ 0.5, G22 ¼ 1, G33 ¼ 3, a0 ¼ 0. 3.8. Thermodynamic constraints The starting point of our considerations is the free energy imbalance of the form

j_  s : ε_

(21)

In addition to an elastic part je, the free energy includes a plastic part jp associated with kinematic hardening,

j ¼ je þ jp ;

(22)

which both must be positive, i.e.

je  0;

(23a)

jp  0:

(23b)

Assuming the elastic strain energy potential

je ¼

   1  C : ε  εp : ε  εp ; 2

(24)

along with the elastic constitutive equation

  vj  ¼ C : ε  εp ; s¼  v ε  εp

(25)

the free energy imbalance may be substituted by the requirement of nonnegative rate of plastic dissipation:

d_ p ¼ s : ε_ p  j_ p  0:

(26)

Next, we assume a plastic free energy of quadratic form

jp ½a1 ; a2  ¼

3 a : ai 4gi Ci i

(27)

which readily satisfies the non-negative requirement (23b) for gi > 0 and Ci > 0. After application of the back stress evolution Equations (13) and (14), the rate of change in the plastic free energy reads

  3 3 j_ p ¼ a1 : ε_ p  a1 : a1 l_ þ b a2 : ε_ p  a2 : a2 l_ 2C1 2C2 Combining Eqs. (26) and (28), yields the rate of plastic dissipation

(28)

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

31

3l_ 3l_ d_ p ¼ ða1 þ a2 þ xÞ : ε_ p  j_ p ¼ fð1  bÞa2 þ xg : ε_ p þ a : a1 þ b a : a2 2C1 1 2C2 2 ¼ ð1  bÞ

a2 : G : x 3l_ 3l_ þ xHill þ a :a þb a :a 2C1 1 1 2C2 2 2 xHill |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

(29)

0

The second term is unconditionally nonnegative. For loading outside the work hardening stagnation regime (b ¼ 1), the thermodynamic constraints are readily satisfied even though a non-associated flow rule is used. A constraint needs to be imposed on the material model parameters if b < 1 and a2:G:x < 0. In that case, the non-zero dissipation condition is satisfied if 2

a2 : G : x  xHill

(30)

To satisfy (30), it is sufficient to impose a constraint on the magnitude of the back stress tensor a2,

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 : G : a2  xHill :

(31)

According to the Frederik-Armstrong law, the evolution of a2 is bound to

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 : G : a2  C2 lmax;G 3

(32)

with lmax,G denoting the largest eigenvalue of G,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lmax;G ¼ max 3; G44 ; 1 þ G12 þ G22 þ 1 þ 4G212 þ 2G12 G22 þ 2G12 þ G222  G22 :

(33)

The free energy imbalance is thus satisfied if the model parameters respect the constraints

gi > 0; Ci > 0; and

2 C l  k0 : 3 2 max;G

(34)

4. Plasticity model identification and validation The plasticity model parameters are identified in a two-step procedure: 1. A first set of parameters is determined using the experimental data of the UTC experiments whose domain of validity is limited by in-plane barreling and out-of-plane buckling at large compressive strains. 2. Subsequently, these parameters are used as seed values for an inverse parameter identification method based on experimental data for the NCT experiments in the pre- and post-necking range. The inverse procedure involves finite element simulations of all experiments on NCT specimens up to the point of fracture and quantifies the difference between the simulation and experimental results. An improved second set of parameters is then obtained from the computational minimization of this difference. 4.1. Identification step I: determination of seed parameters The anisotropy coefficients G12, G22 and G33 of the Hill’48 flow potential are determined from the measured Lankford coefficients r0, r45, r90. Application of Eq. (6) yields G12 ¼ 0.43, G22 ¼ 0.88 and G33 ¼ 2.59. According to the considerations made in Beese and Mohr (2011), the initial back stress for the cold-rolled material is assumed to take the special form

X ¼ X0 ðe1 5e1  e3 5e3 Þ

(35)

where the 1-direction corresponds to the rolling direction, and the 3-direction to the through thickness direction. The initial deformation resistance k0 and the initial back stress X0 are determined from the tension/compression asymmetry of the material response (Fig. 2a). Denoting the absolute values of the initial yield stresses under tension and compression as Yt and Yt, respectively, we have

k0 ¼

1 1 ðYt þ Yc Þ and X0 ¼ ðYt  Yc Þ: 2 3

At 0.2% plastic proof strain, we have Yt y 450 MPa and Yc y 510 MPa, and thus, k0¼480 MPa and X0¼20 MPa.

(36)

32

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

The hardening parameters a ¼ {A,n,ε0,Q,t,C1,g1,C2,g2,h} are identified through optimization. Simulations for pure uniaxial tension followed by compression are performed on a single element. The true stress versus logarithmic strain curve is computed for each experiment and compared with the corresponding experimental result. The cost function is expressed as

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u   !2 uX M SIM εEXP u 4 Xj sj m  1 GI ½a ¼ t sEXP εEXP m j j¼1 m¼1

(37)

where the subscript j ¼ 1,2,..,4 differentiates among the stressestrain curves for different levels of pre-strain; Mj denotes the total of experimental data points used for the computation of the residual for the experiment j. The seed model parameters listed in Table 1 are obtained from minimizing G[a]I. 4.2. Identification step II: full inverse parameter identification 4.2.1. Finite element model A finite element model is built to simulate all NCT experiments. Special attention is paid to the modeling of the postnecking response. Assuming a symmetric mechanical system, one eighth of the entire specimen is modeled using firstorder solid elements with reduced integration (element type C3D8R of the Abaqus element library). The mesh features a total of 34,488 elements with eight elements along the half-thickness of the specimen. The size of the mesh is chosen such that the predicted equivalent plastic strain at fracture initiation converged (less than 2% change upon further mesh refinement). The computed axial displacement is reported for the gage section point that corresponds to the position of the DIC extensometer in the experiments. All simulations are performed under displacement control with the displacement prescribed at the top boundary of the specimen shoulder (of the numerical specimen). The explicit solver of the FE software Abaqus is used to solve the computational boundary value problem using at least 100,000 time steps for each simulation run. 4.2.2. FEA-based optimization of model parameters The set of seed parameters obtained in step I may not accurately capture the behavior of the material at large strains. Buckling limits the range of strains after loading reversal for which the parameters are identified. A strong asymmetry of the material behavior in tension and compression after loading reversal could also affect the relevance of the set of parameters: the seed parameters are obtained from tension followed by compression tests, whereas the fracture experiments feature compression followed by tension. An improved set of parameters is computed through non-linear unconstrained optimization. Each iteration includes the following steps: 1. Definition of the material model input data card for a given set of parameters a. 2. Simulation of all five NCT experiments (pre-strain levels 0.0, 0.026, 0.064, 0.091 and 0.132) and extraction of the corresponding force versus DIC extensometer displacement curves, FjFEA ðuÞ, with the subscript j ¼ 1,2,..,5 differentiating among the results for different pre-strain levels. 3. Evaluation of the global cost function

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X GII ½a ¼ t D2 N j¼1 j

(38)

with N denoting the number of experiments included in the optimization. For each experiment, a cost function Dj is defined

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j2 j2 j2 D1 þ D2 þ D3 Dj ¼ 3

(39)

The cost function is a measure of how well the model performs with respect to three criteria:

Table 1 Plasticity model parameters for the Combined Chaboche-Yoshida (CCY) model and the Chaboche model.

CCY Chaboche

Seed Final Seed Final

k0

X0

A

n

ε0

Q

t

w

C1

g1

C2

g2

MPa

MPa

MPa

e

e

MPa

e

e

MPa

e

MPa

e

e

480. 512. 480. 495.

30. 30. 30. 30.

946.5 997.2 707.2 743.9

0.15 0.19 0.22 0.11

0.023 0.031 0.013 0.019

164.3 143.7 41.8 63.2

22.1 32.2 29.4 40.4

0.5 0.5 0.5 0.5

204.6 235.4 322.7 301.2

51.3 50.8 51.2 48.5

108.9 92.1 497.3 453.3

6.1 5.9 1.1 1.0

0.7 0.6 0. 0.

h

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

33

(1) The first criterion D1 evaluates the relative error in the value of the predicted maximum load.

j D1





 j j max FEXP  max FFEA

 ¼ j max FEXP

(40)

(2) The second criterion D2 evaluates the relative error in the overall predicted drop of the load from its maximum value to its value at fracture.

   

 

j j j max F j  F  F  max F EXP EXP uj ¼uf FEA FEA uj ¼uf EXP EXP EXP FEA   Dj2 ¼

 j j max FEXP  FEXP f

(41)

uEXP ¼uEXP

(3) The third criterion D3 evaluates the relative error in the value of the displacement at predicted maximum load.

Dj3 ¼

j u EXP F j

EXP

j ¼maxðFEXP Þ



ujFEA

j j FFEA ¼maxðFFEA

ujEXP j j FEXP ¼maxðFEXP Þ

Þ

(42)

This set of criteria has been constructed in an attempt to obtain an accurate prediction of the macroscopic specimen response in the post necking range. The optimal set of parameters aopt is then determined through the minimization of the cost function,

aopt ¼ argminGII ½a:

(43)

The minimization is performed using a derivative-free Nelder-Mead algorithm as implemented in the software Matlab (v.7). The final set of parameters as obtained after 84 iterations is given below the seed values in Table 1. In the material chosen for this study, the set of seed parameters already provided satisfying accuracy, and the optimization of the seed parameters provided only minor improvements. However, this method was found necessary for other materials featuring early buckling or pronounced tension-compression asymmetry. 4.3. Model verification for reverse loading The comparisons of the measured and predicted stressestrain curves for uniaxial tension-compression (UTC) are shown in Fig. 6a. The remainder plots in Fig. 6 show the forceedisplacement curves extracted from finite element simulations (solid lines) of the NCT specimens next to the experimental measurements (solid dots). The corresponding surface strain (axial strain at the specimen center as averaged over a distance of 2 mm) evolutions are shown on the blue secondary axis. The proposed model predicts accurately the loadedisplacement relation for all experiments. The magnitude of the maximum load as well as the corresponding displacement are well captured. The post necking behavior is well predicted as the load drop corresponds to the experimental measurement for all levels of pre-compression. A comparison of the computed and measured surface strain evolutions shows that the proposed plasticity model is able to provide reasonable estimates of the surface strains in the area of localized necking for all levels of pre-strain. However, the surface strain history plots are more sensitive to model inaccuracies. For example, the simulation for loading reversal after a pre-strain of 0.09 (Fig. 6e) overestimates the surface strain at the instant of fracture by as much as 14% (0.41 versus 0.36 in the experiment). The predicted forceedisplacement curve on the other hand agrees well with the experiment. With the goal of identifying the loading paths to fracture for monotonic and reverse loading, the model performance has only been assessed for these loading conditions in the context of the present work. The reader is referred to the work by Yoshida and co-workers for a more comprehensive validation of the main model ingredients. For reference, we also repeated the entire calibration and validation procedure with the Chaboche model, i.e. for the CCYmodel with no work hardening stagnation and associated plastic flow. The corresponding results are shown as dashed curves in Fig. 6. The Chaboche model also provides good predictions of the macroscopic forceedisplacement curves, but it tends to overestimate the strains inside the neck. Remark. Thickness profiles can also be used for “local” validation. We had prepared these for specimens extracted from a different batch of sheets that had been tested much earlier (and hence featured slightly different mechanical properties). We also had calibrated both the CCY and Chaboche model for those sheets which resulted in similar agreement of the simulated and measured forceedisplacement curves as for the material above. Selected experimental (solid dots) and simulated

34

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

Fig. 6. Comparison of CCY model (solid lines) and experiments (solid dots) for (a) uniaxial tension-compression (UTC), and (b)e(f) notched compressionetension (NCT); the Chaboche model predictions are depicted as dashed lines.

thickness profiles (solid lines) are shown in Fig. 7. As for the surface strain measurements, the simulation results agree reasonably well with the experimental measurements. In particular, good agreement is observed for the ratio of the strains outside and inside the neck (e.g. 0.55 outside the neck vs 0.35 inside the neck in the case of monotonic tension). Unfortunately, we could not repeat the thickness profiles measurements on the current batch of specimens due to the shortage of material. 4.4. Model verification for monotonic multi-axial loading The stress state sensitivity of the fracture initiation on DP780 steel sheets had been determined in Mohr and Marcadet (2015) using  a flat 10 mm wide specimen with a lateral notch (R ¼ 6.67 mm) for plane strain tension (Dunand and Mohr, 2010);  a flat 20 mm wide specimen with a central hole (R ¼ 20 mm) for uniaxial tension (Dunand and Mohr, 2010);  a butterfly specimen (Dunand and Mohr, 2011) with reduced thickness gage section for shear;

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

35

Here, all experiments are repeated on the specimens extracted from the new batch of sheets closely following the experimental procedures outlined in Mohr and Marcadet (2015). The monotonic experiments are simulated using the CCY plasticity model as calibrated based on the reverse-loading experiments. The comparison of the measured forceedisplacement curves (dotted curve) with the simulation results (solid curves) shown in Fig. 8 confirms the model's ability to provide good predictions of the overall specimen responses. For reference, we also performed these simulations with the Chaboche

Fig. 7. Thickness profiles at the instant of inset of fracture as measured experimentally (solid dots) and extracted from numerical simulation (CCY ¼ solid lines, Chaboche ¼ dashed line) after pre-compression up to a surface strain of (a) 0. (b) 0.08, and (c) 0.12; note that the results are shown for specimens extracted from a different batch of DP780 sheets.

36

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

model which yields less accurate predictions. In particular, the force decrease for notched tension is predicted too early. As a result, the strain at the specimen center at the instant of onset of fracture is overestimated. 5. Effect of loading reversal on ductile fracture initiation The effect of pre-compression is investigated at the specimen level to gain insight into the different ductility of material points that undergo monotonic stretching as compared to those that are subject to compressionetension cycles when a sheet is drawn over a radius or when hinge lines propagate during the folding of box structures. The uniform state of the sheet after production (rolling, etc.) is therefore considered as the reference state of zero strain and the strain introduced thereafter up to the point of failure is reported as fracture strain. 5.1. Characterization of the stress state The plasticity model is anisotropic due to the Hill’48 flow potential and the non-zero back stress. However, for the sake of simplicity and due to the lack of experimental data for different material orientations, it is assumed that the fracture response is not affected by the orientation of the stress tensor with respect to the material coordinate system. The stress state is thus characterized through the stress triaxiality and the Lode angle parameter (which are both isotropic measures). The stress triaxiality h is proportional to the ratio of the first invariant of the Cauchy stress tensor, I1, and the second invariant of the deviatoric stress tensor, J2,

sm h¼ sVM

h i tr s with sm ¼

3

I ¼ 1 3

and sVM

pffiffiffiffiffiffiffi ¼ 3J2 ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 dev½s : dev½s 2

(44)

The dimensionless Lode angle parameter, q, measures the ratio of the second and third invariants of the deviatoric stress tensor, J2 and J3. Its mathematical definition reads

2 3 pffiffiffi 2 63 3 J3 7 qffiffiffiffi5 q ¼ 1  arccos4 p 2 J 32

(45)

J3 :¼ detðdev½sÞ:

(46)

with

5.2. Effect of pre-compression on results for notched tension Fig. 9b shows the computed equivalent plastic strain distribution inside the neck of the notched tensile specimens at the instant of fracture initiation, while Fig. 9c depicts the corresponding thickness profiles. The red dots on the specimen surface in Fig. 9b highlight the position of the 2 mm DIC surface extensometer. In a first approximation, both the average axial strain and the thickness reduction at the instant of fracture are unaffected by the amount of pre-compression. The axial (engineering) strain measured by the DIC surface extensometer is about 0.35 and the thickness reduction is about 30%. However, the shape of the neck and the strain distribution inside the neck are strongly influenced by the amount of pre-compression which results in an increase of the local equivalent plastic strain to fracture as a function of the pre-strain (Fig. 9a), ranging from εf ¼ 0:57 for monotonic tensile loading to εf ¼ 0:77 for fracture after pre-compression to a local equivalent plastic strain of 0.13. In the initial compressive phase of loading, the thickness of the specimen is increased. At the same time, the isotropic hardening capacity in the tensile loading phase decreases as a function of the pre-strain; as a result, the localization of deformation inside the neck is more pronounced which causes a steeper thickness gradient in the vicinity of smallest crosssection (compare the slopes of the thickness profiles in Fig. 9c). A more detailed description of the evolution of the thickness and state variables at the location of fracture initiation in the NCT-13 experiment is given in Fig. 10. In addition to the forceedisplacement curve (Fig. 10a),  Fig. 10b shows the evolutions of the deformation resistance k (solid line), the back stress tensor component {a1}11 (dashed line), and the back stress tensor component {a2}11 (dotted line);  Fig. 10c shows the evolution of the von Mises equivalent stress;  Fig. 10d shows the evolution of the longitudinal thickness profile;  Fig. 10e shows the evolution of the equivalent plastic strain as a function of the stress triaxiality (so-called “loading path to fracture”);

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

37

Fig. 8. Results from monotonic experiments: (a) notched tension (NT6), (b) tension with a central hole (CH), and (c) butterfly specimen (SH); solid dots ¼ experiments, solid lines ¼ CCY model, dashed lines ¼ Chaboche model; the contour plots show the distribution of the equivalent plastic strain at the instant of fracture initiation.

Stagnation (b < 1) takes place between points ① and ③ and reaches its maximum (b ¼ 0) at point ②. At this point, both isotropic hardening and transient softening are inactive, while the kinematic hardening term associated with the Bauschinger effect is the only active hardening mechanism. At the same time, through-thickness necking initiates during the phase of hardening stagnation (compare the thickness profiles ② and ③ in Fig. 10d). Full hardening resumes beyond point ③ (Fig. 10c). It is well known in the literature that the strain hardening capacity affects the geometry of the neck (e.g. Pardoen et al., 2004). Here, the strain hardening capacity is affected in several ways by the pre-compression: the total isotropic hardening capacity decreases as a function of the pre-strain. However, this effect is competing with the hardening capacity in the Bauschinger transition, modelled by the kinematic hardening. During the stagnation phase, the rate of strain hardening is very low (between point ② and ③ in Fig. 10c), leading to a more localized evolution of the neck. The consequence on the evolution of triaxiality is observed in Fig. 10e, where there is a rapid increase of the triaxiality as a function of the plastic strain between point ② and ③. Following the stagnation phase (after point ③), there is a short phase of significant increase of the rate of

38

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

Fig. 9. (a) Strain to fracture after compressionetension loading of NCT specimens as a function of the equivalent plastic strain at the point of loading direction reversal; the star symbols present the hybrid experimental-numerical results, the solid dots care model predictions; (b) Similar plot as (a) replacing the total strain to fracture by net fracture strain (c) Equiv. plastic strain distribution in the longitudinal specimen cross-section, and (d) thickness distribution at the instant of fracture initiation.

hardening (Fig. 10c), leading to a more diffused evolution of the neck. In Fig. 10e, there is temporarily almost no evolution of the triaxiality with respect to the plastic strain after point ③. 5.3. Hosford-Coulomb fracture initiation model Assuming that the onset of fracture is imminent with the onset of localization at the microscale, the Hosford-Coulomb model has been developed based on the results from 3D unit cell computations for proportional loading (Dunand and Mohr, 2014). In particular, it is postulated that ductile fracture initiates after proportional loading when the linear combination of the Hosford equivalent stress and the normal stress acting on the plane of maximum shear exceeds a critical value (Mohr and Marcadet, 2015),

sHF þ cðsI þ sIII Þ ¼ b0

(47)

with

sHF ¼

 1 ðsI  sII Þa þ ðsII  sIII Þa þ ðsI  sIII Þa 2

1 a

:

(48)

After transformation of the Hosford-Coulomb criterion into mixed strain-stress state space, the strain to fracture for proportional loading reads,

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

39

Fig. 10. Detailed analysis of the NCT-13 results: (a) forceedisplacement curve and local surface strain; Evolutions of (b) the CCY hardening terms: isotropic ð1Þ ð2Þ hardening resistance k (solid line), axial component of the back stress tensors aNL (dashed line) and aNL (dotted line), (c) the von Mises stress, (d) the thickness profiles, and (e) the stress triaxiality; the dashed line in (e) depicts the Chaboche model prediction; the labels ①, ②, ③ and ④ indicate the point of load reversal, onset of stagnation (also maximum load in tension), end of stagnation, and onset of fracture.

  εpr h; q f

11p 0

1   a 1 a a a ðf  f2 Þ þ ðf2  f3 Þ þ ðf1  f3 Þ ¼ bð1 þ cÞ @ þ cð2h þ f1 þ f3 ÞA 2 1 1 p

(49)

with the Lode angle parameter dependent trigonometric functions

hp  i   2 f1 q ¼ cos 1q ; 3 6

hp  hp  i i   2   2 f2 q ¼ cos 3 þ q and f3 q ¼  cos 1þq : 3 6 3 6

(50)

with the transformation parameter p ¼ 0.1 (Roth and Mohr, 2014), and the fracture model parameters {a,b,c}. Note that the Hosford exponent a controls the effect of the Lode angle parameter, while the friction coefficient c primarily controls the effect of the stress triaxiality on the strain to fracture. The model parameter b is a multiplier controlling the overall magnitude of the strain to fracture. It is defined such that it is equal to the strain to fracture for uniaxial tension (which is the same as that for equi-biaxial tension).

40

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

To predict the onset of fracture after non-proportional and non-monotonic loading, the above criterion is embedded into a damage indicator model framework. Let D 2 [0,1] denote a scalar damage indicator, with the initial value D ¼ 0 for the undeformed material, and D ¼ 1 for the deformed material at the instant of fracture initiation. The evolution of the damage indicator is then related to the evolution of the equivalent plastic strain using a stress state dependent non-linear damage accumulation rule (Papasidero et al., 2015),

dD ¼ m

εp   εpr h; q f

!m1

dεp  : εpr h; q f

(51)

Irrespective of the choice of the damage accumulation exponent m > 0, the condition D ¼ 1 is fully equivalent to the direct application of (49) for proportional loading. In the case of non-proportional loading, values of m < 1 put more weight on the effect of the stress state at the early stage of loading, whereas values of m > 1 emphasize the effect of the stress state right before fracture initiation. In the case of m ¼ 1, the so-called linear damage accumulation rule is retrieved (e.g. Bai and Wierzbicki, 2010). 5.4. Model identification and verification For each experiment performed, we extract the loading path to fracture, εp ¼ εp ½h; q, at the location within the specimen, where the highest equivalent plastic strain is achieved in the corresponding numerical simulation. The solid lines in Fig. 11a depict the loading paths to fracture in the strain versus stress triaxiality plane for eight different experiments. The end of each path corresponds to the instant of fracture initiation.

Fig. 11. (a) Loading paths to fracture as extracted from finite element simulations of all fracture experiments up to the instant of fracture initiation (end point of solid lines); the Hosford-Coulomb fracture initiation model predictions are shown as solid dots; (b) strain to fracture for proportional loading as a function of the Lode angle parameter and the stress triaxiality; (c) evolution of the damage indicator at the location of fracture initiation in compressionetension experiments.

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

41

The four fracture initiation model parameters {a,b,c,m} are identified based on the loading paths using a gradient free inverse identification algorithm (Nelder-Mead minimization in Matlab). The results for the shear (SH), central hole tension (CH) and the monotonically loaded NCT-specimen (NCT-0) are included in the calibration procedure for {a,b,c} to cover a wide range of stress states; in addition, the reverse loading experiment NCT-9 is included in the data basis to identify the damage accumulation exponent m. After launching the identification procedure with the seed values {1.5,0.7,0.01,0.8}, the “optimal” parameters a ¼ 1.65, b ¼ 0.64, c ¼ 0.05 and m ¼ 0.45 are obtained after 76 iterations. Fig. 11b shows a 3D plot of the identified Hosford-Coulomb criterion εpr ¼ εpr ½h; q for proportional loading. The characteristic strains to fracture of pure shear, uniaxial f f tension and plane strain tension are 0.81, 0.64 and 0.51. The resulting model predictions of the instants of fracture initiation are shown as solid dots in Fig. 11a. The solids dots lie exactly on top of the ends of the loading paths for the calibration experiments (black lines) which indicates that the model has sufficient mathematical flexibility to be fitted to the experimental data. The blue dots predict the instants of fracture for the four experiments that have not been included in the calibration procedure. As for the calibration experiments, the blue dots (model predictions) lie approximately on top of the ends of the blue solid lines (hybrid experimental-numerical data) which is seen as a partial validation of the proposed phenomenological model. Fig. 9a allows for a more detailed comparison of the model predictions (solid dots) with corresponding experimental results (star symbols), revealing a maximum relative error of about 5% for the NCT-13 experiment. The evolution of the damage indicator for the compressionetension experiments is shown in Fig. 11c. The curves all show an abrupt increase in the damage accumulation rate dD=dεp at the instant of loading reversal. This is due to the jump in stress triaxiality (and Lode parameter) from 0.39 (0.8) to þ0.39 (þ0.8). The corresponding strains to fracture for proportional loading at these stress states are 1.43 and 0.60, respectively. According to Eq. (51), this jump implies an increase in the rate of damage accumulation, i.e. at a given equivalent plastic strain, the damage accumulation is much lower under compression than under tension. 5.5. Discussion of the effect of pre-strain on ductile fracture The basic mechanism responsible for the apparent ductility increase is that void like defects (which trigger the shear localization at the microscale) do not nucleate or grow under compression. Hence the main tendency of an increase of the strain to fracture as a function of the pre-compression strain. However, even after computing the net fracture strain (i.e. subtracting the pre-compression strain from the final fracture strain, Table 2 and Fig. 9b), we observe an increase in ductility due to pre-compression. This second order effect is attributed to the local thickening of the specimen in pre-compression phase, thereby delaying the necking related stress triaxiality increases after loading reversal. From a phenomenological point of view, it is worth noting that the present results suggest that the stress state at the beginning of the loading history has an important effect on the final strain to fracture. This becomes apparent when comparing the results obtained with the proposed nonlinear damage accumulation rule (m ¼ 0.45) with those obtained using a linear damage accumulation rule (m ¼ 1). A slight increase of the ductility in terms of equivalent plastic strain at fracture is predicted using the linear rule (open square dots in Fig. 9a), but it is well below the prediction of the non-linear rule (star symbols in Fig. 9a). An alternative modeling approach using the linear rule would be to consider the material after pre-strain as a new virgin material with zero initial value of the damage indicator. Even though the predictions with this approach (open diamond symbols in Fig. 9a) lie above those obtained with the conventional linear rule, there is still a significant gap between the model predictions and the experimental results. 6. Conclusions Large strain compressionetension fracture experiments are performed on uniaxial and notched flat specimens extracted from dual phase steel sheets (DP780). High compressive in-plane strains (of up to 13%) are achieved using a floating antibuckling device. The relative displacement of the specimen boundaries as well as local strains on the specimen surface are measured using digital image correlation. A Combined Chaboche-Yoshida (CCY) model is proposed to account for the observed Bauschinger effect, transient softening and work hardening stagnation. The material parameter identification based on notched compressionetension experiments for very large strains is shown in detail. Subsequently, the model is applied to predict the material response in different monotonic experiments (notched tension, tension with a central hole and pure shear) and for different levels of reverse loading. The plasticity model predictions agree well with the results from all experiments, including the evolution of the surface strains at the specimen center. The extracted loading paths to fracture show a significant increase of the strain to fracture as a monotonic function of the applied pre-strain. For example, applying a pre-strain of 0.13 increases the strain to fracture from 0.57 (for monotonic loading)

Table 2 Net fracture strain (equivalent plastic strain at fracture minus compressive pre-strain) as a function of the compressive pre-strain.

Pre-compression strain Fracture strain Net fracture strain

NCT-0

NCT-3

NCT-6

NCT9

NTC-13

0. 0.56 0.56

0.03 0.60 0.66

0.06 0.73 0.67

0.09 0.75 0.66

0.13 0.77 0.64

42

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

to 0.77. The transient hardening of the material and the local thickening of the sheet during compression delay the formation of a neck and the consequent increase in triaxiality (and the consequent decrease in the Lode parameter). A Hosford-Coulomb damage indicator model with a non-linear damage accumulation rule is calibrated and validated based on the experimental data. The model predictions agree well with all experimental results for the DP780 steel with a maximum relative error of 5% in the strains to fracture. It is worth noting that the same phenomenological damage accumulation rule provided an accurate description of proportional and non-proportional experiments on aluminum 2024-T351 (Papasidero et al., 2015). Acknowledgments The partial financial support through MIT Industrial Fracture Consortium and the French National Research Agency (Grant ANR-11-BS09-0008, LOTERIE) is gratefully acknowledged. The authors are grateful to Professors Tomasz Wierzbicki (MIT) and Allison Beese (Penn State) for valuable discussions. References Abaqus, 2012. Reference Manuals v6.12. Dassault Systemes. Amstrong, P.J., Frederick, C.O., 1966. A Mathematical Representation of the Multiaxial Bauschinger Effect. G.E.G.B. Report RD/B/N 731. Bai, Y., Wierzbicki, T., 2008. A new model of metal plasticity and fracture with pressure and lode dependence. Int. J. Plast. 24 (6), 1071e1096. Bai, Y., Wierzbicki, T., 2010. Application of extended Mohr-Coulomb criterion to ductile fracture. Int. J. Fract. 161, 1e20. Bao, Y., Treitler, R., 2004. Ductile crack formation on notched Al2024-T351 bars under compressionetension loading. Mater. Sci. Eng. A 384 (1), 385e394. Barlat, F., Ferreira Duarte, J.M., Gracio, J.J., Lopes, A.B., Rauch, E.F., 2003. Plastic flow for non-monotonic loading conditions of an aluminum alloy sheet sample. Int. J. Plast. 19 (8), 1215e1244. Beese, A.M., Mohr, D., 2011. Effect of stress triaxiality and lode angle on the kinetics of strain-induced austenite-to-martensite transformation. Acta Mater. 59 (7), 2589e2600. Boger, R.K., Wagoner, R.H., Barlat, F., Lee, M.G., Chung, K., 2005. Continuous, large strain, tension/compression testing of sheet material. Int. J. Plast. 21 (12), 2319e2343. Brünig, M., Chyra, O., Albrecht, D., Driemeier, L., Alves, M., 2008. A ductile damage criterion at various stress triaxialities. Int. J. Plast. 24 (10), 1731e1755. Cao, J., Lee, W., Cheng, H.S., Seniw, M., Wang, H.P., Chung, K., 2009. Experimental and numerical investigation of combined isotropic-kinematic hardening behavior of sheet metals. Int. J. Plast. 25 (5), 942e972. Chaboche, J.L., 2008. A review of some plasticity and viscoplasticity constitutive theories. Int. J. Plast. 24 (10), 1642e1693. Chaboche, J.L., Rousselier, G., 1983. On the plastic and viscoplastic constitutive equationsdpart I: rules developed with internal variable concept. J. Press. Vessel Technol. 105 (2), 153e158. Chaboche, J.L., Dang Van, K., Cordier, G., 1979. Modelization of the strain memory effect on the cyclic hardening of 316 stainless steel. Struct. Mech. React. Technol. Trans. L. Chun, B.K., Jinn, J.T., Lee, J.K., 2002. Modeling the Bauschinger effect for sheet metals, part I: theory and part II: applications. Int. J. Plast. 18 (5), 571e616. Chung, K., Ma, N., Park, T., Kim, D., Yoo, D., Kim, C., 2011. A modified damage model for advanced high strength steel sheets. Int. J. Plast. 27, 1485e1511. Cvitanic, V., Vlak, F., Lozina, Z., 2008. A finite element formulation based on non-associated plasticity for sheet metal forming. Int. J. Plast. 24, 646e687. Dafalias, Y.F., Popov, E.P., 1976. Plastic internal variables formalism of cyclic plasticity. J. Appl. Mech. 43 (4), 645e651. Dietrich, L., Turski, K., 1978. New method for testing thin sheets subjected to compression. Rozpr. Inzynierskie 26 (1), 91e99. Dunand, M., Mohr, D., 2010. Hybrid experimentalenumerical analysis of basic ductile fracture experiments for sheet metals. Int. J. Solids Struct. 47 (9), 1130e1143. Dunand, M., Mohr, D., 2011. Optimized butterfly specimen for the fracture testing of sheet materials under combined normal and shear loading. Eng. Fract. Mech. 78 (17), 2919e2934. Dunand, M., Mohr, D., 2014. Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities. J. Mech. Phys. Solids. 66, 133e153. Eggertsen, P.A., Mattiasson, K., 2009. On the modelling of the bendingeunbending behaviour for accurate springback predictions. Int. J. Mech. Sci. 51 (7), 547e563. Eggertsen, P.A., Mattiasson, K., 2010. On constitutive modeling for springback analysis. Int. J. Mech. Sci. 52 (6), 804e818. Eggertsen, P.A., Mattiasson, K., 2011. On the identification of kinematic hardening material parameters for accurate springback predictions. Int. J. Material Form. 4 (2), 103e120. Geng, L., Wagoner, R.H., 2000. Springback Analysis with a Modified Hardening Model (No. 2000-01-0768). SAE Technical Paper. Gruben, G., Fagerholt, E., Hopperstad, O.S., Børvik, T., 2011. Fracture characteristics of a cold-rolled dual-phase steel. Eur. J. Mech. A/Solids 30 (Issue 3), 204e218.
nchez, P.J., 2012. A finite strain, finite band method for modeling ductile fracture. Int. J. Plast. 28 (1), 53e69. Huespe, A.E., Needleman, A., Oliver, J., Sa Khan, A.S., Liu, H., 2012. A new approach for ductile fracture prediction on Al 2024-T351 alloy. Int. J. Plast. 35, 1e12. Kuwabara, T., Morita, Y., Miyashita, Y., Takahashi, S., 1995. Elastic-plastic behavior of sheet metal subjected to in-plane reverse loading. J. Japan Soc. Technol. Plast. 36, 768e768. Lecarme, L., Tekoglu, C., Pardoen, T., 2011. Void growth and coalescence in ductile solids with stage III and stage IV strain hardening. Int. J. Plast. 27, 1203e1223. Li, H., Fu, M.W., Lu, J., Yang, H., 2011. Ductile fracture: experiments and computations. Int. J. Plast. 27 (2), 147e180. Lou, Y., Yoon, J.W., Huh, H., 2014. Modeling of shear ductile fracture considering a changeable cut-off value for stress triaxiality. Int. J. Plast. 54, 56e80. Luo, M., Dunand, M., Mohr, D., 2012. Experiments and modeling of anisotropic aluminum extrusions under multi-axial loadingepart II: ductile fracture. Int. J. Plast. 32, 36e58.
sar de Sa
, J.M.A., 2012. An assessment of isotropic constitutive models for ductile fracture under high and low stress Malcher, L., Andrade Pires, F.M., Ce triaxiality. Int. J. Plast. 30e31, 81e115. Mohr, D., Marcadet, S.J., 2015. Micromechanically-motivated phenomenological Hosford-Coulomb model for predicting ductile fracture initiation at low stress triaxialites. Int. J. Solids Struct. http://dx.doi.org/10.1016/j.ijsolstr.2015.02.024. Mohr, D., Dunand, M., Kim, K.H., 2010. Evaluation of associated and non-associated quadratic plasticity models for advanced high strength steel sheets under multi-axial loading. Int. J. Plast. 26 (7), 939e956. Mroz, Z., 1967. On the description of anisotropic workhardening. J. Mech. Phys. Solids 15 (3), 163e175. Papasidero, J., Doquet, V., Mohr, D., 2015. Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: BaoWierzbicki results revisited. Int. J. Solids Structures. http://dx.doi.org/10.1016/j.ijsolstr.2015.05.006. Pardoen, T., Hachez, F., Marchioni, B., Blyth, P.H., Atkins, A.G., 2004. Mode I fracture of sheet metal. J. Mech. Phys. Solids 52 (2), 423e452. Prager, W., 1956. A new method of analyzing stresses and strains in work hardening, Plastic Solids. J. Appl. Mech. ASME 23, 493.

S.J. Marcadet, D. Mohr / International Journal of Plasticity 72 (2015) 21e43

43

Roth, C.C., Mohr, D., 2014. Effect of strain rate on ductile fracture initiation in advanced high strength steel sheets: experiments and modeling. Int. J. Plast. 56, 19e44. Stoughton, T.B., 2002. A non-associated flow rule for sheet metal forming. Int. J. Plast. 18, 687e714. Sun, X., Choi, K.S., Liu, W.N., Khaleel, M.A., 2009. Predicting failure modes and ductility of dual phase steels using plastic strain localization. Int. J. Plast. 25 (10), 1888e1909. Sung, J.H., Kim, J.N., Wagoner, R.H., 2010. A plastic constitutive equation incorporating strain, strain-rate, and temperature, Int. J. Plast. 26 (12), 1746e1771. Yoshida, F., Uemori, T., 2002. A model of large-strain cyclic plasticity describing the Bauschinger effect and work hardening stagnation. Int. J. Plast. 18 (5), 661e686. Yoshida, F., Uemori, T., Fujiwara, K., 2002. Elasticeplastic behavior of steel sheets under in-plane cyclic tensionecompression at large strain. Int. J. Plast. 18 (5), 633e659. Zang, S.L., Guo, C., Thuillier, S., Lee, M.G., 2011. A model of one-surface cyclic plasticity and its application to springback prediction. Int. J. Mech. Sci. 53 (6), 425e435.