Crystal plasticity based constitutive model for uniaxial ratchetting of polycrystalline magnesium alloy

Computational Materials Science 84 (2014) 63–73

Contents lists available at ScienceDirect

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

Crystal plasticity based constitutive model for uniaxial ratchetting of polycrystalline magnesium alloy Chao Yu a, Guozheng Kang a,⇑, Qianhua Kan b a b

State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu, Sichuan 610031, PR China School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu, Sichuan 610031, PR China

a r t i c l e

i n f o

Article history: Received 10 September 2013 Received in revised form 4 November 2013 Accepted 19 November 2013 Available online 20 December 2013 Keywords: Magnesium alloy Crystal plasticity Ratchetting Dislocation slipping Twinning

a b s t r a c t In the framework of crystal plasticity, a constitutive model is constructed to describe the uniaxial ratchetting of polycrystalline magnesium alloy at room temperature. At the scale of single crystal, three kinds of slip systems, i.e., basal hai, prismatic hai and pyramidal hc + ai slip systems and one kind of twinning system, i.e., tension twinning system, are introduced into the model. New orientations of slip systems in the twinned region of single crystal are obtained by rotating the original ones, and then a rotation tensor is used in the proposed model. Evolution equations of dislocation slipping and twinning deformation are obtained in the form of power-law. Also, the residual twinning and its accumulation during the cyclic loading are considered in the proposed model. Finally, an explicit scale-transition rule is adopted to extend the proposed single crystal constitutive model to the polycrystalline version. The capability of the proposed model to describe the uniaxial ratchetting of the polycrystalline magnesium alloy is verified by comparing the predictions with corresponding experimental results. More predictions of the ratchetting of magnesium single crystal and polycrystalline alloys are also discussed by the proposed model. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Magnesium and magnesium alloy have been widely used for structural components due to their high specific strength/stiffness, electromagnetic shielding, recyclability and damping capacity [1–4]. It is well-known that magnesium is the lightest structural metal which can be used in automotive and aerospace industries. Its density is about two-thirds of aluminum and less than onequarter of iron. As a low symmetrical hexagonal close packed (HCP) crystal structure, the slip systems in magnesium are limited. Four modes of slip systems had been observed by experiments [5], such as basal hai, prismatic hai, pyramidal hai and pyramidal hc + ai slip systems. However, the critical resolved shear stresses (CRSS) of non-basal slip systems are very high, and the dislocation slipping cannot be easily activated in such systems. Thus, the twinning deformation plays an important role in the plastic deformation of magnesium and magnesium alloys [6,7]. Two twinning systems had been reported for the magnesium  2gh1 0 1  1i and comand magnesium alloys, i.e., the tensile f1 0 1  1gh1 0 1  2i twinning systems [8,9]. Since the CRSS pressive f1 0 1 of the compressive twinning system is much higher than that of the tensile one, the compressive twinning system cannot be ⇑ Corresponding author. Tel.: +86 28 87603794; fax: +86 28 87600797. E-mail addresses: [email protected], [email protected] (G. Kang). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.11.054

activated easily during the loading. At present, the mechanical properties of magnesium single crystals and polycrystalline alloys have been experimentally observed under the monotonic loading conditions [7,10–14]. It is shown that the mechanical responses of the extruded and rolled magnesium polycrystalline alloys are strongly asymmetric under the tensile and compressive loading conditions due to the basal plane textures. The initial compressive yielding stress is smaller than the tensile one. A usual convex curve occurs during the tension, while a concave curve takes place during the compression. In practice, the structural components made from the magnesium single crystal and its polycrystalline alloys often endure to a cyclic loading, and then the cyclic deformation of them is a key issue which should be discussed in advance in order to predict the fatigue life and assess the reliability of such components. Recently, some experimental observations were performed to the cyclic deformation of magnesium single crystals and polycrystalline alloys under the strain-controlled loading conditions [15– 23]. However, the ratchetting of magnesium single crystal and polycrystalline alloys was observed just by [24–27] under the stress-controlled cyclic loading conditions. It is shown that the ratchetting of magnesium single crystal is strongly dependent on the crystallographic orientation [24]; and the uniaxial ratchetting of hot-rolled AZ91D magnesium polycrystalline alloy depends greatly upon the applied mean stress, stress amplitude and stress rate [25–27].

64

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

Based on experimental observations, many constitutive models have been constructed to describe the mechanical responses of polycrystalline alloys. Such models can be classified into two groups: i.e., macroscopic phenomenological models and crystal plasticity based ones. As the typical macroscopic phenomenological models, the cyclic plastic and visco-plastic constitutive models originally provided by Armstrong and Frederick [28], Chaboche [29] and Ohno and Wang [30] and consequently extended by other researchers have been widely used to describe the ratchetting of many metals and alloys. However, such models are only appropriate for the initial isotropic materials without texture, because they adopt an isotropic von-Mises yielding surface. As mentioned above, the initial yield stresses and the strain hardening of magnesium polycrystalline alloys are different under the tensile and compressive loading conditions. In order to describe these features, some phenomenological models [31,32,17] have been constructed from the two-surface model [33] to describe the mechanical responses of magnesium polycrystalline alloys. However, these phenomenological models focus just on the responses of magnesium polycrystalline alloys presented under the monotonic tension and compression, rather than the cyclic deformation of the alloys. Thus, such models cannot used directly to describe the ratchetting of magnesium polycrystalline alloys. To describe the anisotropic mechanical responses of magnesium polycrystalline alloys, the crystal plasticity based constitutive models have been constructed recently, because the plastic mechanisms (slipping, twinning and detwinning) and the initial texture of the alloys can be easily introduced into the proposed models [21,34–39]. These models described the cyclic deformation of the magnesium single crystal and polycrystalline alloy under the strain-controlled cyclic loading conditions well. However, the ratchetting of the alloys presented in the stress-controlled cyclic loading cannot be described reasonably by the proposed models, because they were mainly established from the experimental results obtained under the strain-controlled cyclic loading conditions. Therefore, in this work, a crystal plasticity based constitutive model is constructed to describe the uniaxial ratchetting of polycrystalline magnesium alloys by considering the dislocation slipping and twinning deformation at the single crystal scale. An explicit scale-transition rule considering the inelastic accommodation of grains is adopted to construct the polycrystalline constitutive model from the single crystal version. The capability of the proposed model to describe the uniaxial cyclic deformation of the magnesium single crystal and polycrystalline alloy is verified by comparing the predicted results with the corresponding experimental ones obtained by [19,25]. It is seen that the predictions are in good agreement with the experiments. Furthermore, the predicted ratchetting of magnesium single crystal in different orientations is discussed by using the proposed model.

ntwin X

e_ ptwin ¼

c_ atwin Patwin

ð2-aÞ

a¼1

c_ atwin ¼ g twin f_ atwin Patwin ¼

ð2-bÞ

1 ðmatwin  natwin þ natwin  matwin Þ 2

ð2-cÞ

where c_ atwin and f_ atwin are the rates of shear deformation and volume fraction of ath twinning system, respectively. gtwin is the characteristic shear. ntwin is the number of twinning systems. Since only the tensile twinning systems are considered, the number ntwin equals to 6. Patwin is the orientation tensor of ath twinning system. matwin and natwin are the twinning plane normal and shear directions, respectively. As mentioned above, the crystal structure in the twinned region does not change, and its orientation can be obtained only by a rotation transformation of the parent (un-twinned) one. The rotation tensor caused by the twinning [40] can be written as:

Ratwin ¼ 2matwin  matwin  I

ð3Þ

where I is the second-ordered identity tensor. The plastic strain caused by the dislocation slipping in a single crystal can be obtained from the following formula:

e_ pslip ¼ ð1 

ntwin X

a ftwin Þ

a¼1

Pi0 slip ¼

nslip X

ntwin X

i¼1

a¼1

c_ i0slip Pi0slip þ

a ftwin

nslip X

a ia c_ islip Pslip

ð4-aÞ

i¼1

1 ðmislip  nislip þ nislip  mislip Þ 2 T

a a Pislip ¼ Ratwin Pi0 slip ðR twin Þ

ð4-bÞ

a ¼ 1; 2; . . . ; 6

ð4-cÞ

Pi0 slip

_ i0 slip

where c and are the slipping rate and the orientation tensor a and Pia are of ith slip system in the parent region, respectively. c_ islip slip the slipping rate and the orientation tensor of ith slip system in the twinned region caused by ath twinning system, respectively. mislip and nislip are the slip plane normal and shear directions in the parent region, respectively. nslip is the number of slip systems and equals to 12 here. For simplicity, it is assumed that the stress and strain fields in a single crystal are uniform. The relation between elastic strain and stress tensors in a single crystal can be written as: a ee ¼ Sðftwin Þ:r

2. Crystal plasticity based constitutive model

ð5Þ

Sðf a

where twin Þ is the forth-ordered elastic compliance tensor and is a function of the volume fraction of twinned region. In this work, it is a Þ is a linear function of f a , i.e., assumed that Sðftwin twin

2.1. Definitions of inelastic strains With the assumption of small deformation, the total strain tensor e can be decomposed into three parts, i.e., elastic strain tensor ee , the plastic strain tensor eptwin caused by twinning and epslip by dislocation slipping:

e ¼ ee þ eptwin þ epslip

slipping in the pyramidal hai slip systems and the compressive twinning system to the plastic deformation are neglected, because such slip and twinning systems are hard to be activated due to their high critically resolved shear stresses. The details about the slip and twinning systems considered in this work are listed in Table 1. The evolution equation of plastic strain caused by the twinning can be written as:

ð1Þ

In this work, one twinning system and three modes of slip systems are considered by referring to [9], i.e., tensile twinning system  2gh1 0 1  1i, and basal slipping system f0 0 0 1gh1 1 2  0i (hai), f1 0 1  0gh1 1 2  0i (hai) and pyramidal one prismatic one f1 0 1  2gh1 1 2  3i (ha + ci). The contributions of the dislocation f1 1 2

a

Sðftwin Þ ¼

1

ntwin X

!

ftwin S0 þ a

a¼1

ntwin X

a ftwin Sa

ð6Þ

a¼1

where S0 and Sa are the elastic compliance tensors for the parent and twinned regions, respectively. Similar to the orientation tensor of slip system, Sa can be obtained only by a rotation of S0, which can be written in the component form as

Saijkl ¼ Rai0 i Raj0 j Rak0 k Ral0 l S0i0 j0 k0 l0 S0i0 j0 k0 l0

a

a

a ¼ 1; 2; . . . ; 6

ð7Þ 0

a

a

where Sijkl and Ri0 i are the component forms of S , S and Rtwin , respectively.

65

C. Yu et al. / Computational Materials Science 84 (2014) 63–73 Table 1 Details about the slip/twin systems.

Basal slip Prismatic slip Pyramidal slip Tension twin

Slip/twin plane

Slip/twin direction

Number of systems

Amount of shear

{0 0 0 1}

 0i h1 1 2  0i h1 1 2  3i h1 1 2  1i  h1 0 1

3 3 6 6

– – – 0.129

 0g f1 0 1  2g f1 1 2  2g f1 0 1

Table 2 Material parameters for predicting the cyclic hardening of magnesium single crystal. C11 = 39.6 GPa; C12 = 17.1 GPa; C13 = 14.3 GPa; C33 = 41.1 GPa; C44 = 11 GPa a a a si0c;basal ¼ 0:5 MPa; si0c;pri ¼ 3 MPa; si0c;pyr ¼ 3 MPa

s0;for ¼ 2:4 MPa; ssat;for ¼ 10 MPa; s0;res ¼ 6 MPa;ssat;res ¼ 15 MPa hbasal = 100 MPa; hpri = 1 GPa; hpyr = 8 GPa; k = 2 GPa; q = 0.5 cbasal = 1 GPa; bbasal = 200; cpri = 1 GPa; bpri = 200; cpyr = 1 GPa; bpyr = 200 fsat = 0.02; btwin ¼ 30; m ¼ 100; k0 ¼ 0:001

2.2. Evolution equations of internal variables The rate of dislocation slipping in each slip system can be obtained in the form of power-law as following: a c_ islip

 ia   s  xi a  m  signðsia  xia Þ i ¼ 1; 2; . . . ; 12; a ¼ 0; 1; . . . ; 6 ¼ k0  sia  c

ð8Þ

(a)

20

(b) b

a

0

c

-20

d

1st cycle

-40

Experiment Simulation with RT Simulation without RT

-60 -80 -0.4

-20

d' -40

2nd cycle Experiment Simulation with RT Simulation without RT

-60 -80

e

-0.6

ð9Þ

where c_ afor and c_ ares are the shear rate caused by the twinning during the forward and reverse parts, respectively. It should be noted that the reverse loading part defined here means the change of loading direction (e.g., from tension to compression or vise versa), rather than just the unloading part.

f

0

Stress, MPa

c_ atwin ¼ c_ afor þ c_ ares

Stress, MPa

20

a where k0 is the reference slipping rate; sia ¼ r : Pislip ; xia and sica are the resolved shear stress, back stress and critical resolved shear stress of ith slip system in the twinned region caused by ath twinning system, respectively. It should be noted that a = 0 represents the dislocation slipping occurred in the parent region. i = 1, 2, . . ., 12 represents the twelve different slip systems and a = 0, 1, . . ., 6 represents the parent phase and six different twining systems. m is the factor of rate sensitivity. In order to distinguish the shear deformation caused by the twinning during the forward and reverse loading parts of cyclic loading, the shear rate c_ atwin can be decomposed into two parts, i.e.,

-0.2

0.0

0.2

0.4

0.6

-0.6

-0.4

-0.2

Strain, %

0.0

0.2

0.4

0.6

Strain, %

40

40

(c)

(d)

20 0 0

-40

d'' -60

400th cycle

-80

Experiment Simulation with RT Simulation without RT

-100 -120 -0.6

-0.4

-0.2

0.0

Strain, %

0.2

0.4

0.6

Stress, MPa

Stress, MPa

-20

-40

d''' -80

1610th cycle Experiment Simulation with RT Simulation without RT

-120

-160 -0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

Strain, %

Fig. 1. Experimental (from Yu et al. [19]) and simulated stress–strain curves for the cyclic hardening of magnesium single crystal in [0 0 0 1] under the strain-controlled cyclic loading: (a) 1st cycle; (b) 2nd cycle; (c) 400th cycle; and (d) 1610th cycle.

66

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

Owing to its polar nature, the twinning can be activated only when the resolved shear stress of twinning system is positive. Thus, the evolution equation of twinning rate occurred in the forward loading part can be written in a power-law form too, but it is different from that of dislocation slipping:

c_ afor ¼

8  m1 > < k0 sa xa

whenever

> :

other conditions

sac;for

0

0

ntwin X

a <1 ftwin

a¼1

a ¼ 1; 2; . . . ; 6 ð10Þ

where k00 is the reference twinning rate; sa ¼ r : Patwin ; xa and sac;for are the resolved shear stress, back stress and critical resolved shear stress of ath twinning system, respectively. m1 is the factor of rate sensitivity. hxi is McCauley bracket: when x > 0, hxi = x; when x 6 0, hxi = 0. It is reported that the plastic deformation caused by the twinning can be recovered when the loading direction is reversed, which is named as the detwinning by [21,37]. Zhang et al. [22] also concluded that the work hardening of magnesium alloy occurred during the cyclic deformation was closely related to such twinning–detwinning process. Furthermore, it is reported that the twinning could not be recovered completely by the detwinning, and the cyclic hardening of magnesium single crystal was partly attributed to the accumulation of residual twining [19]. Therefore, the evolution equation of the detwinning is set in this work as:

c_ ares ¼

8 D E < k00 sa þxa m2

a > fa whenever ftwin residual

:

other conditions

0

where k000 is the reference detwinning rate; sa ¼ r : Patwin ; xa and sac;res are the resolved shear stress, back stress and critical resolved shear stress of ath twinning system during the detwinning, respectively. a m2 is the factor of rate sensitivity. fresidual is the volume fraction of residual twining for the ath twinning system. Based on experimental observations [19], the evolution equation of residual twinning during the cyclic deformation is given in an exponential function form, i.e.,

" a

a

fresidual ¼ fsat

fa 1  exp  ac btwin

ð12-aÞ

ð12-bÞ

a ð¼ f Þ is the saturated volume fraction of residual twinwhere fsat sat a ning, btwin ð¼ btwin Þ is a material parameter governing the rate of saturation, and fca is the accumulated volume fraction of twinning for each twinning system.

12

(b) 8

4

4

Stress, MPa

8

0 -4

[0001] -3±7MPa [0001] 3±7MPa

-8 -12 -0.1

0.0

0.1

0.2

0.3

[1014] -3±7MPa [1014] 3±7MPa

0 -4 -8

0.4

-12 -0.4

0.5

-0.2

Strain,% 0.6

[1120] -3±7MPa [1120] 3±7MPa Ratchetting strain,%

0 -4

-0.1

0.0

Strain,%

0.4

0.6

40

50

[0001] 3±7MPa [0001] -3±7MPa [1014] 3±7MPa [1014] -3±7MPa

0.2 0.0 -0.2 -0.4

-8

-0.2

0.2

(d)

0.4

4

-12 -0.3

0.0

Strain,%

(c)

8

Stress, MPa

!#

f_ ac ¼ jf_ atwin j

(a)

12

a ¼ 1; 2; . . . ; 6 ð11Þ

12

Stress, MPa

sac:res

0

0.1

0.2

-0.6 0

[1120] 3±7MPa [1120] -3±7MPa 10

20

30

Number of cycles

Fig. 2. Prediction the uniaxial ratchetting of single crystal under cyclic tension-unloading and compression-unloading in different directions: (a) stress–strain curves for  0 1 4 direction; (c) stress–strain curves for ½1 1 2  0 direction; and (d) evolution of ratchetting strain vs. number of cycles. [0 0 0 1] direction; (b) stress–strain curves for ½1

67

C. Yu et al. / Computational Materials Science 84 (2014) 63–73 Table 3 Material parameters for predicting the rachetting of magnesium polycrystalline alloy. C11 = 58 GPa; C12 = 25 GPa; C13 = 20.8 GPa; C33 = 61.2 GPa; C44 = 16.6 GPa a a a si0c;basal ¼ 30 MPa; si0c;pri ¼ 100 MPa; si0c;pyr ¼ 150 MPa; s0;for ¼ 48 MPa; s0;res ¼ 48 MPa hbasal ¼ 30 MPa; hpri ¼ 30 MPa; hpyr ¼ 30 MPa; k ¼ 100 MPa; q ¼ 0:5 cbasal = 2 GPa; bbasal = 0.5; cpri = 0.4 GPa; bpri = 10; cpyr = 0.4 GPa; bpyr = 10 m ¼ 100; k0 ¼ 0:001 C = 2 GPa; D = 5

400

where Hij is the hardening coefficient matrix containing self- and latent-hardening. h is the hardening modulus. q is a constant determining the magnitude of latent-hardening. From Eq. (14-a), it is seen that the effect of latent hardening on s_ ica is considered only when the slip systems i and j are located in the same twinned region. It can be deduced from the experimental observation by [19] that both the kinematic and isotropic hardening features for the twinning/detwinning play important roles on the cyclic hardening of magnesium single crystal. Thus, a linear kinematic hardening rule is adopted here for twinning/detwinning, i.e.,

Monotonic tension Monotonic compression Tension-compression

300 200

Stress, Mpa

100 0 -100 -200

x_ a ¼ kc_ atwin

-300 -400 -10

-8

-6

-4

-2

0

2

4

6

8

ð15Þ

where k is the kinematic hardening modulus. On the other hand, the isotropic hardening rules are set as following,

10

Strain, %



Fig. 3. Predicted stress–strain curves of polycrystalline alloy under monotonic tension, monotonic compression and tension–compression conditions.



Stress σ22 , MPa

400

p

ε =0.04

p

ε =0.05

ε =0.01

300

ε =0.02

200

ε =0.03



sac;res ¼ s0;res þ ðssat;res  s0;res Þ 1  exp 

p p

fca btwin fca btwin

 ð16-aÞ  ð16-bÞ

where s0,for and ssat,for are the initial and saturated critical resolved shear stresses of each twinning system for the twinning; while s0,res and ssat,res are the initial and saturated critical resolved shear stresses of each twinning system for the detwinning, respectively.

p

100 0

2.4. Explicit scale-transition rule

-100 -200 p

-300

ε =0.08 p

ε =0.09

p

ε =0.10

ε =0.06

-400

ε =0.07

-500 -400

-300

-200

-100

0

100

200

p p

300

400

Stress σ11 , MPa Fig. 4. Predicted yield surfaces of polycrystalline alloy under biaxial loading condition.

2.3. Hardening rules for plastic deformation The Armstrong–Frederick kinematic hardening rule [28] is adopted for all the slip systems, i.e., ia a a  bx jc_ islip j x_ ia ¼ cc_ islip

ð13Þ

where c and b are two material parameters. A linear isotropic hardening rule is adopted for all the slip systems here, i.e.,

s_ ica ¼



sac;for ¼ s0;for þ ðssat;for  s0;for Þ 1  exp 

nslip X

a Hij jc_ jslip j

ð14-aÞ

i¼1

The constitutive equations mentioned above are derived to describe the cyclic deformation of magnesium single crystals. For the polycrystalline magnesium alloys, an effective scale-transition rule is required to deduce the polycrystalline constitutive model from a single crystal version. Hill [41] established a well-known self-consistent model to obtain certain polycrystalline responses from that of a single crystal. However, the self-consistent model requires an iterative algorithm due to its implicit expression, which is very time-consuming. Cailletaud and Pilvin [42] formulated an explicit scale-transition rule denoted as b-rule, with which the stress– strain response of polycrystalline aggregates can be easily obtained from that of a single crystal. The b-rule has been used to describe the cyclic deformation of polycrystalline materials by [21,43,44]. Thus, in this work, the b-rule is used. The b-rule provides that the local stress tensor r in a single crystal can be obtained from the applied uniform macroscopic stress tensor R by using the following formula:

r ¼ R þ Cðb  bg Þ

ð17-aÞ

b_ g ¼ e_ in  Dbg ke_ in k

ð17-bÞ

ein ¼ eptwin þ epslip

ð17-cÞ

g

Hij ¼ h½q þ ð1  qÞdij 

ð14-bÞ

where b = [b ] and the symbol [] denotes the volume average for the polycrystalline aggregations; C and D are material parameters.

68

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

The proposed crystal plasticity based cyclic constitutive model consists of Eqs. (1)–(17) in the intra-granular and inter-granular scales. In the next section, the capability of the proposed model to predict the uniaxial cyclic deformation of the magnesium single crystal and polycrystalline alloy is verified and discussed by comparing the predictions with the corresponding experiments [19,25]. The numerical integration procedure of the proposed model is outlined in the Appendix. 3. Verification and discussion 3.1. Simulations for the cyclic deformation of magnesium single crystal Firstly, the cyclic stress–strain responses of magnesium single crystal presented under the symmetrical strain-controlled cyclic tension–compression loading conditions in [0 0 0 1] direction [19] are simulated by the proposed model. Referring to [19], the maximum and minimum controlled strains are 0.5% and 0.5%, respectively, and the number of cycles is set to be 1610. The material parameters used in the proposed model are determined by the trial-and-error method from the experimental results and are listed in Table 2. For simplicity, the parameters of kinematic hardening for each slip systems are set as the same, i.e., cbasal = cpri = cpyr and bbasal = bpri = bpyr. Fig. 1a shows the experimental and simulated stress–strain curves obtained for the 1st cycle. It is seen that the plastic deformation occurs when the applied stress reaches to 4.8 MPa (Point

a). This is attributed to the activation of tensile twinning. The twinned regions expand with the increasing applied strain and an obvious twinning hardening feature is observed. After the maximum strain (Point b) is reached, the unloading starts. Loading segment b–c shows the elastic unloading. After then, the plastic deformation takes place again as shown in the segment c–d, which is caused by the detwinning deformation occurred in the compression. The detwinning was defined as the shrinking of the twinned regions by [21,37]. It should be noted that the hardening rate significantly increases in the segment d–e, which is explained by the activation of non-basal slip [19]. After the compressive unloading (segment e–f), twinning occurs again when the applied stress reaches to Point f. Such evolution process is reasonably predicted by the proposed model, and the predicted cyclic stress–strain curve is in a good agreement with the experimental one. Fig. 1b–d shows the experimental and predicted stress–strain curves obtained in the 2nd, 400th and 1610th cycles, respectively. It can be seen that all the peak stress, valley stress and the hysteresis loop increase with the increasing number of cycles, and the cyclic hardening of magnesium single crystal is well predicted by the proposed model with the residual twinning and its cyclic accumulation taken into account. As shown in Fig. 1, the predicted results of the model without considering the residual twinning and its cyclic accumulation cannot reflect the gradual movement of the inflection point of plastic hardening modulus with the increasing number of cycles (which is reflected by Points d0 , d00 and d000 in Fig. 1).

350

350

(b) Simulated

(a) Experimental Δσ/2=125 MPa 300 σ =150MPa Stress rate=50MPa/s m Δσ/2=75 MPa

Δσ/2=75 MPa

250

Stress, MPa

Stress, MPa

250

Δσ/2=125 MPa σm=150MPa Stress rate=50MPa/s

300

200 150 100

200 150 100 50

50

Δσ/2=100MPa

Δσ/2=100 MPa 0

0 0

1

2

3

4

5

0

1

2

3

4

5

Strain, %

Strain, % 5

(c)

Ratchetting strain, %

4

Δσ/2=125MPa experimental Δσ/2=125MPa simulated Δσ/2=100MPa experimental Δσ/2=100MPa simulated

3

2

Δσ/2=75MPa experimental Δσ/2=75MPa simulated

1

0 0

200

400

600

800

1000

Number of cycles Fig. 5. Uniaxial ratchetting of polycrystalline alloy with the constant mean stress of 150 MPa and stress rate of 50 MPa/s: (a) experimental stress–strain curve (from Lin et al. [25]); (b) simulated stress–strain curve; and (c) evolution of ratchetting strain vs. number of cycles.

69

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

As mentioned above, the capability of the proposed model to predict the cyclic deformation of magnesium single crystal under the strain-controlled cyclic loading conditions has been verified by comparing the predicted results with the experimental ones. Thus, the uniaxial ratchetting of magnesium single crystal is here predicted by the proposed model in the cyclic tension–compression cases. Also, the effect of different crystallographic orientations on the ratchetting of magnesium single crystal is discussed based on the predictions of the proposed model in the prescribed three  0 1 4 and ½1 1 2  0, respectively. The orientations, i.e., [0 0 0 1], ½1 loading conditions for the prescribed orientations are set as 3 ± 7 MPa and 3 ± 7 MPa (i.e., the applied mean stresses are 3 and 3 MPa, and the stress amplitudes are both 7 MPa). The number of cycles is set as 50. The material parameters are set to be the same as those listed in Table 2. Here, the ratchetting strain er is defined as er = (emax + emin)/2, where emax and emin are the maximum and minimum strains of each cycle, respectively. From the predictions shown in Fig. 2a, it is seen that the yielding stress of single crystal in [0 0 0 1] with the positive mean stress (3 MPa) is about 4.5 MPa, and the plastic deformation is mainly caused by the twinning deformation. However, with the negative mean stress (3 MPa), the yield stress is about 6 MPa. Since all the twinning systems cannot be activated under a compressive stress, the plastic deformation only comes from the dislocation slipping on the pyramidal slip systems. The ratchetting presented under the two loading conditions are mainly caused by the viscosity, because the reverse yielding has not occurred.

As shown in Fig. 2b and c, the predicted ratchetting strains of  0 1 4 and ½1 1 2  0, single crystal in other two orientations, i.e., ½1 are also dependent on the loading modes. For the orientation  0 1 4, the ratchetting strain predicted with the positive mean ½1 stress is higher than that with the negative one. However, for the  0, the ratchetting strain with the negative mean orientation ½1 1 2 stress is higher than that with the positive one. Fig. 2d gives the evolution of ratchetting strain in all cases. It should be noted that the evolutionary rule of predicted ratchetting of the magnesium single crystal and its dependence on the crystallographic orientation are similar to that observed by the experiments [24], but the magnitude of predicted ratchetting is apparently different from that of experimental ratchetting because the material parameters used in the proposed model are obtained by the experimental data from [19]. The magnesium single crystal used in [19] is different from that in [24]. It means that if the material parameters used in the proposed model are determined from the experimental data of the single crystal used in [24], the ratchetting observed there can also be reasonably described. 3.2. Simulations for the ratchetting of polycrystalline magnesium alloy The uniaxial ratchetting of polycrystalline AZ91D magnesium alloy observed in the cyclic tension–compression tests [25] is simulated and predicted by the proposed model, and then the effects of mean stress, stress amplitude and stress rate on the ratchetting are discussed in this subsection.

350

350

(b) Simulated

(a) Experimental

300 Δσ/2=125MPa Stress rate=50MPa/s

Δσ/2=125MPa Stress rate=50MPa/s

250 σm=100MPa

250 σm=100MPa

200

200

Stress, MPa

Stress, MPa

300

150 100

150 100 50

50

σm=125MPa

0

σm=150MPa

0

σm=125MPa

σm=150MPa

-50

-50 0

1

2

3

4

5

0

1

Strain, %

2

3

4

5

Strain, %

5

(c)

Ratchetting strain, %

4

σm=150MPa experimental σm=150MPa simulated

3

σm=125MPa experimental σm=125MPa simulated

2

σm=100MPa experimental σm=100MPa simulated

1

0 0

200

400

600

800

1000

Number of cycles Fig. 6. Uniaxial ratchetting of polycrystalline alloy with the constant stress amplitude of 125 MPa and stress rate of 50 MPa/s: (a) experimental stress–strain curve (from Lin et al. [25]); (b) simulated stress–strain curve; and (c) evolution of ratchetting strain vs. number of cycles.

70

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

same (i.e., k0 ¼ k00 ¼ k000 ), and the factors of rate sensitivity are set as m = m1 = m2 in the prediction of uniaxial ratchetting of the AZ91D magnesium alloy. Fig. 3 shows the predicted stress–strain curves of polycrystalline alloy under the monotonic tension, monotonic compression and tension–compression conditions. It is seen that the initial yield stress under the monotonic tension is about 220 MPa and is much higher than that under the compression (about 100 MPa). The usual convex curve is observed under the monotonic tension, whereas a concave curve is exhibited under the compression. The different hardening behaviors observed in the monotonic tensile and compressive tests ([15,23]) are described by the proposed model. Moreover, the reversed yield stress after the tensile loading is about 100 MPa. Thus, as mentioned above, the repeated twinning–detwinning mechanism can be neglected here. Fig. 4 shows the yield surfaces of polycrystalline magnesium alloy obtained by increasing the equivalent plastic strain from 0.01 to 0.10 with an increment of 0.01. It is seen that the material exhibits a strong anisotropy and a high asymmetry of yielding in different loading directions. Fig. 5 shows the experimental and predicted uniaxial ratchetting of AZ91D magnesium alloy in the loading cases with various stress amplitudes and a constant mean stress of 150 MPa, at the stress rate of 50 MPa/s, respectively. It is seen from Fig. 5 that the predictions are in good agreement with the experimental results. Both the cyclic hysteresis loops and the ratchetting strains are well

The extruded and rolled magnesium alloys are strongly anisotropic in the tensile and compressive cases due to their strong basal plane textures, and the c-axis of each grain is nearly perpendicular to the loading direction. Thus, in this work, 50 single crystal grains are used to represent the polycrystalline aggregates, and each grain is assigned an orientation which makes the direction of its c-axis being located in a scattered band of 80–100° to the loading direction. The elastic constants (C11, C12, C13, C33 and C44) are obtained by referring to [36,37]. It should be noted that, owing to the strong basal plane textures, the twinning deformation of magnesium alloy can be activated only when the compressive stress is large enough (e.g., 100 MPa or higher [15,35,37]). However, the maximum compressive stress in all the prescribed loading cases of [25] is only 25 MPa, and then the repeated twinning–detwinning mechanism considered in the proposed model does not occur (more details can be found in Fig. 3 in the following section). Thus, for simplicity, the isotropic hardening of twinning and the accumulated residual twinning caused by the repeated twinning–detwinning are neglected here, and then set s0,for = ssat,for, s0,res = ssat,res and fsat = 0. The value of s0,res is set as 48 MPa (referring to [34]). The other material parameters used in the proposed model are also determined by a trial-and-error method from the experimental data obtained in the load case of 150 ± 125 MPa and at the stress rate of 50 MPa/s, and are listed in Table 3. Furthermore, the reference slipping, twinning and detwinning rates are also set to be the

350

350

(a) Experimental

300

Δσ/2=125MPa σm=125MPa

100MPa/s

100

200 150

100MPa/s

100 50

50

10MPa/s

10MPa/s

0

0 -50 0.0

50MPa/s

Δσ/2= 125MPa 250 σm=125MPa

200 150

(b) Simulated

300

Stress, MPa

Stress, MPa

250

50MPa/s

0.5

1.0

1.5

2.0

2.5

3.0

-50 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Strain, %

Strain, % 3.0

(c) Ratchetting strain, %

2.5 2.0

100MPa/s experimental 100MPa/s simulated 50MPa/s experimental 50MPa/s simulated 10MPa/s experimental 10MPa/s simulated

1.5 1.0 0.5 0.0 0

200

400

600

800

1000

Number of cycles Fig. 7. Uniaxial ratchetting of polycrystalline alloy with the constant stress amplitude of 125 MPa and mean stress of 125 MPa with different loading rate: (a) experimental stress–strain curve (from Lin et al. [25]); (b) simulated stress–strain curve; and (c) evolution of ratchetting strain vs. number of cycles.

71

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

γ, % -1/2

-1/2

1.0

1.5

(a) Rhombus 3τ

Torsion equivalent strain 3

γ, %

1.5

Torsion equivalent strain 3

th

100 cycle σ

0.5

0.0

-0.5

-1.0 0.0

0.5

1.0

1.5

2.0

2.5

(b) Octagon 3τ

1.0 σ

th

100 cycle

0.5

0.0

-0.5

-1.0 0.0

3.0

0.5

1.0

Axial strain, % 1.5

γ, %

(c) Square -1/2

1.0

Torsion equivalent strain 3

Torsion equivalent strain 3

-1/2

γ, %

1.5

th

100 cycle 0.5

0.0 3τ

-0.5

σ

-1.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Axial strain, %

1.5

2.0

Axial strain, %

(d) Butterfly

1.0

th

100 cycle

0.5

0.0 3τ

-0.5

-1.0 0.0

σ

0.5

1.0

1.5

2.0

Axial strain, %

Fig. 8. Predicted axial strain–shear strain responses of the polycrystalline alloy under the non-proportionally multiaxial loading conditions with different loading paths: (a) rhombus path; (b) octagon path; (c) square path; and (d) butterfly path.

Axial ratchetting strain, %

3.0 2.5 2.0 1.5 1.0

Rhombus Octagon Square Butterfly

0.5 0.0 0

20

40

60

80

100

Number of cycles Fig. 9. Evolution of axial ratchetting strain vs. number of cycles under the nonproportionally multiaxial loading conditions with different loading paths.

predicted by the proposed model, because the reverse dislocation slipping of basal slip system during the unloading part of cyclic loading which is observed and discussed by [23,37] is described by the proposed model reasonably through setting a relatively low critical resolved shear stress in the basal slip system. The

dependence of ratchetting on the applied stress amplitude is also reasonably described. Figs. 6 and 7 show the experimental and predicted uniaxial ratchetting of AZ91D magnesium alloy in the load cases with various mean stresses and at different stress rates, respectively. It is seen from the figures that the predictions are in good agreement with the corresponding experimental results, and the effects of mean stress and stress rate on the ratchetting are also predicted reasonably. The predicted axial strain–shear strain responses of the magnesium polycrystalline alloy are also provided under the nonproportionally multiaxial cyclic loading conditions, as shown in Fig. 8, and the capability of the proposed model to predict the multiaxial ratchetting of the alloy is discussed. Here, four kinds of non-proportionally multiaxial loading paths are used, i.e., rhombus, octagon, square and butterfly paths. The minimum axial stresses are set as 0, 0, 36.6 and 213.4 MPa, respectively, and the maximum ones are set as 250, 250, 213.4 and 213.4 MPa, respectively. The minimum (negative) and maximum (positive) equivalent shear stresses are set to be the same in magnitude, i.e., 125, 125, 88.4 and 88.4 MPa, respectively. Fig. 9 gives the evolution of axial ratchetting strain vs. number of cycles with different non-proportional multiaxial loading paths. It is seen that the axial ratchetting is strongly dependent on the loading path and the stress amplitude.

72

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

4. Conclusions

If

*

ntwin X

a ðftwin Þ < 1 and

sanþ1  xan > 0 Dcatwin ¼ k00

a¼1

A crystal plasticity based cyclic constitutive model is constructed to describe the ratchetting of magnesium single crystal and polycrystalline alloy. In the proposed model, the dislocation slipping of basal, prismatic and pyramidal slip systems and the tensile twinning deformation are introduced as two microscopic mechanisms of plastic deformation. The evolution equation of each internal variable is obtained in the form of power-law. The polycrystalline model is extended from the single crystal version by using the explicit scale-transition b-rule. The capability of the proposed model to describe the uniaxial cyclic deformation of the magnesium single crystal and polycrystalline alloy are verified by comparing with the corresponding experiments. It is seen that the predictions agree well with the experimental results. Surely, much effort will be necessary to improve the proposed model in the future after more uniaxial and multiaxial ratchetting data are obtained.

* a a Else if ftwin > ðfresidual Þn and

sanþ1  xan < 0 Dcatwin ¼ k000

+m1

Dt

sanþ1 þ xanþ1 ðsac;res Þnþ1

Else Dcatwin ¼ 0 xa

nþ1 ,

dx Dca ¼ xan þ kDcatwin dcatwin twin

ðsac;for Þnþ1  ðsac;for Þn þ

dsac;for dcatwin

ðA5-1Þ

Dcatwin

¼ ðsac;for Þn     ðf a Þ ssat;for  s0;for þ exp  c n signðDcatwin ÞDcatwin g twin btwin btwin ðA5-2Þ ðsac;res Þnþ1  ðsac;res Þn þ

dsac;res a Dc dcatwin twin

¼ ðsac;res Þn     ðf a Þ ssat;res  s0;res exp  c n signðDcatwin ÞDcatwin : þ g twin btwin btwin Substituting Eq. (A5) into Eqs. (A4-1) and (A4-2), the nonlinear algebraic equations of Dcatwin can be written as:

* F 1 ðDcatwin Þ ¼ Dcatwin  k00

a

sanþ1  xan  dcdxa Dcatwin twin

dsa

ðsac;for Þn þ dcac;for Dcatwin

+m1

Dt ¼ 0

2

a

a

F ðDctwin Þ ¼ Dctwin þ

k000

a * a +m1 snþ1 þ xan þ dcdxa Dcatwin twin

dsa

ðsac;res Þn þ dcac;res Dcatwin

Dt ¼ 0:

twin

ðA6-2Þ The above mentioned equations can be easily solved by the Newton–Raphson method.

The discretized form of Eq. (17) can be written as:

Dr ¼ DR þ CðDb  Dbg Þ

ðA1Þ

However, an iterative process is needed to solve Eq. (A1). In order to save the computing time, an explicit algorithm is adopt here, which can be written as:

rnþ1 ¼ Rnþ1 þ Cðbn  bgn Þ

ðA2Þ

By Eq. (A2), the local stress for any grain of polycrystalline aggregates can be calculated directly. Step 2: Calculating the resolved shear stress The resolved shear stress for a grain in the polycrystalline aggregates at step n + 1 can be calculated as:

¼ rnþ1 :

ðA6-1Þ

twin

Step 1: Calculating the local stress

s

ðsac;for Þnþ1

ðA5-3Þ

In this appendix, the numerical integration procedure of the proposed model is outlined. Considering the interval from step n to n + 1 with the time increment Dt. For a grain in the polycrystalline aggregates, it is assumed that the stress tensor rn, strain tensor a Þ , ðxia Þ , en , and other internal variables such as xn ¼ fðcatwin Þn ;ðcislip n n a ðsica Þn , ðxa Þn , ðsac;for Þn , ðsac;res Þn , ðfresidual Þn , ðfca Þn , bgn g, ðeptwin Þn and ðepslip Þn in the step n are all known. For a given macroscopic stress tensor increment DR, the macroscopic stress tensor in the step n + 1 can be written as: Rn+1 = Rn + DR, the unknown quantities are: xn+1, enþ1 , ðeptwin Þnþ1 and ðepslip Þnþ1 .

a Pislip

Dt

a

xanþ1  xan þ

A.1. Numerical integration procedure

ia nþ1

+m 2

ðA4-3Þ

Noted that the first-order Taylor series expansions of and ðsac;res Þnþ1 can be respectively written as:

Appendix A

sanþ1 ¼ rnþ1 : Patwin

ðA4-1Þ

ðA4-2Þ

Acknowledgements Financial supports by the National Natural Science Foundation of China (11025210), the special fund for Sichuan Provincial Youth Science and Technology Innovation Team (2013TD0004), the Cultivation Foundation of Excellent Doctoral Dissertation of Southwest Jiaotong University (2013) and the Fundamental Research Funds for the Central Universities (China) are appreciated.

sanþ1  xanþ1 ðsac;for Þnþ1

a Step 4: Calculating Dcislip

The discretized form of Eq. (8) can be written as: ia slip

Dc

  sia  xia m a

 nþ1 nþ1  a ¼ k0   xinþ1 Dt:  sign sinþ1 a  sic;nþ1 

a and sia The first-order Taylor series expansions of xinþ1 c;nþ1 can be written as:

h i dx ia a a a Dcislip ¼ xina þ c  bxnþ1 signðDcislip Þ Dcislip ia dcslip ia

a xinþ1  xina þ

ðA3-1Þ ðA3-2Þ

Step 3: Calculating Dcatwin The discretized forms of Eqs. (9)–(11) can be written as:

ðA7Þ

a a sic;nþ1  sic;n þ

nslip X dsia j¼1

a ¼ sic;n þ

nslip X j¼1

a dcjslip

ðA8-1Þ

a Dcjslip

a a Hij signðDcjslip ÞDcjslip

ðA8-2Þ

C. Yu et al. / Computational Materials Science 84 (2014) 63–73

Substituting Eq. (A8) into Eq. (A7), the nonlinear algebraic equaa can be written as: tions of Dcislip a a a ia f ia ðDc1slip ; Dc2slip ; . . . ; Dc12; slip Þ ¼ Dcslip m      a xina þcDcislip  ia snþ1  1þbjDcia j   a

  slip a  k0   xinþ1 Dt  sign sinþ1  nslip    sia þ X dsia Dcja   c;n ja slip  dcslip   j¼1

¼ 0: ðA9Þ

Similarly, Eq. (A9) can be solved by the Newton–Raphson method. Step 5: Updating the variables a are All the other variables can be updated after Dcatwin and Dcislip known, i.e.,

a a ðftwin Þnþ1 ¼ ðftwin Þn þ

Dcatwin g twin

ðA10-1Þ

 a  Dc  ðfca Þnþ1 ¼ ðfca Þn þ  twin  g twin

ðA10-2Þ

   ðf a Þ a a ðfresidual Þnþ1 ¼ fsat 1  exp  c nþ1 btwin

ðA10-3Þ

ðxia Þnþ1 ¼

a xina þ cDcislip

ðA10-4Þ

a j 1 þ bjDcislip

ðsica Þnþ1 ¼ ðsica Þn þ

nslip X a Hij jDcjslip j

ðA10-5Þ

j¼1

ðxa Þnþ1 ¼ ðxa Þn þ kDcatwin

ðA10-6Þ

   ðf a Þ ðsac;for Þnþ1 ¼ sa0;for þ ðsasat;for  sa0;for Þ 1  exp  c nþ1 btwin

ðA10-7Þ

   ðf a Þ ðsac;res Þnþ1 ¼ sa0;res þ ðsasat;res  sa0;res Þ 1  exp  c nþ1 btwin

ðA10-8Þ

ðeptwin Þnþ1 ¼ ðeptwin Þn þ

ntwin X

Dcatwin Patwin

ðA10-9Þ

a¼1

ðepslip Þnþ1 ¼ ðepslip Þn þ 1  þ

ntwin X

a¼1

ntwin X

! a ðftwin Þnþ1

a¼1

nslip X a a ia ðftwin Þnþ1  Dcislip Pslip i¼1



nslip X i0 Dci0 slip Pslip i¼1

ðA10-10Þ

a enþ1 ¼ Sðftwin Þ : rnþ1 þ ðeptwin Þnþ1 þ ðepslip Þnþ1

bgnþ1 ¼

bgn þ Dein 1 þ DkDein k

73

ðA10-11Þ ðA10-12Þ

References [1] D. Eliezer, E. Aghion, F.H. Froes, Adv. Perform. Mater. 5 (1998) 201–212. [2] B.L. Moedike, T. Ebert, Mater. Sci. Eng. A 302 (2001) 37–45. [3] Y. Chino, T. Furuta, M. Hakamada, M. Mabuchi, J. Mater. Res. 41 (2006) 3229– 3232. [4] H.J. Kim, S.C. Choi, K.T. Lee, H.Y. Kim, Mater. Trans. 49 (5) (2008) 1112–1119. [5] S. Graff, W. Brocks, D. Steglich, Int. J. Plast. 23 (2007) 1957–1978. [6] M. Gharghouri, G. Weatherly, J. Embury, J. Root, Philos. Mag. 79 (1999) 1671– 1695. [7] C.H. Caceres, P. Lukac, Philos. Mag. 88 (7) (2008) 977–989. [8] A. Staroselsky, L. Anand, Int. J. Plast. 19 (2003) 1843–1864. [9] A.L. Oppedal, H.E. Kadiri, C.N. Tomé, G.C. Kaschner, S.C. Vogel, J.C. Baird, M.F. Horstemeyer, Int. J. Plast. 30–31 (2012) 41–61. [10] O. Muránsky, D.G. Carr, P. Šittner, E.C. Oliver, Int. J. Plast. 25 (2009) 1107–1127. [11] T.S. Srivatsan, S. Vasudevan, M. Petraroli, J. Alloys Compd. 461 (2008) 154–159. [12] K. Mueller, S. Mueller, J. Mater. Process. Technol. 187–188 (2007) 775–779. [13] Z. Trojanová, P. Lukácˇ, J. Mater. Process. Technol. 162–163 (2005) 416–421. [14] C.J. Geng, B.L. Wu, X.H. Du, Y.D. Wang, Y.D. Zhang, F. Wagner, C. Esling, Mater. Sci. Eng. A 559 (2013) 307–313. [15] X.Y. Lou, M. Li, R.K. Boger, S.R. Agnew, R.H. Wagoner, Int. J. Plast. 23 (2007) 44– 86. [16] L. Wu, A. Jain, D.W. Brown, G.M. Stoica, S.R. Agnew, B. Clausen, D.E. Fielden, P.K. Liaw, Acta Mater. 56 (2008) 688–695. [17] M. Li, X.Y. Lou, J.H. Kim, R.H. Wagoner, Int. J. Plast. 26 (2010) 820–858. [18] Q. Yu, J.X. Zhang, Y.Y. Jiang, Q.Z. Li, Int. J. Fatigue 44 (2010) 225–233. [19] Q. Yu, J.X. Zhang, Y.Y. Jiang, Philos. Mag. Lett. 91 (2011) 757–765. [20] J. Albinmousa, H. Jahed, S. Lambert, Int. J. Fatigue 33 (2011) 1127–1139. [21] G. Guillemer, M. Clavel, G. Cailletaud, Int. J. Plast. 27 (2011) 2068–2084. [22] J.X. Zhang, Q. Yu, Y.Y. Jiang, Q.Z. Li, Int. J. Plast. 27 (2011) 768–787. [23] T. Hama, Y. Kariyazaki, N. Hosokawa, H. Fujimoto, H. Takuda, Mater. Sci. Eng. A 551 (2012) 209–217. [24] Q.Z. Li, Mater. Sci. Eng. A 556 (2012) 301–308. [25] Y.C. Lin, X.M. Chen, G. Chen, J. Alloys Compd. 509 (2011) 6838–6843. [26] Y.C. Lin, X.M. Chen, Z.H. Liu, J. Chen, Int. J. Fatigue 48 (2013) 122–132. [27] Y.C. Lin, Z.H. Liu, X.M. Chen, J. Chen, Mater. Sci. Eng. A 573 (2013) 234–244. [28] P. Armstrong, C. Frederick, A Mathematical Representation of Multiaxial Bauschinger Effect, CEGB Report rd/b/n731, Berkeley Nuclear Laboratories, Berkeley, UK, 1965. [29] J.L. Chaboche, Int. J. Plast. 2 (1986) 149–188. [30] N. Ohno, J.D. Wang, Int. J. Plast. 9 (1993) 375–390. [31] D.L. McDowell, J. Appl. Mech. 52 (1985) 298–302. [32] M.G. Lee, R.H. Wagoner, J.K. Lee, K. Chung, H.Y. Kim, Int. J. Plast. 24 (2008) 545– 582. [33] J.H. Kim, D. Kim, Y.S. Lee, M.G. Lee, K. Chung, H.Y. Kim, R.H. Wagoner, Int. J. Plast. 50 (2013) 66–93. [34] T. Hama, H. Takuda, Int. J. Plast. 27 (2011) 1072–1092. [35] H. Wang, P.D. Wu, C.N. Tomé, J. Wang, Mater. Sci. Eng. A 555 (2012) 93–98. [36] H. Wang, P.D. Wu, J. Wang, Int. J. Plast. 47 (2013) 49–64. [37] H. Wang, P.D. Wu, J. Wang, C.N. Tomé, A crystal plasticity model for hexagonal close packed (HCP) crystals including twinning and de-twinning mechanisms, Int. J. Plast. 49 (2013) 36–52. [38] T. Hama, T. Hirohiko, Comput. Mater. Sci. 51 (2012) 156–164. [39] A. Izadbakhsh, I. Kaan Inal, R.K. Mishra, N. Marek, Comput. Mater. Sci. 50 (2011) 2185–2202. [40] J.W. Christian, S. Mahajan, Progr. Mater. Sci. 39 (1995) 1–157. [41] R. Hill, J. Mech. Phys. Solids 13 (1965) 89–101. [42] G. Cailletaud, P. Pilvin, Rev. Eur. Elem. Fin 3 (1994) 515–541. [43] G. Cailletaud, K. Sai, Mater. Sci. Eng. A 480 (2008) 24–39. [44] G.Z. Kang, O.T. Bruhns, K. Sai, Comput. Mater. Sci. 50 (2011) 1399–1405.