Constitutive modeling of evolving plasticity in high strength steel sheets

International Journal of Mechanical Sciences 107 (2016) 43–57

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

Constitutive modeling of evolving plasticity in high strength steel sheets Zhengyang Cai a, Keshan Diao b, Xiangdong Wu a,n, Min Wan a a b

Beijing University of Aeronautics and Astronautics, School of Mechanical Engineering and Automation, 704#, Beijing 100191, PR China State Key Laboratory of Development and Application Technology of Automotive Steels (Baosteel Group), Shanghai Baoshan 201900, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 29 January 2015 Received in revised form 22 October 2015 Accepted 5 January 2016 Available online 11 January 2016

In this paper, biaxial tensile tests with cruciform specimens were conducted to obtain the experimental plastic work contours of DP590, DP780 and DP980 sheets. It was found that for high strength dual phase steel sheets, the outlines of the experimental work contours in the first quadrant evolve evidently with the increase of plastic dissipation. In order to calibrate the evolving isotropic and anisotropic plastic behaviors simultaneously, the Yld2000-2d yield criterion with variable exponent and coefficients (Yld2000-2d-Var) was introduced. In this model, the exponent m and coefficients αi of the Yld2000-2d was assumed as the function of the equivalent plastic strain. The identification procedures of the introduced parameters were then proposed based on minimizing the difference between the experimental and the model predicted yield points, using Levenberg–Marquardt optimization method. Afterwards, the proposed model was implemented into ABAQUS as user subroutine VUMAT with semiimplicit back-Euler integration algorithm. The numerical simulations of hydro-bulging and cylindrical deep drawing tests were then carried out to verify the validity of the proposed model. The results showed that the proposed model could enhance the predicting accuracy in both tests due to its flexibility. & 2016 Elsevier Ltd. All rights reserved.

Keywords: High strength steel Evolving plasticity FEM Biaxial tensile test Earing prediction

1. Introduction The application of the advanced high strength steel alloys (AHSS) has been increasing in the automobile industry for the weight reduction and collision-safe performance. However, AHSS sheets are difficult to form in sheet metal forming processes, as they have relatively complex microstructures which lead to a typical plasticity. Compared with tradition trial-and-error method, Finite Element Method (FEM) is a promising way in designing the process of AHSS, since it can lower the cost remarkably with better precision in forming prediction. Nevertheless, the confidence in the FEM formability analysis strongly relies on the accuracy of the constitutive model which describes the plastic behavior of the materials [1]. This is crucial to the material exhibiting anisotropic behavior, as most cold-rolled sheet metals do. Anisotropic behavior of the sheet metals has been studied in great depth and consequently a number of anisotropic yield criteria have been proposed. Generally, the anisotropy is introduced by adding a number of anisotropic coefficients to the original isotropic yield function. For instance, Hill’s quadratic yield function [2], widely used for its simplicity, is derived from von Mises yield n

Corresponding author. Tel.: þ 86 10 82338613; fax: þ 86 10 82338788. E-mail address: [email protected] (X. Wu).

http://dx.doi.org/10.1016/j.ijmecsci.2016.01.006 0020-7403/& 2016 Elsevier Ltd. All rights reserved.

function by introducing six orthotropic parameters. It has been proven that Hill’s quadratic yield function is capable in predicting the formability of some ferrous metal alloys [3]. To precisely describe the deformation behavior of the aluminum alloys, Barlat et al. [4–7] proposed a series of anisotropic yield functions based on Hosford isotropic criterion [8]. Similar approach was then performed by Banabic et al. [9–12]. These Hosford derived criteria improve their precision by adding more coefficients to the functions, which would inevitably complicate the yield function and thus lower the computational efficiency in the simulations [12]. For the propose of introducing more coefficients without complicating the form of the yield function, Karafillis and Boyce [13] proposed the idea of isotropic plastic equivalent (IPE), which introducing the anisotropy by linearly transforming the stress tensor. Inspired by Karafillis and Boyce, Barlat et al. [14] then performed two independent linear transformations on the deviatoric stress tensor and formulated the Yld2000-2d yield function. Afterwards, more flexible yield function were proposed by Barlat et al. [15], Cazacu et al. [16], Aretz et al. [17], etc. using the similar approach. It should be noted that these yield criteria are all phenomenological and need to be verified and tested before applying on the AHSS sheets. Recently, some studies have been carried out on the applicability of present yield criteria for AHSS sheets. Mohr et al. [18] proposed that the Hill quadratic yield criterion combined with

44

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

non-associated flow rules could provide good estimates for dual phase and TRIP-assisted steel under multi-axial loading. Padmanabhan et al. [19] used this yield criterion with Swift hardening law to predict the square cup deep-drawing of the dual phase steel tailor-welded blanks. Nevertheless, Andar et al. [20] performed a series of biaxial tensile tests on DP590 and suggested that the Yld2000-2d yield function could describe the yield behavior of the dual phase steel sheets with better precision. Detailed experimental investigations were then carried out by Kuwabara et al. [21] to reveal the yield behavior of DP980 under tension–compression and pure shear states. It indicated that the Yld2000-2d yield function with exponent of 4 could have better agreement with the experimental data. The identical yield function, combined with modified Chaboche hardening model, was adopted by Ahn et al. [22] in simulating the springback of TWIP-assisted steel sheets with enhanced accuracy. It could be concluded that both Hill quadratic and Yld2000-2d yield criteria are commonly applied in the AHSS sheet, while the Yld2000-2d might exhibit some advantages in predicting with higher precision. Conventionally, the anisotropic coefficients adopted in the yield functions were associated with material properties at the initial yield point. In this paper, however, it would indicate that the plastic behavior of some AHSS sheets evolves obviously during the deformation, thus the constitutive model could lead to an inaccurate result if just using the material properties at one single stage. Recently, the phenomenon of evolving anisotropy was studied mostly for hexagonal-close-packed alloys and aluminum alloys for their complex microscopic deformation mechanism. Plunkett et al. [23,24] observed strong evolving anisotropy during the deformation of the high-purity zirconium and then proposed a modified CPB06ex3 yield criterion with variable anisotropic parameters. The parameters corresponding to an arbitrary equivalent plastic strain are determined by linearly interpolating between the pre-computed parameters associated with several fixed levels of equivalent plastic strain. Based on Plunkett, Ghaffari Tari et al. [25] suggested another CPB06ex3 derived model to calibrate the evolving anisotropy and asymmetry of the AZ31B-O magnesium alloys. In Ghaffari’s model, the anisotropy and asymmetry parameters of the CPB06ex3 yield function are replaced with Voce-type functions expressed in terms of plastic strain. Yoon et al. [26] found the similar phenomenon in the AA5042-H2 aluminum alloy and proposed the modified Yld2000-2d and CPB06ex2 yield functions with dynamically updated anisotropy coefficients. Since the anisotropy coefficients of Yoon’s yield functions should be computed and updated in each incremental step instantly by solving nonlinear equations, the computational efficiency of these two modified yield functions is unsatisfying. Wang et al. [27] also observed the similar evolving anisotropy in the AA5754-O aluminum alloy, and proposed a modified Yld20002d yield criterion which considers the evolving of anisotropy. In Wang’s model, the 8 anisotropy coefficients of the Yld2000-2d are expressed as sixth-order polynomial functions of the equivalent plastic strain. With the modified yield model, Wang successfully enhanced the accuracy of predicting earing profile in deep drawing test. For the above studies, they all focused on the changes of the anisotropic coefficients in the yield functions to calibrate the anisotropic part of the evolving plasticity. The change of the exponent in the yield function, however, has not been taken into consideration, which leads to incapability in describing the isotropic part of the evolving plasticity. The earing behavior of cold-rolled sheet metal was noticed in the early applications of the cylindrical deep drawing process. Detailed analyses have been conducted on the dependency between the earing amplitude and the planar anisotropy. Yoon et al. [28] proposed an analytical approach with arbitrary anisotropic yield function to determine the earing profile in the drawing

and redrawing processes. It indicated that the distribution of earing has strong dependency on the r-value and yield stress directionalities. Further works conducted by Chung et al. [29] proved that the earing profile could be determined analytically only with the anisotropic properties of yield stress and the r-value. With the development of the numerical simulating techniques, the earing prediction could be obtained by finite element simulations. The analyses of the FEM predictions revealed that earing profile has significant dependency on the chosen constitutive model, especially the yield criteria. Therefore, it has been widely accepted that the comparison of FEM predicted earing profile with the experimental one is a convincing way to evaluate the capability of proposed anisotropic yield criteria [12,27,30]. In this paper, the experimental plastic work contours of DP590, DP780 and DP980 sheets were obtained with biaxial tensile tests. Inspired by the evolving trend of the subsequent plastic work contours, the Yld2000-2d yield criterion with variable exponent and coefficients (Yld2000-2d-Var) was proposed to calibrate the anisotropic and isotropic parts of the evolving plasticity simultaneously. The proposed model was then implemented into ABAQUS as user subroutine VUMAT. Afterwards, the finite element analyses of hydro-bulging tests were carried out with different constitutive models (Mises, Hill48, Yld2000-2d and Yld2000-2d-Var). Besides, the deep drawing tests were also carried out and were simulated with different yield criteria. The validity of the proposed model was verified by comparing the experimental and FEM predicted results.

2. Materials and experiments 2.1. Materials The materials used in this paper were the dual phase steel sheets with tensile strength of 590 MPa, 780 MPa and 980 MPa, which produced by BaoSteel, Inc. The thickness of DP590 sheets was 1.5 mm, while DP780 and DP980 were 1.2 mm. The material properties obtained from uniaxial tensile test are listed in Table 1. Fig. 1 shows the uniaxial nominal stress–nominal strain curves of the DP590, DP780 and DP980 from different tensile directions. These tests were performed on the Instron tensile test machine. From Fig. 1 it could be observed that the DP780 and DP980 sheets exhibited more anisotropy than that of DP590. 2.2. Testing apparatus and specimen The biaxial tensile testing system used in this study includes a loading test machine, a control unit and related software. It was improved based on the prototype originally designed by Wu et al. [31]. The loading test machine, shown in Fig. 2, has two pairs of Table 1 Mechanical properties of the materials. Materials

Direction from RD (deg)a

σ0.2 (MPa)

σb (MPa)

r-valueb

DP590

0 45 90 0 45 90 0 45 90

398.4 404.5 387.6 552.5 545.7 578.2 706.7 685.6 724.9

631.7 645.3 630.0 851 833 887 1023 991 1048

0.963 0.815 1.096 0.582 0.943 0.827 0.609 0.914 0.749

DP780

DP980

a b

RD is the abbreviation of Rolling Direction. Measured at the uniaxial nominal strain 0o εN o 0:06

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

45

Fig. 1. Nominal stress–nominal strain curves of the dual phase steel sheets.

Fig. 3. Cruciform specimen used in biaxial tensile tests (dimensions in mm) [27]

2.4. Measurement of plastic work contours

Fig. 2. The biaxial tensile testing apparatus.

hydraulic cylinders, and the hydraulic pressure of these cylinders is servo-controlled independently. With the closed-loop control circuits, this machine could perform biaxial loading test under given loading path. Various kinds of cruciform specimens were designed for the biaxial tensile test [32]. For the purpose of lowering the shear stress in the gage area, cruciform specimen with slots on each arm was adopted, as shown in Fig. 3. The geometry of the specimen was proposed by Kuwabara et al. [33] and Wu et al. [31]. The rolling and transverse directions of the specimen were set to be parallel to the x and y loading axes of the biaxial tensile testing apparatus respectively, while the center of the specimen was set to be the origin of the coordinates. The specimens were fabricated by laser cutting and 11 slots (82 mm long and 0.2 mm wide) were made on each arm, so as to lower the shear stress in the gage area.

2.3. Test procedures The biaxial tensile tests were conducted under linear loading path for the purpose of obtaining the yield loci in the normal stress space. For the specimen of DP590, DP780 and DP980, 9 fixed ratios were selected as the biaxial loading proportions, e.g. Fx:Fy ¼4:0, 4:1, 4:2, 4:3, 4:4, 3:4, 2:4, 1:4, 0:4. For the load ratios of 4:0 and 0:4, the standard uniaxial tensile tests were adopted. The normal strain components, εx and εy, were measured with a pair of extensometers. These extensometers were placed on both side of the 96 mm  96 mm gage area and were parallel to the rolling and transverse directions of the specimen separately. The normal stress components, σx and σy, were defined as dividing the measured tensile loads by the current cross-sectional area of the gage region, which could be calculated with the assumption of constant volume [33].

In order to investigate the initial and subsequent yield behavior of the materials, the concept of plastic work contour was adopted [21,33,34]. Plastic work contours are commonly considered as yield loci for simplicity. According to this concept, two different stress states are considered to be on the same yield locus only if their plastic work per unit volume is equal. Therefore, the relationship between the biaxial and uniaxial stress state with the same equivalent plastic strain could be written as Z Z Z σ 1 dεp1 þ σ 2 dεp2 ¼ σ dεp ð1Þ where the σ 1 , εp1 and σ 2 , εp2 are the stress and plastic strain along the rolling and transverse directions respectively. The σ and εp are the equivalent stress and equivalent plastic strain, normally referring to the true stress and logarithmic strain obtained by uniaxial tensile test from rolling direction. A detailed description of calculating the plastic work contour was given by Kuwabara et al. [33].

3. Anisotropic yield functions In this paper, three representative yield functions, Mises, Hill’s quadratic (Hill’48) and Yld2000-2d, were adopted. Hill’s quadratic yield function under plane stress state could be expressed as ðG þ HÞσ 2x þ ðF þ HÞσ 2y  2H σ x σ y þ2N σ 2xy ¼ σ

ð2Þ

where F, G, H and N are anisotropic coefficients of the material. Hill’s quadratic yield function reduces to Mises yield function if F¼G ¼H ¼0.5 and N ¼1.5. The Yld2000-2d yield function was proposed based on the idea of isotropic plastic equivalent (IPE) [13]. The anisotropy was introduced to this yield function by linear transforming the stress tensor, thus the convexity of this yield function could be proven theoretically [14]. The Yld2000-2d yield function has as many as eight anisotropy coefficients. Therefore, it could predict the yield behavior of the rolled sheets with suitable precision. The Yld2000-2d yield function could be expressed as

ϕ ¼ ϕ0 þ ϕ″ ¼ 2σ m

ð3Þ

46

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

4. Analysis on the yield behavior of the dual phase steel sheets

where m is a material coefficient and (

ϕ0 ¼ j Y 01  Y 02 j m ϕ″ ¼ j 2Y ″2 þ Y ″1 j m þ j 2Y ″1 þ Y ″2 j m

ð4Þ

Here, ϕ is the combination of two isotropic functions, which are symmetric with respect to Y 0i andY ″j . The Y 0i and Y ″i are the principle values of the matrices X 0 and X″:  8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0 0 0 0 0 2 02 > > Y i ¼ 12 X 11 þX 22 7 ðX 11 X 22 Þ þ 4X 12 <  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > Y ″i ¼ 12 X″11 þ X″22 7 ðX ″11  X ″22 Þ2 þ 4X ″2 : 12

ð5Þ

X 0 ¼ L0 σ X″ ¼ L″σ

ð6Þ

where 2

L011

3

2

2=3

0

0

3

2 3 6 L0 7 6  1=3 0 0 7 6 12 7 6 7 α1 6 0 7 6 76 7 6 L21 7 ¼ 6 0  1=3 0 7 6 7 6 74 α 2 5 6 L0 7 6 7 α 0 2=3 0 5 4 22 5 4 7 L066 0 0 1 2

L″11

3

2

2

2

8

6 L″ 7 6 1 4 4 6 12 7 6 6 7 6 6 L″21 7 ¼ 16 4 4 4 6 7 96 6 L″ 7 6 2 4 22 5 4 2 8 L″66 0 0 0

2 4 1 2 0

ð7Þ

32

3

α3 6 7 07 76 α 4 7 7 76 6 7 0 76 α 5 7 7 7 76 0 54 α 6 5 α8 9 0

The experimental plastic work contours of DP590, DP780 and DP980 for different equivalent plastic strains, shown in Fig. 4(a), (c) and (e), were determined according to Eq. (1). Also depicted in the figures are the theoretical yield loci based on the von Mises, Hill’s quadratic (Hill’48) and Yld2000-2d yield criteria with different exponents. For the propose of evaluating the theoretical yield criteria, a deviation function is introduced for each yield locus [27]:

δ¼

The X 0ij and X ″ij in Eq. (5) are linearly transformed stress components and are defined as follows: (

4.1. Plastic work contours

ð8Þ

In Eqs. (6)–(8), σ is Cauchy stress tensor and α1  α8 are 8 anisotropy coefficients. These anisotropy coefficients are determined by 8 material properties of the sheet metal, such as σ 0 , σ 45 , σ 90 , r0 , r45 , r 90 which are the uniaxial tension yield stresses and r value along the rolling, 45° and the transverse directions, as well as the σ B and r B under the balanced biaxial tension. The procedure of obtaining α1  α8 was implemented based on the numerical method proposed by Barlat et al. [14]. Table 2 presents the computed anisotropy coefficients of Yld2000-2d for DP590, DP780 and DP980, respectively. The values of exponent m in Table 2 were determined according to the previous researches of Barlat et al. [14] and Kuwabara et al. [20, 21].

n 1X di qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  100% n i ¼ 1 ðσ i Þ2 þ ðσ i Þ2 1

ð9Þ

2

where σ i1 and σ i2 are the stress components of the experimental yield points measured from the rolling and transverse directions, respectively; di stands for the distance between the experimental yield point and corresponding theoretical one at the same stress ratio; n stands for the number of the experimental points at each yield locus. The deviation values between the theoretical and the experimental plastic work counters of DP590, DP780 and DP980 were calculated and plotted in Fig. 4(b), (d) and (f), separately. Fig. 4(a) and (b) indicate that the Yld2000-2d yield criterion with exponent of 4 has the closest agreement with the experimental work contours of DP590 at all three equivalent plastic strains. Hill’s quadratic yield criterion overestimates the work contours for the stress ratios σ 1 : σ 2 ¼0:4, 1:4 and 2:4, while the Yld2000-2d yield criterion with exponent of 6 underestimates the work contours for the same stress ratios. Fig. 4(c) and (d) indicate that the Yld2000-2d yield criterion with exponent of 6 has the minimal average error for DP780. Both Mises and Hill’s quadratic yield criteria overestimate the work contours markedly. It could also be observed from Fig. 4(c) that the experimental work contours of DP780 do not maintain the same shape. In fact, the outlines of the experimental work contours evolve from roughly ellipse to more Tresca-like yield locus with the increase of plastic dissipation. Fig. 4(e) and (f) demonstrate that the Yld2000-2d yield criterion with exponent of 8 has the minimal average error for DP980. Mises and Hill’s quadratic yield criteria, however, are not capable in describing the yield behavior of DP980 due to significant deviation. It could also be concluded from Fig. 4(f) that the evolution trend of the DP980’s experimental work contours is similar to that of DP780, though less obviously. 4.2. Evolution of plasticity During the plastic deformation of the sheet metal, the texture of the material would evolve due to the rearrangement of

Table 2 Anisotropy coefficients of Yld2000-2d for DP590, DP780 and DP980. Materials

m

α1

α2

α3

α4

α5

α6

α7

α8

DP590

4 6 8 4 6 8 4 6 8

0.8472 0.9179 0.9374 0.8740 0.9644 0.9918 0.8878 0.9586 0.9807

1.1678 1.1005 1.0835 0.9988 0.9414 0.9284 0.9850 0.9555 0.9510

0.9341 0.9332 0.9338 1.0651 1.0777 1.0835 1.0400 1.0570 1.0644

0.9769 0.9948 0.9989 0.9817 0.9907 0.9918 1.0068 1.0080 1.0068

1.0138 0.9977 0.9930 1.0667 1.0461 1.0391 1.0582 1.0405 1.0343

0.9341 0.9332 0.9338 1.0651 1.0777 1.0835 1.0400 1.0570 1.0644

0.9454 0.9611 0.9679 1.0113 1.0116 1.0118 1.0166 1.0223 1.0247

1.0636 1.0370 1.0251 0.9881 0.9877 0.9874 1.0536 1.0428 1.0381

DP780

DP980

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

47

Fig. 4. Comparison between the theoretical and experimental plastic work counters for (a) DP590, (c) DP780 and (e) DP980 (for clarity, only two closest results of Yld20002d for each material were depicted); Deviations between the theoretical and the experimental plastic work counters of (b) DP590, (d) DP780 and (f) DP980.

crystallographic orientation, which could consequently lead to the evolving plastic behaviors of the material in the macroscale. As indicated in Fig. 4, the suitable parameters of the Yld2000-2d for the DP780 and DP980 at the initial yield stage would result in relatively inaccurate description at subsequent yield stage. Therefore, the evolving plasticity should be taken into account to achieve sufficient accuracy during the plastic deformation. For the Yld2000-2d yield function, the plastic behavior could be calibrated by two aspects: the coefficients αi (i¼1, 2, …, 8) describe the anisotropic part of plasticity, while the exponent m mainly calibrates the isotropic part of the plasticity. In the former studies, the evolving anisotropy was paid much attention especially in the HCP alloys, as the flow stresses and r-values of these materials usually change dramatically due to the strong texture evolution. To reveal the evolving anisotropy of the given DP steels, the r-values and homogenized flow stresses under different equivalent plastic strains were measured and plotted in Figs. 5 and 6 separately. Fig. 5 indicates that the r-values of the given dual phase steel sheets change considerably if the equivalent plastic strain falls within the range of 0–0.01, which probably because of the

nonlinear elastic strains at the initial yield stage. For the rest strain range, however, the r-values approximately stay the same. Fig. 6 shows that the homogenized flow stresses of the DP590 vary with clear tendency, which reveals perceptible evolving anisotropy during plastic deformation; while for DP780 and DP980, this phenomenon is less obvious. It could be concluded that these given materials have weak evolving anisotropy which could be taken into consideration in the material modeling for higher precision. As shown in Fig. 4, the relatively optimal exponent of Yld20002d varies with the increase of the plastic dissipation for DP780, which indicates that the isotropic part of the evolving plasticity should be considered to achieve sufficient precision. Nevertheless, the isotropic part of the evolving plasticity has not been sufficiently investigated in the previous studies, since the exponent of the yield function was mostly set to be fixed once it had been identified. Conventionally, the optional value of the exponent is determined by the crystal structure of the material, i.e. m ¼6 for bcc materials and m ¼8 for fcc materials, which is consistent with the original Hershey–Hosford yield function. Yet flexible value of

48

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

1.2

r90

1.2 0.9

0.0 0.00

r0

r45

0.3

r-value

r-value

0.9 0.6

r45

0.04

r0

0.3

Material: DP590 0.02

r90

0.6

0.06

1.2

Material: DP780

0.0 0.00

0.02

0.04

0.06

r45

r-value

0.9 0.6 0.3 0.0 0.00

r0

r90 Material: DP980

0.02

0.04

0.06

Fig. 5. The r-values vs. equivalent plastic strain along different directions for (a) DP590, (b) DP780 and (c) DP980. The equivalent plastic strains of 45° and 90° tensile curves were determined by the principle of plastic work contour (mentioned in Section 2.4).

Fig. 6. The homogenized flow stresses vs. equivalent plastic strain along different directions for (a) DP590, (b) DP780 and (c) DP980. The equivalent plastic strains of 45° and 90° tensile curves were determined by the principle of plastic work contour (mentioned in Section 2.4).

the exponent was suggested by Barlat et al. [14] and Kuwabara [21] to adapt variety kinds of sheet metals. In order to calibrate the isotropic part of the evolving plasticity for the dual phase steel sheets, the regularity of how the exponent changes the shape of the theoretical yield surface should be investigated. Taking the

isotropic case for simplicity, the Yld2000-2d yield function reduces to Hershey–Hosford yield function when all eight anisotropy coefficients equal to one (α1 ¼ α2 ¼…¼ α8 ¼1) [8] jσ 1 jm þ jσ 2 jm þ jσ 1  σ 2 jm ¼ 2σ m

ð10Þ

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

1.25

1.25

1.00

1.00

0.75

0.75

m=2 m=6 m=10 m=20 m→

0.50 0.25 0.00 0.00

0.50

m=2 m=6 m=10 m=100

0.25

0.75

0.50

0.25

Anisotropic

2

2

Isotropic

0.00 0.00

1.25

1.00

49

0.25

0.75

0.50

1

1.25

1.00

1

Fig. 7. Evolution of the theoretical yield loci of Yld2000-2d with the growth of the exponent m: (a) isotropic, (b) adopting the anisotropy coefficients of DP780 for m ¼6.

0.20

0.20 Material: DP590

ε p = 0.2%

ε p = 0.4%

0.10

0.05

0.00

4

5

6

ε p = 0.6%

0.10

0.05

ε p = 1.7% 3

ε p = 0.2%

0.15

δ

δ

0.15

Material: DP780

7

0.00

ε p = 1.5% 3

4

5

Exponent m

6

7

8

9

10

Exponent m

0.20 Material: DP980

ε p = 0.6%

0.15

δ

ε p = 1.0% 0.10

0.05

0.00

ε p = 0.2% 4

5

6

7 8 9 Exponent m

10

11

12

Fig. 8. Deviation of the Yld2000-2d yield function vs. the exponent m of the yield function at different deformation stages for: (a) DP590, (b) DP780, and (c) DP980. Note that the anisotropic coefficients used in each data point of these figures were numerically determined by the given exponent and the material properties at corresponding deformation stage.

where m is the exponent of the yield function. If m ¼2, it then reduces to Mises yield function

σ 21 þ σ 22 þðσ 1  σ 2 Þ2 ¼ 2σ 2

ð11Þ

and the theoretical yield locus turns out to be an ellipse in the σ1– σ2 plane. If m tends to infinity, however, it reduces to Tresca yield function

σ ¼ max fjσ 1 j; jσ 2 j; jσ 1  σ 2 jg

ð12Þ

and the theoretical yield locus becomes two perpendicular line segments in the first quadrant of the σ1–σ2 plane. The yield loci of the rest cases are placed between the Mises and the Tresca yield loci separately, as shown in Fig. 7(a). For the anisotropic case, taking the anisotropy coefficients of DP780 as example, it also shares the same characteristic, i.e. a higher exponent of the

Yld2000-2d yield function could lead to a more Tresca-like theoretical yield locus in the first quadrant, as presented in Fig. 7(b). Considered the role of the exponent in Yld2000-2d mentioned above and the evolution trend of the experimental work contours noted in Section 4.1, it may be suitable to assume that the optimum exponent for DP780 and DP980 sheets should change during the deformation to calibrate the isotropic part of the evolving plasticity. To verify this assumption, the deviations of the Yld2000-2d yield function with a series of exponents were calculated for each experimental work contour, as presented in Fig. 8. It confirms that for DP780 and DP980, the optimal exponent of the Yld2000-2d yield function do increase with the growth of the equivalent plastic strain; while for DP590, this trend is less significant.

50

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

9

5.0 Material: DP590

8

4.8

Material: DP780

mopt

7 mopt

4.6 4.4 4.2 0.0

6 5

0.5

1.0

4 0.0

1.5

0.5

1.0

1.5

ε p (%)

ε (%) p

9.5 9.0

Material: DP980

mopt

8.5 8.0 7.5 7.0 0.2

0.4

0.6

0.8

ε (%) p

Fig. 9. Dependence of the optimal exponent of Yld2000-2d yield function and the equivalent plastic strain for (a) DP590, (b) DP780 and (c) DP980.

4.3. Modeling and identification of the plasticity evolution In this paper, the exponent m and coefficients αi (i¼1, 2, …, 8) of Yld2000-2d were set to be variable to calibrate the isotropic and anisotropic parts of evolving plasticity separately. Since the coefficients αi are numerically determined by the exponent m and the eight material properties [14], the evolving trends of the coefficients αi implicitly rely on the evolving trend of the exponent m. Therefore, it is practicable to model the trend of evolving exponent m in the first place. For the purpose of obtaining the optimal exponent mopt for the given experimental plastic work contour, the following optimization procedure was proposed: min

α A ℝ ;m A ℝ 8

s:t:

N X

g 1 ðα; mÞ ¼ (

‖P itheory ðα; m; εp Þ  P iexp ðεp Þ‖2

i¼1 p

Fðα; m; σ θ ðε ÞÞ ¼ 0

σθ

Δm

p1

p2

R-Square

DP590 DP780 DP980

4.339 6.036 7.407

0.4056 1.659 1.619

285.3 423.2 597.9

2.354E-14 4.100 1.136

0.7837 0.8615 0.7747

Gðα; m; r θ Þ ¼ qx ðr θ Þ

∂ϕðα; mÞ ∂ϕðα; mÞ  qy ðr θ Þ ¼0 ∂sxx ∂syy

α A ℝ ;m A ℝ

ð14Þ

ð15Þ

where the detailed definition of ϕ, qx , qy ,σ θ and r θ are adequately elaborated in the appendix section of Ref. [14]. With classic Courant Penalty function method, the constrained minimization problem (13) could be converted into unconstrained one expressed as follows:

8

where g 1 ðα; mÞis the objective function of the optimization procedure; P iexp stands for the experimental yield point; P itheory stands for the corresponding theoretical yield point predicted by Yld2000-2d yield function; α stands for the set of the anisotropy coefficients; εp stands for the equivalent plastic strain of the given yield locus;σ θ ðεp Þ and r θ ðεp Þ are the material properties (yield stress and r-value along specific direction) at given equivalent plastic strain.N stands for the total number of the experimental yield points at current plastic work contour. The nonlinear equations Fðα; m; σ θ Þ ¼ 0 and Gðα; m; r θ Þ ¼ 0 are the constraints of the anisotropy coefficients, which were firstly proposed by Barlat et al. [14] aiming to determine the anisotropy coefficients using eight material properties: ¼0

m0

min

p

2σ

Materials

ð13Þ

Gðα; m; r θ ðε ÞÞ ¼ 0

Fðα; m; σ θ Þ ¼ ϕðα; mÞ 

Table 3 Results of the nonlinear fitting on the optimal exponents mopt ðεp Þ

g 2 ðα; mÞ ¼

N  X

‖P itheory ðα; m; εp Þ  P iexp ðεp Þ‖2 =‖P iexp ðεp Þ‖2



i¼1

  þ p U ‖Fðα; m; σ θ Þ‖2 þ‖Gðα; m; r θ Þ‖2

ð16Þ

where g 2 ðα; mÞ is the objective function of the optimization procedure; p is the penalty parameter and could be set to one in this case. The minimization (16) has been performed using Levenberg–Marquardt algorithm implemented in MATLAB and the results are then plotted in Fig. 9. Note that the constant r-values (presented in Table 1) was adopted in this approach [26], considered the nonlinear elastic behavior at the early yield stage affecting the precision of the measured r-value. For the purpose of modeling this trend, a bonded function was assumed as the evolution function of the exponent, which could be expressed as m ¼ m0 þ Δm U½expðp1 U εp  p2 Þ 1=½expðp1 U εp  p2 Þ þ 1

ð17Þ

where m0 stands for the initial exponent, Δm stands for the

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

51

Fig. 10. Dependence of the anisotropy coefficients of the Yld2000-2d yield function and the corresponding equivalent plastic strain for (a) DP590, (b) DP780 and (c) DP980.

theoretical fluctuation range of the exponent, p1 and p2 are the fitting coefficients. This evolution function was chosen because it is well-behaved (continuous, differentiable, monotonic and bounded) in the range of plastic strain considered here. The fitting results of obtained ðεp ; mopt Þ are then presented in Table 3 and Fig. 9. It is proven again by Fig. 9 that the optimal exponent of the Yld2000-2d yield function for DP590, DP780 and DP980 increase evidently with the growth of equivalent plastic strain. It can also be confirmed from Fig. 9 that the proposed nonlinear function (17) can fit this trend with convincing precision. With the obtained fitting results mopt ðεp Þ and the material properties (σ θ ðεp Þ and r θ ðεp Þ), the anisotropic coefficients of Yld2000-2d for DP590, DP780 and DP980 at different plastic stages could be determined by Eqs. (14) and (15), which were then plotted in Fig. 10. For the purpose of calibrating the evolving trend of the anisotropic coefficients, a nonlinear fitting function was introduced

αi ¼ ai U expðbi U ε Þ þci Uexpðdi U ε Þ ði ¼ 1; 2; :::; 8Þ p

p

ð18Þ

where αi stands for the anisotropy coefficient of the Yld2000-2d and ai, bi, ci, di are the fitting parameters. This kind of combined exponential function was chosen because it exhibits distinct capability in characterizing the variation trend of α, as shown in Fig. 10. The lengthy results of the fitting parameters are listed in the Appendix section of this paper. In order to prove the benefit of the proposed variable exponent method, the optimum capability of the Yld2000-2d yield function with fixed exponent and coefficients (original Yld2000-2d) should be taken into consideration. Therefore, the deviations of the original Yld2000-2d yield function with commonly adopted exponents (4, 6 and 8) and the Yld2000-2d yield function with variable coefficients and exponent (Yld2000-2d-Var) were computed and depicted in Fig. 11. It could be observed from Fig. 11 that for these three kinds of materials, the Yld2000-2d-Var yield function has a

significantly lower deviation than that of the original Yld2000-2d with most suitable exponent.

5. Implementation of the proposed model In this section, the Yld2000-2d-Var yield criterion mentioned above would be implemented in a general purpose finite element code via VUMAT subroutine of ABAQUS. The semi-implicit backward-Euler algorithm, firstly proposed by Moran et al. [35], was adopted as the integration method for the proposed constitutive model. This algorithm can simplify the computation significantly and would be suitable for explicit dynamics simulations. For each strain increment Δε provided by the previously computed displacement field, the stress tensor σ is updated using the following semi-implicit scheme:

Δεpn þ 1 ¼ Δεpn þ 1 r n

ð19Þ

Δσ n þ 1 ¼ C : ðΔεn þ 1  Δεpn þ 1 Þ Δmn þ 1 ¼ Δεpn þ 1



∂m ∂εp

ð20Þ

 ð21Þ n

f n þ 1 ¼ σ ðσ n þ 1 ; mn þ 1 ; αn þ 1 Þ  σ Y ðεpn þ 1 Þ ¼ 0

ð22Þ

where εp is the plastic strain tensor; εp is the equivalent plastic strain; σ is the equivalent stress associated with the yield criterion; α is a set of the anisotropy coefficients for the Yld2000-2d yield function; C is the fourth-order elasticity tensor; σ Y is the hardening curve of the material; r is the direction of plastic flow and could be expressed as r¼

∂f ∂σ ¼ ∂σ ∂σ

ð23Þ

with the assumption of associated flow rule. The variables in Eqs. (19)–(22) with subscript n are of the previous incremental step,

52

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

Material: DP590

2.5

1.5

2.0

m=8 m=6

δ(%)

δ(%)

2.0

2.5 Material: DP780

1.0 m=4 0.5 0.0 0.2

1.5

m=6

1.0 0.5

0.7

Var-Exp-Coef 1.2

δ(%)

3 2

Var-Exp-Coef

0.0 0.2

1.7

ε p (%)

m=4

m=8

0.7

ε p (%)

1.2

1.7

Material: DP980

m=4

m=6 m=8

1 0 0.2

Var-Exp-Coef

0.7

1.2

1.7

Fig. 11. Comparison between the original Yld2000-2d yield function and the Yld2000-2d yield function with variable exponent and coefficients (Yld2000-2d-Var) on the deviation between theoretical and experimental plastic work counters for (a) DP590, (b) DP780 and (c) DP980. Note that the ‘Var-Exp-Coef’in the figures is the abbreviation of ‘Variable Exponent and Coefficients’.

while variables with subscript n þ1 are of the current step and need to be solved. The return mapping method, consisting of an elastic predictor stage and a corresponding plastic corrector stage, was adopted to solve the integration problem presented in Eqs. (19)–(22). In the predictor stage, the current strain increment is assumed to be pure elastic. Therefore, the initial value of the stress tensor, equivalent plastic strain and exponent for the current step could be expressed as

σ trial n þ 1 ¼ σ n þ C : Δεn þ 1 

ε

p nþ1

mtrial nþ1

trial

¼ε

p n

¼ mn

αtrial n þ 1 ¼ αn

ð24Þ ð25Þ

and the corresponding incremental forms of Eqs. (19), (21) and (29) are presented as follows:

ΔεpðkÞ ¼ ΔεpðkÞ r n ΔmðkÞ ¼ ΔεpðkÞ



ð26Þ ð27Þ

where superscript ‘trial ‘denotes the trial state of the current step. If the following criterion is satisfied within a proper numerical tolerance (Tol) at the trial state; i.e.   p trail trail oTol ð28Þ σ ðσ trail n þ 1 ; mn þ 1 Þ  σ Y ðεn þ 1 Þ the current step would be confirmed as pure elastic. Otherwise, this step is considered to be elasto-plastic and corrector stage is needed to pull the stress tensor back to the new yield surface. The trial elastic stress σ trial n þ 1 computed in the predictor stage will be adopted as the initial value for the corrector stage presented as follows [36]: p p trial σ n þ 1 ¼ σ trial n þ 1  C : Δεn þ 1 ¼ σ n þ 1  Δεn þ 1 C : r n

equations cannot be solved directly since the consistency condition Eq.(22) is a nonlinear implicit function for σ , m, and εp . Thus the Newton–Raphson algorithm is introduced to solve the equations numerically. The iteration form of Eq. (22) could be written as  ðkÞ  ðkÞ  ðkÞ ∂f ∂f ∂f ðkÞ pðkÞ f þ : Δσ ðkÞ þ Δε þ ΔmðkÞ ¼ 0 ð30Þ ∂σ ∂m ∂ εp

ð29Þ

where σ n þ 1 and Δεpn þ 1 need to meet the consistency condition Eq. (22) to ensure the new stress tensor stays on the new yield surface. The corrected stress tensor σ n þ 1 could be obtained by solving Eqs. (19), (21), (22) and (29) concurrently. Nonetheless, this set of

∂m ∂ εp

ð31Þ  ð32Þ n

Δσ ðkÞ ¼  ΔεpðkÞ C : rn

ð33Þ

where the variables with superscript (k) are of the kth iteration for the current incremental step. Then, the increment of equivalent plastic strain ΔεpðkÞ could be solved by substituting Eqs. (31) and (32) into Eq. (30):

ΔεpðkÞ ¼   ∂σ ðkÞ ∂σ

f

ðkÞ

 ðkÞ  ∂σ ðkÞ  ∂m  : C : r n  ∂m þ ∂∂σεpY p ∂ε

ð34Þ

n

where the partial derivative ∂σ =∂m is presented as follows: m m P1 P2 P3 ∂σ σ Pm 1 ln σ þ P 2 ln σ þ P 3 ln σ ¼ m m m ∂m m P1 þ P2 þ P3

ð35Þ

P 1 ¼ Y 01  Y 02

ð36Þ

P 2 ¼ 2Y ″2 þ Y ″1

ð37Þ

P 3 ¼ 2Y ″1 þ Y ″2

ð38Þ

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

53

Fig. 12. Results of experimental hydro-bulgingtests: (a) DP590, (b) DP780, and (c) DP980.

Here, the Y 0i and Y ″i are defined in Eq. (5) and the detailed deduction of ∂σ =∂σ was adequately elaborated by Barlat et al. [14]. Afterwards, The increment of σ , m, and εp could be easily obtained by substituting Eq. (34) in Eqs. (31)–(33) and the anisotropy coefficients could be updated by Eq. (18). The iteration procedure then continues with the newly updated stress tensor and other state variables until the consistency condition is satisfied within a proper numerical tolerance (Tol): f

ðkÞ

¼ σ ðσ ðkÞ ; mðkÞ Þ  σ Y ðεpðkÞ Þ o Tol

ð39Þ

As shown in Eqs. (31)–(34), the iteration scheme of the semiimplicit algorithm would not update the plastic flow r during the local iteration and consequently the partial derivative of the plastic flow ∂r=∂σ is no longer needed. This feature would simplify the computation and shorten the simulation time significantly compared with the classic fully-implicit algorithm [36]. It should be noted that the anisotropy coefficients would be fixed as constants if the current equivalent plastic strain is out of the fitting range shown in Fig. 10, as Eq. (18) is not a bounded fitting function which could consequently lead to failing or meaningless extended fitting results.

6. Results 6.1. Hydro-bulging tests and simulations In this subsection, the hydro-bulging tests and corresponding numerical simulations were proceeded to verify the constitutive models of the dual phase steel sheets. The experimental strain at the top of the bulging area could be obtained by measuring the etched gird on the specimens digitally. The balanced bulging stress, however, could be calculated as follows:

σ ¼ p U ρ=2t

ð40Þ

where p stands for the hydro-forming pressure; ρ stands for the curvature radius at the top of the bulging area and could be

measured by the curvature meter; t stands for the current thickness at the top position. The experimental results of the hydro-bulging tests are shown in Fig. 12 and the corresponding stress–strain points are presented in Fig. 13. Four yield criteria, i.e. Mises, Hill’s quadratic, Yld2000-2d and Yld2000-2d with variable exponent and coefficients (Yld2000-2dVar), were adopted as the constitutive models in the simulations. Mises and Hill’s quadratic yield functions are already included in ABAQUS; while the Yld2000-2d and Yld2000-2d-Var were implemented in a finite element code via VUMAT subroutine of ABAQUS. The values of the coefficients used in the simulations are listed in Table 2 (for original Yld2000-2d yield function with most suitable exponent), Table 3 and Table A.1 (for Yld2000-2d-Var yield function). The numerical predicted stress–strain curves of hydro-bulging tests are then presented in Fig. 13. As shown in Fig. 13, it could be observed that the Yld2000-2dVar yield criterion shows distinct advantage in describing the deformation behaviors of the dual phase steel sheets under hydrobulging condition, partially because the Yld2000-2d-Var model implicitly contains the subsequent material properties under balanced tension condition. The original Yld2000-2d (with most suitable exponent) and von Mises yield criteria also have relatively good accuracy, though less precise than the Yld2000-2d-Var. The predictions made by Hill’s quadratic yield criterion, however, is relatively rough and imprecise. 6.2. Cylindrical deep drawing tests and simulations In this subsection, the cylindrical deep drawing tests and corresponding numerical simulations were carried out to verify the constitutive models for the dual phase steel sheets.The geometry dimensions of the tools for the deep drawing tests are specified as follows: The punch diameter is 50.00 mm and the punch profile radius is 5.00 mm. As the thicknesses of these sheets are not identical, the dies with different dimensions were adopted. For DP590, the opening diameter of the die is 55.20 mm and the die profile radius is 18.6 mm; while for DP780 and DP980, the

54

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

Material: DP590

900 Experimental Mises Hill48 Yld2000-2d(m=4) Yld2000-2d-Var

800

700 0.0

0.1

Material: DP590

1000

Rolling Direction

σ2(MPa)

σ1(MPa)

1000

0.2

0.3

0.4

Transverse Direction

900 Experimental Mises Hill48 Yld2000-2d(m=4) Yld2000-2d-Var

800

700 0.0

0.5

0.1

0.2

0.3

ε1 1300

Material: DP780 Transverse Direction

Rolling Direction

1200

1100

Experimental Mises Hill48 Yld2000-2d(m=6) Yld2000-2d-Var

1000 0.0

0.1

0.2

0.3

0.4

σ2(MPa)

1200 σ1(MPa)

0.5

1300

Material: DP780

1100

Experimental Mises Hill48 Yld2000-2d(m=6) Yld2000-2d-Var

1000 0.0

0.5

0.1

0.2 ε2

ε1 1500

0.3

0.4

1500 Material: DP980

Material: DP980

Rolling Direction

Transverse Direction

1400

1300

Experimental Mises Hill48 Yld2000-2d(m=8) Yld2000-2d-Var

1200 0.0

0.1

0.2

0.3

σ2(MPa)

1400 σ1(MPa)

0.4

ε2

1300

Experimental Mises Hill48 Yld2000-2d(m=8) Yld2000-2d-Var

1200 0.0

ε1

0.1

ε2

0.2

0.3

Fig. 13. Comparison between the experimental and simulated true stress–logarithmic strain curves under hydro-bulging tests: (a), (b) DP590; (c), (d) DP780; (e), (f) DP980. Note that the ‘Yld2000-2d-Var’ is the abbreviation of ‘Yld2000-2d with Variable exponent and coefficients’.

corresponding values are 53.64 mm and 13.00 mm respectively. The initial thicknesses of the blanks used in the tests are 1.5 mm for DP590 and 1.2 mm for DP780 and DP980.All these blanks share the same diameter of 90.00 mm. The blank holding forces (BHF) used in the deep drawing tests were 15 kN equally, and the Coulomb coefficient of friction was set to be 0.2 in the simulations. The 4-node shell element (ABAQUS S4R) was adopted in the finite element modeling. Same as before, four yield criteria, i.e. Mises, Hill’s quadratic, Yld2000-2d and Yld2000-2d-Var, were adopted as the material models in the simulations. The values of the coefficients used in the simulations are listed in Table 2 (for original Yld2000-2d yield function with most suitable exponent), Table 3 and Table A.1 (for Yld2000-2d-Var yield function). The experimental results of the deep drawing tests are shown in Fig. 14, and the comparison between the experimental and numerical predicted earing profiles are presented in Fig. 15. As shown in Fig. 15, the simulated results based on the Yld2000-2d-Var yield function are in better agreement with the experimental results than other yield functions. Although the Hill’s quadratic and Yld2000-2d yield function could predict the locations of the peaks and valleys, the amplitudes of the earing profiles

(maximum and minimum height difference), however, are overestimated by Hill’s quadratic yield function and underestimated by original Yld2000-2d yield function for the most cases. What’s more, the Hill’s quadratic and Yld2000-2d yield function are not capable in predicting the heights of the peaks and valleys accurately; while the Yld2000-2d-Var yield function could provide precise predictions on the peaks’ heights, and relatively reasonable estimates on the valleys’ heights. Since no anisotropy is considered, the Mises yield function could not predict the earing phenomenon. The fluctuation predicted by Mises might be related with the numerical errors or the gap between die and punch.

7. Discussion The yield behaviors of the high strength dual phase steel sheets have been investigated with proportional biaxial tensile tests. It revealed that the given dual phase steel sheets share distinct evolving plastic behavior during the early plastic deformation, i.e. the experimental work contours in the first quadrant of σ1–σ2 plane evolve from roughly elliptic to relatively more Tresca-like yield loci with the increase of plastic dissipation, which means

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

55

Fig. 14. Results of experimental cylindrical deep drawing tests: (a) DP590, (b) DP780, and (c) DP980.

Earing Profile (mm)

27

26 Material: DP590 Experimental Yld2000-2d(m=4) Yld2000-2d-Var

25 0

45

90

Hill48 Mises 135

180

θ(Deg.) 28

27

Earing Profile (mm)

Earing Profile (mm)

28

26 Material: DP780 Experimental Yld2000-2d(m=6) Yld2000-2d-Var

25 0

45

90 θ(Deg.)

Hill48 Mises 135

180

27 26 Material: DP980 Experimental Yld2000-2d(m=8) Yld2000-2d-Var

25 24

0

45

90

Hill48 Mises 135

180

θ(Deg.)

Fig. 15. (a) FE analysis result of completely drawn cup for DP590 (showing the distribution of the equivalent plastic strain); (b)–(d) comparison between the experimental and numerical predicted earing profiles for DP590, DP780 and DP980 respectively. Note that the ‘Yld2000-2d-Var’ is the abbreviation of ‘Yld2000-2d with Variable exponent and coefficients’.

56

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

that the isotropic part of the plasticity evolution should not be ignored. Therefore, the Yld2000-2d yield function with variable exponent and coefficients (Yld2000-2d-Var) was proposed. In this model, the isotropic plasticity evolution of the given dual phase steel sheets could be well characterized by the increase of the exponent in the yield function. Due to the variable exponent of the yield function, the proposed model shows distinct flexibility in modeling the sheet metals. If the exponent was fixed to a certain integer, the proposed model would reduce to original Yld2000-2d or even Mises yield function (exponent equals to 2 in the isotropic case). Thus the sheet metals which have been previously verified to be suitable for these two yield functions could be modeled adequately by the Yld2000-2d-Var. Moreover, the variable exponent might provide the capability in enhancing the modeling precision for these sheet metals, as the given dual phase steel sheets do. It should be noted that the proposed model was developed and identified under small range of strain conditions (usually εp r 1:7%), because the cruciform specimens under biaxial tension conditions mostly fractured at the slit area of the arm before the gauge area could have reached a larger strain. Other scholars using the same specimen shape also met similar problems [3,20,27]. Even using the materials with better formability, like IF steel or aluminum alloys, the maximum strain in the gauge area could only reached around 4%. Further approaches has been made to enlarge the plastic deformation of the biaxial tension test by eliminate the slits or redesign the shape of the cruciform specimen [32], but it might gain more shear stress in the gauge area and consequently leads to imprecise experimental yield locus, which would compromise the precision in the identification procedures of the proposed model. With limited strain range in the identification procedures, the proposed model seems not to be convincing if the strain do not fall within the calibrated range. Therefore, the hydro-bulging tests and cylindrical deep drawing tests were conducted to verify the accuracy of the model under large deformation conditions. The results of the earing prediction proved that the proposed model still exhibit higher accuracy than the original Yld2000-2d at this circumstance. Nevertheless, the precision advantage of the proposed model still needs to be further tested at arbitrary large plastic deformation conditions, as this model is pure phenomenological. It should also be noted that the proposed model was developed and verified in the proportional loading condition only. The validity of the proposed model under complex loading paths, in particular reversed loading, is still needed to be verified in the further studies. Besides, the texture evolution mechanism of the given AHSS sheets has not been investigated in this work, thus the description of the evolving plasticity presented in this paper is restricted to be phenomenological.

8. Conclusion In this paper, biaxial tensile tests under the proportional loading condition were conducted to obtain the experimental plastic work contours of DP590, DP780 and DP980 sheets. It was found that for high strength dual phase steel sheets, the outlines of the experimental work contours in the first quadrant of σ1–σ2 plane do not maintain the same: it usually evolves from roughly elliptic to relatively more Tresca-like yield loci with the increase of the plastic dissipation. In order to describe the evolution trend of the subsequent work contours, the Yld2000-2d yield criterion with variable coefficients and exponent (Yld2000-2d-Var) was proposed. In this model, the exponent m and the coefficients αi of the Yld2000-2d was assumed as the function of the equivalent plastic strain. The identification procedures of the introduced parameters were then proposed based on minimizing

the difference between the experimental and the model predicted yield points, using Levenberg–Marquardt optimization method. Afterwards, the Yld2000-2d-Var was implemented into ABAQUS as user subroutine VUMAT with semi-implicit back-Euler integration algorithm. The finite element analyses of hydro-bulging tests were carried out with different constitutive models (Mises, Hill’s quadratic, Yld2000-2d and Yld2000-2d-Var). The results showed that the stress–strain curves predicted by the Yld2000-2d-Var were in better agreement with the experimental data than those of the Mises, Hill’s quadratic and original Yld2000-2d yield criteria. Besides, the deep drawing tests were also carried out and were simulated with different yield criteria. It was found again that the Yld2000-2d-Var performed better in predicting the earing profile of the dual phase steel sheets. It could be concluded that the proposed method, though with defects on the identification procedures with small strain range, might be a promising way to calibrate the isotropic and anisotropic parts of the evolving plasticity for the sheet metals with higher precision.

Acknowledgments The author would like to thank the National Natural Science Foundation of China (No. 51275026) and State Key Laboratory of Development and Application Technology of Automotive Steels (Baosteel Group) for the support given to this research. The helpful revise comments of the reviewers are gratefully appreciated.

Appendix See Table A.1. Table A.1 Nonlinear fitting parameters of the anisotropy coefficients αi. Fitting parameters

DP590

DP780

DP980

a1 b1 c1 d1 a2 b2 c2 d2 a3 b3 c3 d3 a4 b4 c4 d4 a5 b5 c5 d5 a6 b6 c6 d6 a7 b7 c7 d7 a8 b8 c8 d8

0.9727  0.3854  0.1114  178.8 0.1247  218.1 1.034 0.0909 1.035  1.456  0.106  61.25 0.9904  0.2023  0.009162  22.69 1.015  0.2011  0.01041  14.65 1.035  1.456  0.106  61.25 0.9631  0.2051  0.01699  26.30 0.07746  136.5 0.9816 0.9345

0.985  0.3423  0.08229  110.1 0.03546  140.4 0.9420 0.02253 1.084  1.508  0.019  52.39 0.9959  0.2982  0.01038  70.83 0.0226  153.2 1.039  0.2656 1.084  1.508  0.019  52.39 1.026  0.1814  0.01556  30.77  69.56 3.504 70.56 3.472

0.9966  0.6061  0.02789  107.7 7.923E  3  112.9 0.9465 7.341E  3 1.193  4.051  0.1406  20.00 1.005  0.5596 9.273E  4  67.44 1.031  0.563 7.566E  3  472.1 1.193  4.051  0.1406  20.00 1.070  0.5225  0.04783  11.75  2.350E 05 61.53 1.039 1.544

Z. Cai et al. / International Journal of Mechanical Sciences 107 (2016) 43–57

References [1] Abedrabbo N, Pourboghrat F, Carsley J. Forming of AA5182-O and AA5754-O at elevated temperatures using coupled thermo-mechanical finite element models. Int J Plast 2007;23:841–75. [2] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc R Soc Lond Ser A. Math Phys Sci 1948;193:281–97. [3] Kuwabara T, Van Bael A, Iizuka E. Measurement and analysis of yield locus and work hardening characteristics of steel sheets wtih different r-values. Acta Mater 2002;50:3717–29. [4] Barlat F, Lian K. Plastic behavior and stretchability of sheet metals. Part I: a yield function for orthotropic sheets under plane stress conditions. Int J Plast 1989;5:51–66. [5] Barlat F, Lege DJ, Brem JC. A six-component yield function for anisotropic materials. Int J Plast 1991;7:693–712. [6] Barlat F, Maeda Y, Chung K, Yanagawa M, Brem JC, Hayashida Y, Lege DJ, Matsui K, Murtha SJ, Hattori S, Becker RC, Makosey S. Yield function development for aluminum alloy sheets. J Mech Phys Solids 1997;45:1727–63. [7] Barlat F, Becker RC, Hayashida Y, Maeda Y, Yanagawa M, Chung K, Brem JC, Lege DJ, Matsui K, Murtha SJ, Hattori S. Yielding description for solution strengthened aluminum alloys. Int J Plast 1997;13:385–401. [8] Hosford W. A generalized isotropic yield criterion. J Appl Mech 1972;39:607–9. [9] Banabic D, Balan T, Comsa D. A new yield criterion for orthotropic sheet metals under plane-stress conditions. In: Proceedings of the 7th conference ‘TPR2000’. Cluj Napoca, Romania; 2000. p. 217–24. [10] Banabic D, Kuwabara T, Balan T, Comsa DS, Julean D. Non-quadratic yield criterion for orthotropic sheet metals under plane-stress conditions. Int J Mech Sci 2003;45:797–811. [11] Banabic D, Aretz H, Comsa DS, Paraianu L. An improved analytical description of orthotropy in metallic sheets. Int J Plast 2005;21:493–512. [12] D.-S.Comsa, D.Banabic, Plane-stress yield criterion for highly-anisotropic sheet metals. In: Proceedings of the 7th international conference and workshop on numerical simulation of 3D sheet metal forming processes. NUMISHEET; 2008. p. 43–8. [13] Karafillis AP, Boyce MC. A general anisotropic yield criterion using bounds and a transformation weighting tensor. J Mech Phys Solids 1993;41:1859–86. [14] Barlat F, Brem JC, Yoon JW, Chung K, Dick RE, Lege DJ, Pourboghrat F, Choi SH, Chu E. Plane stress yield function for aluminum alloy sheets—part 1: theory. Int J Plast 2003;19:1297–319. [15] Barlat F, Aretz H, Yoon JW, Karabin ME, Brem JC, Dick RE. Linear transfomation-based anisotropic yield functions. Int J Plast 2005;21:1009–39. [16] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. Int J Plast 2006;22:1171–94. [17] Aretz H, Barlat F. New convex yield functions for orthotropic metal plasticity. Int J Non-Linear Mech 2013;51:97–111. [18] Mohr D, Dunand M, Kim K-H. Evaluation of associated and non-associated quadratic plasticity models for advanced high strength steel sheets under multi-axial loading. Int J Plast 2010;26:939–56.

57

[19] Padmanabhan R, Baptista AJ, Oliveira MC, Menezes LF. Effect of anisotropy on the deep-drawing of mild steel and dual-phase steel tailor-welded blanks. J Mater Process Technol 2007;184:288–93. [20] Andar MO, Kuwabara T, Yonemura S, Uenishi A. Elastic–plastic and inelastic characteristics of high strength steel sheets under biaxial loading and unloading. ISIJ Int 2010;50:613–9. [21] Kuwabara T, Nakajima T. Material modeling of 980 MPa dual phase steel sheet based on biaxial tensile test and in-plane stress reversal test. J Solid Mech Mater Eng 2011;5:709–20. [22] Ahn K, Yoo D, Seo M, Park S-H, Chung K. Springback prediction of TWIP automotive sheets. Met Mater Int 2009;15:637–47. [23] Plunkett B, Lebensohn RA, Cazacu O, Barlat F. Anisotropic yield function of hexagonal materials taking into account texture development and anisotropic hardening. Acta Mater 2006;54:4159–69. [24] Plunkett B, Cazacu O, Lebensohn RA, Barlat F. Elastic–viscoplastic anisotropic modeling of textured metals and validation using the Taylor cylinder impact test. Int J Plast 2007;23:1001–21. [25] Ghaffari Tari D, Worswick MJ, Ali U, Gharghouri MA. Mechanical response of AZ31B magnesium alloy: experimental characterization and material modeling considering proportional loading at room temperature. Int J Plast 2014;55:247–67. [26] Yoon J-H, Cazacu O, Whan Yoon J, Dick RE. Earing predictions for strongly textured aluminum sheets. Int J Mech Sci 2010;52:1563–78. [27] Wang H, Wan M, Wu X, Yan Y. The equivalent plastic strain-dependent Yld2000-2d yield function and the experimental verification. Comput Mater Sci 2009;47:12–22. [28] Yoon JW, Dick RE, Barlat F. A new analytical theory for earing generated from anisotropic plasticity. Int J Plast 2011;27:1165–84. [29] Chung K, Kim D, Park T. Analytical derivation of earing in circular cup drawing based on simple tension properties. Eur J Mech A—Solid 2011;30:275–80. [30] Yoon JW, Barlat F, Dick RE, Karabin ME. Prediction of six or eight ears in a drawn cup based on a new anisotropic yield function. Int J Plast 2006;22:174–93. [31] Xiang-dong W, Min W, Fei H, Hai-bo W. Stress–strain curves for different loading paths and yield loci of aluminum alloy sheets. Trans Nonferr Metal Soc 2006;16:1460–4. [32] Bruschi S, Altan T, Banabic D, Bariani PF, Brosius A, Cao J, Ghiotti A, Khraisheh M, Merklein M, Tekkaya AE. Testing and modelling of material behaviour and formability in sheet metal forming. CIRP Ann-Manuf Technol 2014;63:727–49. [33] Kuwabara T, Ikeda S, Kuroda K. Measurement and analysis of differential work hardening in cold-rolled steel sheet under biaxial tension. J Mater Process Technol 1998;80–81:517–23. [34] Wang H-B, Wan M, Wu X-D, Yan Y. Subsequent yield loci of 5754O aluminum alloy sheet. Trans Nonferr Metal Soc 2009;19:1076–80. [35] Moran B, Ortiz M, Shih CF. Formulation of implicit finite element methods for multiplicative finite deformation plasticity. Int J Numer Methods Eng 1990;29:483–514. [36] Belytschko T, Liu WK, Moran B, Elkhodary K. Nonlinear finite elements for continua and structures. John Wiley & Sons; 2013.