A reduced Yld2004 function for modeling of anisotropic plastic deformation of metals under triaxial loading

Accepted Manuscript

A Reduced Yld2004 Function for Modeling of Anisotropic Plastic Deformation of Metals under triaxial Loading Yanshan Lou , Saijun Zhang , Jeong Whan Yoon PII: DOI: Article Number: Reference:

S0020-7403(19)30314-5 https://doi.org/10.1016/j.ijmecsci.2019.105027 105027 MS 105027

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

25 January 2019 2 July 2019 15 July 2019

Please cite this article as: Yanshan Lou , Saijun Zhang , Jeong Whan Yoon , A Reduced Yld2004 Function for Modeling of Anisotropic Plastic Deformation of Metals under triaxial Loading, International Journal of Mechanical Sciences (2019), doi: https://doi.org/10.1016/j.ijmecsci.2019.105027

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

ACCEPTED MANUSCRIPT

Highlights 

Yld2004-18p function (Barlat et al., 2005) is modified to provide satisfactory predictability of anisotropic behavior of BCC and FCC metals with reduced experimental costs for the calibration of anisotropic parameters.



The reduced Yld2004 function is also implemented into numerical analysis to assess its accuracy and computation efficiency compared with the Hill48 yield



CR IP T

function, the Yld2004-18p function and the anisotropic Drucker yield function. The reduced Yld2004 function provides competitive modeling of anisotropic plastic deformation for BCC and FCC metals with moderate directionality under triaxial stress states, but the cost for the calibration of anisotropic parameters is cut down by about half compared to the Yld2004-18p function.

The reduced Yld2004 function is expected to reliably reproduce actual stamping processes of the majority of advanced high strength steel and aluminum alloys with virtual models.

AC

CE

PT

ED

M

AN US



-1-

ACCEPTED MANUSCRIPT

A Reduced Yld2004 Function for Modeling of Anisotropic Plastic Deformation of Metals under triaxial Loading

CR IP T

Yanshan Lou (a), Saijun Zhang (b), Jeong Whan Yoon (c,d,*)

(a) School of Mechanical Engineering, Xi’an Jiao Tong University, 28 Xianning West Road, Xi’an, Shaanxi 710049, China

AN US

(b) Guangdong Provincial Key Laboratory of Precision Equipment and Manufacturing Technology, School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou 510640, China

M

(c) Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon, 305-701, Republic of Korea (d) School of Engineering, Deakin University, Waurn Ponds, VIC 3220, Australia

AC

CE

PT

ED

* Email: [email protected] / [email protected]

-2-

ACCEPTED MANUSCRIPT

Abstract Solid elements are increasingly utilized in finite element analysis of sheet metal stamping processes with the rising application of advanced high strength steels (AHSSs) and aluminum alloys in auto industries as

CR IP T

these advanced materials fail in a form of rupture with negligible necking during sheet metal stamping. However, there are rare orthotropic yield functions for the description of directional dependence of plastic deformation of body-centered cubic (BCC) materials and face-centered cubic (FCC) metals under spatial loading, especially for metals with moderate anisotropy. For this purpose, the Yld2004-18p function is

AN US

revisited and modified to provide satisfactory predictability of orthotropic behavior of BCC and FCC materials under spatial loading with reduced experimental costs for the calibration of orthotropic coefficients. The reduced Yld2004 equation is utilized to characterize the directional dependence of plastic behavior of AHSSs and aluminum alloys to verify its promising accuracy and cost efficiency. The

M

reduced Yld2004 function is also implemented into finite element analysis to assess its accuracy and

ED

computation efficiency compared with the Hill48, Yld2004-18p and the anisotropic Drucker equations. The application in both analytical prediction and numerical analysis proves that the reduced Yld2004

PT

function provides competitive predictability of anisotropic plastic deformation for BCC and FCC metals with moderate directionality under triaxial stress states, but the cost for the parameter calibration is cut

CE

down by about half compared to the Yld2004-18p function. Considering the satisfactory accuracy, reduced cost for the parameter calibration as well as the fact that most metals are of gentle dependence on

AC

loading direction, the reduced Yld2004 function is expected to reliably reproduce actual stamping processes of the majority of AHSSs and aluminum alloys with virtual models thereby improving the credibility of both numerical design of tools and process optimization. Key Words: Plastic anisotropy, advanced high strength steel, aluminum alloy, body-centered cubic, facecentered cubic, metal forming

-3-

ACCEPTED MANUSCRIPT

1. Introduction As the rising employment of advanced metals for the lightweight design of auto bodies, ductile fracture

is generally accepted as the primary cause of failure during stamping of advanced high strength steel

CR IP T

(AHSS) sheets, aluminum alloys, magnesium alloys, etc. During the last decade, numerous stamping and

plastic deformation processes of advanced metals are numerically modeled by solid elements instead of

shell elements to better estimate rupture initiation during shaping of metal sheets, such as failure

AN US

estimation of three-point bending test of martensitic hat assembly [1], rupture estimation between shear

and plane strain stretching of DP980 [2], failure modeling in distinct loading conditions from shear to the

M

balanced biaxial tension [3], rupture in torsion of sheet metals [4], fracture prediction of bulk metal

forming of AA6061-T6 [5], edge fracture prediction [6], independent validation tests of fracture models

ED

including tension of a hole specimen, a hole expansion test, a hemispherical punch test, etc [7], and

PT

fracture during tube spinning of titanium alloy [8]. Numerical simulation of plastic deformation also

CE

requires the use of solid elements (i.e. fracture in high velocity perforation [9, 10], shear crack prediction

AC

of AA2024-T351 rods in the Taylor test [11], fracture under non-proportional loading for AA2024-T351

[12], and bending of titanium tubes [13]).

Increasing employment of solid elements during numerical simulation for shaping of sheet metals and

bulk metals requires accurate modeling of the directional dependence of plastic behavior for metals under

spatial loading. Among the orthotropic yield functions proposed in the last 70 years since the pioneering

-4-

ACCEPTED MANUSCRIPT

Hill48 yield function [14], more than half of the proposed yield equations are only capable of depict

plastic directionality under plane stress conditions, such as Barlat and Lian [15], Barlat et al. [16],

Banabic et al. [17] and Lou et al. [18]. The rest yield functions [14, 19-29] can model anisotropic

CR IP T

deformation of metals under triaxial loading. The popular yield equations are tabulated in Table 1 for

spatial loading conditions. Most of the functions in Table 1 are proposed for special microstructures of

AN US

metals. However, the proper microstructures are not specified for the Hill48 function, but it is generally

accepted that the Hill48 equation can model anisotropic behavior of steel with high accuracy. The

applicable microstructures are also not specified for the Cazacu-Barlat’2001 functions [23] and that

M

proposed by Yoshida et al. [27] which were both modified from the Drucker equation. These two

ED

anisotropic yield equations are suggested for the modeling of orthotropic plastic deformation of body-

centered cubic (BCC) metals and face-centered cubic (FCC) materials, but not for hexagonal close packed

PT

(HCP) metals. Lou and Yoon [29] differentiated the Drucker-based anisotropic yield equations for BCC

CE

and FCC metals and proved that the stress invariant-based yield equations can reduce 60% of

AC

computation costs compared with non-quadratic equations in a form of principal stresses in case that solid

elements are utilized in simulations. It is very clear that the non-quadratic Yld91, Yld96 and Yld2004-18p

functions as well as that proposed by Aretz and Baralat [30] are for BCC metals when the exponent is six,

and are for FCC materials in case that the exponent is set to eight. The function developed by Karafillis

and Boyce [22] can differentiate the difference in yield envelopes between BCC and FCC materials but

-5-

ACCEPTED MANUSCRIPT

the calibration of its parameters need further investigation. The rest functions of Cazacu and Barlat [24],

Cazacu et al. [25] and Yoon et al. [26] are generally utilized to model the strength differential effect of

orthotropic behavior of FCC metals.

Table 1 Summary of anisotropic yield functions for spatial loading

CR IP T

HCP metals even though it can be found sometimes that these functions are utilized to illustrate

anisotropic parameter # (number in Type of microstructure

blanket is # of anisotropic

AN US

yield function

parameter under plane stress)

Not specified

6 (4)

Yld91 [19]

BCC & FCC

6 (4)

Karafillis and Boyce [22]

BCC & FCC

7 (5)

Yld96 [20]

BCC & FCC

9 (7)

Cazacu and Barlat [23]

M

Hill48 [14]

Not specified

18 (10)

HCP

18 (10)

BCC & FCC

18 (14)

HCP

7 (5)

Not specified

6×n (4×n)

Aretz and Barlat [30]

BBC & FCC

27 (21)

Yoon et al. [26]

HCP

12 (8)

Lou and Yoon [29]

BCC & FCC

6×n (4×n)

Yld2004-18p [21] Cazacu et al. [25]

AC

CE

PT

Yoshida et al. [27]

ED

Cazacu and Barlat [24]

The other issue compared in Table 1 is the number of anisotropic parameters introduced in each yield

functions. It should be stressed that the number of anisotropic parameters in Table 1 is for spatial loading

cases while the number inside the blanket is the number of anisotropic parameters when the metals are

loaded under plane stress. Besides, it must be noticed that the number of anisotropic parameters is flexible

-6-

ACCEPTED MANUSCRIPT

for the functions in Yoshida et al. [27] and Lou and Yoon [29] because these two functions utilized a well-

known approach to improve the predictability of the yield functions by adding n anisotropic Drucker

functions. It is obvious that the accuracy is improved when the number of anisotropic parameters rises,

CR IP T

but the calibration cost for anisotropic parameters also increases with the number of anisotropic

parameters. For example, the calibration of the Yld91 function requires uniaxial tensions along rolling

AN US

direction (RD), diagonal directions (DD) and transverse direction (TD) as well as the equibiaxial tensile

yield stress using bulging tests or the biaxial tensile test [31], but the calibration of the Yld2004-18p

function needs seven uniaxial tensile tests along every 15° from RD, the equibiaxial tensile yield stress

M

and the r-value under equibiaxial tension measured by disk compression [16]. More experiments for

ED

parameter calibration indicate the big increase in calibration costs due to more specimens to be prepared,

increasing time for the conduct of tests and handling of experimental data, etc. Also it should be noted

PT

that the function of Aretz and Barlat [30] has high opportunity to provide higher accuracy for the

CE

modeling of anisotropic plastic deformation of metals since the equation introduces 27 anisotropic

AC

parameters under triaxial loading.

Then the question is whether the improvement of predictability is worth raising experimental costs for

the calibration of anisotropic parameters. There is no doubt that it is worthy of the increased cost for

strong anisotropic metals, such as AA2090-T3, AA5042-H2 and AA3104-H19, because six or eight ears

can only be predicted in cup drawing by high flexible yield functions with improved accuracy [32, 33].

-7-

ACCEPTED MANUSCRIPT

However, anisotropy is not as strong and complicated as that of AA2090-T3 and AA3104-H19 for most

metals. The question above is then questioned for metals with moderate anisotropy. For the moderately

anisotropic metals, yield functions with intermediate accuracy are expected to provide satisfactory

CR IP T

accuracy in numerical analysis of metal stamping and shaping processes, but experimental cost is

significantly reduced for parameter calibration. Anisotropic plastic deformation of metals with moderate

AN US

anisotropy can be represented by eight experimental data: three yield stresses and three r-values under uniaxial tension along RD, DD, and TD (denoted as 𝜎0 , 𝜎45 , 𝜎90 , 𝑟0 , 𝑟45 𝑎𝑛𝑑 𝑟90 ), the equibiaxial tensile

yield stress (𝜎𝑏 ) and the r-value under the balanced biaxial tension (𝑟𝑏 ). For accurate modeling of most or

M

all of these eight experimental results for BCC and FCC materials, it is obvious that the Yld96 equation

ED

and that based on addition of two anisotropic Drucker functions [29] are two good choices if associate

flow rule (AFR) is adopted. However, the numerical application of the Yld96 function is very

PT

complicated and the adding approach is somewhat weird and sophistic to enhance the accuracy of the

CE

orthotropic Drucker equation [29]. With the help of non-associate flow rule (non-AFR), the Yld91

AC

function [19] and the anisotropic Drucker function [29] are recommended for description of anisotropic

plastic deformation of BCC and FCC materials with gentle anisotropy. Nevertheless, the numerical

application of non-AFR needs special care to erase the error caused by the gap between the yield function

and the potential [34]. It needs experiences to guarantee correct implementation of yield functions under

non-AFR.

-8-

ACCEPTED MANUSCRIPT

In this study, a non-quadratic yield equation under AFR is proposed to illustrate anisotropic behavior of

metals with moderate anisotropy with reduced costs for the parameter calibration based on the Yld2004-

18p equation. The cost reduction in parameter calibration is realized by substituting two nine-parameter

CR IP T

tensors of the Yld2004-18p function by two linear transformation tensors with six anisotropic parameters.

The reduced Yld2004 equation is then used to describe orthotropic characteristics of some BCC and FCC

AN US

materials of DP980, AA6022-T4E32 and AA2090-T3 for the verification of its predictability under spatial

loading. The reduced Yld2004 function is further implemented into ABAQUS/Explicit to investigate its

accuracy and computation efficiency compared with the Hill48, Yld2004-18p and anisotropic Drucker

M

functions in numerical simulation of tension of notched specimens and numerical prediction of earing

AC

CE

PT

ED

profiles during cup deep drawing of AA2008.

-9-

ACCEPTED MANUSCRIPT

2. Reduced Yld2004 function The most popular yield functions used now for BCC and FCC metals are proposed in a form of the

non-quadratic Hershey equation as below:

]

2

= 𝜎̅𝐻

(1)

CR IP T

⁄ |𝜎1 −𝜎2 |𝑚+|𝜎2 −𝜎3 |𝑚+|𝜎3 −𝜎1 |𝑚 1 𝑚

𝑓(𝜎𝑖𝑗 ) = [

where 𝜎𝑖𝑗 is the Cauchy stress tensor, 𝜎1 , 𝜎2 and 𝜎3 are three principal stresses, and 𝜎̅𝐻 is the

AN US

Hershey equivalent stress. The exponent 𝑚 is suggested to be six for BCC metals and eight for FCC

materials [35]. Thereafter, many non-quadratic orthotropic equations were proposed to reproduce the

anisotropic characteristics of both BCC and FCC materials in finite element simulation of sheet metal

M

forming.

ED

Another family of yield functions is developed as functions of stress invariants, but these yield

functions do not differentiate distinct plastic behavior between BCC and FCC materials until Lou and

PT

Yoon [29] calibrated the influence of the third stress invariant for BCC and FCC materials in the Drucker

CE

equation [36] as below:

(2)

AC

𝑓(𝜎𝑖𝑗 ) = 𝑎(𝐽23 + 𝑐𝐽32 )1⁄6 = 𝜎̅𝐷

where 𝑎 and 𝑐 are two material coefficients, and 𝐽2 and 𝐽3 are the second invariant and the third invariant of the deviatoric stress tensor. Effect of the third invariant 𝐽3 is adjusted by the material constant 𝑐 which is calibrated to be about 1.226 for BCC materials and 2.0 for FCC metals.

Deng et al. [37] and Kuwabara et al. [38] analyzed yield surfaces of several metals including AA2090-

- 10 -

ACCEPTED MANUSCRIPT

T3, AA6061-O and AA6061-T4 and concluded that the suggested exponent [35] does not provide the best

description for these alloys. Higher exponents were found to give better prediction of biaxial tension yield

stress for these aluminum alloys. Therefore, the upper bound and the lower bound of the non-quadratic

CR IP T

Hershey function in Eq. (1) are illustrated as compared in Figure 1 with the lower and upper bounds of

the Drucker function of Eq. (2). The upper bound is reached for the Hershey yield function when the

AN US

exponent is about 2.767, while its lower bound is the Tresca yield surface at 𝑚 = 0(𝑜𝑟 ∞). For the Drucker equation, the material constant 𝑐 ranges from − 27⁄8 to 9⁄4 in order to guarantee that the function is strictly semi-convex. The Drucker yield surface reaches its maximum when 𝑐 = − 27⁄8, and

M

the smallest yield envelope is observed by the Drucker function at 𝑐 = 9⁄4. The comparison of the upper

ED

and lower bounds of these two yield functions indicates that the non-quadratic Hershey yield function can illustrate lower strength around plane strain tension, while the 𝐽2 -𝐽3 based Drucker function can

PT

reproduce stronger strength of materials under plane strain tension.

CE

Plenty of non-quadratic orthotropic functions were introduced motivated by the Hershey function and

AC

the recommendation of Hosford [35] for the exponent of BCC and FCC materials. One of the most

popular non-quadratic yield equations with high accuracy is the Yld2004-18p function [21] to describe

orthotropic characteristics of metals with strong anisotropy. The Yld2004-18p function is formulated in a

form of

𝑓(𝜎𝑖𝑗 ) = [

𝑚

𝑚

𝑚

|𝑆1̅ ′ −𝑆1̅ ′′ | +|𝑆1̅ ′ −𝑆2̅ ′′ | +|𝑆1̅ ′ −𝑆3̅ ′′ |

𝑚

𝑚

𝑚

𝑚

𝑚

𝑚 1⁄𝑚

+|𝑆2̅ ′ −𝑆1̅ ′′ | +|𝑆2̅ ′ −𝑆2̅ ′′ | +|𝑆2̅ ′ −𝑆3̅ ′′ | +|𝑆3̅ ′ −𝑆1̅ ′′ | +|𝑆3̅ ′ −𝑆2̅ ′′ | +|𝑆3̅ ′ −𝑆3̅ ′′ | 4

]

(3)

- 11 -

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 1 Comparison of the non-quadratic Hershey yield function and the 𝐽2-𝐽3 based Drucker function:

- 12 -

ACCEPTED MANUSCRIPT

(a) on π-plane; and (b) under biaxial loading where 𝑆1̅ ′ , 𝑆2̅ ′ , 𝑆3̅ ′ , 𝑆1̅ ′′ , 𝑆2̅ ′′ and 𝑆3̅ ′′ are the principal values of 𝐒̅′ and 𝐒̅′′ computed by two linear transformations as follows:

𝐒̅′ = 𝐂 ′ 𝐬 = 𝐂 ′ 𝐓𝛔 = 𝐋′ 𝛔

CR IP T

(4)

𝐒̅′′ = 𝐂 ′′ 𝐬 = 𝐂 ′′ 𝐓𝛔 = 𝐋′′ 𝛔

(5)

AN US

with 𝛔 as the Cauchy stress tensor. The matrices of 𝐂 ′ and 𝑪′′ introduce anisotropic coefficients and 𝐓 is the transformation tensor to transform the Cauchy stress tensor to the deviatoric stress tensor, which

are given by 0 0 0 0 ′ 𝑐55 0

0 0 0 ′′ 𝑐44 0 0

0 0 0 0 ′′ 𝑐55 0

0 0 0 0 0 ′ 𝑐66 ]

M

0 0 0 ′ 𝑐44 0 0

′′ −𝑐13 ′′ −𝑐23 0 0 0 0

0 0 0 0 0 ′′ 𝑐66 ]

2 −1 −1 0 0 0 −1 2 −1 0 0 0 1 −1 −1 2 0 0 0 𝐓=3 0 0 0 3 0 0 0 0 0 0 3 0 [0 0 0 0 0 3]

AC

(6)

ED

′′ −𝑐12 0 ′′ −𝑐32 0 0 0

′ −𝑐13 ′ −𝑐23 0 0 0 0

CE

0 ′′ −𝑐21 ′′ −𝑐31 𝐂 ′′ = 0 0 [ 0

′ −𝑐12 0 ′ −𝑐32 0 0 0

PT

0 ′ −𝑐21 ′ −𝑐31 𝐂′ = 0 0 [ 0

(7)

(8)

Therefore, the tensors of 𝐋′ and 𝐋′′ transforming 𝛔 to 𝐒̅′ and 𝐒̅′′ are computed in forms of

- 13 -

ACCEPTED MANUSCRIPT

′′ ′′ 𝑐12 + 𝑐13 ′′ ′′ −2𝑐21 + 𝑐23 ′′ ′′ 1 −2𝑐31 + 𝑐32 𝐋′′ = 𝐂 ′′ 𝐓 = 3 0 0 [ 0

′ ′ −2𝑐12 + 𝑐13 ′ ′ 𝑐21 + 𝑐23 ′ ′ 𝑐31 − 2𝑐32 0 0 0 ′′ ′′ −2𝑐12 + 𝑐13 ′′ ′′ 𝑐21 + 𝑐23 ′′ ′′ 𝑐31 − 2𝑐32 0 0 0

′ ′ 𝑐12 − 2𝑐13 ′ ′ 𝑐21 − 2𝑐23 ′ ′ 𝑐31 + 𝑐32 0 0 0 ′′ ′′ 𝑐12 − 2𝑐13 ′′ ′′ 𝑐21 − 2𝑐23 ′′ ′′ 𝑐31 + 𝑐32 0 0 0

0 0 0 ′ 3𝑐44 0 0 0 0 0 ′′ 3𝑐44 0 0

0 0 0 0 ′ 3𝑐55 0 0 0 0 0 ′′ 3𝑐55 0

0 0 0 0 0 ′ 3𝑐66 ] 0 0 0 0 0 ′′ 3𝑐66 ]

(9)

CR IP T

′ ′ 𝑐12 + 𝑐13 ′ ′ −2𝑐21 + 𝑐23 ′ ′ 1 −2𝑐31 + 𝑐32 𝐋′ = 𝐂 ′ 𝐓 = 3 0 0 [ 0

(10)

AN US

It is proved by Barlat et al. [21] that the function is convex with regard to its arguments for 𝑚 ∈ [0, ∞).

The isotropic Hershey function [39] can be recovered in case that unity is set to all the parameters in the transformation tensors of 𝑪′ and 𝑪′′ . The Yld2004-18p equation reduces to the von Mises function

M

when 𝑚 = 2(𝑜𝑟 4) and with 𝑐𝑖𝑗′ = 𝑐𝑖𝑗′′ = 1, and to the Tresca function when 𝑚 = 0(𝑜𝑟 ∞) and with

ED

𝑐𝑖𝑗′ = 𝑐𝑖𝑗′′ = 1. In spite of the exponent, there are totally 18 anisotropic coefficients introduced in two tensors of 𝑪′ and 𝑪′′ to illustrate the anisotropic characteristics of materials. 14 of these orthotropic

PT

coefficients are related with in-plane orthotropic characteristics of sheet metals and suggested to be

CE

calibrated by using seven uniaxial tensile tests of dogbone specimens cut along every 15° from RD, the

AC

bulging test or the equibiaxial tensile test with the cruciform specimen, and the disk compression test. The

Yld2004-18p function was widely used to model anisotropic characteristics of strong anisotropic

materials of AA2090 T3 [32, 40] and AA5042-H2 [41].

However, anisotropy of most sheet metals is not as strong as AA2090-T3 and AA5054-H2. For metals

with moderate anisotropy, a yield function with less anisotropic coefficients is expected to numerically

- 14 -

ACCEPTED MANUSCRIPT

reproduce the sheet metal forming processes with satisfactory accuracy but is economically effective with

reduced cost and time in experiments, handling of experimental data, calibration of anisotropic

parameters and numerical simulations. Therefore, the Yld2004-18p function in Eq. (3) is reduced to a

0 0 0 0 𝑐5′ 0

0 0 0 0 0 𝑐6′ ]

(11)

AN US

(𝑐2′ + 𝑐3′ )⁄3 −𝑐3′ ⁄3 −𝑐2′⁄3 0 ′⁄ ′ ′ (𝑐3 + 𝑐1 )⁄3 −𝑐3 3 −𝑐1′ ⁄3 0 ′⁄ ′⁄ ′ ′ )⁄ (𝑐 −𝑐 3 −𝑐 3 + 𝑐 3 0 2 1 1 2 𝐋′ = 0 0 0 𝑐4′ 0 0 0 0 [ 0 0 0 0

CR IP T

finance-friendly form by substituting two fourth-order tensors: 𝐋′ and 𝐋′′ in Eqs. (9) and (10) with

0 0 0 0 𝑐5′′ 0

0 0 0 0 0 𝑐6′′ ]

(12)

M

(𝑐2′′ + 𝑐3′′ )⁄3 −𝑐3′′ ⁄3 −𝑐2′′ ⁄3 0 ′′ ⁄ ′′ ′′ (𝑐3 + 𝑐1 )⁄3 −𝑐3 3 −𝑐1′′ ⁄3 0 ′′ ⁄ ′′ ⁄ ′′ ′′ )⁄ (𝑐 −𝑐 3 −𝑐 3 + 𝑐 3 0 2 1 1 2 𝐋′′ = 0 0 0 𝑐4′′ 0 0 0 0 [ 0 0 0 0

PT

𝐒̅′ = 𝐋′ 𝛔, 𝐒̅′′ = 𝐋′′ 𝛔

ED

Then the transformed tensors denoted by 𝐒̅′ and 𝐒̅′′ are computed in forms of

(13)

with the Cauchy stress tensor as

CE

𝛔 = [𝝈𝒙𝒙 , 𝝈𝒚𝒚 , 𝝈𝒛𝒛 , 𝝈𝒚𝒛 , 𝝈𝒛𝒙 , 𝝈𝒙𝒚 ]

𝑻

(14)

AC

In the reduced Yld2004 function in Eq. (3), 𝑆1̅ ′ , 𝑆2̅ ′ , 𝑆3̅ ′ , 𝑆1̅ ′′ , 𝑆2̅ ′′ and 𝑆3̅ ′′ are the principal values of 𝐒̅′ and 𝐒̅′′ calculated in Eq. (13) with the reduced linear transformation tensors of Eqs. (11) and (12). Since

the linear transformation fully considers six components of the Cauchy stress tensor, the reduced Yld2004

function is capable of modeling plastic deformation under triaxial loading conditions. With the reduced

Yld2004 function, 12 anisotropic parameters are introduced in two linear transformation tensors of Eqs.

- 15 -

ACCEPTED MANUSCRIPT

(11) and (12). Among these parameters, four (𝑐4′ , 𝑐5′ , 𝑐4′′ and 𝑐5′′ ) are related with the through-thickness shear properties. These four parameters are suggested to be unity or to be set as 𝑐4′ = 𝑐5′ = 𝑐6′ and 𝑐4′′ = 𝑐5′′ = 𝑐6′′ considering that the through-thickness shear test is a big challenge for sheet metals. The

CR IP T

rest eight anisotropic parameters are used to describe orthotropic plastic deformation in the plane of sheet

metals and calibrated with tensile tests of dogbone specimens along RD, DD and TD, the equibiaxial

AN US

tensile yield stress and the r-value under the balanced biaxial tension from disk compression test. The

calibration of anisotropic coefficients is generally realized through minimization of Eq. (14) below:

𝑒𝑟𝑟 = ∑𝑙𝑖=1 𝑤𝑖 (

𝑝𝑟𝑒𝑑.

𝜎𝑖

𝑒𝑥𝑝. 𝜎𝑖

2

2

𝑝𝑟𝑒𝑑.

𝑟𝑗

− 1) + ∑𝑛𝑗=1 𝑤𝑗 (

𝑒𝑥𝑝. 𝑟𝑗

− 1)

(14)

M

where 𝑤𝑖 and 𝑤𝑗 are the weights which need adjustment for better prediction of anisotropic plastic

ED

behavior of metals. With the minimization of the error function in Eq. (14), the reduced Yld2004 function

AC

CE

PT

is applied to several BCC and FCC materials for its predictability assessment.

- 16 -

ACCEPTED MANUSCRIPT

3. Application of the reduced Yld2004 function to DP980 The reduced Yld2004 function is utilized to illustrate the anisotropic plastic deformation of a dual

CR IP T

phase steel of DP980 (t1.2). Experiments of the steel sheet were conducted by Kuwabara et al. [42].

Tensile tests of dogbone specimens were conducted every 15° from RD at a strain rate of about

0.0005/sec to measure anisotropic behavior of uniaxial tension including the yield stresses and r-values.

AN US

The biaxial stretching experiments were performed using cruciform specimens with stress ratios of 𝜎𝑥𝑥 : 𝜎𝑦𝑦 = 4: 1, 2: 1, 4: 3, 1: 1, 3: 4, 1: 2 and 1: 4 . Cruciform specimens were stretched by a servocontrolled machine [43]. The biaxial tensile stress pairs of 𝜎𝑥𝑥 and 𝜎𝑦𝑦 were measured with respect to

M

plastic deformation during the biaxial tensile test. The flow directions were also recorded for each stress

ED

paths. The uniaxial and biaxial tensile yield stresses were measured when the density of plastic work is

identical with that at 0.002 plastic strain along RD. The r-values were computed when strain is about 0.1

PT

for uniaxial tension, while the ratio of the plastic strain rates of the balanced biaxial tension, 𝑟𝑏 , was

CE

obtained when the plastic work was identical to that of uniaxial tension along RD at plastic strain of 0.002.

AC

All the experimental results are tabulated below in Table 2.

Table 2 Characteristics of DP980 measured from experiments (unit: MPa) [42] 𝜎0

𝜎15

𝜎30

𝜎45

𝜎60

𝜎75

𝜎90

𝜎𝑏

684.2

676.9

669.9

676.6

685.2

691.4

701.5

684.2

𝑟0

𝑟15

𝑟30

𝑟45

𝑟60

𝑟75

𝑟90

𝑟𝑏

0.691

0.760

0.888

1.050

1.080

1.020

0.959

0.703

- 17 -

ACCEPTED MANUSCRIPT

The reduced Yld2004 function is calibrated by experimental data points of tension of dogbone

specimens along three directions of RD, DD and TD, and equibiaxial stretching by minimizing prediction error calculated by Eq. (14). The anisotropic parameters are calibrated with 𝑚 = 6 as summarized in

CR IP T

Table 3 since the dual phase steel is a BCC metal. The calibrated anisotropic coefficients are then

employed to draw the yield surface of the reduced Yld2004 function as compared with experimental

AN US

measurement in Figure 2. In Figure 2 (a), the predicted yield surface is compared with the yield stress

pairs under biaxial tension along RD and TD measured by Kuwabara et al. [42]. The load ratios of biaxial tension are set to be 𝜎𝑥𝑥 : 𝜎𝑦𝑦 = 4: 1, 2: 1, 4: 3, 1: 1, 3: 4, 1: 2 and 1: 4. Uniaxial tensile yield stresses

M

along RD and TD are also presented in Figure 2 (a) for the comparison purpose. Figure 2 (a) also

ED

illustrates the effect of the shear stress on the yield surface of biaxial tension. Effect of shear on the yield

envelope is illustrated in 3D space in Figure 2 (b). The biaxial tensile yield stress pairs and the uniaxial

PT

tensile yield stresses along different loading directions are illustrated in the 3D space of (𝜎𝑥𝑥 , 𝜎𝑦𝑦 , 𝜎𝑥𝑦 ). It

CE

is clearly observed from the comparison between the predicted yield surface and the experimental data

AC

points that biaxial tensile yield stresses are modeled with high accuracy compared with experimental

results when the exponent is six for the metal. The 3D yield envelope is also compared with experimental

data points in Figure 2 (b) with the effect of shear stress. The comparison not only proves the accurate

prediction of the equibiaxial tensile yield stress points but also demonstrates that the reduced Yld2004

equation reasonably describes the directionality of the uniaxial tensile yield stresses.

- 18 -

ACCEPTED MANUSCRIPT

Table 3 Orthotropic coefficients of the reduced Yld2004 function for DP980 𝑐2′

𝑐3′

𝑐6′

0.7093

1.1238

1.3727

1.1486

𝑐1′′

𝑐2′′

𝑐3′′

𝑐6′′

1.1909

0.8949

0.4877

0.8688

𝑚

6

AC

CE

PT

ED

M

AN US

CR IP T

𝑐1′

(a)

- 19 -

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

M

Figure 2 Comparison of the reduced Yld2004 surfaces with experimental results for DP980: (a) contour

considering effect of shear.

ED

of yield envelope with influence of shear; (b) 3D yield surface of the reduced Yld2004 function

PT

Flow directions under uniaxial tension (denoted by r-values) are estimated by the reduced Yld2004

CE

function with anisotropic parameters in Table 3 for comparison with r-values measured from tests in

AC

Figure 3. It is also observed that r-values under uniaxial tension are correctly described by the reduced

Yld2004 equation along RD, DD and TD since these experimental flow directions are utilized during the

calibration of anisotropic coefficients of the reduced Yld2004 equation. The comparison also shows that

the r-values are reasonably modeled by the reduced Yld2004 equation for uniaxial tension along 15°, 30°,

60° and 75° away from RD because the anisotropy of DP980 steel sheet is mild.

- 20 -

ACCEPTED MANUSCRIPT

(b) 1.1 1.0 0.9 0.8

0.6 0.5

Reduced Yld2004 (m = 6) Exp. data

0.4 0.3 0.2 0.1 0.0 0

15

30

45

60

75

90

AN US

Loading angle [0]

CR IP T

r-value

0.7

Figure 3 Comparison of r-values predicted by the reduced Yld2004 equation with experiments for DP980

M

The reduced Yld2004 function is also applied to illustrate the flow directions of plastic strain rates for

the material stretched along RD and TD with anisotropic parameters in Table 3 as compared to

ED

experiments in Figure 4. It is also observed that the reduced Yld2004 function accurately describes the

PT

flow directions at 𝜑 = 00 , 450, and 900 . It is because 𝜑 = 00 corresponds to the uniaxial tension along

CE

RD, 𝜑 = 450 denotes the balanced biaxial tension with 𝜎𝑥𝑥 = 𝜎𝑦𝑦 = 𝜎𝑏 , and 𝜑 = 900 is homologous

AC

with the uniaxial tension along TD. That is, the flow directions of these three loading conditions are corresponding with 𝑟0 , 𝑟𝑏 and 𝑟90 . Accordingly, it is natural for the good prediction of these three

experimental data points considering that these experimental results are utilized in the parameter

calibration of the reduced Yld2004 function. It is also noted that the reduced Yld2004 function accurately

describes the other flow directions of plastic strain rates for materials stretched biaxially.

- 21 -

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Figure 4 Comparison of the flow directions of plastic deformation under biaxial tension predicted with

AC

CE

the reduced Yld2004 function with experimental results for DP980

- 22 -

ACCEPTED MANUSCRIPT

4. Application of the reduced Yld2004 function to AA6022-T4E32 The reduced Yld2004 equation is employed to model orthotropic mechanical properties of AA6022-

CR IP T

T4E32 with experimental results in Stoughton and Yoon [44]. The reduced Yld2004 function is calibrated

by experimental results and the optimized parameters are tabulated in Table 4. The reduced Yld2004

envelope is compared with experimental results of uniaxial tension along every 15° and the balanced

AN US

biaxial tension in Figure 5. The isolines in Figure 5 (b) indicate the equivalent yield envelope under the

identical effect of the shear stress. Take the uniaxial tension along 30° as an example, the ratio of the shear stress to the uniaxial tensile yield stress of 𝜎30 is 𝜏⁄𝜎30 = sin(30°) ∗ cos(30°) = √3⁄4 ≈ 0.433.

M

It is obvious that 𝜎60 is on the same isoline of 𝜎30, and 𝜎15 and 𝜎75 are on the same isoline with

ED

𝜏⁄𝜎15 = 0.25. Here 𝜎𝜃 represents uniaxial tensile yield stress along 𝜃-direction. The uniaxial tensile

PT

yield stress of 𝜎45 is on the isoline of 𝜏⁄𝜎45 = 0.5.

CE

Obviously the reduced Yld2004 function accurately depicts the strength dependence on loading

directions of the AA6022-T4E32 alloy under uniaxial tension along distinct stretching orientations and

AC

equibiaxial tension. The anisotropy of AA6022-T4E32 is also modeled by the Hill48, Yld2004-18p and

orientation-dependent Drucker functions for the purpose of comparison. The predicted yield envelopes

are compared to the reduced Yld2004 yield surface and experimental results in Figure 6. The comparison

reveals that all the equations reasonably depict yield stresses of uniaxial tension along two directions of

RD and TD and equibiaxial tensile yield stress. The difference among different yield functions lies in the

- 23 -

ACCEPTED MANUSCRIPT

strength between uniaxial and equibiaxial tension. However, the highest strength of plane strain stretching

is predicted by the Hill48 yield equation, while the Yld2004-18p and anisotropic Drucker functions model

poor strength of plane strain tension. The strength under plane strain stretching illustrated with the

CR IP T

reduced Yld2004 equation is between that modeled by the Hill48 yield equation and those of the

Yld2004-18p and anisotropic Drucker functions. It is obvious that the curvature between uniaxial and

AN US

equibiaxial tension can be modified by changing the exponent of the non-quadratic yield equations

including the Yld2004-18p, Yld2000-2d and reduced Yld2004 equations. Therefore, the influence of the

exponent on the reduced Yld2004 yield envelope is investigated by changing the value of the exponent

M

during the calibration. Two sets of anisotropic parameters are calibrated for two exponents of 20 and 50 as

ED

summarized in Table 4. The calibration of the eight anisotropic parameters in the reduced Yld2004

function is conducted by minimization of the error function of the predicted yield stresses with the

PT

downhill Simplex approach, which is identical to the calibration of the anisotropic parameters when the

CE

exponent is six or eight. With these calibrated anisotropic parameters, the reduced Yld2004 yield surfaces

AC

for different exponents are compared with measurement from experiments in Figure 7. It is noted that the

reduced Yld2004 function predicts the uniaxial tensile yield stresses along RD and TD and equibiaxial

tension yield stress with good agreement no matter whether the exponent is 8, 20 or 50, but the predicted

yield stresses under plane strain tension slightly drops with the rise of the exponent for calibration of

orthotropic coefficients of the reduced Yld2004 equation. It is widely accepted that the exponent in non-

- 24 -

ACCEPTED MANUSCRIPT

quadratic yield function is six for BCC and eight for FCC metals. However, the strength under plane

strain tension is also affected by textures formed during manufacturing processes. Biaxial tensile test

results of Kuwabara et al. [38] for AA6016-O and AA6016-T4 suggest that the exponent should be higher

CR IP T

than eight for accurately modeling of the strength of these alloys between uniaxial tension and balanced

biaxial tension. Accordingly, plane strain tensile tests or biaxial tension of cruciform specimens are

AN US

recommended to be conducted for accurate calibration of the exponent in the reduced Yld2004 function.

For the purpose of simplicity and cost control, the widely accepted value of the exponents can still be

used for the non-quadratic function.

M

The uniaxial tensile yield stresses are predicted and plotted together with experimental results in Figure

ED

8 (a) while comparison between the predicted r-values and experimental measurement is illustrated in

PT

Figure 8 (b). The comparison validates that all the equations predict the Lankford parameter with high

accuracy and the difference is negligible for the r-values modeled by different yield functions.

CE

Concerning the predictability of the uniaxial tensile yield stresses compared in Figure 8 (a), the

AC

anisotropic Drucker yield function (𝑛 = 3) provides the best prediction and accuracy of the other three

equations is similar. Improvement of the anisotropic Drucker yield equation is attributed to the increased

number of anisotropic coefficients calibrated by tension of dogbone specimens along every 15° from RD,

which causes dramatic rise in experimental cost for parameter calibration of the anisotropic Drucker

function. Besides, improvement of the predictability is negligible for the anisotropic uniaxial tensile yield

- 25 -

ACCEPTED MANUSCRIPT

stresses compared with the Hill48, Yld2004-18p and reduced Yld2004 equations.

Table 4 Anisotropic coefficients of the reduced Yld2004 equation for AA6022-T4E32 𝑐1′

𝑐2′

𝑐3′

𝑐6′

-0.0561

-0.1425

1.7688

1.3016

𝑐1′′

𝑐2′′

𝑐3′′

1.5365

1.3932

0.2691

0.2852

𝑐1′

𝑐2′

𝑐3′

𝑐6′

-0.0402

-0.0852

1.7067

1.1572

𝑐1′′

𝑐2′′

𝑐3′′

𝑐6′′

1.5086

1.3588

0.1103

0.4148

𝑐1′

𝑐2′

𝑐3′

𝑐6′

-0.0211

-0.0424

1.6812

1.1010

𝑐1′′

𝑐2′′

𝑐3′′

𝑐6′′

1.5023

1.3499

0.0391

0.4502

CR IP T

𝑚

𝑐6′′

8

𝑚

AN US

20

𝑚

M

50

90

PT

Reduced Yld2004 yield surface

-1.00

1.00

75 60

0.75

0.50

CE AC

-1.25

b

yy

ED

1.25

30 45

0.25

15 0

-0.75

-0.50

-0.25

0.25

0.50

0.75

1.00

1.25

xx -0.25

-0.50

-0.75

-1.00

-1.25

- 26 -

ACCEPTED MANUSCRIPT

AN US

CR IP T

(a)

(b)

M

Figure 5 Comparison of the reduced Yld2004 surfaces with experimental data points for AA6022-T4E32:

ED

(a) contour of yield envelope with effect of shear; (b) 3D yield surface of the reduced Yld2004 function

AC

CE

PT

considering the influence of the shear stress

- 27 -

1.25 Hill48 plastic potential Yld2004-18p Anisotropic Drucker (n = 3) Reduced Yld2004 (m = 8) Exp. data

yy

ACCEPTED MANUSCRIPT

1.00

0.75

0.50

0.25

-1.25

-1.00

-0.75

-0.50

-0.25

0.25

0.50

-0.25

-0.50

0.75

1.00

1.25

AN US

-0.75

CR IP T

xx

-1.00

-1.25

Figure 6 Comparison of yield envelopes predicted by various yield equations with experimental results of

1.25

Reduced Yld2004 yield surface

ED

1.00

m=8 m = 20 m = 50

-1.00

90 75 60

0.75

PT CE AC

-1.25

b

yy

M

AA6022-T4E32

0.50

30 45

0.25

15 0

-0.75

-0.50

-0.25

0.25

0.50

0.75

1.00

1.25

xx

-0.25

-0.50

-0.75

-1.00

-1.25

Figure 7 Effect of the exponent on the reduced Yld2004 yield surfaces calibrated for AA6022-T4E32

- 28 -

ACCEPTED MANUSCRIPT

1.01

(a)

0.99

0.98

CR IP T

Normalized uniaxial yield stress

1.00

Hill48 plastic potential Yld2004-18p Anisotropic Drucker (n = 3) Reduced Yld2004 (m = 8) Exp. data

0.97

0.96

0.94 0

15

30

AN US

0.95

45

60

75

90

60

75

90

0

Loading angle [ ]

0.9

M

(b)

0.8

ED

0.7 0.6

PT

0.4

Hill48 plastic potential Yld2004-18p Anisotropic Drucker (n = 3) Reduced Yld2004 (m = 8) Exp. data

CE

r-value

0.5

0.3

AC

0.2 0.1 0.0

0

15

30

45 Loading angle [0]

Figure 8 Comparison of the directional dependence computed by various yield equations with

experimental data points for AA6022-T4E32: (a) yield stress at uniaxial tensile; and (b) r-values

- 29 -

ACCEPTED MANUSCRIPT

5. Application of the reduced Yld2004 Function to AA2090-T3

The reduced Yld2004 function is applied to illustrate orthotropic plastic deformation of a strong

CR IP T

anisotropic aluminum alloy of AA2090-T3. The yield stresses and r-values of uniaxial tension are

referred from Yoon et al. [32, 45] as tabulated in Table 5, while the yield stresses under biaxial tension

were measured by Deng et al. [2015]. The reduced Yld2004 function is calibrated with the experimental

AN US

data points in Table 5 with calibrated orthotropic coefficients in Table 6. The exponent of the reduced

Yld2004 function is also taken as a variable during calibration. With the calibrated parameters, the

reduced Yld2004 yield surface is constructed with the influence of shear stress on the yield surfaces as

M

compared with the biaxial tensile data points in Figure 9 (a) and (b) for 2D and 3D representations. It is

ED

obvious that the reduced Yld2004 function illustrates the balanced biaxial tensile yield stresses with high

PT

accuracy compared with testing results. From the observation in Figure 9 (b), it is easy to conclude that

the reduced Yld2004 function slightly underestimates yield stress of dogbone specimens of 30° from RD,

CE

but overestimates 𝜎60. Overestimate of 𝜎60 is the reason why the data point of 𝜎60 is inside of the

AC

reduced Yld2004 surface in Figure 9 (a). The uniaxial tensile yield stresses along other stretching

directions are all predicted with acceptable error by the reduced Yld2004 function compared with

experimental results.

- 30 -

ACCEPTED MANUSCRIPT

Table 5 Experimental results of AA2090-T3 [32, 45] 𝜎15 ⁄𝜎0

𝜎30⁄𝜎0

𝜎45⁄𝜎0

𝜎60⁄𝜎0

𝜎75⁄𝜎0

𝜎90⁄𝜎0

𝜎𝑏 ⁄𝜎0

1.0000

0.9605

0.9102

0.8114

0.8096

0.8815

0.9102

1.0350

𝑟0

𝑟15

𝑟30

𝑟45

𝑟60

𝑟75

𝑟90

𝑟𝑏

0.2115

0.3269

0.6923

1.5769

1.0385

0.5384

0.6923

0.6700

CR IP T

𝜎0⁄𝜎0

Table 6 Anisotropic parameters of the reduced Yld2004 function for AA2090-T3 () 𝑐2′

𝑐3′

𝑐4′

𝑐5′

1.3239

0.6896

-2.5993

1.6951

1.6951

1.6951

𝑐1′′

𝑐2′′

𝑐3′′

𝑐4′′

𝑐5′′

𝑐6′′

1.0478

0.6616

0.1340

0.7375

0.7375

0.7375

𝑐6′

𝑚

14.2115

AC

CE

PT

ED

M

AN US

𝑐1′

- 31 -

Reduced Yld2004 yield surface 1.00

b 90 75 60

0.75

0.50

0.25

-1.25

-1.00

-0.75

-0.50

-0.25

30 15

45

0.25

0

0.50

CR IP T

1.25

yy

ACCEPTED MANUSCRIPT

0.75

1.00

1.25

xx

-0.25

-0.50

AN US

-0.75

-1.00

-1.25

M

Figure 9 Comparison of the reduced Yld2004 surface with experimental results for AA2090-T3: (a) 3D

ED

representation; and (b) 2D representation

The predictability of the reduced Yld2004 function is compared to the Hill48, Yld2004-18p and

PT

anisotropic Drucker equations for AA2090-T3. It should be noted that the Hill48 yield function is

CE

calibrated by the yield stresses of uniaxial tension along RD, DD and TD, and the equibiaxial tensile yield

AC

stress, while the anisotropic Drucker equation is used under the scheme of non-AFR. Details of the

anisotropic Drucker equation under non-AFR are referred from Lou and Yoon [29]. The anisotropic

parameters of the Yld2004-18p equation is referred from Barlat et al. [21] as summarized in Table 7,

while the anisotropic Drucker equation is calibrated again with the calibrated parameters in Table 8. The

calibrated coefficients of the anisotropic Drucker equation are different from those in [29] for the better

- 32 -

ACCEPTED MANUSCRIPT

modeling of yield surfaces under balanced biaxial tension.

The predicted strength under biaxial stretching along RD and TD is compared with experimental

results in Figure 10 for these four yield functions. It is evident that all the equations predict strength of

CR IP T

uniaxial stretching of RD and TD and equibiaxial tensile strength perfectly. However, the Hill48 yield

function describes higher strength from uniaxial tension to the equibiaxial tension for RD and TD. The

AN US

Yld2004-18p and anisotropic Drucker functions model the strength between uniaxial and equibiaxial

tension with the best agreement compared to experimental results for both directions. The reduced

Yld2004 function perfectly depicts the strength around plane strain tension along RD, but slightly

ED

M

overestimates the yield stress between uniaxial and equibiaxial tension stretched along TD.

Table 7 Anisotropic coefficients of Yld2004-18p for AA2090-T3 with 𝑚 = 8 ′ 𝑐21

′ 𝑐23

′ 𝑐31

′ 𝑐32

′ 𝑐44

′ 𝑐55

′ 𝑐66

0.9364

0.0791

1.0030

0.5247

1.3631

1.0237

1.0690

0.9543

′′ 𝑐13

′′ 𝑐21

′′ 𝑐23

′′ 𝑐31

′′ 𝑐32

′′ 𝑐44

′′ 𝑐55

′′ 𝑐66

0.5753

0.8668

1.1450

-0.0792

1.0516

1.1471

1.4046

′ 𝑐13

-0.0698 ′′ 𝑐12

CE

PT

′ 𝑐12

0.4767

AC

0.9811

- 33 -

ACCEPTED MANUSCRIPT

Table 8 Anisotropic coefficients of anisotropic Drucker equation (𝑛 = 2) and plastic potential equation (𝑛 = 2) for AA 2090 T3 (𝑐 ′ = 𝑐̂ = 2) ′(1)

′(1)

𝑐2

𝑐3

′(1)

𝑐4

′(1)

′(1)

𝑐5

𝑐6

′(2)

𝑐1

′(2)

𝑐2

′(2)

′(2)

′(2)

′(2)

CR IP T

′(1)

𝑐1

𝑐3

𝑐4

𝑐5

𝑐6

0.0730 0.4181 3.0191 2.5942 2.5942 2.5942 3.9787 2.4289 -3.0597 1.4032 1. 4032 1. 4032 (1)

𝑐̂1

(1)

(1)

𝑐̂2

𝑐̂3

(1)

𝑐̂4

(1)

(1)

𝑐̂5

𝑐̂6

(2)

𝑐̂1

(2)

𝑐̂2

(2)

𝑐̂3

(2)

(2)

𝑐̂4

𝑐̂5

(2)

𝑐̂6

1.25

Hill48 Yld2004-18p ADrucker yield surface ADrucker potential reduced Yld2004 Exp. data

yy

AN US

-4.5160 2.1525 3.1762 2.9128 2.9128 2.9128 1.4151 1.8376 -1.3083 1.4843 1. 4843 1. 4843

1.00

M

0.75

0.50

-1.00

-0.75

AC

CE

PT

-1.25

ED

0.25

-0.50

xx

-0.25

0.25

0.50

0.75

1.00

1.25

-0.25

-0.50

-0.75

-1.00

-1.25

Figure 10 Predictability assessment of various yield functions on the modeling of yielding under biaxial loading along RD and TD (ADrucker: anisotropic Drucker)

These four yield functions are evaluated for their modeling of anisotropic behavior of uniaxial tension

- 34 -

ACCEPTED MANUSCRIPT

in Figure 11 (a) for uniaxial tensile yield stresses and in Figure 11 (b) for r-values of uniaxial tension.

Obviously, the Hill48, anisotropic Drucker and reduced Yld2004 equations predict the orthotropic

uniaxial stretching yield stresses with satisfactory accuracy. The Yld2004-18p function provides the most

CR IP T

accurate description for the anisotropic yielding of AA2090-T3, but it requires more experiments for the

calibration of its anisotropic coefficients. The estimated Lankford parameters are compared with testing

AN US

data points in Figure 11 (b). Obviously the Hill48 equation cannot model the flow directions of plastic

deformation of uniaxial stretching when the Hill48 equation is calibrated using uniaxial yield stresses

along different orientations, and equibiaxial tensile yield stress. Therefore, the Hill48 equation is

M

recommended to be used together with non-AFR for proper modeling of the anisotropic behavior of both

ED

strength and flow directions of a metal. The Yld2004-18p equation and the anisotropic Drucker potential

predict the r-values of uniaxial stretching along different tensile directions with best agreement compared

PT

with experimental results. The reduced Yld2004 function is noted to be capable of accurate illustration of

CE

the anisotropy of the flow directions under uniaxial tension with less anisotropic parameters to be

AC

calibrated by experiments than the Yld2004-18p function and the anisotropic Drucker potential.

- 35 -

ACCEPTED MANUSCRIPT

1.02

(a)

1.00

Hill48 Yld2004-18p ADrucker reduced Yld2004 Exp. data

0.96 0.94 0.92 0.90 0.88 0.86 0.84 0.82 0.80 0.78 0

15

30

45

60

2.4

(b)

2.2 2.0 1.8

90

M

1.6 1.4

r-value

75

AN US

Loading angle [0]

CR IP T

Normalized uniaxial yield stress

0.98

1.2

ED

1.0 0.8 0.6

PT

0.4 0.2

Hill48 Yld2004-18p ADrucker reduced Yld2004 Exp. data

0.0

CE

0

15

30

45

60

75

90

Loading angle [0]

AC

Figure 11 Evaluation of the Hill48, Yld2004-18p and reduced Yld2004 equations for uniaxial tensile tests: (a) uniaxial tensile yield stresses; and (b) Lankford parameters

- 36 -

ACCEPTED MANUSCRIPT

7. Implementation into numerical simulation

A user subroutine is written for a commercial software of ABAQUS/Explicit to implement the reduced

CR IP T

Yld2004 equation for numerical analysis of metal forming. For the comparison purpose, numerical

simulations are also conducted by the Hill48, Yld2004-18p and anisotropic Drucker equations. Plastic increments are calculated using the Euler’s backward approach and return mapping method [34]. Tension

AN US

of the notched specimen in Lou and Yoon [46] is simulated with these four yield functions. Taking the

advantage of symmetry, FE model is developed for one eighth of notched specimens, and symmetric

boundary condition is thereby employed for nodes on the surface of x, y and z in Figure 12 during

M

numerical simulation. The numerical test stretches the specimen at a fixed velocity of 0.01 mm/sec along

ED

the longitudinal direction by setting the velocity boundary condition for nodes on the top surface in

PT

Figure 12 along x-direction. The specimen is assumed to be manufactured from AA2090-T3 sheet whose

thickness is 1.0 mm. The simulations are conducted for the Hill48, Yld2004-18p, anisotropic Drucker and

CE

reduced Yld2004 functions with anisotropic parameters in Section 6. Strain hardening of AA2090-T3 is

AC

assumed to follow the Swift rule of 𝜎̅ = 646 × (0.025 + 𝜀̅)0.227 MPa. Essential conditions for numerical

analysis are tabulated in Table 9.

The computation time is compared in Table 10 for simulations with these four equations. The

simulation with the Hill48 yield function is observed to take the shortest time due to its simplicity. The

computation cost of the anisotropic Drucker equation with non-AFR is about two times of the Hill48

- 37 -

ACCEPTED MANUSCRIPT

yield function because of the high order of the anisotropic Drucker equation. Simulations with the

Yld2004-18p and reduced Yld2004 equations take the longest time with little difference because of both

the high order of these two equations and the complication in the calculation of principal stresses under

CR IP T

spatial loading.

The distribution of the von Mises effective stress is extracted from the numerical simulation of the

AN US

stretching of notched specimens at a stroke of 1.2 mm as shown in Figure 13. The comparison

demonstrates that the distribution of the effective stress estimated by the reduced Yld2004 equation is

very close to that estimated by the Yld2004-18p equation. Besides, the results show that the highest

M

effective stress is predicted by the anisotropic Drucker equation with non-AFR represented in gray at the

ED

outer corner of the notch in Figure 13 (c), but the Hill48 yield function estimates the largest area of the

von Mises equivalent stress between 0.45 GPa and 0.5 GPa denoted in orange in Figure 13 (a). The

PT

highest von Mises equivalent stress estimated by the anisotropic Drucker equation is attributed to the fact

CE

that the anisotropic Drucker equation is used together with non-AFR in the numerical simulation, while

AC

simulations with the other three yield functions is coupled with AFR.

The load-stroke curves estimated by different equations are compared in Figure 14. The comparison

shows that the load-stroke curve estimated by the Hill48 equation is about 7% larger than those predicted

by the other three equations. The predicted load capability of the notched specimen is negligible

difference for the other three yield functions. This is due to the fact that Lode parameter is approximately

- 38 -

ACCEPTED MANUSCRIPT

about zero for most part of the notch as depicted by green for the distribution of the Lode parameter in

Figure 14. The zero Lode parameter indicates that the loading condition in the notch of the specimen is

very close to the normalized plane strain [47, 48]. Strength under plane strain tension along RD estimated

CR IP T

by the Hill48 equation is 15% higher than that characterized by the other three yield functions, as

observed in the comparison of four yield surfaces with experimental results in Figure 10. Thus the 7%

AN US

overall difference in load capability of the notched specimen is the result of the complex stress state of the

notched specimens from uniaxial tension at the edge of outer notch to plane strain tension for central part

of notch.

M

The numerical simulation with the reduced Yld2004 function and the other three yield functions proves

ED

that the reduced Yld2004 yield function provides reasonable prediction of the strength around plane strain

tension compared with the Yld2004-18p equation and the anisotropic Drucker equation with non-AFR.

PT

However, there is no obvious difference in computation cost for numerical simulation between the

CE

reduced Yld2004 function and Yld2004-18p equation, both of which are higher than simulation with the

AC

anisotropic Drucker equation under non-AFR. The numerical simulation also indicates that tension of the

notched specimen is a simple but efficient method to calibrate the strength of plane strain tension when

the complicated biaxial stretching test of the cruciform specimen is not available.

- 39 -

ACCEPTED MANUSCRIPT

Table 9 Summary of numerical simulations of tension of notched specimens Element number

Elastic modulus

Poisson ratio

6705

5504

0.69 GPa

0.33

Memory

software

System

Processor

8.0 GB

ABAQUS/Explicit

64-bit Windows 7

Single Precision

CR IP T

Node number

Intel (R) Core (TM) i5-6500 CPU @ 3.2GHz 3.2GHz

tension of notched specimens

AN US

Table 10 Comparison of the computation time between different yield functions for the simulations of

Anisotropic Drucker

Yld2004-18p

10 min 13 sec

31 min 39 sec

under non-AFR 20 min 27 sec

Reduced Yld2004 30 min 44 sec

CE

PT

ED

M

Hill48

AC

Figure 12 Boundary and loading conditions for tensile tests of one eighth FE model of the notched specimen

- 40 -

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

- 41 -

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 13 von Mises effective stress distribution in numerical simulation of tension of notched specimens at a stretch stroke of 1.2 mm predicted by different equations: (a) Hill48; (b) Yld2004-18p; (c) anisotropic Drucker; and (d) reduced Yld2004.

- 42 -

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 14 Load-displacement curve comparison of notched specimens predicted by different yield

AC

CE

𝜎1 ≥ 𝜎2 ≥ 𝜎3)

PT

functions (Lode parameter denoted by 𝐿 is calculated as 𝐿 = (2𝜎2 − 𝜎1 − 𝜎3 )⁄(𝜎1 − 𝜎3 ) with

- 43 -

ACCEPTED MANUSCRIPT

8. Earing prediction for AA2008-T4

For the further assessment of predictability of the reduced Yld2004 equation, the earing profile is

CR IP T

predicted numerically for an aluminum alloy of 2008-T4 during cup deep drawing by the reduced

Yld2004 function implemented into the commercial software of ABAQUS/Explicit above. Cup deep

drawing of AA2008-T4 is also simulated by using the Yld91 equation with non-AFR, the anisotropic

AN US

Drucker equation with non-AFR (n=1 and n=2), and the Yld2004-18p function for the comparison

purpose. The setting up of the tools and their dimensions are given in Yoon et al. [49]. The thickness of

AA2008-T4 before deformation is 1.24 mm and original diameter is 162 mm. The friction is set to 0.05

M

between blank and die and between blank and holder, and 0.2 between punch and blank. The total holding

ED

force is 22.2kN and therefore 5.55kN is used as blank holding force in numerical simulation since only

PT

quarter model is utilized in numerical analysis taking the advantage of the symmetry of blank properties

with respect to RD and TD. The detailed geometry and tests conditions are suggested to refer to Yoon et

CE

al. [49].

AC

The material used in the simulation is AA2008-T4, which is an FCC metal with moderate anisotropy. Strain hardening is modeled by the Voce equation of 𝜎̅ = 408 − 223exp(−6.14𝜀̅)𝑀𝑃𝑎, and isotropic

hardening is assumed for the alloy. Its anisotropy in plastic deformation and strength is measured by

experiments as summarized in Table 11. The experimental data points are then used to calibrate the Yld91 equation with non-AFR, anisotropic Drucker equation with non-AFR (𝑛 = 1 and 𝑛 = 2), the Yld2004-

- 44 -

ACCEPTED MANUSCRIPT

18p equation and the reduced Yld2004 equation. Table 12 summarizes the orthotropic coefficients for the Yld91 equation with non-AFR, Table 13 for the anisotropic Drucker equation with non-AFR with 𝑛 = 1, Table 14 for the anisotropic Drucker equation with non-AFR with 𝑛 = 2, Table 15 for the Yld2004-18p

CR IP T

equation, and Table 16 for reduced Yld2004 equation.

Table 11 Experimental results of AA2008-T4 [49] 𝜎15 ⁄𝜎0

𝜎30⁄𝜎0

𝜎45⁄𝜎0

𝜎60⁄𝜎0

𝜎75⁄𝜎0

𝜎90⁄𝜎0

𝜎𝑏 ⁄𝜎0

1.0000

0.9963

0.9835

0.9459

0.9303

0.9171

0.9044

0.9010

𝑟0

𝑟15

𝑟30

𝑟45

0.8674

0.8077

0.6188

0.4915

AN US

𝜎0⁄𝜎0

𝑟60

𝑟75

𝑟90

𝑟𝑏

0.4955

0.5114

0.5313

1.0000

M

Table 12 Anisotropic coefficients of the Yld91 equation with non-AFR for AA2008-T4 with 𝑚 = 8 𝑦

𝑦

𝐶1

yield function

1.2080

Parameters for

𝐶1

potential

1.0890

1.0045 𝑝

𝑦

𝑦

𝑦

𝑦

𝐶3

𝐶4

𝐶5

𝐶6

0.9955

1.0391

1.0391

1.0391

𝑝

𝑝

𝑝

𝑝

𝐶2

𝐶3

𝐶4

𝐶5

𝐶6

1.0194

0.9803

0.9711

0.9711

0.9711

CE

PT

𝑝

𝐶2

ED

Parameters for

AC

Table 13 Anisotropic coefficients of the anisotropic Drucker equation (𝑛 = 1) and plastic potential (𝑛 = 1) for AA2008-T4 (𝑐 ′ = 𝑐̂ = 2) 𝑦

𝑦

𝑦

𝑦

𝑦

𝑦

Parameters for

𝐶1

𝐶2

𝐶3

𝐶4

𝐶5

𝐶6

yield function

2.2190

1.8448

1.8282

1.9082

1.9082

1.9082

Parameters for

𝑝 𝐶1

𝑝 𝐶2

𝑝 𝐶3

𝑝 𝐶4

𝑝 𝐶5

𝐶6

potential

2.1913

1.8729

1.7995

1.7829

1.7829

1.7829

𝑝

- 45 -

ACCEPTED MANUSCRIPT

Table 14 Anisotropic coefficients of the anisotropic Drucker equation with 𝑛 = 2 and plastic potential (𝑛 = 2) of AA2008-T4 (𝑐 ′ = 𝑐̂ = 2) 𝑦(1)

𝑐1

𝑦(1)

𝑦(1)

𝑐2

𝑐3

𝑦(1)

𝑐4

𝑦(1)

𝑦(1)

𝑐5

𝑦(2)

𝑐6

𝑦(2)

𝑐1

𝑐2

𝑦(2)

𝑐3

𝑦(2)

𝑦(2)

𝑦(2)

𝑐5

𝑐4

𝑐6

2.3007 2.5455 -3.7739 1.8589 1.8589 1.8589 2.0831 1.1176 2.1676 -1.8779 -1.8779 -1.8779 𝑝(1)

𝑝(1)

𝑐2

𝑐3

𝑝(1)

𝑐4

𝑝(1)

𝑝(1)

𝑐5

𝑝(2)

𝑐6

𝑝(2)

𝑐1

𝑐2

𝑝(2)

𝑐3

𝑝(2)

𝑝(2)

𝑐4

𝑝(2)

𝑐5

𝑐6

CR IP T

𝑝(1)

𝑐1

-0.0985 0.5202 0.2988 0.1809 0.1809 0.1809 4.4777 3.2780 3.2277 3.4380 3.4380 3.4380

′ 𝑐13

′ 𝑐21

′ 𝑐23

1.0601

1.0407

1.1240

1.1968

′′ 𝑐12

′′ 𝑐13

′′ 𝑐21

′′ 𝑐23

0.9218

1.0694

0.9894

1.1843

′ 𝑐31

′ 𝑐32

′ 𝑐44

′ 𝑐55

′ 𝑐66

1.0225

1.0270

1.3407

-0.2526

1.2443

′′ 𝑐31

′′ 𝑐32

′′ 𝑐44

′′ 𝑐55

′′ 𝑐66

1.0237

1.1199

0.9354

0.8848

0.7336

ED

M

′ 𝑐12

AN US

Table 15 Anisotropic coefficients of the Yld2004-18p equation for AA2090-T3 with 𝑚 = 8

Table 16 Anisotropic coefficients of reduced Yld2004 𝑐2′

𝑐3′

𝑐4′

PT

𝑐1′

𝑐5′

𝑐6′

𝑐1′′

𝑐2′′

𝑐3′′

𝑐4′′

𝑐5′′

𝑐6′′

8

CE

1.7643 1.2895 0.6796 0.8267 0.8267 0.8267 0.6689 0.3646 -1.7634 1.1057 1.1057 1.1057

𝑚

AC

With the anisotropic parameters calibrated for AA2008-T4 above, directionality in plastic flow and

uniaxial tensile yield stress is predicted and compared with the measurements from experiments in Figure

15(a) and (b). Predicted yield loci are compared in Figure 16 and the potential loci are plotted in Figure

17 for the three yield functions under non-AFR. It is observed that the difference between the Yld91 function and the anisotropic Drucker equation (𝑛 = 1) with non-AFR is barely noticeable in predicted

- 46 -

ACCEPTED MANUSCRIPT

uniaxial tensile yield stress, r-values, yield envelope and potential loci. It means that the predictability of the Yld91 and the anisotropic Drucker equation (𝑛 = 1) is almost equal. The Yld2004 equation and the anisotropic Drucker equation (𝑛 = 2) with non-AFR provide slightly better depiction of the anisotropic

CR IP T

behavior of uniaxial tensile strength and r-values as observed in Figure 15, but the improvement is not as

obvious as that for AA2090-T3 in Figure 11. This is because the anisotropy of AA2008-T4 is not as

AN US

strong as AA2090-T3, so the advantage of the Yld2004-18p function cannot be observed for metals with

moderate anisotropy (i.e. AA2008-T4).

(a)

M

0.98

ED

0.96

0.94

PT

Yld91 Drucker nonAFR (n=1) Drucker nonAFR (n=2) Yld2004-18p Reduced Yld2004 Exp.

0.92

CE

Normalized uniaxial yield stress

1.00

0.90

AC

0

15

30

45

60

75

90

Loading angle [0]

- 47 -

ACCEPTED MANUSCRIPT

1.0

(b)

0.9 0.8 0.7

0.5

CR IP T

r-value

0.6 Yld91 Drucker nonAFR (n=1) Drucker nonAFR (n=2) Yld2004-18p Reduced Yld2004 Exp.

0.4 0.3 0.2

AN US

0.1 0.0 0

15

30

45

60

75

90

0

Loading angle [ ]

Figure 15 Assessment of various yield functions for AA2008-T4: (a) yield stresses of uniaxial stretching;

AC

CE

PT

ED

M

and (b) Lankford parameters

- 48 -

0.50

0.25

CR IP T

1.25 Yield surface Yld91 Drucker nonAFR (n=1) 1.00 Drucker nonAFR (n=2) Yld2004-18p Reduced Yld2004 0.75 Exp.

yy

ACCEPTED MANUSCRIPT

xx

-1.00

-0.75

-0.50

-0.25

0.25

0.50

0.75

1.00

1.25

AN US

-1.25

-0.25

-0.50

M

-0.75

-1.25

PT

ED

-1.00

AC

CE

Figure 16 Assessment of the yield envelopes predicted by various equations for AA2008-T4

- 49 -

1.25 Potential surface Yld91 Drucker nonAFR (n=1) Drucker nonAFR (n=2)

yy

ACCEPTED MANUSCRIPT

1.00

0.50

0.25

CR IP T

0.75

xx

-1.00

-0.75

-0.50

-0.25

0.25

0.50

0.75

1.00

1.25

AN US

-1.25

-0.25

-0.50

PT

ED

M

-0.75

-1.00

-1.25

AA2008-T4

AC

CE

Figure 17 Comparison of the potential loci predicted by three yield functions under non-AFR for

The calibrated anisotropic yield functions are then implemented into numerical analysis of deep

drawing of circular cup for AA2008-T4. The cut height is measured after the cup is fully drawn into the

die and plotted together with experimental results in Figure 18. It demonstrates that the Yld91 equation

with non-AFR and the anisotropic Drucker equation (n=1) with non-AFR predict the earing profile of the

- 50 -

ACCEPTED MANUSCRIPT

alloy with little difference. The difference of predicted earing profile is also hardly observed between the

Yld2004-18p and the anisotropic Drucker equation (n=2) with non-AFR. These observations above are

natural considering that the difference in the analytical comparison in Figures 15-17 is not apparent for

CR IP T

these two pairs of yield equations. The earing profile of AA2008-T4 is observed to be correctly predicted

the reduced Yld2004 function as well as the other four functions even though the largest difference in

M

AN US

45.0 44.5 44.0 43.5 43.0 42.5 42.0 41.5 41.0 40.5 40.0 39.5 39.0 30

60

90

PT

0

ED

Cup height [mm]

height is about 1.5 mm located at about 45º for all the yield functions.

120

150

Yld91 nonAFR Drucker (n=1) nonAFR Drucker (n=2) nonAFR Yld2004-18p reduced Yld2004 Exp. data 180

210

240

270

300

330

360

Angle from rolling direction [0]

AC

CE

Figure 18 Earing profiles of numerical prediction and testing measurements

In addition, the computation cost is compared in Table 17 for numerical simulations of circular cup

deep drawing of AA2008-T4 with these five yield functions. All the simulations are conducted with a

personal computer with Inter(R) Xeon(R) CPU E3-1225 v6 @ 3.3GHz and 16GB memory. Two CUPs are

used in the simulation with double precision. The comparison in Table 17 shows that the anisotropic

- 51 -

ACCEPTED MANUSCRIPT

Drucker functions under non-AFR with n=1 and n=2 take the shortest time among these five yield

functions since the calculation of derivatives is simple for stress-invariant-based functions under spatial

loading. It takes a little longer for the anisotropic Drucker equation under non-AFR with n=2 than that

CR IP T

with n=1 considering that the anisotropic Drucker equation and its derivatives are computed for two times

when n=2. The simulations with non-quadratic equations of Yld91 under non-AFR, Yld2004-18p and

AN US

reduced Yld2004 take longer time because the computation of the principal stresses and derivatives are

extremely complicated and time consuming for triaxial loading. The Yld91 function under non-AFR takes

shortest time among these three non-quadratic yield functions considering that the Yld91 function is

M

relatively simpler than the Yld2004-18p equation and the reduced Yld2004 equation. The reduced

ED

Yld2004 equation takes slightly shorter time than the Yld2004-18p yield function in the numerical

simulation considering that the six-parameter linear transformation tensors are comparatively easy

PT

compared to the nine-parameter linear transformation tensors in the Yld2004-18p function. It takes the

CE

longest time for the Yld2004-18p equation in simulation of circular cup deep drawing, but this is the only

AC

yield function under AFR which can predict six or more ears in circular cup deep drawing for strong

anisotropic metals, such as AA3104, AA5042-H2 and AA2090-T3.

- 52 -

ACCEPTED MANUSCRIPT

Table 17 Comparison of computation efficient of different yield functions

Yield function

Yld91 under Anisotropic Drucker non-AFR

Time (second) Time compared

Anisotropic Drucker

Yld2004-18p

(n=1) under non-AFR (n=2) under non-AFR

Reduced Yld2004

556

414

479

649

600

0.8567

0.6379

0.7180

1.0000

0.9245

CR IP T

to Yld2004-18p

In summary, the anisotropic Drucker function under non-AFR with n=1 has similar accuracy compared

AN US

with the Yld91 function under non-AFR, but can improve the computation efficiency by about 25% for

the circular cup deep drawing simulation. There is barely difference in term of accuracy between the

anisotropic Drucker equation with non-AFR (n=2) and the Yld2004 function for earing profile prediction,

M

but the computation cost is reduced by about 28% for the simulation using the anisotropic Drucker

ED

equation with non-AFR (n=2) compared to the Yld2004-18p equation. However, the implementation of

PT

the anisotropic Drucker equation with non-AFR (n=2) needs special experiences and cares considering

that non-AFR is coupled. At last, there is no improvement in computation cost for the reduced Yld2004

CE

function, but the reduced Yld2004 function is based on AFR. Accordingly, its implementation in

AC

numerical simulation is with ease. Besides, it reduces the calibration cost by about half compared with the

Yld2004-18p function.

- 53 -

ACCEPTED MANUSCRIPT

9. Conclusions The Yld2004-18p function is modified into a form with less anisotropic parameters to characterize

anisotropic plastic deformation of moderately anisotropic metals. The reduction of the Yld2004-18p

CR IP T

equation is realized by substituting two nine-parameter tensors with two six-parameter linear

transformation tensors in Barlat et al. [19]. The reduced Yld2004 equation is successfully employed to

AN US

illustrate yielding and anisotropic deformation of an AHSS of DP980 and two aluminum alloys: AA6022-

T4E32 and AA2090-T3. The reduced Yld2004 functions is also used in FEA of stretching of notched

specimens and deep drawing of circular cup for AA2008-T4 alloy. Both analytical and numerical

M

applications prove that the reduced Yld2004 function can describe orthotropic characteristics of materials

ED

under spatial loading with satisfactory accuracy, but the calibration cost is cut down by about half for the

reduced Yld2004 function compared with the Yld2004-18p function. The numerical application also

PT

shows that there is little difference in computation cost between the Yld2004-18p and reduced Yld2004

CE

equations, while the anisotropic Drucker function enhances computation efficiency of numerical

AC

simulation between 25% and 60% compared with the Yld2004-18p and reduced Yld2004 equations. All

the non-quadratic functions and anisotropic Drucker equations model the strength of plane strain tension

with better accuracy than the Hill48 quadratic function even though the Hill48 equation takes the shortest

time for simulation. Accordingly, the reduced Yld2004 equation is recommended for modeling of

anisotropic mechanical properties of metals with moderate anisotropy especially for triaxial loading

- 54 -

ACCEPTED MANUSCRIPT

considering its satisfactory accuracy, reduced cost in parameter calibration and easy implementation in

numerical analysis under AFR.

CR IP T

Acknowledgement

The authors acknowledge the ARC linkage project: LP150101027 and the financial support by the “Young Talent Support Plan” of Xi’an Jiaotong University under the grant number of JX6J006, for their

AN US

supports of this study. The authors would like to thank the Guangdong National Natural Science

AC

CE

PT

ED

M

Foundation under Grant Nos. 2016A030313453 and 2016A030313519, for their supports of this study.

- 55 -

ACCEPTED MANUSCRIPT

References [1] Pack K, Marcadet SJ. Numerical failure analysis of three-point bending on martensitic hat assembly using advanced plasticity and fracture models for complex loading. Int J Solids Struct 2016;85/86:144-59. [2] Lou Y, Huh H. Prediction of ductile fracture for advanced high strength steel with a new criterion:

CR IP T

experiments and simulation. J Mater Process Technol 2013;213:1284-302.

[3] Lou Y, Chen L, Clausmeyer T, Tekkaya AE, Yoon JW. Modeling of ductile fracture from shear to balanced biaxial tension for sheet metals. Int J Solids Struct 2017;112:169-84.

[4] Yin Q, Soyarslan C, Isik K, Tekkaya AE. A grooved in-plane torsion test for the investigation of shear fracture in sheet metals. Int J Solids Struct 2015;66:121-32.

AN US

[5] Li H, Fu MW, Lu J, Yang H. Ductile fracture: experiments and computations. Int J Plasticity 2011;27:147-80.

[6] Mu L, Wang Y, Zang Y, Stemler PMA. Edge fracture prediction using uncoupled ductile fracture models for DP780 sheet. J Fail Anal Prev 2017;17:321-9.

[7] Anderson D, Butcher C, Pathak N, Worswick MJ. Failure parameter identification and validation for a dual-phase 780 steel sheet. Int J Solids Struct 2017;124:89-107.

M

[8] Xu W, Wu H, Ma H, Shan D. Damage evolution and ductile fracture prediction during tube spinning of titanium alloy. Int J Mech Sci 2018;135:226-39.

ED

[9] Vershinin VV. Validation of metal plasticity and fracture models through numerical simulation of high velocity perforation. Int J Solids Struct 2015;67/68:127-38. [10] Xiao X, Pan H, Bai Y, Lou Y, Chen L. Application of the modified Mohr–Coulomb fracture

PT

criterion in predicting the ballistic resistance of 2024-T351 aluminum alloy plates impacted by blunt projectiles. Int J Impact Eng 2019;123:26-37.

CE

[11] Xiao X, Mu Z, Pan H, Lou. Effect of the Lode parameter in predicting shear cracking of 2024-T351 aluminum alloy Taylor rods. Int J Impact Eng 2018;120:185-201. [12] Zhuang X, Wang T, Zhu X, Zhao Z. Calibration and application of ductile fracture criterion under

AC

non-proportional loading condition. Eng Fract Mech 2016;165:39-56.

[13] Li H, Hu X, Yang H, Li L. Anisotropic and asymmetric yielding and its distorted evolution: modeling and applications. Int J Plasticity 2016;82:127-58.

[14] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc Roy Soc Lond A 1948;193:281-97. [15] 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 Plasticity 1989;5:51-66. [16] Barlat F, Brem JC, Yoon JW, Chung K, Dick RE, Lege DJ, Pourboghrat F, Choi S-H, Chu E. Plane stress yield function for aluminum alloy sheets – Part I: Theory. Int J Plasticity 2003;19:1297-319.

- 56 -

ACCEPTED MANUSCRIPT

[17] Banabic D, Aretz H, Comsa DS, Paraianu L. An improved analytical description of orthotropy in metallic sheets. Int J Plasticity 2005;21:493-512. [18] Lou Y, Huh H, Yoon JW. Consideration of strength differential effect in sheet metals with symmetric yield functions. Int J Mech Sci 2013;66:214-23. [19] Barlat F, Lege DJ, Brem JC. A six-component yield function for anisotropic materials. Int J Plasticity 1991;7:693-712.

CR IP T

[20] 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.

[21] Barlat F, Aretz H, Yoon JW, Karabin ME, Brem JC, Dick RE. Linear transformation-based anisotropic yield functions. Int J Plasticity 2005;21:1009-39.

[22] Karafillis AP, Boyce MC. A general anisotropic yield criterion using bounds bad a transformation

AN US

weighting tensor. J Mech Phys Solids 1993;41:1859-86.

[23] Cazacu O, Barlat F. Generalization of Drucker’s yield criterion to orthotropy. Math Mech Solids 2001;6:613-30.

[24] Cazacu O, Barlat F. A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals. Int J Plasticity 2004;20:2027-45.

Plasticity 2006;22:1171-94.

M

[25] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. Int J

[26] Yoon JW, Lou Y, Yoon JH, Glazoff MV. Asymmetric yield function based on the stress invariants

ED

for pressure sensitive metals. Int J Plasticity 2014;56:184-202. [27] Yoshida F, Hamasaki H, Uemori T. A user-friendly 3D yield function to describe anisotropy of steel

PT

sheet. Int J Plasticity 2013;45:119-39.

[28] Hu Q, Li X, Han X, Li H, Chen J. A normalized stress invariant-based yield criterion: modeling and validation. Int J Plasticity 2017;99:248-73.

CE

[29] Lou Y, Yoon JW. Anisotropic yield function based on stress invariants for BCC and FCC metals and its extension to ductile fracture criterion. Int J Plasticity 2018;101:125-55.

AC

[30] Aretz H, Barlat F. New convex yield functions for orthotropic metal plasticity. Int J Non Linear Mech 2013;51:97-111.

[31] Kuwabara T, Yoshida K, Narihara K, Takahashi S. Anisotropic plastic deformation of extruded aluminum alloy tube under axial forces and internal pressure. Int J Plasticity 2005;21:101-17.

[32] 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 Plasticity 2006;22:174-93. [33] Dick RE, Yoon JW. Plastic anisotropy and failure in thin metal: Material characterization and fracture prediction with an advanced constitutive model and polar EPS (effective plastic strain) fracture diagram for AA3014-H19. Int J Solids Struct 2018;151:195-213.

- 57 -

ACCEPTED MANUSCRIPT

[34] Zhang SJ, Lou Y, Yoon JW. Earing prediction of AA 2008-T4 with anisotropic Drucker yield function based on the second and third stress invariants. J Phys NUMISHEET’2018 2018;1063:12113. [35] Hosford WF. A general isotropic yield criterion. J Appl Mech 1972;39:607-09. [36] Drucker DC. Relations of experiments to mathematical theories of plasticity. J Appl Mech 1949;16:349-57.

identification of anisotropic sheets. Exp Mech 2015;55:1005-22.

CR IP T

[37] Deng N, Kuwabara T, Korkolis YP. Cruciform specimen design and verification for constitutive

[38] Kuwabara T, Mori T, Asano M, Hakoyama T, Barlat F. Material modeling of 6016-O and 6016-T4 aluminum alloy sheets and application to hole expansion forming simulation. Int J Plasticity 2017;93:164-86.

[39] Hershey AV. The plasticity of an isotropic aggregate of anisotropic face centred cubic crystals. J

AN US

Appl Mech 1954;21:241-49.

[40] Yoon JW, Dick RE, Barlat F. A new analytical theory for earing generated from anisotropic plasticity. Int J Plasticity 2011;27:1165-84.

[41] Yoon J-H, Cazacu O, Yoon JW, Dick RE. Earing prediction for strongly textured aluminum sheets. Int J Mech Sci 2010;52:1563-78.

M

[42] Kuwabara T, Hakoyama T, Maeda T, Sekiguchi C. Benchmark-1: hole expansion of a high strength steel sheet. NUMISHEET’2018, Tokyo, Japan 2018. [43] Kuwabara T, Ikeda S, Kuroda T. Measurement and analysis of differential work-hardening in cold-

ED

rolled steel sheet under biaxial tension. J Mater Process Technol 1998;80/81:517-23. [44] Stoughton TB, Yoon JW. Anisotropic hardening and non-associated flow in proportional loading of

PT

sheet metals. Int J Plasticity 2009;25:1777-817. [45] Yoon JW, Barlat F, Chung K, Pourboghrat F, Yang DY. Earing predictions based on asymmetric nonquadratic yield function. Int J Plasticity 2000;16:1075-104.

CE

[46] Lou YS, Yoon JW. Anisotropic ductile fracture criterion based on linear transformation. Int J Plasticity 2017;93:3-25.

AC

[47] Lou Y, Huh H. Extension of a shear controlled ductile fracture model considering the stress triaxiality and the Lode parameter. Int J Solids Struct 2013;50:447-55.

[48] Lou YS, Yoon JW, Huh H. Modeling of shear ductile fracture considering a changeable cut-off value for the stress triaxiality. Int J Plasticity 2014;54:56-80.

[49] Yoon JW, Song IS, Yang DY, Chung K, Barlat F. Finite element method for sheet forming based on an anisotropic strain-rate potential and the convected coordinate system. Int J Mech Sci 1995;37:733-52.

- 58 -