asymmetric hardening behavior of magnesium alloy sheets: Application to sheet springback

Available online at www.sciencedirect.com

International Journal of Plasticity 25 (2009) 70–104 www.elsevier.com/locate/ijplas

Constitutive modeling for anisotropic/asymmetric hardening behavior of magnesium alloy sheets: Application to sheet springback M.G. Lee a,*, S.J. Kim a, R.H. Wagoner b, K. Chung c, H.Y. Kim d a

Ferros Alloys Research Group, 531 Changwondaero, Korea Institute of Materials Science, Changwon 641-010, Republic of Korea b Department of Materials Science and Engineering, 2041 College Road, Ohio State University, Columbus, OH 43210, USA c School of Materials Science and Engineering, San 56-1, Shinlim-dong, Kwanak-gu, Seoul National University, Seoul 151-742, Republic of Korea d Division of Mechanical Engineering and Mechatronics, 192-1 Hyoja 2-Dong, Kangwon National University, Chunchon 200-701, Republic of Korea Received 14 September 2007; received in final revised form 7 December 2007 Available online 23 December 2007

Abstract The constitutive model for the unusual asymmetric hardening behavior of magnesium alloy sheet presented in a companion paper (Lee, M.G., Wagoner, R.H., Lee, J.K., Chung, K., Kim, H.Y., 2008. Constitutive modeling for anisotropic/asymmetric hardening behavior of magnesium alloy sheet, Int. J. Plasticity 24(4), 545–582) was applied to the springback prediction in sheet metal forming. The implicit finite element program ABAQUS was utilized to implement the developed constitutive equations via user material subroutine. For the verification purpose, the springback of AZ31B magnesium alloy sheet was measured using the unconstrained cylindrical bending test of Numisheet (Numisheet ’2002 Benchmark Problem, 2002. In: Yang, D.Y., Oh, S.I., Huh, H., Kim, Y.H. (Eds.), Proceedings of 5th International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, Jeju, Korea) and 2D draw bend test. With the specially designed draw bend test the direct restraining force and long drawn distance were attainable, thus the measurement of the springback could be made with improved accuracy comparable with conventional U channel draw

*

Corresponding author. Tel.: +82 55 280 3432; fax: +82 55 280 3599. E-mail address: [email protected] (M.G. Lee).

0749-6419/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2007.12.003

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

71

bend test. Besides the developed constitutive models, other models based on isotropic constitutive equations and the Chaboche type kinematic hardening model were also considered. Comparisons were made between simulated results by the finite element analysis and corresponding experiments and the newly proposed model showed enhanced prediction capability, which was also supported by the simple bending analysis adopting asymmetric stress–strain response. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: A. Springback; B. Constitutive behavior; Asymmetric material; C. Finite elements

1. Introduction Magnesium alloys have been reported to have unusual stress–strain responses compared with other conventional cubic crystalline metals used for the automotive and electronic industries (Bettles and Gibson, 2005). For instance, magnesium alloy sheet has three deformation modes during the in-plane uni-axial cyclic tests. Uni-axial tensile deformation is dominated by the slip mode and the flow curve is nominal concave-down shape. However, the deformations during the in-plane compression (or compression following tension) and tension following compression are dominated by the twinning and untwinning modes which show unusual concave-up or S-shape curves. In addition to the different types of hardening behavior, the size of linear (elastic) region is also unique. The size of elastic region increases monotonously from the undeformed state during the slip deformation, while instant reduction of the linear region and subsequent (re)increase are observed for the twinning and untwinning deformations. In terms of the directional difference in the initial yield and difference between yield stresses in tension and compression are also unique properties of magnesium alloy sheet. For more details of the experimental behavior of magnesium alloy sheet in room temperature are well documented in elsewhere (Lou et al., 2007; Duygulu and Agnew, 2003; Agnew and Duygulu, 2005). Therefore, to increase prediction capability in the sheet metal forming with magnesium alloys, proper constitutive equations which reflect the aforementioned unique mechanical responses might be inevitable. The companion paper of the present study by Lee et al. (2008) proposed constitutive modeling for the magnesium alloy sheets to represent their unusual mechanical behavior in a practical way. The previous paper was devoted to the models for anisotropic yield and asymmetric hardening. First of these is the modified Drucker–Prager yield criterion (Drucker, 1949) with pressure dependent term and the later is practical two-surface model (Lee et al., 2007) based on classical Dafalias and Popov (1976) and Krieg models (1975). The plasticity theory based on continuum mechanics was formulated mathematically to be implemented into the finite element method. The efficient numerical algorithm based on Newton–Raphson iteration was proposed. In addition to the phenomenological plasticity theory and its numerical implementation, the paper also introduced the characterization procedure of material parameters necessary for the developed constitutive equations. From the cyclic stress–strain curves which were measured by specially designed tension/compression test (Boger et al., 2005), empirical hardening equations were proposed. Finally, the developed theory and implementation were validated by correlation between measured stress–strain responses and corresponding calculations for the two magnesium alloy sheets. In this study, as a continuing work of the companion paper, the developed constitutive models for AZ31B magnesium alloy sheet have been applied to the finite element simula-

72

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

tions for the sheet springback. The springback is elastically driven response of sheet by releasing the forming tools after stamping (Wagoner et al., 2006). Lots of previous researches have been carried out to accurately predict material responses during the forming and springback by proposing new or modified constitutive equations. Regarding the deformation behavior in sheet metals, most of works have been focusing on proper description of real mechanical behavior such as transient behavior (Doucet and Wagoner, 1987; Chung and Wagoner, 1986), Bauschinger effect (Bauschinger, 1886; Prager, 1956; Ziegler, 1959; Hodge, 1957) and both by single surface (Amstrong and Frederick, 1966; Chaboche, 1991; Chung et al., 2005; Ohno and Wang, 1993; Ohno and Kachi, 1986) or multi-surface model (Mroz, 1967; Mroz and Niemunis, 1987; Krieg, 1975; Dafalias and Popov, 1976; Hashiguchi, 1981, 1997; McDowell, 1985; Khoei and Jamali, 2005). More recently, permanent softening induced from changing deformation paths has also been considered on top of Bauschnger and transient behavior (Geng and Wagoner, 2002; Chun et al., 2002; Lee et al., 2007). The detailed descriptions of single and multi-surface plasticity theories are well documented in other articles (Khan and Huang, 1995; Flores et al., 2007). Although these models have been able to represent the complex mechanical behavior observed in conventional automotive sheets, they cannot reasonably represent the asymmetry in yield and hardening of magnesium alloys. Though several models based on crystal plasticity have been proposed for the constitutive behavior of hcp metals (Lebenson and Tome, 1993; Brown et al., 2005; Kalidindi, 1998; Staroselsky and Anand, 2003), there are limited number of phenomenological models which are based on continuum plasticity and can be readily implemented in the finite element simulations for metal forming simulations (Li, 2006; Lee et al., 2008; Kim et al., submitted). Especially, to the best of our knowledge, there are few tries for the application of finite element simulations to the prediction of springback in a systematic ways. Therefore, in the following sections, prediction capability of the proposed constitutive model considering asymmetry and anisotropy of magnesium alloy sheet has been validated by comparing with experimentally measured springback. Two frequently used numerical examples are considered: the unconstrained cylindrical bending test proposed in Numisheet (2002) benchmark problem (Numisheet, 2002) and recently developed 2D the draw bend test. In addition to the numerical works, intensive experimental procedure to measure the springback of AZ31B magnesium sheet has also been performed with a specially designed 2D draw bend test with direct controllable restraining back force (Carden et al., 2002). 2. Summary of constitutive modeling The constitutive equations for the asymmetric/anisotropic hardening behavior of the magnesium alloy sheets were developed in the companion paper of the present work and summarized in the following sections. For more detailed theoretical derivations, refer to Lee et al. (2007). 2.1. Flow rule and yield function The loading surface in general form is iso ¼ 0: Uðr  aÞ  r

ð1Þ

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

73

iso is the effective stress measuring where r, a are Cauchy and back-stress, respectively and r the size of the yield surface. The evolution of back stress adopted in this paper is   d a d a m de da ¼ m¼ : ð2Þ iso ðmÞ iso ðmÞ r de r where the translational direction of the loading surface follows Ziegler type law, m  (r  a). Using the generalized linear elasticity with the additive decomposition of total strain increment (de) and normality rule for the plastic strain increment (dep) under the plastic work equivalence principle, the equivalent plastic strain increment can be updated every loading step as follows: de ¼

o riso  C  de  d a oðraÞ o riso o riso  C  oðraÞ þ hiso oðraÞ

;

ð3Þ

Þ is the slope of the isotropic hardening curve. where C is the elasticity matrix, hiso ð drdiso e The bounding surface is similarly defined as PðR  AÞ  Riso ¼ 0;

ð4Þ

where R, A are the stress and back-stress of the bounding surface, respectively, and Riso represents the size of the bounding surface. The general assumptions for the relative motion of the inner (loading) and bounding surfaces are Riso ðr  aÞ; iso r dA  da ¼ dlðR  rÞ

RA¼

ð5Þ or dA ¼ da  dlðR  rÞ ¼ dA1  dA2 ;

ð6Þ

2 iso ðdA2 Þ and r iso ¼ r iso ðR  rÞ. where dA2 ¼ dA  ðR  rÞ with dA2 ¼ r r  iso The hardening can be expressed as a nonlinear mixture rule with proper functions f(ml) and g(mb).

r þ ð1  f ðml ÞÞd r ¼ d riso þ d a; d r ¼ f ðml Þd

ð7Þ

dR ¼ gðmb ÞdR þ ð1  gðmb ÞÞdR ¼ dRiso þ dA;

ð8Þ

In this paper, a simple linear mixture rule with constant f = ml and g = mb is adopted.  þ d (or With the relationship among the two stresses and gap function ðdÞ, R ¼ r dR d r d d ¼ de þ de Þ, the hardening of loading surface can be obtained. de Under the current two-surface scheme, the gap distance between bounding and loading curves should be properly defined. From the reference 1-D stress–strain curves, the following form is proposed.   ¼ R  ð t ÞÞ ¼ ð1  #ÞðR  r t Þ; d ¼ UðRxx  rxx Þ ¼ Rxx  rxx ¼ R  r rt þ #ðR  r ð9Þ where Rxx and rxx are uni-axial tensile stresses on the loading and bounding curves, t is the reference stress as indicated in Fig. 1, and # is a sigmoid function respectively, the r with S-shape, which gives a value between 0 and 1. Here, for simplicity without deteriorating the meaning of bounding surface, static bounding curve is assumed as R ¼ R0 .

74

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

1 0

1

0 t

Fig. 1. Schematic illustration of S-shaped loading curve and corresponding bounding curve.

In addition to the hardening law, the yield function should be defined to represent the anisotropy and asymmetry of the initial yield stress. In this paper, the Drucker–Prager model is slightly modified by adding the anisotropic coefficients U ¼ pðr2xx  b2 rxx ryy þ b22 r2yy þ 3b23 r2xy Þ

1=2

iso ¼ 0; þ qðrxx þ b4 ryy Þ  r

ð10Þ

iso denotes the size of the yield surface and the five parameters p, q, b2, b3 and b4 where r are the material constants to be determined experimentally. 2.2. Characterization of material properties One-dimensional continuous in-plane tension–compression and compression–tension tests were performed to verify the implementation of developed theory (Boger et al., 2005). Since the test device uses the clamping side force to prevent buckling during compression, two sources which lead to over (or under)-estimate the real uni-axial stress–strain response are needed to be corrected: friction and bi-axial effects. For the friction effect, the previous work by Balakrishnam (1999) was adopted, in which several stress data were obtained using various clamping forces near the optimum force and then the correct stress was obtained by extrapolating the data to the value at the zero clamping force. In terms of iso in Eq. (10) is evaluated by the the effect of bi-axial stress state, the effective stress r known reference stress rxx and measured stress by the clamping force and contacting area in the direction of thickness. Because the thickness stress was much smaller than the stress in the testing direction, the bi-axial effect was not critical, thus ignored. Three important features were observed from the experiments: (1) Strong asymmetry in tensile and compressive yield stresses, (2) unusual concave-up shape during the compression and S-shape curves during tension following compression, and (3) reduction in the size of elastic range during the compression and compression following tension. Considering these features, hardening relations by defining empirical gap functions for each case were proposed. For gap function of the slip mode which corresponds to the initial tensile mode, the following exponential function is used  d ¼ as ð din Þ þ bs ð din Þ expðcs ð din Þel Þ; ð11Þ

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

75

where as, bs, cs are material parameters which depend on the initial gap distance din and el is the plastic strain whose value is re-initialized for each reverse loading. Therefore, the total equivalent plastic strain is the sum of el from each load reversal. In terms of the gap functions for twin mode (which corresponds to initial compression or compression following tension in the uni-axial test) and untwin mode (which corresponds to tension folt and # in Eq. lowing compression in the uni-axial test), the following empirical forms for r (9) were adopted in Þel ; t ¼ aT or UT ðdin Þ þ bT or UT ðd ð12Þ r # ¼ y 0;T or UT ðdin Þ þ

cT or UT ð d Þ  l in ; e x0;T or UT ð d Þ 1 þ exp  d T or UT ðd Þin

ð13Þ

in

where aT or UT, bT or UT, y0, T or UT, cT or UT, x0, T or UT, and dT or UT are material parameters obtained from measured gap distance with different initial gap distance. Note that the subscripts ‘T’ and ‘UT’ denote the twinning and untwinning modes, respectively. Thus, total 15 material parameters are needed to properly represent the unusual stress–strain response of asymmetric materials. To consider the reduction of the yield surface (which corresponds to the reduction of elastic linear region in the uni-axial test) during the twinning mode, a constant value k as a shrinkage ratio is assumed. Since the active stress shares the same point on the loading surfaces before and after instant softening, the back stress of loading surface can be calculated as following: a0 ¼ r 

0iso r ðR  AÞ; Riso

ð14Þ

0

where a and a are updated back-stress and previous back-stress before instant softening, 0iso is the size of the yield surface after instant softening considering the respectively. And, r constant reduction ratio: r0iso ¼ k riso . Considering the instant softening, the different isotropic hardening ratio ml of Eq. (7) during the twinning mode with the one during slip mode is introduced. This is based on the experimental observations. That is, the isotropic hardening during the slop deformation is different with that by twinning deformation after instant softening. The two parameters ml,T and ml,s are considered to obtain proper expansion of yield surface during the twinning and slip modes, respectively. In addition to the hardening functions, the material parameters for the (initial) asymmetric/anisotropic yield function needs to be characterized. Five material parameters in Eq. (10) can be determined from the two tensile yield stresses rTxx , rTyy , two compressive yield stresses rCxx , rCyy in the x (0° from rolling direction), y (90° from rolling direction) directions, and the shear yield stress rxy or tensile yield stress rTxx45 in the 45 degree direction from rolling direction. More complex yield functions have been proposed to include the effect of r-values in addition to the yield stress, but this issue is beyond the scope of the current work (Cazacu and Barlat, 2001, 2004). For more details, refer to the companion paper by Lee et al. (2007). 3. Application and discussion To validate the constitutive equations developed in the companion paper (and summarized in the previous section), the AZ31B magnesium alloy sheet was employed for the measurement and prediction of springback for two cases: (1) the unconstrained cylindrical

76

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

bending test and (2) recently developed novel draw bend test. The simulated results by FE method were compared with those by experiments. Also, the simple bending theory incorporating asymmetry in tension and compression was developed to explain a unique springback phenomenon in magnesium alloy sheets. 3.1. Material properties Although two model materials, O-tempered AZ31B and as-received AZ31B, were characterized in the companion paper, only as-received AZ31B magnesium alloy sheet was considered in the present paper due to the availability of sample size for the springback tests. AZ31B is one of the most widely available of the magnesium grades and is commonly used as an alternative to aluminum alloys due to its high strength to weight ratio. The application area of AZ31B covers from the aerospace to general commercial applications. The chemical composition of the current AZ31B alloy sheet is Mn = 0.2, Al = 3.0%, Zn = 1.0%, Si = 0.1% and Mg = balance in weight percent and the density is 1700 kg/m3. The modulus of elasticity and Poisson’s ratio are 45 GPa and 0.33, respectively, and the initial yield stresses in the rolling direction are 220 MPa in tension and 120 MPa in compression. Based on the developed constitutive theory characterizing asymmetric yielding and hardening, each material parameter is obtained by the directional tension and compression tests and continuous tension tension/compression tests. The basic mechanical properties and constitutive parameters used in Eqs. (10)–(14) of the previous section are summarized in Tables 1 and 2, respectively. Note that the values between measuring points with different initial gap distances are obtained by linear interpolations and for the exterior range the linear extrapolations are used. For the isotropic hardening rate, two different values 0.88 and 0.81 are used for the slip mode and twin mode, respectively, and the instant shrinkage ratio k = 0.54 is used. Fig. 2a shows the correlations between calculated stress–strain curves and their measurement counterparts for the uni-axial tension–compression–tension (T–C–T) tests with various pre-strains. Also, in Fig. 2b the initial yield function calculated from the five parameters is shown. In general, the developed model could reproduce the experimental behavior with good accuracy and therefore the parameters will be adopted for the subsequent numerical applications. In addition to the developed constitutive model, the (conventional) isotropic hardening model and the Chaboche type combined isotropic-kinematic hardening model (hereafter, it is denoted as ‘‘the Chaboche IK model”) are also compared. For simplicity, isotropic von Mises yield function is considered. In terms of the Chaboche IK model, the detailed description of the model can be referred from Chaboche (1991) and Chung et al. (2005). Since the isotropic hardening model is described by the single hardening curve for the reference deformation state, two stress–stress curves are considered: the reference state with

Table 1 Mechanical properties of AZ31B Direction from RD (°)

Young’s modulus (E)

Poisson’s ratio

Yield stress in tension (MPa)

Yield stress in compression (MPa)

0 45 90

45 GPa

0.33

220 210 250

120 – 140

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

77

Table 2 Parameters for the hardening and yield function

din (MPa) 120 130 138 143 208

p 1.32 aT (MPa) 190 200 185 240 480

q 0.32 bT (MPa) 130 130 123 125 110

din (MPa)

aUT (MPa)

bUT (MPa)

350 380 392 400 440

980 1100 1500 1720 4500

147 130 114 94 90

aS (MPa)

bS (MPa) 114

186

b2 1.14 cT

b3 0.99 dT

b4 0.94 y0,T

x0,T

2.4 2.5 2.5 2.3 2.5

0.012 0.012 0.012 0.012 0.012

0 0 0 0 0

0.065 0.06 0.06 0.063 0.045

cUT

dUT

y0,UT

x0,UT

0.006 0.005 0.002 0.002 0.003

0 0 0.03 0.03 0.03

0.08 0.07 0.053 0.055 0.015

cS

ml,T

ml,S

k

R0 (MPa)

15.4

0.81

0.88

0.54

520

0.04 0.04 0.041 0.038 0.043

(1) stress–strain curve in tension (isotropic-tensile model) and (2) stress–strain curve in compression (isotropic-compressive model) in the rolling direction. The following two equations for the two isotropic assumptions were successful to fit the stress–strain curves. For the isotropic-tensile model, iso ¼ 191:5 þ 114:8ð1  e16:14e Þ ðMPaÞ r

ð15Þ

and for the isotropic-compressive model, iso ¼ 98:8 þ 10:9e32:0e ðMPaÞ: r

ð16Þ

Regarding the Chaboche IK model, the procedure is somewhat complex because of the intrinsic exponential nature of the back stress (Chaboche, 1991; Chung et al., 2005). The material parameters for the back stress increment are obtained considering the slopes of forward and reverse loading curves at the instance of unloading. However, the slope (or hardening) of reverse loading in real stress–strain curves of magnesium alloy is (almost) negligible without fast transient behavior. When the isotropic and kinematic hardening parts are separated from the measured stress–strain curves as shown in Fig. 3, the back stress shows saturation during early strain. However, the hardening parameters obtained from the stress–strain curves of magnesium alloys cannot properly fit the behavior of the back stress. Therefore, in the present paper, the hardening parameters for the Chaboche IK model are assumed constant so that the exponential type back stress evolution can be analytically fitted (Khan and Huang, 1995). The resulting hardening equations are as following. For the isotropic part, iso ¼ 148 þ 88:4ð1  e13:3e Þ ðMPaÞ r

ð17Þ

and for the kinematic hardening  a ¼ 72:2ð1  e113:17e Þ ðMPaÞ:

ð18Þ

78

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

a

400 300

Stress (MPa)

200 100 0 -100 -200 -300 -0.10

-0.05

0.00

0.05

0.10

0.15

0.20

True strain 2.0

Normalized Yield Stress at TD

b

1.5

1.0

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0

0.5

1.0

1.5

Normalized Yield Stress at RD Fig. 2. Examples of (a) correlation between experimentally measured consecutive tension–compression tests and corresponding simulations with various prestrains, (b) asymmetric/anisotropic yield stress calculated by Eq. (12).

Note that only tensile stress–strain curve is used for the characterization of the Chaboche IK model due to the unusual shape of compressive stress–strain response. 3.2. Limitations of conventional hardening models Most of automotive sheets show lower yield stress when the loading is reversed, which is so-called as the ‘‘Bauschinger effect”. The Bauschinger effect has been successfully represented by assuming the initial yield surface to translate: the kinematic hardening model. By combining translation and expansion of the yield surface during the plastic deformation, the Chaboche type models could account for the Bauschinger effect and the transient behavior of commonly used sheet metals. Here, in this paper, the limitations of the

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

79

300

Stress (MPa)

250

Isotropic hardening (Measured)

iso

148 88.4 1 e

200

13.3 p

Isotropic hardening (Fitted)

150 100 50 0 0.00

72.2 1 e Kinematic hardening (Measured) 0.05

113.17

p

Kinematic hardening (Fitted)

0.10

0.15

Plastic Strain Fig. 3. Decomposition of tensile stress-plastic strain curve into the kinematic hardening and isotropic hardening parts.

Chaboche type kinematic hardening model for the prediction of stress–strain response of magnesium alloy sheets are explored. Besides the intrinsic exponential nature of the hardening of back stress discussed in the previous section, the Chaboche IK model has another drawback for describing the stress– strain curve of magnesium alloy sheet: the reverse hardening curve converges to the forward curve in the tension–compression test as shown in Fig. 4. This is briefly explained as following. Fig. 4 shows usual hardening curve of tension–compression test where the reverse loading curve is re-plotted into the first quadrant. For the forwarding curve shown in this

Fig. 4. Schematic view of the loading–unloading-reverse loading curve calculated by the Chaboche type isotropic-kinematic hardening model.

80

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

figure, the limit value of the curve is the sum of limit values of isotropic and kinematic parts. The limit value of back stress becomes applying the initial values e0 ¼ 0,  aðe0 Þ ¼ 0 and the L’Hospital theorem. R e a ðeÞ da1 e2 de ea2 ðeÞ  h1 ðeÞ h1 ðeÞ : ð19Þ lim  ¼ lim ¼ lim aðeÞ ¼ lim 0 a ðeÞde e!1 e!1 e!1 h2 ð e2 eÞ  ea2 ðeÞ e!1 h2 ðeÞ where h1 ¼ d a1 =de and h2 ¼ da2 =de. When the loading direction is reversed, the direction of the back stress is also reversed. If the effective strain is denoted as eO when the back stress passes the initial position ð~ a ¼ 0Þ, then the equation h1  h2 ~a ¼ d~a=de is also valid for the back stress from this position and Z e R e  R e  O h2 ðeÞde h2 ðeÞde O O ~ e e h1 ðeÞde þ ~ aðe Þ : ð20Þ aðeÞ ¼ e e eO eÞ . Since the value ~ aðeO Þ vanishes, the limit value of Eq. (20) also becomes lim aðeÞ ¼ lim hh12 ð e!1 e!1 ðeÞ Therefore, the reverse loading predicted by the Chaboche IK model eventually converges to the forward hardening curve which is frequently expressed by the pure isotropic hardening model.

3.3. Unconstrained cylindrical bending 3.3.1. Problem description and characterization of material properties The unconstrained cylindrical bending test (UCB test) proposed for the benchmark of Numisheet (2002) was chosen to evaluate the springback of magnesium alloy sheet with different constitutive equations. The test has been frequently taken by many researches on the forming and springback of sheet metals because it involves simple die geometry but complex contact condition, thus the numerical stability and accuracy of the employed finite element analysis can be conveniently evaluated. The geometry with dimensions and the corresponding finite element mesh of the test are shown in Fig. 5a and b, respectively. The initial rectangular blank sheet is 120 mm (length)  30 mm (width) and the sheet thickness is 2 mm. The rolling direction coincides with the length direction. Four node shell element with reduced integration (S4R) and nine integration points through the thickness is utilized for the descretization of blank sheet, while three-dimensional rigid-body element is used for the tools. The implicit finite element code, ABAQUS/ Standard with the user-defined material subroutine UMAT is used in the simulations for both forming and springback (ABAQUS, 2006). The friction coefficient is assumed constant with 0.1. 3.3.2. Finite element analysis The deformed shapes after forming and springback are captured from the experiment in Fig. 6 and the simulated results were compared in Fig. 7 and Table 3, respectively. In Fig. 7a, the forming behavior is similar for all the hardening models although magnified figure shows the smallest angle (equivalent to the largest forming angle) with the isotropic-compressive model and the largest with the isotropic-tensile model. The result with the Chaboche IK model is almost the same as the one by the isotropic-tensile model, which is expected by the fact that the two models share the same hardening curve during the bending (The slight difference comes from the fitting error for the same experimental

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

81

a Punch 23.5

Sheet .0

4.0

25.5

50

Die

b

Fig. 5. (a) Initial set-up with dimensions and (b) finite element mesh of the unconstrained cylindrical bending test proposed in Numisheet (2002).

curve). Note that the formed shape by the experiment is omitted in Fig. 7a, while averaged formed angle of the four simulations is used because of the difficulty in measurement. The magnitude of springback is defined by the angle difference before and after springback as schematically illustrated in Fig. 8. The simulated results in Fig. 7b and c show that the largest springback is attainable by the isotropic-tensile model, while the lowest with the isotropic-compressive model. Also, like the forming result, the result by the Chaboche IK model is almost the same as the one by the isotropic-tensile model by the same reason explained above. The prediction with the present constitutive equations for AZ31B alloy sheet is in-between the two isotropic cases and shows the best agreement with the measurement. The difference may come from the asymmetry of the yield and hardening behavior. The longitudinal stresses at the integration points close to the top and bottom surfaces of sheet at the end of forming are schematically compared in Fig. 9. The material point near the bottom surface of the blank experiences tension during the bending, and vice versa. Therefore, the stress states at the beginning of springback by removing tool surface become squared points A and B in Fig. 9. The stress values of the two points are (almost)

82

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Fig. 6. Deformed blank shape after (a) forming and (b) springback.

symmetric for the isotropic models and the Chaboche IK model except small deviations which are due to the small tensile and/or shear force by the friction between blank and tool die (if the pure bending is attainable, the isotropic models should result in perfect symmetric stress distribution through the thickness). In the tensile region (or bottom surface), the stress by the present asymmetric model is the same as that by the isotropic-tensile model and the Chaboche IK model, while the latter two models quite overestimate in the compressive region (or top surface). Thus, the isotropic-tensile model and the Chaboche IK model predict the largest springback. Similar explanation can be made for the lowest springback with the isotropic-compressive model by this simple one-dimensional analogy. On the other hand, the stress (therefore, the bending moment) with the developed asymmetric constitutive equations (present model) is averaged value of tensile and compressive flow stresses through the thickness. Note that the springback amount is proportional to the magnitude of bending moment before springback for the elastic unloading. To investigate detailed deformation state by the developed asymmetric constitutive laws during the forming and springback, three deformation modes – slip, twin and untwin – at two outer integration points of the blank are shown in Fig. 10. Here, the deformation

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

a

50

45 Chaboche IK model Isotropic-Tensile

40

Y-Coord. (mm)

83

40

35 30

30

Present Model

25 Isotropic-Compressive 20 26

20

27

28

29

30

10 0

-40

-20

0

20

40

X-Coord. (mm)

b

50 Present Model Elastic/Plastic Unloading

Y-coord. (mm)

40 IsotropicCompressive

30

Present Model Elastic Unloading EXP.

20 IsotropicTensile

Chaboche IK model

10 0

-40

-20

0

20

40

X-coord. (mm)

c

50 Present Model Elastic Unloading

Y-coord. (mm)

40

EXP.

Isotropic-Compressive

30

Present Model Elastic/Plastic Unloading

20

Chaboche IK model Isotropic-Tensile

10 0 20

25

30

35

40

45

X-coord. (mm) Fig. 7. Deformed shapes of blank sheets with different constitutive models after (a) forming, (b) springback and (c) enlargement showing detailed springback profiles shown in (b).

mode ‘imode-1, 2 and 3’ denotes slip, twin and untwin modes, respectively. The figure shows that the elements on the bottom side undergo only slip mode during the whole process. This means that the initially originated ‘tensile’ deformation continues until the forming is finished and the unloading during the springback is ‘elastic’. On the other hand, most of the elements on the top side experience ‘twin’ mode during the whole forming stage but the deformation mode is changed to ‘untwin’ at a certain point during the springback, which means ‘elastic–plastic’ unloading occurs during the springback. The ‘asymmetric’ unloading behavior is also originated from the asymmetric stress–strain response of AZ31B magnesium alloys. Especially, the instant shrinkage of elastic region during

84

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Table 3 Deformed angles after drawing and springback with different models

Isotropic-tensile model Chaboche IK model Isotropic-compressive model Present model Present model w/elastic unloading Exp. a

hf

hs

Dh(=hs  hf) (°)

37.64 37.46 35.0 36.78 36.78 36.73a

66.78 66.16 48.66 61.44 57.94 59.34

29.14 28.70 13.66 24.66 21.16 22.61

Averaged value of the five simulation models.

Fig. 8. Definition of springback angle.

the twin deformation plays an important role on initiating early plastic deformation during the unloading deformation. From the deformation process observed in Fig. 10, it is inferred that the plastic unloading might influence the magnitude of springback although most of the conventional analyses have adopted elastic unloading scheme. Therefore, two different unloading schemes are simulated with the developed constitute equations; elastic unloading and elasto-plastic unloading. Fig. 7b and Table 3 show that the springback decreases when only elastic unloading is allowed. The early re-yielding during the equilibration of moment causes larger springback, which was previously reported in the U-channel draw bend test with pure kinematic hardening (Pouboghrat et al., 1998). Besides the deformed shape, Fig. 11 shows the simulated punch forces of the UCB test and their comparison with the measurement. The prediction using the proposed (asymmetric) constitutive equations falls within the two bounds of simulated results with the two isotropic cases, and shows good agreement except high oscillatory behavior. Note that high oscillations that initiate around 15 mm punch displacement may come from the severe change of contact status with somewhat coarse mesh used in the current finite element code due to the unconstraint boundary condition (without blank holder).

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

B

85

A

Asymmetric model B

A Isotropic model Chaboche s Kinematic model

Isotropic model

Fig. 9. Schematic view of stress–strain curves in tension–compression and compression–tension under the present model, isotropic-tensile model and the Chaboche type isotropic-kinematic hardening model.

3.4. Draw bend test 3.4.1. Description of the draw bend test The draw bend test considered as a second example is schematically shown in Fig. 12, which has advantages over the conventional U-shaped forming type test being widely used for the study of springback. For example, the direct control of constraint force and longer draw distance are attainable. Also, the test does not have conventional die and punch so that the analysis does not involve complex process effects such as friction. Detailed descriptions of the draw bend test are well reported elsewhere (Carden et al., 2002) and only brief explanation is followed. As shown in Fig. 12a, a restraining back force is prescribed to the upper grip and the sheet is drawn over the cylindrical tool surface by the lower grip with constant speed. Therefore, the sheet blank undergoes tensile force by the prescribed back force, bending over the tool radius and reverse bending (straightening) at the end of tool travel. The springback occurs by releasing the lower grip. Fig. 12b shows a typical schematic view of the blank sheet before and after springback. In the figure, the springback in the tool region represents the change of angle from the formed angle over the cylindrical tool (in this case, p2) to that of the relaxed sheet (hb). The angle change Dhb is  p p r Dhb ¼ hb  ¼  1 ; ð21Þ 2 2 r0 where r and r0 are radii of curvature of the tool region before and after springback, respectively. Conversely, the springback in the sidewall region represents the change

86

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Fig. 10. Deformation modes after forming and springback at the integration points (a) at the bottom and (b) at the top surface. Note that imode = 1, 2 and 3 denote slip, twin and untwin modes, respectively. 1000 Isotropic-Tensile Model

Punch Force (N)

800

Chaboche IK model Present Model

EXP.

600

400 Isotropic-Compressive Model

200

0 0

5

10

15

20

25

30

Punch Dislpacement (mm) Fig. 11. Variations of punch force vs. punch displacement in the unconstrained cylindrical bending tests.

of curvature from zero to a finite value. By assuming a constant curvature, the angle change becomes d Dhs ¼ hs ¼ 00 ; ð22Þ r where r00 is the radius of curvature of the region after springback and d is drawing distance. For simplicity, following the suggestion made previously (Carden et al., 2002; Li,

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

87

Fig. 12. (a) Geometry of the draw bend test with (b) deformed shape before and after springback.

2006), the sum of the two angle changes Dh = Dhb + Dhs is adopted as a total springback amount. Among three major parameters – back force, the ratio of the tool radius over the sheet thickness (R/t ratio) and friction coefficient, the restraint back force normalized by the tensile yield stress and R/t are known to be the most sensitive to the springback amount. Therefore, to exclude the effect of friction on the springback, the friction coefficient is minimized by rotating the cylindrical tool with the same speed of drawing.

88

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

3.4.2. Experimental results A rectangular sheet strip with 635 mm (length)  50.8 mm (width) is installed in the draw bend test machine and drawn d = 127 mm with constant speed 50 mm/s. The springback angles for three different thicknesses are measured with increasing normalized back forces, shown in Figs. 13–15. Note that the restraining force is normalized by the initial tensile yield stress. The measurements stop when the restraining force leads to fracture. The fracture points are indicated in each figure with arrows. The geometry of springbacked samples and magnitude of springback angle are measured manually in approximately 10 s after releasing from the grips. The figures show that four features are noted in the behavior of forming and springback of AZ31B magnesium alloy sheet in room temperature.  (Total) springback angle and angle of sidewall curl decrease as the restraint back force increases. The springback angle of the sidewall region dominates over the springback in the cylindrical die, thus the pattern of total springback decrease with respect to the back force is similar to that of the springback of the sidewall. Also, the angle changes at the roller and sidewall regions have different signs.  For the same sheet thickness, the springback increases as the R/t ratio decreases. However, the springback in the roller region is less sensitive to the R/t ratio than those at the sidewall region.  The critical restraint force at the initiation of fracture decreases as the R/t ratio decreases. In the figures, the critical normalized back forces are Fb = 0.8 for R/t = 9.5, Fb = 0.4 for R/t = 5.64 with 2 mm thickness, Fb = 1.0 for R/t = 19.1, Fb = 0.8 for R/ t = 11.1 with 1 mm thickness, and Fb = 1.2 for R/t = 38.1, Fb = 1.0 for R/t = 22.2 with 0.5 mm thickness.  Total springback angle in the lower ranges of back force slightly increases and then decreases rapidly after some point. This point was explained that the most section through the thickness becomes plastic and sometimes anticlastic curvature occurs (Wang et al., 2005).

Fig. 13. Springback angles with varying normalized back force in the draw bend test with 2 mm AZ31B magnesium alloys sheets: (a) R/t = 9.5, (b) R/t = 5.6.

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

89

Fig. 14. Springback angles with varying normalized back force in the draw bend test with 1 mm AZ31B magnesium alloys sheets: (a) R/t = 19.1, (b) R/t = 11.1.

Fig. 15. Springback angles with varying normalized back force in the draw bend test with 0.5 mm AZ31B magnesium alloys sheets: (a) R/t = 38.1, (b) R/t = 22.2.

Among the above features observed in the present experimental study, the first three were commonly observed in the springback of other metal sheets such as aluminum alloys and conventional high strength steels (Geng and Wagoner, 2002). However, there are two unique features observed only in the springback of AZ31B magnesium alloys. The first is that the slope of springback reduction after the critical back force is relatively sharper than that shown in aluminum alloys. Secondly, more important feature is that AZ31B magnesium alloy has increasing or (almost) constant springback in the region of lower back force and then decreases again, while aluminum alloy sheet has monotonically decreasing springback in the whole region until the samples are fractured. In Fig. 16, the two main features observed in the springback of magnesium alloys sheet in the draw bend test are compared with the springback of aluminum alloy 6022-T4 which was obtained from the previous report with similar R/t ratio (R/t = 10.5) by one of the present authors (Carden et al., 2002).

90

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Fig. 16. Comparison of springback with varying normalized back forces: (a) Al 6022-T4, R/t = 10.5 (Carden et al., 2002), (b) Magnesium alloy AZ31B, R/t = 11.1.

3.4.3. Finite element simulations of the draw bend test Similar to the UCB simulation, ABAQUS/Standard with the user-defined subroutine UMAT is used. Taking into account the geometric symmetry of the process, only half side of the blank is simulated. Based on pre-performed sensitivity results, 2400 shell elements and nine integration points are selected as optimum numbers for the simulations. Also, for the accuracy of springback, five times smaller elements are used in the roller and sidewall regions as shown in Fig. 17. Although it is very difficult to measure precise friction coefficient between AZ31B sheet and cylindrical tool, it is reasonably assumed constant with 0.05, which is small enough value by considering rotating roller in accordance to the drawing velocity. The simulated springback results with the developed model are compared with the two isotropic models, the Chaboche IK model and experiments obtained for two different roller radii. The simulated springback profiles and corresponding experiments for three different normalized back forces are shown in Fig. 18 in case of R/t = 19.05 as an example. The figures show that the springback profiles with the developed constitutive model reasonably well predict the deformed shapes by experiments. Conversely, the isotropic-tensile model and the Chaboche IK model quite overestimate the deformed shapes especially when the restraining back force is low. When the developed model is compared with the isotropic-compressive model, the deformed shapes seem similar at first glance. However, the simulated springback profiles at the roller region reveal that the present model has much better prediction capability than that with isotropic-compressive model (this will be shown in the later figures). Especially, the present model results in high springback angle in the roller region, which is one of unique feature of magnesium alloy sheet shown in Fig. 16. Note that the simulation with the isotropic-compressive model for Fb = 1.0 could not be completed because the force is normalized by the tensile yield stress and such too large force leads the material element to be numerically instable. The springback angle defined in Fig. 12b is calculated and compared for the four constitutive models with varying normalized back force in Fig. 19. The prediction presented in the figures shows that all four constitutive equations can reproduce the main features of springback observed in conventional automotive sheet metals as explained in the previous

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

91

Fig. 17. Initial and deformed finite element meshes of the 2D draw bend test.

Section 3.4.2. In other words, the springback decreases as the restraining back force and the R/t ratio increase. Also, the results show the slope change when the normalized back force is around 0.6, where the sheet might become nearly plastic through the thickness and the anticlastic deformation remains after springback, thus the springback is prohibited significantly (Carden et al., 2002). When the calculated results are compared in more detail, the isotropic-tensile model shows the largest springback compared with the other three constitutive models. Especially, in the lower range of the back force, the relative errors from the measured angles become over 100% for two R/t ratios. However, as the normalized back force becomes larger, the errors are diminished and converge to those of the other three cases. As for the Chaboche IK model, the predicted springback angles are lower than those by the isotropic-tensile model but still have large scatter from the measured springback angles, especially at the lower range of back force. The similarity as the isotropic-tensile model in predicting the springback angle is well explained by Fig. 9, which shows the convergence of the Chaboche IK model to the isotropic model as the deformation proceeds. In case of isotropic-compressive model, the prediction of springback angle is good enough in the lower back force region compared with experiments. But, as the back force increase, the calculated springback angles deviate from those by experiments. The slopes of angle decrease are smaller than those by the other two models. Contrary to the three conventional models, the developed constitutive model considering asymmetry and anisotropy of AZ31B sheet reasonably well predicts the springback with respect to the restraining back force for two R/t ratios. In addition to the capability in predicting common features of springback shown in the conventional symmetric materials, the new constitutive model can represent almost constant springback angles before the sharp decrease

92

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Fig. 18. Comparison of deformed shape after springback with various normalized back force and R/t = 19.05. (a) Fb = 0.2, (b) Fb = 0.65 and (c) Fb = 1.0. *Converged solution for isotropic-compression with Fb = 1.0 was not obtained.

at certain back forces. Summing up the above results and comparisons between four different constitutive models and experiments, the consideration of asymmetry and anisotropy by developing proper constitutive equations may be important in predicting springback of magnesium alloy sheet. Finally, to validate whether the developed constitutive model is properly implemented in the finite element program, the deformation status of two representative elements are

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

degrees)

a

100

AZ31B R/t=19.05 =0.05

Chaboche IK model 80 Iso-Ten.

60

Springback Angle

93

Measured 40 20 Iso-Comp. 0

Present Model

0.2

0.4

0.6

0.8

1.0

Normalized Back Force (Fb)

Springback Angle

degrees)

b

120 AZ31B R/t=11 =0.05

Chaboche IK model 100

Chaboche IK model 80 60 40 Measured

Iso-Comp.

20 0 0.2

Present Model

0.3

0.4

0.5

0.6

0.7

0.8

Normalized Back Force (Fb) Fig. 19. Comparison of springback angles predicted by four different models with various normalized back force for (a) R/t = 19.05, and (b) R/t = 11.1.

investigated. The two elements 4255 and 4200, which represents the elements in the roller region and sidewall region at the end of drawing, respectively, are selected and shown in Fig. 20. The element 4255 in the roller region experiences applied tensile force by the upper grip and bending strain by the constant curvature of the roller. During the traveling along the roller, the element is further strained by the additional tensile forces by the friction and bending energy. The schematic view of strain history in the element is shown in Fig. 21a. On the other hand, the deformation of element 4200 in the sidewall region further experiences unbending (or straightening) at the exit of the roller. During the straightening, the element is further stretched by the unbending energy and the final strain state is shown in Fig. 21b. Therefore, the material points at the top and bottom surfaces experience different tangential strains. For example, the material point A located near the top surface is first stretched by bending at the roller region and then compressed again as the same material point travels along the sidewall region. Fig. 22 shows the trace of deformation modes

94

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

A T0

T0 Mb

B

Mb Tb

Tb Element 4255

Mub

Mub Tub

Tub Element 4200

Fig. 20. Deformation patterns at three typical elements representing flange, roller radius and sidewall regions.

a +

Strain by bending and applied load

=

Additional tensile strain

Total strain of roller region

b +

+

Strain at the exit of roller region

Strain by unbending (or straightening)

=

Additional strain by unbending

Total strain of sidewall region

Fig. 21. Schematic view of tangential strain distributions: (a) tool radius region and (b) sidewall region.

of the two representative elements during the drawing. In Fig. 22a, the deformation modes of the material points A and B which correspond to the top and bottom surface of the element in the roller region (element 4255) are obtained from the drawing simulations. The deformation mode of the point A is ‘imode = 1 (or slip mode)’ following the elastic deformation, while ‘imode = 2 (or twin mode)’ at point B. Therefore, the tangential stress values of the material points at the end of drawing are squared A and B on the solid line

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

95

a B

Imode

2

A 1

0 0

Element 4255 imode def. mode 1 slip 2 twin 3 untwin 20

40

60

80

100

120

Drawing distance (mm)

b

B 3

Imode

A 2

Element 4200 imode def. mode 1 slip 2 twin 3 untwin

1

0 0

20

40

60

80

100

120

Drawing distance (mm) Fig. 22. History of deformation modes at two typical elements during the drawing: (a) element 4255 on the roller region, and (b) element 4200 on the sidewall region.

shown in Fig. 9. When the stress values are compared with those by isotropic-tensile model and the Chaboche IK model, the stress state of the point A has (almost) the same value.1 However, the stresses at the point B become quite different because of the asymmetry. Therefore, the stress predicted by isotropic-tensile model (squared B on the dashed line in Fig. 9) overestimates the real stress–strain response of AZ31B sheet especially under the neutral surface. Note that the stress-state by the isotropic-tensile model and Chaboche IK model are almost the same as shown in Fig. 9. Similarly, the deformation modes of the material point A and B of element in the sidewall region (element 4200) are calculated in Fig. 22b. The material point A near the top surface experiences slip (‘imode = 1’) followed by twin (‘imode = 2’) at the drawing distances around 20 mm and 40 mm, respectively. On the other hand, the material point B near the bottom surface takes the twin mode first and 1 The distribution of tangential strain with the present asymmetric model is a little different from the isotropic models even for the pure bending deformation. This is because force and moment equilibrations lead to shift of neutral line considering asymmetric stress–strain response.

96

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

then untwin (‘imode = 3’) during the drawing. In this case, as illustrated in Fig. 9, the stress values are located circled A and B just before the unloading. Therefore, the tangential stresses by the two models are considerably different for both top and bottom surfaces, which leads significant overestimation of springback angles in the sidewall region predicted by isotropic-tensile model. Contrary to the roller region, the stress state by the Chaboche IK model is also different with the one by the isotropic-tensile model. This is because the Chaboche model can represent the transient and Bauschinger effect of the reversed loading curve, which explains lower stress state than the isotropic-tensile model in Fig. 9. However, even though the Chaboche IK model can introduce more realistic unloading behavior, the constitutive equation cannot represent the unusual asymmetry in tension–compression test. To back up the intuitive interpretation made so far, the equivalent stresses are plotted for four different constitutive models in Fig. 23. Only the integration point which has the largest compressive tangential strain is considered. As explicitly shown in Fig. 23, by the isotropic models and the Chaboche IK model, the stresses always increase as the material travels along the roller. In other words, the stress value at the sidewall region is much larger than that in the roller region because of the expansion of yield surface and symmetric nature. On the contrary to the conventional models, the novel asymmetric model based on the two-surface plasticity theory shows decrease in stress value at the sidewall region after traveling tool radius. This is because the present model can consider the asymmetric stress–strain response as well as instant shrinkage of elastic region observed in real AZ31B magnesium alloy sheets. 3.5. Bending theory with tension/compression asymmetry To understand the unusual phenomenon of springback in magnesium alloy sheets, especially the reverse effect of the constraint tensile force on springback, the simple but analytical bending theory with applied tensile force is extended to the asymmetric materials. The common assumptions for the classical bending theory are employed. For example, the plane sections remain planar and no change occurs in sheet thickness. The model problem is schematically illustrated in Fig. 24a, in which the part is bent to a starting bend radius, R under the prescribed moment, M and superimposed tensile force T during bending. Elastoperfect plasticity with different tensile and compressive yield stresses is considered (Fig. 24b). The stress distribution for the elastic and plastic regions becomes 8 T r ðe > eTy Þ; > < y r ¼ Ee ðeCy 6 e 6 eTy Þ; ð23Þ > : C C ry ðe < ey Þ; where rTy and eTy are stress and strain at tensile yield, respectively, while rCy and eCy are stress and strain at compressive yield. From the stress–strain curve for AZ31B magnesium alloy, the condition jrCy j < rTy , C jey j < eTy can be assumed. The total strain generated through the cross-section is e ¼ eb þ et ¼ jðy  y 0b Þ þ et ¼ jy þ et  jy 0b ¼ jy þ e0t ;

ð24Þ

where eb is the strain by bending with curvature, j = 1/R, and et is the membrane strain by the prescribed tensile force. Since the present model considers high asymmetry in tension and compression, the neutral surface during the pure bending may not be preserved and

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

97

Fig. 23. Distribution of equivalent stresses for three constitutive models: (a) present model, (b) isotropic-tension model, and (c) isotropic-compression model, and (d) Chaboche type isotropic-kinematic hardening model.

the position of neutral surface itself is not a convenient parameter to calculate stress and strain distributions. However, note that in the elastic range the neutral surface is preserved at the centerline. Therefore, the pseudo-neutral surface, y0b, where the strain equals the membrane strain, et is defined as shown in Fig. 25. When the strain distribution is not sufficient to initiate yield from both tensile and compressive regions, the stress–strain through the thickness is purely elastic as shown in Fig. 26a. From the prescribed tensile force normalized by sheet width (w), T , the tensile T strain, et ¼ Eh and the shifted neutral surface, y0 is obtained with vanishing strain and static pseudo-neutral surface (y0b = 0) as following:

98

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Fig. 24. Schematic views of (a) simple bending geometry and (b) elasto-perfect plasticity with tension– compression asymmetry.

Fig. 25. Assumed strain decomposition into the bending strain and uniform tensile strain.

y0 ¼ 

T : jEh

Corresponding bending moment becomes Z h=2 Z h=2 1 M¼ rðyÞy dA ¼ E  ðky þ et Þ  y  b dy ¼ Ebh3  j: 12 h=2 h=2

ð25Þ

ð26Þ

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

99

Fig. 26. Strain and stress distributions through the thickness for various deformation cases: (a) elastic case and (b) plastic in both tensile and compressive region.

The initial yield initiates either in the compressive or tensile regions depending on the rT applied tensile force and the ratio jrCy j. When yielding in the compressive region is initiated y first at y =  h/2, the moment–curvature curve meets elastic limit and the strain, curvature and moment at this instant are !   h h et C C ey ¼ e y ¼  ¼ jy y 0b þ  C : ð27Þ 2 2 jy ! et  eCy rCy 1 2 C C and M y ¼ Ebh et  jy ¼ : ð28Þ 6 y 0b þ h=2 E Similarly for the first yield in tensile region at y ¼ h2 when the tensile force is sufficiently large, the strain and curvature at the initiation of tensile yielding becomes !   eTy  et h h e t  y 0b þ T ; and jTy ¼ eTy ¼ e y ¼ : ð29Þ ¼ jTy 2 2 jy h=2  y 0b The critical moment is obtained by substituting Eq. (29) into Eq. (26). When the curvature keeps increasing up to jCy < j 6 jTy , the deformation is sufficient to make the compressive region plastic, while the tensile region is still elastic. The posi-

100

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

tion of updated pseudo-neutral surface is calculated from the force equilibrium with prescribed tensile force, which is in general nonlinear equation with respect to. Similar procedure can be used for the case jTy < j < jCy , in which the plastic stress occurs from the tensile region, while the compressive region is still elastic. Finally, when both the tensile and compressive regions experience plastic deformations for a given curvature as shown in Fig. 26b, the position of new pseudo-neutral surface is obtained by the equilibrium equation and the resultant bending moment can be calculated. For example, the equations of the force equilibrium and the bending moment for the final stage are T ¼

Z

yC h2

rCy dy þ

Z

yT yC

Eðjy þ e0t Þ dy þ

Z

h 2

yT

rTy dy

   h h 1 T T C C  y ry þ ry y þ þ Eðy T  y C Þðjðy T þ y C Þ þ 2e0t Þ; ¼ 2 2 2 Z yC Z yT Z h2 M¼ rCy  y dy þ Eðjy þ e0t Þ  y dy þ rTy  y dy 

h2

"

yC

ð30Þ

yT

 2 # 1 C C 2 h 1 1 3 3 2 2 ¼ ry ðy Þ  þ Ej½ðy T Þ  ðy C Þ  þ Ee0t ½ðy T Þ  ðy C Þ  2 2 3 2 "  # 2 1 T h 2  ðy T Þ ; þ ry 2 2

ð31Þ

where yT and yC are the coordinate of yield points in tensile and compressive regions, respectively. For the unloading by releasing the bending moment after initial bending, which is equivalent to the springback, the elastic slope in the moment–curvature curve is obtained by Eq. (26) Dj=w ¼

Mu 1 Eh3 12

ð32Þ

where M u is moment at the instance of unloading for a given curvature. Therefore, under this elastic unloading scheme the amount of springback (or curvature change) is proportional to the magnitude of M u . The developed simple analytical theory was applied to the AZ31B magnesium alloy sheet considered in the paper. The unit thickness and width are assumed, while the asymmetry in tension and Tcompression was focused on by setting rTy ¼ 220 MPa and rTy ¼ 120 MPa r (therefore, jrCy j ¼ 1:83) with constant Young’s modulus E = 45 GPa. The pseudo-neutral y surface by the equilibrium equation is solved by adopting bi-section method. Fig. 27a shows the moment–curvature curves with various tensile forces normalized by the initial tensile yield stress. It is notable that the moment curve increases as the tensile force increases up to T/Ty = 0.2, which is not known results for the conventional symmetric materials. When the force further increase to T/Ty = 0.5 the moment significantly drops, which means abrupt springback reduction. Therefore, these results reasonably support the experimental and simulation results in the previous sections where the reverse effect of constraint force is presented. For comparison, two other cases with symmetric material

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

101

60

Moment (N.mm)

50 40

Asymmetric Pefect Plasticity ( yT > | yC|)

30

T/Ty=0 T/Ty=0.05 T/Ty=0.1 T/Ty=0.2 T/Ty=0.5

20 10 0 0.00

0.02

0.04

0.06

0.08

0.10

Curvature (mm-1)

Moment (N.mm)

80

60

Symmetric Perfect Plasticity ( yT=| yC|)

40

T/Ty=0 T/Ty=0.05 T/Ty=0.1 T/Ty=0.2 T/Ty=0.5

20

0 0.00

0.02

0.04

0.06

0.08

0.10

Curvature (mm-1) 60

Asymmetric Perfect Plasticity ( yT<| yC|)

Moment (N.mm)

50 40 30

T/Ty=0 T/Ty=0.05 T/Ty=0.1 T/Ty=0.2 T/Ty=0.5

20 10 0 0.00

0.02

0.04

0.06

0.08

0.10

Curvature (mm-1 ) Fig. 27. Moment–curvature curves with various tensile constraint forces: (a) asymmetric perfect plasticity with higher tensile yield stress, (b) symmetric perfect plasticity, and (c) asymmetric perfect plasticity with higher compressive yield stress.

102



rT y

jrC yj

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104



¼ 1 and the asymmetric material with larger yield stress in compression



rT y jrC yj

¼ 0:55



are also considered and their results are shown in Fig. 27b and c, respectively. Contrary to the asymmetric case with higher yield stress in tension, the moment monotonously decreases as the tensile force increases. Therefore, from the results by simple analytical model, the slow or reverse change of draw bend springback in lower back force and rapid decrease after certain value can be attributed to the high asymmetry in tension and compression along with higher tensile yield stresses. Of course, the current simple analysis cannot give quantitative evaluations of springback but the method might be effectively used for a rapid estimation of springback in asymmetric materials. 4. Concluding remarks The constitutive model developed in the companion paper (Lee et al., 2007) to consider the asymmetry and anisotropy of magnesium alloy sheet AZ31B was applied to the finite element simulations for the sheet springback and the results were compared with other constitutive equations based on the isotropic models and the Chaboche type isotropickinematic hardening model. Also, the predictions were validated by the experimentally measured springback. Two examples were considered: the unconstrained cylindrical bending test and the novel draw bend test. Besides the numerical works, experimental results were obtained with the efficient draw bend test which corresponds to the drawing of sheet blank over a tool radius and significant springback by the sidewall curl reproducing a typical forming operation. The followings were concluded regarding the springback of AZ31B magnesium alloy sheet and the prediction capability of the proposed constitutive models.  The developed constitutive equations in the companion paper were successfully implemented in the implicit finite element program ABAQUS/Standard via user material subroutine UMAT. Shell elements were used for the sheet blank considering the plane stress formulations.  In the unconstrained cylindrical bending test the magnitudes of springback angles predicted by the isotropic-tensile and isotropic-compressive models showed two extremes, while the present model predicted intermediate springback angle. The results were explained by the 1D analogy of stress states at the end of bending for four different constitutive models. The new constitutive equations could consider both high and low flow stresses in tensile and compressive regions, respectively. The punch force during the bending also showed that the present model predicted intermediate between the values by isotropic-tensile and isotropic-compressive models.  Unloading scheme was deemed to affect the amount of springback especially for the present asymmetric materials. At the integration point near the top blank surface where twin dominates during the bending, untwinning mode by early re-yielding was observed during the unloading because of the reduction of elastic region.  The measurement of springback in the draw bend test showed several features typically observed in the conventional springback tests: reduction in springback angle with higher restraining force and larger R/t ratio. However, in addition to the typical results, AZ31B showed increase or constant springback at lower back force and rapid decrease after normalized back force around Fb = 0.6  0.8.

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

103

 The simulations employing the developed constitutive model agreed reasonably well with experiments. Especially, the results could predict the slow (almost constant) change of springback in lower back force and rapid decrease after certain points.  A simple analytical bending/elastic unloading model adopting asymmetry in tension and compression was developed to calculate moment–curvature relationship with different amount of curvature. From the results, the reverse effect of the constraint tensile force on the springback (that is, increase of springback with higher tensile force) could be explained for only asymmetric material with higher tensile yield stress. Conversely, symmetric material and asymmetric material with lower tensile yield stress than compressive yield stress could not show the unusual springback phenomenon of the magnesium alloy sheets. Acknowledgements This work was partially supported by the SRC/ERC program of MOST/KOSEF (R11– 2005–065) and by the Center for Advanced Materials Processing (CAMP) of the 21st Century Frontier R&D Program in Korea, which are greatly appreciated. References ABAQUS, 2006. User’s Manual for Version 6.5.1, Hibbit, Karlson & Sorensen Inc. Agnew, S.R., Duygulu, O., 2005. Plastic anisotropy and the role of non-basal slip in magnesium alloy AZ31B. Int. J. Plasticity 21, 1161–1193. Amstrong, P.J., Frederick, C.O., 1966. A mathematical representation of the multiaxial Bauschinger effect. G.E.G.B. Report RD/B/N 731. Balakrishnam, V., 1999. Measurement of in-plane Bauschinger effect in metal sheet. MS Thesis. The Ohio State University, Columbus, OH. Bauschinger, J., 1886. On the change of the elastic limit and the strength of iron and steel, by drawing out, by heating and cooling, and by repetition of loading. In: Minutes of Proceedings of the Institution of Civil Engineers with Other Selected and Abstracted Papers LXXXVII, p. 463. Bettles, C.J., Gibson, M.A., 2005. Current wrought magnesium alloys: strengths and weakness. JOM 57, 46. 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. Plasticity 21, 2319. Brown, D.W., Agnew, S.R., Bourke, M.A.M., Holden, T.M., Vogel, S.C., Tome, C.N., 2005. Internal strain and texture evolution during deformation twinning in magnesium. Mater. Sci. Eng. A 399, 1–12. Carden, W.D., Geng, L.M., Matlock, D.K., Wagoner, R.H., 2002. Measurement of springback. Int. J. Mech. Sci. 44, 79–101. Cazacu, O., Barlat, F., 2001. Generalization of Drucker’s yield criterion to orthotropy. Math. Mech. Solids 6, 613–630. Cazacu, O., Barlat, F., 2004. A criterion for description of anisotropy and yield differential effects in pressureinsensitive metal. Int. J. Plasticity 20, 2027–2045. Chaboche, J.L., 1991. Time independent constitutive theories for cyclic plasticity. Int. J. Plasticity 2, 149. Chun, B.K., Jinn, J.T., Lee, J.K., 2002. Modeling the Bauschinger effect for sheet metals-part I: theory. Int. J. Plasticity 18, 569. Chung, K., Wagoner, R.H., 1986. Effect of stress–strain-law transients on formability. Metall. Trans. A 17, 1001. Chung, K., Lee, M.G., Kim, D., Kim, C., Wenner, M.L., Barlat, F., 2005. Spring-back evaluation of automotive sheets based on isotropic-kinematic hardening laws and non-quadratic anisotropic yield functions, part I: theory and formulation. Int. J. Plasticity 21, 861. Dafalias, Y.F., Popov, E.P., 1976. Plastic internal variables formalism of cyclic plasticity. J. Appl. Mech. ASME 98, 645. Doucet, A.B., Wagoner, R.H., 1987. Plane-strain work hardening and transient behavior of interstitial-free steel. Metall. Trans. A. 18, 2129.

104

M.G. Lee et al. / International Journal of Plasticity 25 (2009) 70–104

Drucker, D.C., 1949. Relation of experiments to mathematical theories of plasticity. J. Appl. Mech. 16, 349–357. Duygulu, O., Agnew, S.R., 2003. The effect of temperature and strain rate on the tensile properties of texture magnesium alloy AZ31B sheet. In: Kaplan, H. (Ed.), Magnesium Technology. TMS, San Diego, CA, USA, pp. 237–242. Flores, P., Duchene, L., Bouffioux, C., Lelotte, T., Henrard, C., Pernin, N., Van Bael, A., He, S., Duflou, J., Habraken, A.M., 2007. Model identification and FE simulations: effect of different yield loci and hardening laws in sheet forming. Int. J. Plasticity 23, 420–449. Geng, L., Wagoner, R.H., 2002. Role of plastic anisotropy and its evolution on springback. Int. J. Mech. Sci. 44 (1), 123. Hashiguchi, K., 1981. Constitutive equations of elastoplastic materials with anisotropic hardening and elasticplastic transition. J. Appl. Mech. ASME 48, 297. Hashiguchi, K., 1997. The extended flow rule in plasticity. Int. J. Plasticity 13, 37. Hodge, P., 1957. A new method of analyzing stresses and strains in work hardening plastic solids. J. Appl. Mech. ASME 24, 482. Kalidindi, S.R., 1998. Incorporation of deformation twinning in crystal plasticity models. J. Mech. Phys. Solids 46, 267–290. Khan, A.S., Huang, S., 1995. Continuum Theory of Plasticity. John Wiely & Sons, Inc., NY, USA. Kim, J., Ryou, H., Kim, D.G., Kim, D., Chung, K., Hong, S.-H., submitted. Constitutive law of AZ31 Mg alloy sheet and finite element simulation of three point bending. Int. J. Mech. Sci. Khoei, A.R., Jamali, N., 2005. On the implementation of a multi-surface kinematic hardening plasticity and its applications. Int. J. Plasticity 21, 1741. Krieg, R.D., 1975. A practical two surface plasticity theory. J. Appl. Mech. ASME 42, 641. Lebenson, R.A., Tome, C.N., 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development in polycrystals-application to zirconium alloys. Acta Metall. 41, 2611–2624. Lee, M.G., Kim, D., Kim, C., Wenner, M.L., Wagoner, R.H., Chung, K., 2007. A practical two-surface plasticity model and its application to springback prediction. Int. J. Plasticity 23, 1189–1212. Lee, M.G., Wagoner, R.H., Lee, J.K., Chung, K., Kim, H.Y., 2008. Constitutive modeling for anisotropic/ asymmetric hardening behavior of magnesium alloy sheets. Int. J. Plasticity 24 (4), 545–582. Li, M., 2006. Constitutive modeling of slip, twinning, and untwinning in AZ31B magnesium. Ph.D Thesis, Ohio State University, USA. Lou, X.Y., Li, Min., Boger, R.K., Agnew, S.R., Wagoner, R.H., 2007. Hardening evolution of AZ31B Mg sheet. Int. J. Plasticity 23, 44. McDowell, D.L., 1985. A two surface model for transient nonproportional cyclic plasticity: part 1: development of appropriate equations. J. Appl. Mech. ASME 52, 298. Mroz, Z., 1967. On the description of anisotropic work-hardening. J. Mech. Phys. Mech. 15, 163. Mroz, Z., Niemunis, A., 1987. On the description of deformation anisotropy of materials. In: Boehler, J.P. (Ed.), Proceedings of the IUTAM Symposium: Plasticity and Damage of Anisotropic Solids. Mech. Eng. Press, Grenoble. Numisheet ’2002 Benchmark Problem, 2002. In: Yang, D.Y., Oh, S.I., Huh, H., Kim, Y.H. (Eds.), Proceedings of 5th International Conference and Workshop on Numerical Simulation of 3D Sheet Forming Processes, Jeju, Korea. Ohno, N., Kachi, Y., 1986. A constitutive model of cyclic plasticity for nonlinear hardening materials. J. Appl. Mech. 53, 395. Ohno, N., Wang, J.D., 1993. Nonlinear hardening rules with critical state of dynamic recovery. Part I: Formulation and basic features for ratcheting behavior. Int. J. Plasticity 9, 375. Prager, W., 1956. A new method of analyzing stresses and strains in work hardening. Plastic solids. J. Appl. Mech. ASME 23, 493. Pouboghrat, F., Chung, K., Richmond, O., 1998. A hybrid membrane/shell method for rapid estimation of springback in anisotropic sheet metals. ASME J. Appl. Mech. 65, 671–684. Staroselsky, A., Anand, L., 2003. A constitutive model for hcp materials deforming by slip and twinning: application to magnesium alloy AZ31B. Int. J. Plasticity 19, 1843–1864. Wagoner, R.H., Wang, J.F., Li, M., 2006. Springback. In: Chapter in ASM Metals Handbook. Metalworking: Sheet Forming, vol. 14B. ASM, Materials Park, OH, p. 733. Wang, J.F., Wagoner, R.H., Matlock, D.K., Barlat, F., 2005. Anticlastic curvature in draw-bend springback. Int. J. Solids Struct. 42, 1287–1307. Ziegler, H., 1959. A modification of Prager’s hardening rule. Quart. Appl. Math. 17, 55.