Numerical procedures for predicting localization in sheet metals using crystal plasticity

Computational Materials Science 72 (2013) 107–115

Contents lists available at SciVerse ScienceDirect

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

Numerical procedures for predicting localization in sheet metals using crystal plasticity Ji Hoon Kim a, Myoung-Gyu Lee b,⇑, Daeyong Kim a, Frédéric Barlat b a b

Materials Deformation Department, Korea Institute of Materials Science, 797 Changwondaero, Changwon, Gyeongnam 642-831, Republic of Korea Graduate Institute of Ferrous Technology, Pohang University of Science and Technology, San 31 Hyoja-dong, Nam-gu, Pohang, Gyeongbuk 790-784, Republic of Korea

a r t i c l e

i n f o

Article history: Received 14 October 2012 Received in revised form 18 January 2013 Accepted 6 February 2013

Keywords: Formability Forming limit Crystal plasticity Simplex method Direct search

a b s t r a c t In order to obtain reliable solutions for the initiation of localized instability in sheet metal forming, numerical procedures are developed for predicting forming limits using crystal plasticity. The Marciniak–Kuczynski (M–K) method is employed to calculate the initiation of the instability with rate-dependent elastic-crystal plasticity for cubic metals. The Nelder–Mead (N–M) direct search and Newton– Raphson (N–R) methods are used to solve the highly nonlinear equilibrium and associated compatibility equations under the M–K scheme. The N–M method eliminates the need of the first-order derivatives of the constitutive equations that may not be readily available for complex material models. The numerical results show that the calculation time using the N–M method is comparable to that using the N–R for the highly nonlinear M–K equations if the simplex size parameters are properly optimized. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction In the sheet metal forming, it is helpful to know the forming limits of a sheet metal so that sheet failures can be avoided by designing the forming process properly. The forming limit curve (FLC) of a sheet metal represents the limits of major and minor principal strains at the onset of localized necking. Various approaches have been developed for predicting forming limits [1– 5]. Among them, a theory proposed by Marciniak and Kuczynski (M–K) [1] has been widely used in which calculate the localized necking occurring at a locally imperfect band. The calculation procedure seeks the deformation and stress in the band and the outer homogeneous region for a given biaxial boundary condition imposed on the homogeneous region until necking occurs. The overall numerical scheme consists of both the integration of stresses by the constitutive relationship and the force equilibration under a given compatibility along the imperfection band. Conventionally phenomenological elastic–plastic constitutive models have been used for the forming limit calculations based on the M–K theory [6–8]. Recently, more attentions have been paid to the forming limit calculations using the crystal plasticity constitutive models [9–11] because the crystal plasticity can capture the evolution of anisotropy caused by texture evolution during plastic deformation on a more physical basis.

⇑ Corresponding author. Tel.: +82 54 279 9034; fax: +82 54 279 9299. E-mail address: [email protected] (M.-G. Lee). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.02.008

In the M–K forming limit calculation, a numerical procedure is inevitably involved in calculating equilibrium equations unless the material behavior is simple and analytical solutions are available. A popular choice for solving the equilibrium equations is the Newton–Raphson (N–R) method in which the equations are linearized and the solutions are iteratively updated [12–14]. These methods require the first-order derivatives of the constitutive equations, or the Jacobian matrix. However, it is not always easy to obtain Jacobian matrix in crystal plasticity constitutive models. Even for simple elasto-viscoplastic crystal plasticity models, the expressions for the Jacobian matrix are lengthy and the calculation procedure is cumbersome (see Appendix). For complex material models, the first-order derivatives may not be readily available in the analytical form. In addition, the derivative-based methods are known to have convergence issues when the equations are highly nonlinear. The Nelder–Mead (N–M) simplex algorithm [15] is one of the most widely used direct search methods for multi-dimensional unconstrained minimization. The N–M method attempts to minimize a scalar-valued nonlinear function of multiple real variables using only function evaluations without any derivative information. During iterations, a simplex – a convex hull of n + 1 vertices – evolves according to given rules until the function values satisfy specified criteria. This method is useful for problems where derivatives are too expensive to obtain or cannot be calculated in closed forms. In this work, numerical procedures are presented for predicting the sheet forming limits using the crystal plasticity based M–K method without the use of the material Jacobian matrix. The

108

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

details on the constitutive equations and the M–K theory utilized in this study are overviewed in the following sections. The numerical procedures for applying the N–R method and the N–M simplex method to the M–K method are presented. The accuracy and the robustness of the numerical procedures employing the N–R and N–M methods are compared. 2. Rate-dependent crystal plasticity A classical Taylor-type isothermal crystal plasticity constitutive equation based on dislocation glide is employed as originally proposed by Asaro [16] and Rice [17]. The total deformation gradient is expressed in terms of the deformation gradient Fp associated with dislocation glide on specific crystal slip planes and the deformation gradient Fe associated with elastic distortion and rigid body rotation. e p

F¼F F

ð1Þ

Gliding of dislocations on active slip plane a with slip plane norðaÞ ðaÞ mal n0 and slip direction s0 , induces a rate of change plastic part p of deformation gradient, F_

F_ p ¼ Lp Fp ;

Lp ¼

NS X

NS X

a¼1

a¼1

c_ ðaÞ sð0aÞ  n0ðaÞ ¼

c_ ðaÞ Sð0aÞ

ð2Þ

where  denotes a tensor product. The rate of shear on slip system a is expressed according to Peirce et al. [18]

 ðaÞ ð1=mÞ s   signðsðaÞ Þ g ðaÞ 

c_ ðaÞ ¼ c_ 0 

ð3Þ

where s(a) is the resolved shear stress onto slip plane a in slip direcðaÞ tion s0 ; c_ 0 is a reference shear rate, and m is the strain rate sensitivity. The slip resistance g(a) has an initial value g0 and evolves with plastic strain on all slip systems according to

g_ ðaÞ ¼

NS X hab jc_ ðbÞ j

ð4Þ

b¼1

where NS is the number of total slip systems. The hardening coefficients are given by

hab ¼ hb ðq þ ð1  qÞdab Þ;

 a g ðbÞ hb ¼ h0 1  gs

ð5Þ

where hab represents interactions of dislocations among different slip planes. Therefore, q is the ratio of latent-to-self hardening coefficients, h0 is a reference self-hardening coefficient, a is the strain hardening exponent, and gs is the slip resistance at which the self and latent hardening coefficients approach zero. Anisotropic elastic-viscoplasticity is used for the single crystal constitutive model. As for the stress/deformation measure, the

100

TD

RD

110

Fig. 1. Initial configuration of the imperfection.

second Piola–Kirchhoff stress (PK2) based on the intermediate configuration and its conjugate Lagrangian strain tensor are used [19], given respectively by,

P ¼ detðFe ÞFe1 rFeT

ð6Þ

and

Ee ¼ ðFeT Fe  IÞ=2

ð7Þ

where P is the intermediate configuration-based PK2 stress tensor, r is the Cauchy stress, and Ee is the elastic Lagrangian strain tensor. The constitutive relation is given by

P ¼ Ce : Ee

ð8Þ

e

where C is the fourth order tensor containing the anisotropic elastic moduli. Implicit integration schemes usually require the tangential response of the constitutive relationship. Differentiating Eq. (8) with respect to the rate of deformation tensor, D, gives the ‘‘conceptual’’ expression for the fourth-order Jacobian matrix, CT, as

r_ ¼ CT : D

ð9Þ T

The analytical expression of C is not always available. For the ratedependent crystal plasticity model used in this work, the analytical form of CT is given elsewhere [20], and a brief summary of the Jacobian matrix is given in Appendix. For application to cubic metals, 12 {1 1 1}h1 1 0i slip systems are considered for FCC structures while 12 {1 1 2}h1 1 1i and 12 {1 1 0}h1 1 1i are usually adopted for BCC [21]. In order to obtain the response of a polycrystal, the Taylor assumption was used [18]. A numerical scheme to update stress tensor is similar to that for solving nonlinear boundary value problems using the finite element method. In the stress integration procedure, the variables ðaÞ p in the previous time increment  n; ðrn ; Fn ; g n Þ, are known along ðaÞ ðaÞ with the slip systems, s0 ; n0 . The stress, plastic deformation gradient, slip resistance, and orientation of slip systems are then updated for the given deformation, Fn+1. The viscoplastic shear rate can be rewritten in a general form:

  aÞ ðaÞ DcðaÞ ¼ c_ ðaÞ Dt ¼ / sðnþ1 ; g nþ1

RD

TD

Fig. 2. Pole figures of the model steel sheet.

111

TD

ð10Þ

RD

109

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

A in Fig. 1) to the thickness of the groove region (region B in Fig. 1), B A f ¼ h0 =h0 , taking a value of 0.99 in this work. With no shear straining, the deformation in Region A is described in terms of the velocity gradient tensor LA as follows:

700

Fitted

True Stress (MPa)

600 500

0

1 0 DA 0 DA11 0 B C B 11 A C B ¼ L A ¼ DA ¼ B @ 0 D22 0 A @ 0 0 0 DA33 0

Target 400 300

RD Strain rate = 0.001/s

0 0.00

0.05

0.10

0.15

0.20

qDA11 0

1

0 A e_ C B 11 0 C A¼@ 0 0 DA33 0

0 qe_ A11 0

1 0 C 0 A e_ A33

ð12Þ

where DA is the rate of deformation tensor, q specifies the deformation paths, taking a value between 0.5 and 1.0, so that the 1-direction is taken as the major principal strain direction, and e_ is the strain rate tensor. The deformation gradient tensor, F, at the end of the current increment (n ? n + 1) is given by

200 100

0

0.25

FAnþ1 ¼ FAn þ LA DtFAnþ1 or FAnþ1 ¼ ðI  LA DtÞ1 FAn

True Strain Fig. 3. Measured and fitted uniaxial stress–strain curves of the steel sheet.

ðaÞ ðaÞ Since the resolved shear stress snþ1 and g nþ1 are functions of Dc(a) the above equation is a system of nonlinear equations in terms of Dc(a). Therefore, this set of equations can be solved by the N–R method. Detailed numerical procedure for rate-dependent crystal plasticity is well documented elsewhere [22,23]. One numerical problem frequently met during the solution procedure is the convergence issue. Eq. (3) becomes very stiff as strain rate exponent m approaches 0. A line search algorithm is used to avoid the convergence issue [24]. The strain increment is updated as

ð13Þ

where Dt is the time increment. In general, the total deformation does not automatically preserve the incompressibility. In such cases, the rate-of-deformation in the thickness direction, DA33 , should be calculated using the plane stress condition,

rA33 ¼ 0

ð14Þ

Eq. (14) can be solved for

DA;mþ1 ¼ DA;m 33 33 

r

DA33

by successive N–R corrections:

A 33

ð15Þ

C T3333 Dt

where the superscripts m and m + 1 are the iteration numbers. The thickness of the region A is updated as A

aÞ ðaÞ ða Þ Dcðkþ1 ¼ Dck þ b  dDckþ1

ð11Þ

A

A

A

A

hnþ1 ¼ hn þ DA33  Dt  hnþ1 or hnþ1 ¼

ðaÞ

where dDckþ1 is the strain increment for the k + 1th iteration obtained by the N–R method, b is a parameter dependent on a norm, Nf = jdDc(a)/Dc(a)j. After selecting b according to the fraction norm, the residual of linearized Eq. (8) is re-estimated and the same procedure is repeated with the b divided by a factor 2 until the residual decreases. When no b is found after seven divisions, the time step is automatically reduced to half of the current time step.

hn

1  DA33 Dt

ð16Þ

As Region A deforms, the angle of the groove rotates as follows [25]:

/ ¼ tan1 ðexpðð1  qÞe11 Þ tan /0 Þ

ð17Þ

In addition, the band region B deforms, satisfying the equilibrium and compatibility between the regions A and B [12,25]. The compatibility equations are given by

3. Forming limit prediction

LBab ¼ LAab þ c_ a nb ;

The numerical procedure using the Newton–Raphson method is first introduced for calculating the Marciniak–Kuczynski forming limits and the application of the Nelder–Mead simplex method is explained later.

where n is a unit vector perpendicular to the groove line, T n ¼ ð cos / sin / Þ , and c_ is an unknown vector. The force equilibrium equations are given by

3.1. The Newton–Raphson method

A

B

na rAab h  na rBab h ¼ 0 A

ð18Þ

ð19Þ

B

where h and h are the current thicknesses of Regions A and B, respectively. In addition, the plane stress condition is imposed:

The M–K theory assumes that a groove-type imperfection exists across the sheet with an angle /0 from the reference direction as shown in Fig. 1. The degree of imperfection, f, is represented by the ratio of the initial thickness of the homogenous region (region

Table 1 Single crystal parameters for the model steel. C11 C12 C44 c_ 0 m s0 h0 ss ahard q

ða; b ¼ 1; 2Þ

237 GPa 141 GPa 116 GPa 0.001/s 0.01 95.33 MPa 888.51 MPa 219.44 MPa 1.5772 1.4

rB33 ¼ 0

ð20Þ

To obtain the deformation in Region B, the system of equations, Eqs. (19) and (20), should be solved for the unknowns c_ 1 ; c_ 2 , and DB33 . The three equations can be expressed in the matrix form as 0 1 0   A1 0 1 0 1 2 3 rB11 B B 0 n1 rA11 þ n2 rA21 h f1 0 n2 h 0 B B C n1 h B r   AC B C 6 7B 22 C C C¼B A A 0A @ f2 A ¼ 4 0 n2 hB n1 hB 0 5B B C  B @ n r þ n r h @ A 1 2 12 22 @ r12 A 0 f3 0 0 0 1 0 B

r33

ð21Þ

or

f ¼ ArB þ b ¼ 0 where

ð22Þ

110

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

2

(a) 0.6



@ rB @DB

0.5

Major Strain

C T1112

C T2222

C T2212

C T1222

C T1212

0

0

0

3

7 0 7 7 7 0 5

ð26Þ

C T3333

and NewtonRaphson

ρ=-0.5 0.3

NelderMead

0.2

0.1

0.0 -0.2

-0.1

0.0

0.1

3

2 n1 " # B 60 @D 6 ¼ 6 n2 4 2 @x

0

0

n2 n1 2

07 7 7 05

0

0

1

ð27Þ

Here, rB and DB have only four non-zero components in the order of 11, 22, 12, and 33. The non-zero terms of Eq. (25) can be expressed as follows:

ρ=1.0

0.2

0.3

0.4

Minor Strain

(b) 30 Initial Angle of Imperfection (degree)

6 T 6 C 2211 ¼ Dt 6 6 T 4 C 1211

C T1122

0

0.4

h   @f1 n2  n2 i B ¼ h Dt n1 C T1111 n1 þ C T1112 þ n2 C T1211 n1 þ C T1212 @ c_ 1 2 2

ð28Þ

h   @f1 n1  n1 i B ¼ h Dt n1 C T1122 n2 þ C T1112 þ n2 C T1222 n2 þ C T1212 @ c_ 2 2 2

ð29Þ

h   @f2 n2  n2 i B ¼ h Dt n2 C T2211 n1 þ C T2212 þ n1 C T1211 n1 þ C T1212 @ c_ 1 2 2

ð30Þ

h   @f2 n1  n1 i B ¼ h Dt n2 C T2222 n2 þ C T2212 þ n1 C T1222 n2 þ C T1212 @ c_ 2 2 2

ð31Þ

Newton-Raphson

20

@f3

10

@DB33

Nelder-Mead

¼ C T3333 Dt

ð32Þ

The thickness of the band region is updated as B

B

B

B

B

h ¼ hold þ DB33 Dth or h ¼

0

-0.4

-0. 2

0.0

0.2

0.4

0.6

0.8

Fig. 4. Comparison of the predicted (a) forming limit strains and (b) critical imperfection angle obtained using the Newton–Raphson and the Nelder–Mead simplex method.

B

n1 h

0

n2 h

B B

6 A¼4 0

n2 h

n1 h

0

0

0

B

0

3

7 0 5 and b ¼  1

A

ðn1 rA11 þ n2 rA21 Þh

!

A

ðn1 rA12 þ n2 rA22 Þh 0 ð23Þ

Eq. (22) can be solved using the N–R method [6,7]. In the N–R method, the solution is iteratively corrected by

xmþ1 ¼ xm 



1 @f f @x

ð24Þ

 where x is the unknown vector, x ¼ c_ 1 matrix @f/@x is given by

c_ 2

DB33

T

. The Jacobian

@f @ rB @DB ¼A B @x @D @x

ð25Þ B

B

B

hold 1  DB33 Dt

ð33Þ

where the subscript ‘old’ denotes a value from the previous increment.

ρ

2



C T1111

@h @r Here, the derivative @D and @D are giB is ignored. The derivatives @x @DB 33 ven from Eqs. (9) and (18), respectively, by

3.2. The Nelder–Mead simplex method The N–R method is only applicable when the Jacobian matrix of the constitutive equations (CT in this work) is available. For more complex material models such as the crystal plasticity model incorporated with the dislocation pileup [26,27], the exact Jacobian matrix may not be readily available in an analytical form. It is possible to obtain the Jacobian matrix numerically but it is computationally expensive. The Nelder–Mead (N–M) simplex method was proposed for function minimization with n variables without seeking the gradient of function [15]. The N–M algorithm consists of four operations – reflection, expansion, contraction, and shrinkage, with the corresponding scalar constants d, v, c, and r, respectively. These four parameters should satisfy d > 0, v > 1, v > d, 0 < c < 1, and 0 < r < 1. The following set of values d = 1, v = 2, c = 0.5, and r = 0.5 are used in this work, as suggested in the Nelder–Mead paper. The iteration begins with an initial non-degenerate simplex with the n + 1 vertices (x1, . . . , xn+1). The iteration results in a new simplex for the next iteration. The iteration procedure for minimizing a function f is illustrated in Box 1. Note that the vertex indices are ordered such that f(x1) 6 f(x2) 6    6 f(xn+1).

111

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uX u n u ðf ðxi ÞÞ2 t i¼1

Iteration procedure for the Nelder–Mead simplex method.

n

1. Reflect Compute the reflection point xr as

 þ dðx   xnþ1 Þ xr ¼ x

ð1Þ

 is the centroid of the vertices except the worst point where x xn+1. If f(x1) 6 f(xr) < f(xn), replace xn+1 with xr (or accept xr) and stop the iteration. 2. Expand If f(xr) < f(x1), compute the expansion point xe as

 þ vðxr  x Þ xe ¼ x

ð2Þ

If f(xe) < f(xr), accept xe and stop the iteration. Otherwise, accept the reflected value and stop the iteration. 3. Contract 3.a. If f(xn) 6 f(xr) < f(xn+1), compute the contraction point xc as

 þ cðxr  x Þ xc ¼ x

ð3Þ

If f(xc) 6 f(xr), accept xc and stop the iteration. Otherwise, move to the Shrink operation. 3.b. If f(xr) P f(xn+1), compute the contraction point xcc as

  cðx   xnþ1 Þ xcc ¼ x

ð4Þ

If f(xcc) < f(xn+1), accept xcc and stop the iteration. Otherwise, move to the Shrink operation. 4. Shrink Replace x2, . . . ,xn+1 with x0 i = x1 + r (xi  x1) ,i = 2, . . . ,n + 1.

To find roots of a non-negative function f, the iteration continues until the following convergence criterion is met: 100

RD

TD

110

< Tol

ð34Þ

where Tol is a tolerance (taking a value of 0.01 in this work). The N–M simplex method can be used to solve the equilibrium equations (Eqs. (14) and (22)) without the need for the Jacobian matrix. The N–M solution procedure for the M–K method is given as follows:

3.2.1. Calculation in Region A (Homogeneous region) For a new time increment ðtnþ1 ¼ tn þ Dt; DeA11 ¼ e_ A11 DtÞ, the velocity gradient is given by 0

1 0 DA 0 DA11 0 B C B 11 A;i A B C B L ¼ @ 0 D22 0 A ¼ @ 0 A;i 0 0 D33 0

0

qDA11 0

1

0 A e_ C B 11 0 C A¼@ 0 0 DA;i 33 0

0 qe_ A11 0

1 0 0 C A ði ¼ 1; 2Þ _eA;i 33 ð35Þ

where the superscript i denotes the ith vertex, Eq. (12). The plane stress constraint, Eq. (14), should be solved for the unknown DA33 . The following constraint scalar function is minimized using the N–M simplex method.

g A ¼ jr33 j

ð36Þ

The numbers of dimension and vertices in the N–M simplex method are 1 and 2, respectively. The values ð1 þ qÞDA11 =2 and 0 are used for the two initial simplex vertices. The deformation gradient tensor in Region A is obtained with the given velocity gradient during iterations, and the constitutive routine is sought by integrating stress. The simplex evolves according to the rules in Box 1 with the constraint function gA. Once converged, the thickness is calculated using Eq. (16). RD

TD

111

RD

TD

(a) 100

TD

RD

110

RD

TD

111

RD

TD

(b) Fig. 5. Pole figures at localization obtained using the Newton–Raphson method for the case with q = 0.8 and h = 0°, in Regions (a) A and (b) B.

112

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

RD

100

RD

110

TD

TD

RD

111

TD

(a) RD

100

RD

110

TD

TD

RD

111

TD

(b) Fig. 6. Pole figures at localization obtained using the Nelder–Mead method for the case with q = 0.8 and h = 0°, in Regions (a) A and (b) B.

6

0

c_ j n1 B 1j B c_ n1 @ 2

LB;j ¼ LA þ

c_ j1 n2 c_ j2 n2

0

0

0

Principal Strain Rate Ratio of B to A

3.2.2. Calculation in Region B (imperfection band) The unknowns of Region B, or a vertex are x ¼ ðc_ 1 ; c_ 2 ; DB33 Þ. Using the compatibility equation, Eq. (18), the velocity gradient is assigned for each vertex within the N–M simplex scheme:

1

C 0 C A ðj ¼ 1; 2; 3; 4Þ B;j D33

ð37Þ

where the superscript j denotes the jth vertex. For the first iteration of the first time increment when no previous solution is available, the vertices of the initial simplex are set as

a1 ;

a2 ;

a3 and ð0; 0; 0Þ:

ð38Þ

where a1 ¼ ðe_ 1 ; 0; 0Þ; a2 ¼ ð0; e_ 1 ; 0Þ, and a3 ¼ ð0; 0; e_ 1 Þ. For the first iteration of the subsequent time increments, the following set of vertices is used for the initial simplex:

 prev þ pa1 ; x

 prev þ pa2 ; x

prev þ pa3 ; and x prev x

variable : p ¼ p0

ð40Þ

e_ BM e_ A11

ð41Þ

where p0 is a constant initial simplex size parameter and e_ BM is the major principal strain rate in Region B. The variable case is to set p proportional to the principal strain rate ratios such that a larger simplex is used when Region B deforms faster than Region A. The effect of the choice of p on the convergence rate is analyzed in Section 4.

5

4

3

NewtonRaphson

2

ρ = 0.2, θ = 5 Δε = 0.001

1

ð39Þ

0.00

0.05

0.10

0.15

0.20

o

0.25

Major Strain

where xprev is the centroid of the vertices except the worst point in the previous time increment and p is a simplex size parameter. Two approaches have been tried for the parameter p: fixed and variable p:

fixed : p ¼ p0

Nelder-Mead

Fig. 7. Principal strain ratio of the regions B–A.

The deformation gradient in Region B is obtained with the given velocity gradient during iterations and the stress integration procedure is sought to update the stress tensor and thickness of the groove region. Then, a constraint scalar function is defined considering the equilibrium equation and the plane stress condition, Eqs. (19) and (20), as follows gB ¼

n



 A   B o2 n1 rA11 þ n2 rA21 h  n1 rB11 þ n2 rB21 h

1 n  B   B o2  B 2 2 þ r33 þ n1 rA12 þ n2 rA22 h  n1 rB21 þ n2 rB22 h

ð42Þ

113

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

(a)

Nelder-Mead

Cumulative # of Material Routine Calls

# of Material Routine Calls / Strain Increment 0.001

1400

Newton-Raphson ρ = 0.2, θ = 5 Δε = 0.001

1200

o

Region B

1000

800

600

400

Region A

200

Newton-Raphson 10000

ρ = 0.2, θ = 5o Δε = 0.001 0.00

0.00

0.05

0.10

0.15

0.20

0.25

Major Strain

1400

Nelder-Mead ρ = 0.2, θ = 5o Δε = 0.001

1000

800

600

Region B 400

Region A

200

0

0.00

0.05

0.10

0.15

0.20

0.25

Major Strain Fig. 8. Number of material routine calls for each strain increment 0.001 for Region A and B obtained using the (a) Newton–Raphson and (b) Nelder–Mead methods for the case with q = 0.2 and h = 5°.

Once the solutions are obtained with a proper tolerance, the forming limit criterion is checked. The criterion for the occurrence of necking is given as follows:

e_ BM P ae_ A11

0.05

0.10

0.15

0.20

0.25

Major Strain Fig. 9. Cumulative number of material routine calls.

(b) # of Material Routine Calls / Strain Increment 0.001

20000

0

0

1200

30000

ð43Þ

where a is a parameter taking a value of 5 in this work. The procedures are repeated until the necking criterion is met. 4. Numerical results The developed numerical procedures are tested for a model BCC material. The pole figures of the undeformed sample are shown in Fig. 2. The Euler angles of 300 grains are sampled and used for the initial texture for the crystal plasticity calculations for both Regions A and B. The single crystal model parameters (s0, h0, ss, a) were obtained

by fitting the predicted uniaxial stress–strain response to the target one in the rolling direction, as shown in Fig. 3. The single crystal model parameters for the model material are listed in Table 1. The forming limit curves are predicted following the M–K procedures using the N–R and N–M simplex methods, as shown in Fig. 4a. The strain ratio in the region A, q, increased from 0.5 to 0.9 by 0.1, for a constant principal strain rate e_ A11 ¼ 0:001=s in Region A. For each q, the initial imperfection angles, /0, from 0° to 90° by 5° were tried, and the lowest critical strains in Region A were chosen as the forming limit strains. The results of the N–R and N–M methods compared well with each other (Fig. 4b). The forming limits increase as q approaches 1.0 and similar results have been reported for BCC sheets [28]. This increase of forming limit is related to texture evolution. Crystallographic texture becomes axisymmetric and leads to a very flat yield surface in the balanced biaxial state. As a result, the strain increment is very stable even with stress fluctuations around this state. Thus, the evolution of the strain increment towards plane strain, the plastic flow localization mode, is inhibited, which explains the high forming limit for this strain path. Additional details on this issue are the object of a subsequent article. The pole figures at localization in Regions A and B for the case with q = 0.8 and h = 0° are compared in Figs. 5 and 6. It is shown that the texture has changed from that of the as-received material (Fig. 2) as the sheet has been stretched considerably (a major strain of 0.38) but the texture of Region B is similar to that of Region A. The texture prediction using the Nelder–Mead method agreed with the prediction using the Newton–Raphson method. The evolution of the principal strain ratios of Region B–A are compared for the case with q = 0.2 and h = 5° in Fig. 7. The curves obtained using the N–R and N–M methods compared well with each other, with a slight difference in the necking strain when the principal strain ratio reaches the critical value, 5 in this work. The numbers of material routine calls in Regions A and B for each strain increment of 0.001 during the forming limit calculation are compared for the case where q = 0.2 and h = 5° in Fig. 8. The number of material routine calls is used instead of the usual iteration numbers because the material routine takes most of the calculation time as it involves nonlinear equations for the 24 slip systems in 300 grains. In addition, the number of material routine

114

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

80000

(a) 80000

o

Total # of Material Routine Calls

ρ = 0.2, θ = 5

Total # of Material Routine Calls

70000

60000

Fixed 50000

40000

Variable

Nelder-Mead 60000

40000

20000

Newton-Raphson

30000 o

ρ = 0.2, θ = 5 Δε = 0.001 20000 0.1

1

10

100

Initial Simplex Size Parameter, p0

(b)

0.19

Major Limit Strain

1e-4

1e-3

1e-2

1e-1

Initial Strain Increment Size, Δε11A Fig. 11. Effect of the size of the initial strain increment on the total material routine calls.

0.20

0.18

0.17

0 1e-5

Fixed

Variable

0.16

ρ = 0.2, θ = 5 Δ ε = 0.001

o

0.15 0.1

1

10

100

Initial Simplex Size Parameter, p0 Fig. 10. Effect of the initial simplex size parameter on (a) the total number of material routine calls and (b) the major limit strains.

calls for each increment is different depending on the equation solving method: for the N–R method, the number of material routine calls equals the number of iterations, whereas the N–M method calls the material routine one to five times more than the N–R method does for each increment, depending on the active simplex operations. In both the N–R and N–M methods, the calculation in Region A takes smaller number of material routine calls than Region B since the number of unknown is only one in Region A. In the Newton–Raphson method (Fig. 8a), the number of material routine calls in Region B is large in the early stage (e1 < 0.01) where the material behavior evolves from the elasticity-dominant to the plasticity-dominant behavior and therefore the tangential modulus changes rapidly. The number of material routine calls then decreases and stays small in the middle stage (0.01 < e1 < 0.16) where plastic deformation is dominant. Finally the number increases rapidly as the principal strain ratio of Region B–A increases and the sheet approaches the necking point. This

behavior shows that the number of material routine calls or the iteration number is closely related to the nonlinearity of the constitutive and equilibrium equations. In the Nelder–Mead method (Fig. 8b), on the other hand, the number of material routine calls in Region B stays at a comparatively constant value from the early to the final stage of the deformation. It implies that the Nelder–Mead method is relatively insensitive to the tangential modulus change and can accommodate large change in the solution. The cumulative numbers of material routine calls for the N–R and N–M methods are compared in Fig. 9. As can be expected from Fig. 8, the number of material routine calls of the N–R method shows an S-shaped curve, whereas that of the N–M method is linear. Interestingly, the total numbers of material routine calls until the necking are similar to each other, resulting in a similar calculation time for the two methods. The effect of the simplex size parameter, p, (Eq. (43)) on the total number of material routine calls are compared for the fixed (Eq. (44)) and variable (Eq. (45)) simplex size parameter in Fig. 10a. In general, too large initial simplex needs many shrink operations, whereas too small initial simplex requires reflection and/or expansion operations. Therefore, there can be an optimum size for the initial simplex size. Overall, the total number of material routine calls is minimum at p0 = 1 for the variable p case, whereas, p0 = 5 resulted in the least number of material routine calls among the fixed p cases. This result is reasonable in that the larger simplex is required as the deformation rate of Region B increases compared to that of Region A. The simplex size parameter, however, does not affect the necking or limit strains as shown in Fig. 10b. The effect of the initial strain increment size, DeA11 , on the total number of material routine calls are compared for the N–R and N–M methods in Fig. 11. The total number of material routine calls increases rapidly as the initial strain increment size decrease for the N–M method, whereas the number stays relatively constant for the N–R method. The large initial strain increment size may be advantageous for the N–M method, but it may not always result in accurate results as the equations are highly nonlinear.

J.H. Kim et al. / Computational Materials Science 72 (2013) 107–115

jijkl

5. Conclusions

X ðaÞ ðaÞ ¼ Iijkl þ C ij Bkl

115

ðA6Þ

a

In this work, numerical procedures are presented for predicting the sheet forming limits using the Marciniak–Kuczynski (M–K) method with crystal plasticity. The numerical procedures for applying the Newton–Raphson (N–R) method and the Nelder– Mead (N–M) simplex method to the crystal plasticity-based M–K method are developed. From the numerical results, the following conclusions are drawn:  Both the N–R and N–M methods can be used to solve the crystal plasticity-based M–K equations and result in similar results (forming limits and texture evolution).  The number of material routine calls changes considerably with the N–R method depending on the nonlinearity of the constitutive and equilibrium equations during the calculation procedure whereas it is relatively insensitive to this factor with the N–M method.  Overall, the total number of material routine calls for the N–M method is comparable to that of the N–R method, even though the N–M method does not use the tangent modulus in the constitutive equations.  The initial simplex and strain increment sizes can be optimized to speed up the calculation in the N–M method.

Acknowledgements This research was performed under programs on Fusion Core Technology for Industry (No. 10040078) funded by Ministry of Knowledge Economy, Korea, the ERC program (NRF2012R1A5A1048294), and Fundamental Research Program of the Korea Institute of Materials Science (KIMS). Appendix A The derivation of the analytical form of the Jacobian matrix of the elasto-viscoplastic crystal plasticity model presented in [20] is briefly summarized here. The Cauchy stress tensor is given by

r¼

1 ½Fe PFeT  det Fe

ðA1Þ

By differentiating Eq. (A1), the Jacobian matrix may be expressed as C Tijkl ¼

 i @ rij 1 h e eT e e eT e1 ¼ Simkl Pmn F eT nj þ F im Hmnkl F nj þ F im Pmn Sjnkl  F im Pmn F nj Spqkl F qp @ ekl detFe

ðA2Þ

e

where e is the logarithmic strain tensor, S ¼ @F , and H ¼ @P . The @e @e equations for calculating Eq. (A2) are as follows:

Sijkl ¼ Rik F eij  Rik F eip

X X ðaÞ ðaÞ aÞ DcðaÞ Sð0;pj  Rim U mn F enp Rkl S0;pj a

ðaÞ

ðaÞ

Rij ¼ Bkl Hklij 1 ijmn Smnkl

Hijkl ¼ j

X aÞ  DcðaÞ Jðmnkl a

ðA3Þ

a

ðA4Þ ðA5Þ

BðaÞ ¼ ðaÞ

J ijkl ¼

 @ DcðaÞ ðaÞ 1  ðaÞ ðaÞT s S þ S0 2 0 @

1 ða Þ C ijkl Cmnkl 2

aÞ ða Þ ðaÞT Cðmnkl ¼ Lmpkl S0;pn þ S0;mp Lpnkl

ðA7Þ ðA8Þ ðA9Þ

1 C ijmn Lmnkl 2

ðA10Þ

C ijkl ¼ Q im Q jn Q ko Q lp C emnop

ðA11Þ

e eT e Lijkl ¼ F eT ik U lm F mj þ F im U mk F lj

ðA12Þ

Sijkl ¼

References [1] Z. Marciniak, K. Kuczynski, Int. J. Mech. Sci. 9 (1967) 609–620. [2] J.L. Bassani, J.W. Hutchinson, K.W. Neale, On the prediction of necking in anisotropic sheets, in: H. Lippman (Ed.), Metal Forming Plasticity, Springer, Berlin, 1979, pp. 1–13. [3] F. Barlat, Mater. Sci. Technol. 91 (1987) 55–72. [4] L.S. Toth, D. Dudzinski, A. Molinari, Int. J. Mech. Sci. 38 (1995) 805–824. [5] N. Boudeau, J.-C. Gelin, S. Sahli, Comp. Mat. Sci. 11 (1998) 45–64. [6] L. Paraianu, D.S. Comsa, D. Banabic, Advan. Sci. Lett. 19 (2013) 1011–1015. [7] A. Assempour, R. Hashemi, K. Abrinia, M. Ganjiani, E. Masoumi, Comput. Mater. Sci. 45 (2009) 195–204. [8] M. Nurcheshmeh, D.E. Green, Int. J. Mech. Sci. 53 (2011) 145–153. [9] C.J. Neil, S.R. Agnew, Int. J. Plasticity 25 (2009) 379–398. [10] M. Yang, X. Dong, R. Zhou, J. Cao, Mat. Sci. Eng. A 527 (2010) 6607– 6613. [11] H. Wang, P.D. Wu, K.P. Boyle, K.W. Neale, Int. J. Solids Struct. 48 (2011) 1000– 1010. [12] P.D. Wu, K.W. Neale, E. Van der Giessen, Proc. R. Soc. Lond. A 453 (1997) 1848– 1931. [13] R. Knockaert, Y. Chastel, E. Massoni, Int. J. Plasticity 18 (2002) 231–247. [14] J.W. Signorelli, M.A. Bertinetti, P.A. Turner, Int. J. Plasticity 25 (2009) 1–25. [15] J.A. Nelder, R. Mead, Comput. J. 7 (1965) 308–313. [16] R.J. Asaro, Acta Met. 27 (1979) 445–453. [17] J.R. Rice, J. Mech. Phys. Solids 19 (1971) 433–455. [18] D. Peirce, R.J. Asaro, A. Needleman, Acta Metall. 31 (1983) 1951–1976. [19] L. Anand, Int. J. Plasticity 1 (1985) 213–231. [20] S. Balasubramanian, Polycrystalline Plasticity: Application to Deformation Processing of Lightweight Metals, Ph.D. Thesis, Massachusetts Institute of Technology, 1998. [21] M.G. Lee, Wang J, Anderson PM, Mater. Sci. Eng. A 463 (2007) 263–270. [22] S.R. Kalidindi, C.A. Bronkshorst, L.J. Anand, J. Mech. Phys. Solids 40 (1992) 537– 569. [23] M. Kothari, L. Anand, J. Mech. Phys. Solids 46 (1998) 51–83. [24] M.G. Lee, Int. J. Mater. Product Technol. 33 (2008) 375–661. [25] J.W. Hutchinson, K.W. Neale, Sheet necking-II. Time-independent behavior, in: K.P. Koistinen, N.-M. Wang (Eds.), Mechanics of Sheet Metal Forming, Plenum Press, New York, NY, 1978, pp. 269–285. [26] M.G. Lee, H.J. Lim, B.L. Adams, J.P. Hirth, R.H. Wagoner, Int. J. Plasticity 26 (2010) 925–938. [27] H. Lim, M.G. Lee, J.H. Kim, B.L. Adams, R.H. Wagoner, Int. J. Plasticity 27 (2011) 1328–1354. [28] K. Inal, K.W. Neale, A. Aboutajeddine, Forming limit comparisons for FCC and BCC sheets, Int. J. Plasticity 21 (2005) 1255–1266.