138
456

 

THESIS

l

. ,. .L‘
a j\,

(2/7/ff‘

LIBRARIES
MICHIGAN STATE UNIVERSITY
EAST LANSING, MICH 48824-1048

This is to certify that the
thesis entitled

FORMING LIMIT DIAGRAMS

presented by

JIMMY ISSA

has been accepted towards fulﬁllment
of the requirements for the

MS. degree in MECHANICAL ENGINEERING

 

dQ Wk

Major Professor’s Signature

12/8/0%

Date

MSU is an Aﬁirmative Action/Equal Opportunity Institution

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/CIRC/DateDue.p65-p.15

FORMING LIMIT DIAGRAMS

BY

Jimmy Issa

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

2004

ABSTRACT

FORMING LIMIT DIAGRAMS

BY

Jimmy Issa

In this work, the Marciniak and Kuczynski model (MK) will be
reviewed and will be generalized to a more general approach for the
calculation of forming limit diagrams. In the new model we will
account for all directions of the imperfection. This new approach has
been used to calculate the theoretical FLDs (forming limit diagrams)
for Vonmises, Barlat Yield96 and Barlat YieleOOO for aluminum alloy
Al2008-T4 and compared with an experimental one. All the derivations
needed for the implementation of the code are derived in the appendix
for all yield functions. Finally the stretch stamping test will be modeled
using LS-DYNA. A very good agreement has been achieved between

the experimental and theoretical results.

To my family

iii

 

 

 

‘ sém-dﬁdldnmnrtu o! ....L..... DDT). .Ih ”a, W}

ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. Farhang Pourboghrat, who
supported me throughout my study, Dr. Cloud and Dr. Liu for being on
my committee. I would also like to thank Nader Abedrabo and the
faculty members of the department of Mechanical Engineering at
Michigan State University who taught me through my masters or
helped me with my research. Finally thanks to the mechanical
engineering department at Michigan state university for giving me

admission and ﬁnancial support.

iv

TABLE OF CONTENTS

LIST OF TABLES ......................................................................... vii
LIST OF FIGURES ...................................................................... viii
LIST OF ABBREVIATIONS ............................................................. x
INTRODUCTION .......................................................................... 1
CHAPTER 1: Experimental forming limits ........................................ 3
CHAPTER 2: MK model ................................................................. 6
2.1 Introduction .................................................................... 6
2.2 Right hand side of forming limit diagrams ............................ 7
2.2.1 Algorithm ................................................................... 8
2.2.2 Vonmises ................................................................. 11

2. 2. 3 Hill’ 5 quadratic 1948 .................................................. 14

2. 2. 4 Hill’ 5 non- quadratic 1979 ........................................... 16

2. 3 Left hand side of forming limit diagrams ............................ 19
2.3.1 Algorithm ................................................................. 20

2. 3. 2 Thichness limit ......................................................... 21
2.3.3 Results .................................................................... 23
CHAPTER 3: General MK model .................................................... 25
3.1 Introduction ................................................................... 25
3.2 Modiﬁed MK model ......................................................... 27
3.2.1 Algorithm ................................................................. 29
3.2.1.1 Region a .............................................................. 29
3.2.1.2 Region b .............................................................. 32

3.3 Yield functions ................................................................ 37
3.4 Hardening rules .............................................................. 37
3.5 Material constants .......................................................... 37
3.6 Results and discussions ................................................... 39
CHAPTER 4: Simulation of Stretch Stamping ................................. 40
4.1 Introduction ................................................................... 40
4.2 Description of the model ................................................. 41
4.2.1 Part dimensions ........................................................ 41

4.2.2 Element type and material model ................................ 41
4.2.3 Material properties .................................................... 42
4.2.4 Constrained and prescribed motion ............................. 44
4.3 LS-DYNA ....................................................................... 46
4.3.1 Swift and Voce FLDs .................................................. 46
4.3.2 DYNA FLDs ............................................................... 47
4.4 Experimental findings ...................................................... 49
4.5 Numerical results ........................................................... 50
CHAPTER 5: Conclusions ............................................................. 53
Appendix: Numerical implementations .......................................... 55
A Vonmises ...................................................................... 55
A1 Yield function-2d ....................................................... 55
A2 First derivative ......................................................... 55
A3 Second derivative ..................................................... 55
B Barlat Y|d96 ................................................................... 56
Bl Yield function-2d ....................................................... 56
32 First derivative ......................................................... 58
83 Second derivative ..................................................... 59
C Barlat Yld2000 ............................................................... 63
C1 Yield function-2d ....................................................... 63
C2 First derivative ......................................................... 64
C3 Second derivative ..................................................... 66

vi

Table 1 :
Table 2 :
Table 3 :
Table 4 :
Table 5 :
Table 6 :
Table 7 :
Table 8 :
Table 9 :

LIST OF TABLES

Constants for voce and swift - Al2008T4 ......................... 38
Constants for Barlat Yield96 - Al2008T4 ......................... 38
Constants for Barlat Yield2000 — Al2008T4 ...................... 38
Properties of Steel ........................................................ 42
Properties of AA3003-H111 ........................................... 42
Constants for voce equation .......................................... 43
Constants for swift equation .......................................... 43
Constants for Barlat Yield2000 ....................................... 43
FLD user ..................................................................... 47

vii

Figure 1 :
Figure 2 :
Figure 3 :
Figure 4 :
Figure 5 :
Figure 6 :
Figure 7 :
Figure 8 :
Figure 9 :

Figure 10 :

Figure 11
Figure 12
Figure 13

Figure 14 :

Figure 15

Figure 16 :
Figure 17 :
Figure 18 :

Figure 19

Figure 20 :

Figure 21

LIST OF FIGURES

Fractured specimen ....................................................... 4
Failure under biaxial tension ........................................... 5
RHS model ................................................................... 7
Different n values ....................................................... 12
Different m values ...................................................... 12
Different f values ........................................................ 13
Different R values ....................................................... 15
Different R values ....................................................... 17
Different a values ....................................................... 18

LHS model ............................................................... 19
: Thickness limit Different n values ................................ 21
: Thickness limit Different m values ............................... 22
: Thickness limit Different f values ................................. 22

LHS Different n values ............................................... 23
: LHS Different m values .............................................. 24

LHS Different f values ................................................ 24

Modiﬁed MK model .................................................... 27

Flds Al2008-T4-Swift ................................................. 38
: Flds Al2008-T4-Voce .................................................. 39

Model stretch stamping .............................................. 44
: Blank holder force ..................................................... 45

viii

Figure 22 :
Figure 23 :
Figure 24 :
Figure 25 :
Figure 26 :
Figure 27 :
Figure 28 :
Figure 29 :
Figure 30 :

Punch prescribed motion ............................................ 45
FLD AA3003-H111 ..................................................... 46
LS-DYNA FLDs .......................................................... 48
Failed sheet .............................................................. 49
FLD AA3003-H111 (t: 15 ms) ..................................... 50
FLD AA3003-H111 (t: 23 ms) ..................................... 51
Middle Element displacement ...................................... 51
FLD contour at failure ................................................ 52
Experimental and numerical ﬁndings ........................... 54

ix

LIST OF ABBREVIATIONS

a— Cauchy stress

3,3” - Yield stress from yield function

5”,. - Yield stress from hardening rule

5,? - Cauchy and effective strain

e- Engineering Strain

n,k,£0 - Swift constant

A,B,C- Voce constants

f - Imperfection parameter

p- Strain ratio

a - Stress ratio

cb- Yield ratio

t - Thickness

R — Ratio of plastic strain in the width and thickness direction
it” - Average of Rvalues

a - Exponent to account for crystallographically textured fcc and bcc metal
J - Jacobien

0- Imperfection orientation

C,,C2,C3,C6‘ax,ay,azo,azl- Barlat yield96 material constants

a,,a2,a3,a4,a5,a6,a7,a8- Barlat yie|d2000 material constants

Introduction

Forming limit diagrams (FLDs) have been valuable tools for analyzing
sheet metal forming after the introduction of the idea in 1965. Stretching
operations of sheet metals will normally lead to a development of a sharp
localized neck on the surface where the plastic ﬂow will localize and lead to
failure. Forming limit diagrams have been ﬁrst introduced by Keeler [1-2]
and are still used as a failure criterion to assess the formability of sheet
metals. Many tried to model the localized necking analytically after the
introduction of the idea, starting with the pioneer work of Marciniak and
Kuczynski [3] who first presented a mathematical model to help understand
sheet metal deformation during forming. They modeled the imperfection as a
long groove, perpendicular to the direction of the largest stress for positive
major strains (RHS) or along Hill’s direction of zero-extension (the strain
increment parallel to the groove inside and outside of it is equal to zero) for
the negative major strains (LHS). Using this approach theoretical FLDs have
been calculated for linear strain path using different yield criteria, von Mises,
Hill’s 1948 [4-5-6], Hill’s 1979 [5-6] and many other yield functions. Graf
and Hosford [6] used an anisotropic yield function (Hill's 1979) and
investigated the effect of the hardening exponent n, the strain rate sensitivity
m, the imperfection f, the R-value and the yield function exponent a. F.
Barlat [7] calculated the FLD using a yield function based on von Mises but

with a high exponent. The effect of nonlinear strain paths has been

investigated on forming limits where experimental and numerical FLD have
been predicted for changing strain paths for different materials [8-9-10-11].
M.C. Butuc [12-13] used a more general MK model and calculated the FLDs

for von Mises, Hill’s 1948, Hill’s 1979 and Barlat Yld96 for Al6061.

In this work the MK model and a more general approach for the calculation of
forming limit diagrams are explained in detail. The general approach is based
on the MK model but is more general since there is no constraint on the
groove direction and the limiting strain is calculated for all possible angels
and minimized versus the angel of the groove. This new approach uses the
Newton-Raphson method to solve for the transverse stress and effective
strain in the groove. For this reason the second derivative of the yield
function is needed for the evaluation of the Jacobian. I calculated the forming
limit diagram for Al2008-T4 aluminum alloy using the modiﬁed MK model
where von Mises, Yield96 and Yield2000 were used to describe the yield locus
and modeled the hardening using both Voce and Swift models. The ﬁrst and
second derivatives for the three yield functions are evaluated in the appendix
to make future implementation easier. F. Barlat [14] calculated the FLD by
incorporating microstructural parameters, internal damage (voids) and
crystallographic textures in his model and it has given very accurate results.
Using a macroscopic approach, a very good agreement has been achieved
between the calculated and the experimental FLD when the yield locus is

described with Yld96 and the material hardening with Swift model.

CHAPTER 1: Experimental Forming Limits

Over the years experimental forming limit diagrams have been
determined [4-15] for different linear strain paths and even for material pre-
strained in different directions (biaxial tension, uniaxial tension and plane

strain).

The most widely used technique involves printing or etching a grid of
small circles of diameter do on the sheet metal before deforming it. After
stretching or forming the sheet, the circles will deform into a shape which
resembles an ellipse. The forming process will stop when the ﬁrst perceptible
neck is observed. To ﬁnd the limiting strains for different strain paths,
lubrication and specimen widths are varied. To obtain the limit for negative
minor strains, very narrow specimens are used and by increasing the width
to full width specimen, we can ﬁnd the limits up to balanced biaxial. After
deformation the ellipses wholly or partially in the neck are considered failed,
while the circles one or more diameters away form the neck are considered
safe. The principal diameters of the failed ellipses are measured (major d1
and minor d2). The limiting strain can then be represented as engineering or

true strains using the two formulas written below.

True Limiting Strains Engineering Limiting Strains

do d0

d2 ‘12 ‘do

= — e :—
82 111((10) 2 do

Finally, Forming limit diagrams determined in different laboratories
may be somewhat different. This is because the determination of the ﬁrst
perceptible neck is subjective. Many techniques have been used to determine
the onset of the neck; visual inspection and feeling with the ﬁngers. Also,
there is no restriction placed on how far one should go from the neck before
considering a circle to be safe. It is recommended that safe readings should
be made at a distance at least equal to 1.5 circle diameter from the center of

the neck. Below are two pictures of two different failed specimens.

    

- an». a lawn-
Figl: Fractured specimen

 

FigZ: Failure under biaxial tension

CHAPTER 2: MK Model

2.1 Introduction

The MK model is the ﬁrst mathematical model presented by Marciniak
and Kuczynski [3] to help understand localized necking in sheet metals. They
assumed a long groove, a locally thinned region within which the strain will
concentrate and will lead to failure due to the onset of localized necking. The
origin of this imperfection is a variation in the sheet thickness [3], difference
of the material strength [16] or a combination of both. The failure is
assumed to occur when the ratio of the strain increment within the groove
and outside of it reaches a critical value. There are some constraints on the
orientation of the imperfection where it has been taken to be perpendicular
to the axis of the largest principal stress and strain for positive major strains
(RHS of the FLD), and in the direction of zero-extension for negative major
strains (LHS of the FLD). In this section we will explain the algorithm for the
calculation of FLDs based on the MK model. We will investigate the effect of
the hardening exponent n, the strain rate sensitivity exponent m, the
imperfection parameter f, the R-value and the yield function exponent a,
where the yield locus will be described by Hill’s 1979 yield function. We will
investigate yield function differences in predicting localized necking. The

corresponding algorithms are all written in Matlab.

2.2 Right hand side of forming limit diagrams

For positive minor strain the groove is assumed to be perpendicular to
the axis of the largest principal stress and strain. The strain increments
parallel to the groove within and outside of it are equal. The algorithm can be
summarized by assuming, for a known strain path, the largest principal strain
increment in the imperfection (region b) and then calculating the one in the
homogeneous region (region a). The failure is assumed to occur when their
ratio reaches a critical value (in the calculation below the critical value is
assumed to be equal to 10). The algorithm below is written for von Mises
yield function. The changes that will need to be implemented for other yield
functions are described in the deﬁnition of each yield function. The Matlab
programs are written for von Mises (Appendix A), Hill’s quadratic (Appendix
B), and non-quadratic yield functions (Appendix C). In the algorithm below
indices a and b refer to the homogeneous and the defected region
respectively. Indices 1 and 2 refer to the principal values. The ﬁgure below

describes well the model.

 

 

 

 

 

FigB: RHS model

2.2.1: Algorithm

Under plane stress condition (03 = O), for a known strain ratio in region a:

The stress ratio is calculated:

0

a =54=(2p.+1)/(p.+2>
0'

al

The yield function is then evaluated:

(1)“ = E“- =‘/a:-aa +1

Impose a value for the strain increment in the groove:

 

dab, = 0.001

Assume a value for the strain increment in region a:

algal

Calculate the rest of the variables needed for region a:

(1802 = padgal
dga3 = —d£al ’dgaz
5;“ : d8al(1+ aapa )/(Da

3;.
ﬂa — d5

 

al

The equation below holds due to the compatibility requirement:

dab2 = (1802

The rest of the variables in region b are evaluated:

dgb3 = _dgbl -—d£b2

E0 = d8bl(1+ 0‘be )/q)b

 

 

35.
ﬂ. — dab]
dabz

b 2:
dab,

ab = (2,01; +1)(pb + 2)

(I),J =,/a: —ab +1

From the force balance equation between region a and b the largest strain
increment in the region a is recalculated. Then it is correspondingly adjusted
and the process is repeated until the force balance is satisﬁed. The force
balance equation is derived below:
F

:Fbl

For a unit width and thickness t, the force in the 1-direction is:

Fl 2 (at = 5(g} = [2}
0' (D

The hardening rule with strain rate dependence can be written:

n-'-m

EzkE 5

Using the variables deﬁned above the force can be expressed by:
F, = kt(§ + d5)" a" /<1> = kt(§ + d3)” (ﬂ/p)§,"' /q>

By applying the above force expression to regions a and b, and after some

adjustments, the force balance can be reduced to:

éta + dga Y[§—:] 55:2 = éexpﬁgb " £30 Xg” + dgbr(—E:L] big;

Where the imperfection parameter f and the current thickness are:

f : kbtbO /kataO
t: to exp£3

To solve the force balance equation, the Newton-Raphson or the Mid-point
method can be used. Both methods were programmed in Matlab, but the
Newton-Raphson proved to be faster. The same algorithm can be used for
any yield function provided that the equations that involve the use of the
yield function are changed correspondingly. Below, the algorithm will be

applied to several different yield functions.

10

2.2.2 von Mises

The von Mises yield function is the most basic yield function describing
the plastic behavior of isotropic materials. An isotropic material responds in
the same manner when strained or stressed in any direction. All the
derivations needed for the implementation of the algorithm for von Mises
yield function are shown above. In plane stress applications the yield function

is:

232 =(0'I —0'2)2 +032 +0“:

The calculated forming limit diagrams below are for von Mises yield function
where the imperfection parameter f, the hardening exponent n and the strain
rate exponent m are variables. The FLD curves shown below have been
plotted using the MK model. The groove does not change direction when the
deformation occurs. So it is assumed to be perpendicular to the direction of
the largest stress and strain and does not rotate during the straining of the
sheet. The plots shown below indicate that formability increases and failure
occurs at higher strains, with increased strain rate sensitivity and strain-

hardening exponent, and decreased imperfection parameter.

11

Major Strain

Major Strain

0.9

0.8

07

0.6

OJ

0.9

0.1

von Mises - different n values

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Minor Strain
FigS: Different m values

12

0.28 I
* $336
(—(622
t (—0.2 -
F / I
f = 0.998
F m = 0.0001 R
l l i L l l l l
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Minor Strain
Fig4: Differentnvalues
von Mises — different m values
n 0.2
f f = 0.998 "
L l l L l l l I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Major Strain

0.8

von Mises — different f values

 

 

 

l T r I I

(—(1998

 

 

 

 

 

 

.1
— 0.2 j
m = 0.0001
I I I L I J I
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Minor Strain

Fig6: Different f values

13

2.2.3: Hill’s Quadratic 1948

The Hill’s quadratic yield function was one of the ﬁrst yield functions to
describe the plastic anisotropy of metals. Hill has formulated a quantitative
treatment of plastic anisotropy without regard to its crystallographic origin. A
new parameter R has been used which is the ratio of plastic strain in the
width and thickness direction for a tensile specimen cut from a sheet. The R-
values usually vary with the test direction and it is common to characterize a

material by an average R-value written below.

 

Where the indices represent the angles (with the rolling direction) to which
the tensile specimen has been cut from the sheet. For R equal to 1 the
quadratic hill's yield function will simply be reduced to von Mises. For plane

stress loading the yield function simplify to the equation below:

_ 2 2R
0 = a, +0, ——0',az
1+R

The only two expressions in the algorithm that should be changed in order to
get the forming limit diagrams using Hill’s 1948 are the stress ratio and the
yield ratio. The stress ratio is a function of the strain ratio and the yield ratio

is a function of the stress ratio. Those two functions are written below.

14

09

08

07

06

05

Major Strain

04

03

02

01

a = (R + p(l +R))/(1 + R(p+1))

 

(DzJﬁsz +1+R(1—a)2)

+

Quadratic Hill's - different R values

 

 

 

 

 

 

 

 

I I T r I I I I T
(—0.75 “
+—1 I
(—1.5 T
+—2 1
q

n = 0.2

m = 0.0001

f = 0.998 "

i 1 1 1 1 1 l I I
01 02 03 04 05 06 07 03 09 I

Minor Strain

Fig7: Different R values

The ﬁgure above shows that the forming limit is very sensitive to the value of

R and it decrease with increasing R. The forming limit for this case is

calculated with the largest strain assumed to occur in the RD direction.

15

2.2.4: Hill’s Non-Quadratic 1979

Hill’s quadratic yield function overestimated the effect of the R-value on
the shape of the yield loci. A new yield function has been developed to
overcome this weakness where a higher exponent “a" has been used to
account for crystallographically textured fcc and bcc metals. With a equal to
2 this criterion simpliﬁes to Hill’s quadratic. However, crystallographic
calculations are much better approximated by a being equal to 6 for bcc and

8 for fcc metals. For plane stress loading this criterion simpliﬁes to:
E“(1+R) = a,“ +0; +R(c7l —02)“

The yield ratio for this yield function can be expressed explicitly as a function
of the stress ratio but the stress ratio cannot be expressed explicitly as a
function of the strain ratio. A nonlinear equation should be solved in order to
ﬁrst get the stress ratio and then evaluate the yield ratio. The Newton
Raphson method or the Midpoint method can be used to solve this equation.

Below the non-linear equation and the yield ratio are written:

p(1+ R(1— a)“l )- a“ + R(l — a)"" = 0

<1) = [ﬁv +1+R(1—a)“)jz

16

Hill's 1979 - different R value

0-5 ﬂ I I I I I ﬂ I ‘I

 

 

 

 

 

 

 

 

 

045 a
04 ~
035 ~
g 0.3 -
m
u
U
U) 025 _
H
o
.3 02 a
z .
015— _
0.1 — n = 0.2 ‘
m = 0.0001
0.05F f = 0.998 2
a = 8
0 I 1 I I I I 1 1 1
0 0.05 0.1 0.]5 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Minor Strain

Fig8: Different R values

The plot above proves that the forming limit diagram is insensitive to the
level of R, where it was very dependant when the quadratic hill's yield
function was used. This independence of the FLD of the R-value gave the
non-quadratic function a great importance because it matched with the

experimental FLD where it has been shown to not depend on R values.

17

Hill's 1979 - different a values

0-5 I T T If ﬂ I I T I

 

045

04

035

03

025

Major Strain

02

 

015~ c

 

OJ ~

0.2 a
0.0001
- 0.998

005—

w m B S
I l
N

 

 

 

 

 

0 1 I l 1 1 1 1 I 1
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Minor Strain

Fig9: Different a values

 

As mentioned before the exponent a, connected to the crystal structure, is
equal to 6 or 8 for bcc and fcc alloys respectively. For the special case where
R is equal to 1 and a equal to 2, Hill’s 1979 will simply reduce to von Mises.

The plot above shows a dependence of the FLD on the value of a.

18

2.3: Left hand side of forming limit diagrams

For the negative minor strain regime the groove direction is assumed to
line up with Hill’s direction of zero-extension such that the plane 'strain
condition is satisﬁed. The existence of a critical thickness strain criterion for
localized necking of sheets in the negative minor strain regime has been
demonstrated. The effective limit strain in region a increases as the strain
ratio becomes more negative and tends towards uniaxial tension. The limit
thickness strain is a constant regardless of the imposed strain ratio, thus
indicating the existence of a critical thickness strain for the onset of localized
necking. The limiting thickness strain is dependant only on n, m and f. It is
independent of the plastic anisotropy, the yield function and the plasticity

theories (ﬂow or deformation theory).

 

 

 

 

 

 

 

 

FiglO: LHS model

19

2.3.1: Algorithm

The algorithm is simple here because it involves the solution of a
nonlinear equation with the thickness strain. A thickness strain increment is
assumed in the imperfection and from the nonlinear equation the thickness
strain in the homogeneous region is calculated until the thickness strain in
the homogeneous region reaches a limiting value. This equation is solved
using either Newton Raphson method or the Mid-Point method. After lengthy
mathematical manipulation where the force balance equation coupled with
the condition of Hill’s zero extension have been used, the nonlinear equation
simpliﬁes to the equation below. The full derivation of this equation is in

[17].

(d8: 7" (8:. + a: )" exp<e§> = flde: )"' (sf. + £5)f exp<a§>

20

2.3.2: Thickness limit

The limiting thickness strain in the homogeneous region is independent
from the strain ratio. For each value of n, m and fthe limiting strain in region
a is constant. The variation of the limiting strain with n, m and f is shown in
the plots below. By increasing the values of n and m, the limiting thickness
strain increases by keeping a constant ratio with the thickness strain in the
imperfection. By increasing the f values, the limiting thickness strain will
increase with a small increase in its ratio with the limiting strain in the

imperfection.

Thickness limit — different n values

 

I I I I I I

0.25 — (— 0.28 .

e-IIZO
ﬂ4
0.2 ~ ﬂ-ZZ .

/+—_0.2

OJ ~ a

— Thickness strain in region a

005~ a

.0001
.998

 

m S
n n
o o

 

 

 

 

 

O l l 1 l l l
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
- Thickness strain in region b

 

Figll: Thickness limit Different n values

21

- Thickness strain in region a

-Thickness strain in region a

025

02

015

01

005

02

018

016

012

O1

008

Thickness limit — different m values

 

 

I T I I I

(—(101

 

(— 0.005

 

f = 0.998

 

 

 

l I l I 4

 

 

0] 02 03 04 05
- Thickness strain in region b

Fing: Thickness limit Different m values

Thickness limit — different f values

 

 

 

0.2
0.0001

n
m

 

 

 

1 I l l

 

 

0.05 0.1 0.15 0.2
—Thickness strain in region b
Fig 13: Thickness limit Different f values

22

025

2.3.3: Results

The plots below show the effect of the hardening exponent n, the strain
rate sensitivity parameter m and the imperfect parameter f. The FLD for the
LHS increases with the increase in the values of n, m, and decrease with the
increasing value of f. This is reasonable because an increase in the

imperfection parameter means a deeper defect, which will lead to an earlier

 

  
  
 

 

 

 

 

 

failure.
LHS — different n values
0-5 f I . f f ._—_
028—)
0.45— 026 _; a
0.24
04» —+
022—»
c:
‘3 0.351
I: 0.2..
U)
H
3. 0.3
In
21
025» \
0.2 e r— 7 ~ a
m = 0.000;j
f = 0.998 I
L.
‘ l l l l
-0.25 -0.2 -O.15 -0.1 -0.05 0

Minor Strain

Fig 14: LHS Different n values

23

Major Strain

Major Strain

05

LBS — different m values

 

045~

9

w

VI
1

9
I»
I

025—

F I I I I

 

02»—

 

 

0.998

 

 

 

 

l 4 l I l

 

-0.25 -0.2 -0 15 -0.1 -0.05 0
Minor Strain

Fing: LHS Different m values
LHS - different f values

 

035

03»—

53
N
u
T

9’
N
I

 

 

I I I I I I I I I

 

 

0.2
0.0001

 

 

 

I I 1 L L I l l I

 

01
-02

-018 -O.l6 -0.14 -0.12 -0.l -0.08 -0.06 -0.04 -0.02 0
Minor Strain

Fig 16: LHS Different f values

24

CHAPTER 3: General MK Model

3.1: Introduction

In this chapter, we will explain in detail a general approach for the
calculation of forming limit diagrams. This present model is an extension of a
previous analysis by Marciniak and Kaczynski (MK-model), where it is
generalized to account for all possible directions of the imperfection. The
limiting strain is calculated for all possible angels and minimized versus the
angel of the groove. We will use this approach to compare different yield
functions starting with the most basic one for isotropic materials von Mises
and two anisotropic yield functions Barlat Yld96 and Barlat Yld2000 proven to
describe well the plastic behavior of aluminum sheet metals. Two work-
hardening laws, swift and voce, will be used to model the effective stress-
strain behavior of the material. This new approach uses the Newton-Raphson
method to solve for the transverse stress and effective strain in the groove.
For this reason, the second derivative of the yield function is needed for the
evaluation of the Jacobian. The ﬁrst and second derivatives of the three yield
functions are evaluated in the appendix for implementation purposes.
Forming limit diagrams will be calculated for aluminum alloy Al2008-T4 and
will be compared with experimental data. Barlat calculated the FLD for
Al2008-T4 by incorporating micro structural parameters, internal damage
(voids) and crystallographic textures in his model and showed very accurate

results. All computed FLDs assume that the largest strain is in the y-direction

25

(TD direction) in order to be on the safe side. A very good agreement
between the calculated and the experimental FLD was obtained, when the
Yld2000 and the voce hardening law were used. M.C. Butuc [12-13] already
developed FLD for the Al6061 aluminum alloy based on von Mises, Hill’s

1948, Hill’s 1979 and Barlat Yld96yield function.

26

3.2: Modified M-K model

This model uses the same assumption proposed by Marciniak and
Kuczynski [3] who assumed a long groove, a locally thinned region within
which the strain will concentrate and will lead to failure due to the onset of
localized necking. The origin of this imperfection is a variation in the sheet
thickness [3], difference of the material strength [16], or a combination of
both. The failure is assumed to occur when the ratio of the strain increment
within the groove and outside of it reaches a critical value. As we mentioned
earlier there is some constraints on the orientation of the imperfection. The
modiﬁed M-K model will consider all possible orientations for the groove and
consider the failure to occur when the ratio of the total strain increment

inside and outside of the imperfection reach the value of 10.

 

 

 

 

 

 

 

 

Figl7: Modified MK model

27

The model above consist of two regions, the homogeneous region “a” where
the axis 1, 2 and 3 represent the principal directions of stress and strain and
the groove region “b” where the axis n, t and 3 are the normal to the
direction of the groove, the longitudinal direction, and the thickness
direction. The M-K imperfection parameter is deﬁned as the initial thickness
ratio of the defected and the homogeneous area .The deformation within the
groove occurs at a faster rate than the one in the homogeneous region and
the plastic ﬂow will localize in the imperfection area leading to eventual
failure. Under deformation, the strains parallel to the imperfection direction
within and outside of the groove are equal and forces perpendicular to the
groove in both regions are the same. Those two constraints are called
compatibility and equilibrium conditions. The groove is assumed to rotate
while the deformation occurs. Let e and e0 be the current and the initial sheet
thickness. The value of the imperfection parameter, the compatibility and the

equilibrium conditions can then be written as

Imperfection parameter: f = b/e"
Compatibility condition: d8; = def;
Equilibrium condition : 0'“ e“ = 0” e” , 0"e“ = abeb

28

3.2.1 Algorithm

Given the stress state in region “a”, the plastic strain is evaluated in
region “b” for a known orientation of the groove. In the ﬁgure above we will
assume the x-, y- and z-axis to be the orthotropic axes of the material
(rolling, transverse and normal) that line up with axes 1, 2 and 3
respectively. For one increment we will explain clearly how to calculate the
plastic strain in region “b”. First, assume a small increment of the effective
strain in region “a” and then follow the steps below. The details below show
how this algorithm should be used in order to calculate the FLD for
proportional loading, by assuming that no shear stress exists in region “a”.
This algorithm could be used to ﬁnd the critical strain for an unknown strain

path where the strain state will be an input from an FEM analysis.

3.2.1.1 Region “a”

. Strain ratio
Usually in the calculation of forming limit diagrams for linear strain path the

strain ratio in region “a” will be an input.

a
dag

dab

II

p:

29

. Stress ratio
The stress ratio is then evaluated by assuming that it is a variable equal to
the stress in the y-direction, if the stress in the x-direction is set equal to 1

and shear stress equal to zero.

 

For most yield functions, to get the stress ratio one will need to solve the

following nonlinear equation:

1
a . 6 ,-
—J’F—p—0” =0 For 0': a
60' 60'

yy xx 0

d5; — pair):r = 0 :>

The Newton Raphson or Mid-Point method could be used for this purpose. For
the Mid-Point method the initial interval where the solution is located can be

taken as [-1; 2].

. Yield ratio

The yield ratio is the ratio of the yield function to the stress in the x-
direction. The yield ratio becomes the yield stress by assuming the stress in
the x-direction to be equal to 1, the stress in the y-direction to be equal to

the stress ratio, and the shear stress to be equal to zero.

QI

 

, ¢=0yr For a: a

C“
O

30

. Stress state

Having assumed an increment for the effective strain, the yield stress can be
evaluated using the hardening rule (swift or voce). Coupled with the yield
ratio, the stress in the x-direction can be calculated, knowing the stress in
the y-direction, and by setting the shear stress equal to zero (no shear stress

in region “a”).

3”,, = hardening?“ + did)

00 =EHR/(D

U
a _ a
a”, —a0ju
a —
0'0 —0

. Strain state
The strain state is calculated from the Flow Rule knowing the stress state in

the orthotropic referential frame.

 

 

0'
a _ — YF a _ a a
dart-d8a¥:81x_£u+d8u
XX
60'
a _ - YF a _ a 0
d5” — den 0'“ :> a”, — a” +d5y)
W
160
a _ - _ YF a _ a a
dew—d :>£ —5xy+d£xy

31

. Output needed

To solve for the effective strain increment in region “b" four values are
needed; the strain increment parallel to the groove, the strain in the
thickness direction, the stress normal to the groove and the shear stress in

the groove coordinate axis. Those four values are calculated as following:

a _ a - 2 a 2 a -
d8" —d£M sm 6+d£w cos 6—day sm 26
a __ a 2 0 ° 2 a '
0",," —O-xx cos 0+0” sm (9+0; sm26
of, = 0'1‘,’y(cos2 6— sin2 (9) + (va — 0;)sin6cos6’

a __ a 0

3.2.1.2 Region “b”

. Nonlinear Equations

For this region to complete our study we have to compute the effective strain
increment and compare it with the one in the homogeneous region. Two
nonlinear equations in the effective strain increment and the stress in the
longitudinal direction must be satisﬁed. The ﬁrst state that the yield stress
computed from the yield function must equal the one computed form the
hardening rule. The second is the compatibility equation. The two equations

are written below.

32

GIIdEb’ab)= 35R ’50:? = 0

fl

G,(d§b,a{: )= d5: — d5: = 0
To solve these equations the Newton-Raphson method is used. This method
involves the computation of the Jacobian needed in the iterative evaluation of

the solution. The solution is reached by iteration using the formula below:

dgil-il _ d5? _J—l Gl(dEib’U:i)
Grim — 6:; 02(dgib’O-3i)

Where J is the Jacobian of the nonlinear equations. The Jacobian is the
matrix of the derivatives of G1 and G2 with respect to the effective strain
increment and the stress in the longitudinal direction. For this case the

Jacobian is symmetric and its components are described below:

_

 

 

' 00, at},

J = 6d?” 60':

at}, at},
Lads—‘1’ 60'": -

 

 

The ﬁrst component of the Jacobian is straightforward and is obtained by

deriving either hardening rules (Swift or Voce)

6Gl nK(§0 + 5” + (15” )""l Swift hardening

6d?” ‘— C x 3 exp[— C(g” + (15” )] Voce hardening

 

33

As we mentioned above the Jacobian is symmetric and the non-diagonal

component involve the ﬁrst derivative of the yield function, which is derived

in the Appendix.

 

 

 

sin26
661 = 602 = _60,,. 63:, 60W cosz 6
60': 6d?” 60': 60'?v 6 0': ,
~ y —sml9cost9

The derivative of G; with respect to on involves the second derivative of the

yield function, which is derived in the appendix.

 

 

 

a__20'rr 625:)? 623W "
6203; 60'2 60' 603x610}, ' sin20
26—i— = —alE”[sin2 6 cos2 6 — sin 26?] ii— 6 0,", -a—UlF—— cos2 (9
60'“ 60'” 60'“ 60'; 60'” 60'”, ,
_1_ 5223;; 1 623,5} 1:3}; L—smdcosl?
L2 60560” 2 6030,60”. 2 60'”, _

 

. System Solution

To solve the system of the two nonlinear equations we will use the Newton-
Raphson method. The initial values are assumed to be the one of the
homogeneous region. Using the equilibrium condition we are able to compute
the normal stress to the groove and the shear stress in the groove. We
assume the initial imperfection parameter to be the current for the ﬁrst
iteration so we can proceed with the calculations and provided that the

variation of the imperfection parameter is too small at the beginning of the

34

deformation. Using the rotation rule below and the flow rule the stresses and
strains in the groove in the orthotropic referential frame are determined and
ﬁnally we calculate the strain increment in the groove in the longitudinal

direction using the rotation formula:

 

... . INITIAL [ ‘15:] [ ‘15.]
d5 ,0'" -) b = a
VALUES Ono or:
of. = as. /f
a; ,o.‘.'. ,f EQUILIBRIUM-) b
a... = a: //
0'2" 0:, = of" 0052 6 + 0': sin2 6 — 0:, sin 26
0': ROTATION-9 0;, = 0:, cos2 6 + 0': sin2 6 + 0:, sin 26
0'5: 0;; = of, (cos2 6 — sin2 6) — (0': — 0:")sin 6cos6
0'; 5:.
of; FLOW RULE-) sf;
0'; 8:;
81’.
a; ROTATION-9 def; = def; sin2 6 + dag, cos2 6 — defy sin 20
b
a”

35

Now having all the values needed to solve the nonlinear equation, the strain
increment in the groove and the stress in the longitudinal direction of the
groove are computed by iterating until a convergence value is reached. After
that the current imperfection parameter and the angle of the band are
updated. The process continue by assuming an effective strain increment in
region a and calculating the one in the groove until the ratio of these two
approaches a critical value corresponding to local instability. It is usually
assumed that failure occurs when the strain increment in the homogeneous

region reaches 10 percent or less than the one in the groove.

 

Imperfection parameter: f = fo exp(g§’ —£;’)
. . 1+ d6“
Groove orIentatIon : tan(6 + d6) = tan 6 "‘
1+ d6;

The study above is done assuming the major strain to occur along the rolling
direction (x-axis), but to obtain the full FLD we have to assume the case
where the major strain occur along the transverse direction where it will give
the lowest values for the strain and will be used as the FLD for material. To
calculate the FLD where the major strain occurs along the (y-axis), the only
thing to do is to switch the x and y material coefﬁcient and use the same
algorithm. For Vonmises both FLDs are the same due to the symmetry of the
yield function with the rolling and transverse directions, but for Barlat Yield96

and Barlat 2000 the full FLD is not symmetric with the ﬁrst bisector.

36

 

 

3.3: Yield Functions

In this work three yield functions are used to calculate the theoretical
FLDs for Al2008-T4, von Mises, Barlat Yield96 and Barlat Yield2000. A
description of the yield functions and the derivation needed for the

implementation of the codes are presented in the appendix.

3.4: Hardening rules

The coefﬁcients of an effective stress-strain curve that deﬁnes the strain
hardening of a material is found by ﬁtting the experimental data obtained
from a uni-axial tensile test to one of the following equations:.

Swift 5 = KG, + a)"
Voce E = A — Bexp(—CE)
Where 0’ and E are the effective stress and effective strain respectively, and

the rest of the variables are constants that are obtained from theby ﬁtting

the above equations of to experimental results from a tensile test.

3.5: Material constants

The material constants, needed for the implementation of the code for
aluminum alloy Al2008-T4, describes for both two hardening rules and yield
functions are presented below. The imperfection parameter is assumed to be

0.996.

37

 

 

 

 

 

Swift Voce
50 0.002 A 447 Mpa
n 0.1882 B 248 Mpa
K 456.82 Mpa C 4.3 Mpa

 

 

 

Table 1: Constants for voce and swift equation — Al2008T4

 

 

 

 

 

 

 

 

 

Yield96
a C1 C2 C3 C5 (1x (1y Orzo all
8 1.1227 0.9154 1.0137 1.0578 1.55 1.75 1.0 0.6

 

 

 

Table 2: Constants for Barlat Yield96 — Al2008T4

Yield 2000

 

a (II (12 (13 04 (15 (16 (17 as

 

 

 

 

 

8 0.91777 1.14937 1.17360 1.14556 1.05262 1.19593 1.01494 1.06188

 

 

 

 

 

 

 

Table 3: Constants for Barlat Yield2000 - Al2008T4

 

0.6 ~

TD Strain

  
 

— von Mises

—o— Yld96

+ Yld2000
0 Experimental

 

 

T

0.2 0.4 0.6

 

 

-O.4 -0.2 0
RD Strain

Fig 18: Flds Al2008-T4-Swift

38

 

 

   

 

 

 

 

0.6 J

.S

«I

t:

(I)

Q

I—
—— von Mises
-——o— Yld96
—+— Yld2000

0 Experimental
{1
~04 -0.2 O 0.2 0.4 0.6

RD Strain

Fig 19: Flds Al2008-T4-Voce

3.6: Results and discussions

The experimental and the computed FLDs for von Mises, Barlat Yld96 and
Barlat Yld2000 yield function for Al2008-T4 using Swift hardening law are
shown in Fig-18. It can be seen that von Mises over predicted the
experimentally measured failure points due to the anisotropy of the metal
and Yld96 and Yl2000 have given good result with Yld96 being more
conservative. In Fig-18 the FLDs are calculated using Swift hardening rule
where it is shown that failure points are well predicted when Yield96
describes the yield locus. For Al2008-T4 it is best to use Yield96 to describe
the yield locus and Swift equation to model the hardening of the material to

be on the safe side.

39

CHAPTER 4: Simulation of Stretch Stamping

4.1: Introduction

In many sheet metal forming simulations a forming limit diagram is used
to predict failure. In sheet metal stamping simulations, the minor and major
strains for all elements are plotted against each other and compared with
with the FLD curve to check whether or not strains at any element exceed
the limiting curve. Usually the FLD is generated using both Voce and Swift
hardening rule, and the most conservative one is used to check for the
possibility of the development of a localized neck in the sheet. In this chapter
stretch stamping of a sheet metal with a hemispherical punch is simulated
using LS-DYNA and the major and minor strains of all the elements are
tracked through the simulation to see how and when they exceed its the
FLDlimit. When one element reaches the limiting curve the simulation is
stopped and the limiting punch depth is recorded. Although strain paths in
this type of simulations are not linear, FLDs are still used to check for
necking and has been shown to provide qualitatively accurate and failure
validprediction. The pure-stretch stamping test of a Al2008-T4 sheet using a
hemispherical punch has been done experimentally and initial failure
occurredwas recorded for at a punch depth of 27mm at room temperature.
The ﬁnite element simulation result with LS-DYNA alsos matched the
experimental results, where it was showed predicted that the ﬁrst time the
strains in a ﬁrst element exceeded reached the FLD, occurred when the

punch depth was at 27mm.

40

 

4.2: Description of the model

In pure-stretch stamping, a sheet is clamped between a die and the
blank holder, and a punch is used to deform the sheet. A full model has been
used in this simulation to capture the anisotropic behavior of the AA3003-
H111 sheet. The dimensions of the part, the material properties, and all the

data needed to model this experiment is described below.

4.2.1: Part dimensions

The punch diameter is 4 in. The diameter of the die and the blank holder
is 7.5 in. The diameter of the sheet is 3.54 in. The height of the die is 1.81
in, and 0.4 in for the blank holder. The radius of the surface of the die and
blank holder parallel to the sheet is 2.3 in. The sheet is assumed to have a
uniform thickness of 1mm. The initial distance between the die and the sheet
was taken to be zero, and the one between the blank holder and the sheet
was taken to be 0.5 mm to prevent any initial penetration. The initial
distance between the punch and the sheet was assumed to be 3 mm. The
radial distance between the punch and the blank holder was assumed to be

0.7mm.

4.2.2: Element type and material models

For all parts the “Belytschko-Tsay” (Shell element number 2 in LS-DYNA)

was used, using 5 integration points for the sheet and 2 for the rest of the

41

 

parts. The punch and the blank holder were modeled as rigid materials

 

“Steel” (material number 33 in LS-DYNA) and constrained to move only
normal to the initial plane of the sheet. The die was modeled as rigid “Steel”
but constrained to remain ﬁxed during the forming process. The sheet was
aluminum AA3003-H111 and modeled to follow BARLAT-YLD2000 (user
material option has been used) and assumed to follow Voce and Swift

hardening rule.

 

4.2.3: Material properties

The material constants for the aluminum and steel are shown below.

 

 

 

 

 

 

 

 

 

 

 

 

Steel

Young’s modulus 210 GPA

Poisson ratio 0.3

Density 7850 kg/m3
Table 4: Properties of Steel

AA3003-H111

Young’s modulus 70 GPA

Poisson ratio 0.33

Density 2710 kg/m3

 

 

 

 

Table 5: Properties of AA3003-H111

42

 

 

 

 

Voce

A 144.015 MPA
B 76.217 MPA
C 13.058

 

 

 

 

Table 6: Constants for voce equation

 

 

 

 

 

 

Swift

Ea 5.7E-4

n 0.2157

k 199.82 MPa

 

 

Table 7: Constants for swift equation

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Barlat YLD2000
a 8

L111 0.620325624
L112 -0.310162812
L121 -0.365965106
L122 0.731930213
L155 1.069533131
L211 0.681653114
L212 -0.331290171
L221 -0.350798498
L222 0.699980075
L255 1.06096938

 

 

Table 8: Constants for Barlat Yield2000

See Appendix C for the deﬁnition of the parameters above.

43

4.2.4: Constraints and prescribed motion

The die was assumed to follow a prescribed motion (displacement), with
the corresponding curve shown below. A force was applied on the blank
holder to clamp the sheet to the die. This force started from zero and
increased until it reaches 450 KN after nine seconds. A plot that describes

this motion is shown below. A picture of the full model is also shown below.

M3003-H111. RT

 

FigZO: Model stretch stamping

The plots show the prescribed motion for the punch and the blank holding

force.

Blank holder - force cuve

 

 

Force KN

2001

‘u

 

 

0 r T I I I I
0 5 10 1 5 20 25 30 35
Time ms

Fing: Blank holder force

 

 

Punch - Prescribed tmtion

 

3 8 8 8

Displacement mrr

a
O

 

 

 

0 5 1 0 15 20 25 30 35
Time ms

Fi922: Punch prescribed motion

45

4.3: LS-DYNA

After calculating the theoretical FLDs of the material using the program
described in the previous chapter for both voce and swift hardening rules, the

more conservative curve was used to check for failure.

4.3.1: Swift and Voce FLDs

For this case the FLD curve based on voce hardening rule turned out to
give the more conservative FLD. Both FLDs were calculated with the larger
principal strain assumed to occur in the TD direction to be on the safe side.
The material coefﬁcients for swift and voce hardening rules are written in the

previous paragraph. Both voce and swift FLDs are plotted below.

 

0.5 ~

 

 

 

 

.S
(G
t:
m
E 0.25 4
— Swiﬂ
— Voce
. c . ,
-O.3 ~01 0.1 0.3
RD Strain

Fi923: FLD AA3003-H 1 1 1

46

 

4.3.2: LS-DYNA FLDS

In LS-DYNA, to open a FLD curve, the points should be c0pied to a txt ﬁle
starting with “fld_user” and followed by double the number of points. For
AA3003-H111 the input FLD ﬁle is written below. In LS-DYNA the curve is
read and plotted. The FLD curve used in this simulation is plotted below.
Usually many curves are generated from one data input ﬁle. The most useful
ones are the original FLD (Cracks curve) and a more conservative (Risk of
Cracks curve) curve to keep one on the safe side. Materials are considered to
be safe if they do not cross the crack curve. But it is better to be on the safer

side where one should not cross the risk of crack curve.

 

 

ﬂd_user

36

0.178426 -0.08921
0.168248 -0.0673
0.159279 -0.04778
0.150106 —0.03002
0.141445 -0.01414
0.132489 0.000002
0.137678 0.011016
0.147395 0.023585
0.15628 0.03751
0.1647 0.052713
0.17261 1 0.069068
0.180267 0.086552
0.187824 0.105236
0.195239 0.124963
0.202804 0.146084
0.210229 0.168262
0.217641 0.191588
0.221089 0.212251

 

Table9: FLD user

47

 

Major Engineering Strain (9B)

LS—DYNA FLD VOCE

 

 

 

 

 

 

 

 

 

 

 

F1924: LS-DYNA FLDS

48

 

 

25
l.
20 j/
//
.mi .7 -

15 \\\// _[
10

5

0 l l l I

-5 0 5 10 15 20
Minor Engineering Strain (96)

 

 

4.4: Experimental Findings

The pure-stretch stamping test has been performed in the lab at room
temperature. The picture below shows the fractured sheet at the depth of 27

mm.

 

FigZS: Failed sheet

49

4.5: Numerical results

The computer simulation was stopped when the strains in an element
crossed the forming limit diagram. The plots below show the strain
distribution of all elements after 10 ms and at the failure time (23 ms). The
strain paths range from plane strain for the element on the sides of the
sheets to biaxial tension for the elements in the middle of the sheet. The

displacement of the middle element is plotted too.

LS-DYNA FLD (1:15 ms)

 

 

 

15 \ ./
_ \/
10

 

Major Engineering Strain (‘5)

Ch

 

 

 

 

 

 

 

 

 

 

 

-5 0 5 10 15 20

Minor Engineering Strain (96)
Fi926: FLD AA3003-H111 (t: 15 ms)

50

Major Engineering Strain (96)

Displacement - mm

LS-DYNA FLD (1:23 ms)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_ fwd...—
20
15 _\\\ ,/
10

0 J l I l
-5 0 5 10 15 20
Minor Engineering Strain (96)
FigZ7: FLD AA3003-H111 (t: 23 ms)
Midde Element
0
_ \
\
-5 \\
-10 \
\\
-15 \\
-20
\.

'25 \\
_30 1 1 I I

0 5 10 15 20

Tune - ms

Fi928: Middle Element Displacement

51

 

 

The failure of the sheet has been found to take place in the cracks region
in the picture below. In this region the strains in some of the elements have
crossed the forming limit diagram and the sheet was considered to have
failed due to the onset of localized necking. A comparison of the actual and
predicted failed sheets indicates that a good qualitative match exists between

experimental and numerical results.

Fonnabll It
Contours of Formablllty: Mid. Surface "y ”y

Cracks

Risk
of cracks

 

Severe
thlnnlng

 

Good

 

 

Inadequate
stretch

ernldlng
tendenqr

Wrinkles

 

   

FigZ9: FLD contour at failure

52

CHAPTER 5: Conclusions

The experimental and numerical results were very close. First the punch
depth at failure was 27 mm for the experiment and 29 mm for the numerical
simulation, which corresponds to a discrepancy of only 7.5%. More
importantly, the failure occurred in the same location in both experimental
and numerical simulation. The results below show that forming limits are a
very useful tool to check for the failure of stamped sheet metals in numerical
simulations. Experimental and numerical ﬁndings are compared in the picture
below. The discrepancy between experimental and numerical ﬁndings can be
a cause of some accuracy issues or even because the strain paths in the
sheet are not linear for all elements. For elements on the boundary of the
deformed sheet the deformation could be assumed to be linear, however it
becomes non-linear as one moves toward the center of the sheet. For the
pure-stretch stamping it has been shown [10,11] that the FLD lowers, in
comparison with linear path deformations. This might explain the difference
between the actual and predicted punch depth at the time of failure. And
ﬁnally in all the numerical simulations we have assumed failure occurs when
the strain increment in the imperfection region is ten times the value of the
strain in the homogeneous region, and by assuming an imperfection
parameter. All these assumptions could be partially responsible for the

discrepancy between the experimental and numerical failure punch height.

53

 

Fig30: Experimental and numerical ﬁndings

54

Appendix A - von Mises

It is the yield function that describes the yield locus of isotropic material.

A1 - Yield function-2d

_—2_ 2 2 2_
¢—0' —0'U+0'W+30'xy 0'

0'

)9’

Where 5 is the yield stress and the variables are the Cauchy stresses in the

x-direction, y direction and the shear stress respectively.

 

 

 

 

 

 

A2 - First derivative
25 05 = 005
60',I 60,,
Where n refers to xx, yy and xy.
ﬂ¢—-~20' — jig—=20 ‘0'” i=6
60'” "‘ y’ 00'”, yy 00”, ’3
A3 - Second derivative
626" _i 1 62¢ _ 66 ® 63
603.60,, E 2 603.60], 603, 60",
Where n and p refer to xx, yy and xy.
2 2 2
60'; 60'», 60'”,
52¢ 62¢ _ 62¢ 62¢ 62¢ 62¢ _ 0

 

 

603M603y = 60W60xx 60360” — 6030,60“ _

55

 

60'”60'Xy _ 6030,60”

Appendix B - Barlat Yld96

Many yield functions have been developed to describe more accurately
the plastic anisotropy of sheet metals. Barlat Yld96 [18] has proven to
describe well the plastic behavior of aluminum alloys. Although its convexity
cannot be proven and has complex derivatives, which make it a difﬁcult yield
function to implement,ation it describes the behavior of aluminum alloys very

accurately.

B1 - Yield function-2d

0

¢ = 23“ = a,|s2 —S,|" +a2|33 —S,|" +a3|S, —S2

 

Where 51, $2 and 53 are the eigen values of a stress tensor modiﬁed by a
linear operator to obtain anisotropic properties, c? is the measure of the yield
surface size, the coefﬁcients a], (12 and 03 are three functions of the angle [3
deﬁned below and the exponent a, connected to the crystal structure, is

equal to 6 or 8 for BCC and FCC alloys respectively. For plane stress the

eight material coefﬁcients needed to describe the anisotropy are

C1 ’C2 ’C3’C6.ax’ay ’a20,azl

56

The stress tensor modiﬁed by a linear transformation is:

Sn (C2 +C3)/3 —C3/3 0 0”
SW = —C3/3 (C1 +C3)/3 O 0'”,
SD. 0 0 C6 0'”.

The eigen values of the modiﬁed stress tensor are:

 

 

The coefﬁcients on. are found using the transformation rule:

a, = 0, cos2 6+0), sin2 ,6
a2 = a, sin2 ,B+ay cos2 ,6

a3 = 0:20 cos2 2,6 +az, sin2 2,6

Where [3 is the angle between the x direction (rolling direction) and the

direction associated with principal stress 5; and is found easily by:

ZS
2,6 = tan’l —”
S —S»'

I!

57

B2 - First derivative

6¢ _6¢ x,6Sx65p+6¢X6a,x6,6x
60,, —,6S 68p 60' 60. 6,6 6Sp 60,,

n I

 

2030-1—63 = 6¢
60' 60'

 

Where n and p denote xx, yy and xy and i for 1, 2 and 3.

66—35: —aa 23IS —S,aI sign(S3 —S,) )+aa 3I,—S —S ,a—I sign(S, -S 2)
a¢ a—I .

6S, —=aa,IS, —S ,aI s,ign(S —S) —aa,IS, -SZI Slgn(SI ’52)
g: = ——aa IS2 —S ,‘HI sign( (,S -S,) )+aa,IS3 —,HSI sign (S —S ,)

If A is not equal to zero the derivative below can be calculated:

A =(S,,,r -SW)2 +4Sfy

 

 

BileIJLSa.‘ 652-1(1312‘1) ass _IﬁiéiIﬂl

 

 

 

 

 

 

 

as” 2 J3 , as“ 2 J; ) as” as,“ as“

as, 1 S... —S,,) as2 1 S... -S,.‘ as f as as I
_:__1___ __1+____ 3 :_ 1+ 2 =__1
68,, 2 JZ ) 6S” 2 JZ , as” (as, as” ,

6S1 = 250 6S, ___ _Eﬂ 6S3 _ ( 65, 619,)— 0
a8, JX as, «E as” ’ — as”, + as”, , "

 

 

58

 

6
Sn =L
60'

P

 

Where L is the linear transformation matrix used to ﬁnd the modiﬁed stress

 

 

 

 

 

 

 

 

tensor.
(C, +C3)/3 —C,/3 0
L = —C,/3 (C, +C3)/3 0
0 0 C6
2g=IS,—S,a %-=(a,—ax)sin2,6 Biz—i1
60, 6,6 65” A
6¢ =IS3—S,a (702 :(ax—a))sin2ﬂ _a_ﬂ__S_n
a , (M as”, A
611’ =ISI ’52 0 6a} — 2(a I 'azo)sm4ﬂ aﬂ SIX —Svy
aa, 5,5 as = A
xy

B3 - Second derivative

62¢ 62¢ 65,- 68. a¢ 625.- 65. 62¢ 6a.- an 65.
-—-= X X 'I' X X + X X X
60,60, 6S ,60', 6S p 60,, 65, 6S p60, 60,, 60,60, 6,6 6S p 60','

2 2
+6¢X6a, xaﬂxasp+a¢xaaix 6’6 anp
60,. 6,660, 6S 60' 60. 6,3 6Sp60'q 60,,

p n I

 

 

p

2— 2 _ _.
6 0' : El-“ _1_ 6 ¢ _(a _1) 80' ® 60' {-7.072
60,60], 20 60,60}, 60" 60'

Indices n, p and q denote for xx, yy and xy and i for 1, 2 and 3.

59

 

Where:

a2 a2 as. as a2 60. a as
as,a:, zasgxasraap +6S,6¢:z,x a/i xaéxaap
azs, = a2s, xas,

aspaa, aspas, 60',

62¢ z 62¢ XBS.XBS.

605,60, 60,65, 65,, 60',

620, :aza,xa,ax68,,

6,660,, apz as, 60'
am = am xas,
aspao, aspas, 60

 

 

q q

 

q

q

oThe second derivative of function Ip with respect to S. is a symmetric

second order tensor (3-3). The elements of this tensor are written below.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

a2 .- a- 62 62 (P

as? ”(a-11015341 2.....5, —s, 2) ﬁg'ags‘mla-Uais. —s2 2

a2 .- a- 62 62 "-

6S? =a(a—1Xa,IS, —S3 2 +a3IS, -—S, 2) 6S 6‘: = 6S 6: = -0(a—1)a,IS3 ’SI 2
2 l 3 3 1

a2 ._ ._ 62¢ 62¢ 0'

as? =a(a-1)(al|52 ’53 2 +a,|s, -51 2) 6S26S3 = as3as, =’a(a_1)a‘lsi —S’ 2

60

. The second derivative of the stress vector 5. with respect to the stress
tensor Sn is a third order tensor (3-3-3). The non-zero elements of this

matrix are derived below.

 

 

 

azs,:_ 525, =_ 525, zﬂzl JZ-(S—nﬁ A
as f, asuasw as,y as,“ as; 2 JZ

azS1 __ aZSI _ 6251 __ 62SI __2 (Sn —Sy.v 0 A
asnasn, asyyasxy aswasn asxyas, JX

azs, _ _ 4S3,
.7511“ IV]

azs,=_ 525, :_ 525, :azS,=_l .fA—_i§£_—_§J A
as; asnasyy aswasn as;y 2 JZ
azs2 __ azs2 _ azs, __ azs2 _2 (Sn—5,, ,, A
asnasxy asyyas, asgasu asxyasyy JZ

als, _ 453,
as; "ZII‘F' JZ VA]

625, __ a2s, + azs2
as,as, as,as, as,as,
azs, _

as,as,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

61

. The second derivative of function (p with respect to a.- and S. are second

order tensor (3-3). The elements of this tensor are written below.

 

 

 

 

 

 

 

 

 

 

 

a2¢ a2¢ a2¢ ..-.
=0 =— = s -s ' s —s
6a,6S, aa,as2 aa,as, ”I 2 3 “gm 2 3)
62¢ a2¢ a2¢ .4
= =— :05 —S si nS —S
6a,6S, 6a,6S, aa,as, l 3 ' g ( 3 ‘)
62¢ a2¢ a2¢ ..
062305. aa,as, 60,65, al ‘ 2 S’g"( ‘ 2)

. The second derivative of function a,- with respect to B is a vector of three

elements derived below.

 

 

 

2
(2;; = 2(0), — ar )cos 26
620
6622 = 2(0I,I —a, )c0526
620
6623 = 8(az, - 0,0)cos 46

o The second derivative of the angle 6 with respect to 5p is a second order

tensor (3-3).

azﬂ azﬂ 2 Sxy(Sxx -80)

 

 

 

 

 

 

 

 

 

 

as; = as; A2

aZﬂ = azﬂ =—2 50(5)“ —Syy)
019305”, 6Syy6Sx, A2

62:6 _ _ azﬂ _ — (Sn —Sy.v)2 + 4Sxy
6511,65,), aSyyaSW A2

6% _ _ am _ -(S... -S....)’ +48.
65,965,“ 659619”, A2
626 = —8 Sxy(Sxx —S_Iy)
as; A2

62

Appendix C - Barlat Yld2000

This is a new yield function that has been proposed by F. Barlat [19] to
overcome the weakness of yield96 where it is proven to be convex and its
mathematical formulation is less complicated. This yield function was
developed for aluminum alloys and has been shown to describe well the

anisotropic behavior of aluminum sheets.

C1 - yield function-2d

a

¢=2a°=|X,—X,

 

“+px+nr+pn+x

 

Where X1, X2, Y1 and Y; are the eigenvalues of a stress tensor X modiﬁed by
a linear operator L to obtain anisotropic properties, the exponent a,
connected to the crystal structure, is equal to 6 or 8 for BCC and FCC alloys
respectively, 6 is the measure of the yield surface size. For plane stress the
eight material coefﬁcients needed to describe the anisotropy and the linear

operator L are written below.

ai’azravauas’asaawas

 

 

‘—2 2

L,,, 2/3 0 'L,,,1 8 —2 a,
L,,, —1/3 0 o'a, L,,, 1 1 —4 —4 4 0 a,
L,,, 0 -1/3 0 a, 1.2,, 5 4 —4 —4 1 0 a,
L,,, 0 2/3 0 _a, L,,, —2 8 2 —2 0 a,
L,66 0 0 1_ _L,,,_ _ 0 0 0 0 9__a,_

 

 

 

63

 

 

 

 

 

The linear operator L is found by a linear transformation form the anisotropy

coefﬁcients. The modiﬁed stress tensor X can be found by a linear

transformation with L:

 

 

”Xx. _
X}? 0
X”, n
X = Lu Y- = L 0,,
xx 0x ’
Y}? )
_ Y... _

 

 

 

 

 

 

 

X +X X +X 2 Y +Y Y +Y 2
X1: ”2 W+J£ ”2 W] “I'Xxyz 1: n2 yy-I-J[ ”2 W] +ny2
X +X X +X 2 Y +Y Y +Y 2
X,= ”2 W—I/I ”2 ”I +sz zzuzw—IIHZWI +13,2

C2 - First derivative

60 z a¢ x 6X.- an.
6X, 6Xp 60,,

 

 

60',
2670-1 60- : 6¢
60,, 60",

Where i denote indices for 1, 2 and n, p for xx, yy, xy.

64

 

ﬂzaIX, —X2
6X,

 

 

6X,

1

 

2.?— : 2aI2),2 + I],
61’

2

If Ax and Av are not equal to zero the derivative below can be calculated:

Ax = (Xxx —ny)2 +4X3)’

 

 

 

 

 

 

i=lK1+————X’“-X’W

aXn 2( JE

6X. =1,__X_n'_Xw

6X. 21 W37

_6_XL—2X”

aXU E

ayl 1[ Yu-Yw\
=_ 1+

61’,“ 2 A, )

ElileI—Y“—Y”\

61’” 2 Ay )

61’,_2ny

6Y 7.];

0431.3"(XI —X2)

065 =—aIX1‘X2Ia_15ig"(XI—X2)

I

\

 

 

% = 2al2Y. + Y,|"“sign(2Y. + Y. )+ aIZYz + Y. 1051318425 + Y1)

“"‘sign(2Y, + Y. )+ aIZYI + BIO-131842"! + Y2)

AY : (Yr: —YY)')2 +412:

6X2 1[1_Xxx—XYY\
,IA,r )

\
8X2 1[1+_/Y_X_";X:£

 

 

 

 

any-Z J3 /
6X, =_2X,,
5X... J4:

 

 

 

 

 

 

 

airy, - 2 JAT ,
6Y, =_ 2Y0
er. i4?

65

C3 - Second derivative

 

 

 

 

 

 

62¢ _ 62¢ anianq+a¢x 62X.- anq
60,60, 6X,.60'p 6X4 60'" 6X, 6Xq60'p 60”,,
2— 2 _. _
60' :E,_,, _1_ 6¢ _(a_1)60'®60' 6'“
60,60], 20 603,60}, 60",, 60,,
Where:
62¢ 62¢ 6X; 6X,

 

 

: X X
6X,60'p 6X,6Xj 6X" 60“,,

62X, ex,

= X
axqaap max, 60’

62X.

1

P

oThe second derivative of function (p with respect to X. is a symmetric

second order tensor (4-4). The non- zero elements are written below.

 

 

 

 

2 2
6 (1:: 6 f =a(a—1]X, —X2 “.2
6X, ax,
2 2
6¢ = 6¢ =—a(a—1)X,—X2 “-2
aX,aX, aX,aX,

 

 

 

 

 

 

 

2
(if: = a(a -1)(4|2Y, + Y2 “‘2 +|2Y2 + Y, ““2)
l
2
:1}: = a(a -—1)(4|2Y2 + Y, “‘2 +|2Y, + Y2 “‘2)
2
62¢ 62¢

 

“—2)

= 2a(a —1)(sz, + Y2 “‘2 +|2Y2 + Y,

 

 

aY,aY2 ___ 6Y26Y,

66

 

. The second derivative of the stress vector X. with respect to Xn is a third

order tensor (4-6-6). The non-zero elements are written below.

 

 

 

62X,__ 62X, __ 62X, _aX-X,__1_ A _(X, —X,,
an, aXUaX, aXWaX, an, 2 X ,/A,
_6_X.___ XX. z XX: ___X_X____2 (X X )X
aXuaX, aX,,aX,, aXWaX, aXWaXW ,/A, AX
:2 A, ———,_A. A,
6X; [ AX
52X, 62X, 62X, 62X, 1 ,— (X,,,-—X,,.)2
———3—:——:— = 2 :-—— AX— AX
6X, aXuaX, aXWaX, 6X”, 2 ,/A,
.122--sz z 222, ___a_£__2 (X,_X,,.)x, 2
aXnaX, anyaX, 6XU6XU 6XW6XW A, X
62X, 4X3,
.—.—2 ,/A ———- A
6X; [l 2 FAX ]/ J
2 2 2 2 _ 2
a);,:_ aY, =_ 6Y2 =a§=l ,[AT—(YX‘ Y“) A,
aY, aYnaY, anaY, aYW 2 ,/A,

62),] :— 62),] = 62),] :_ 62"! :_2 (YXJ —Y)'.V)yxy A
aYuaY, aYyyaY, aYngH aYnan A, ’

é

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

62y,:_ azy, _ azy2 __22_6Y «(Elk—Y”?

aY; aY,aY,, aY, aY, "'26ng 2
622,2 _— 622,2 _ aZYz __ azYz _2 (YU_Y)O')YX)’ A
6Y,,6Y,y anaY, aYUaY, aYnaY, ,/A, ’

: —2 A), — —4— A}:
an: [ FA,

 

 

 

 

 

67

 

 

References

[1] S. P. Keeler and W. A. Backofen, "Plastic Instability and Fracture in
Sheets Stretched Over Rigid Punches" Trans ASM 56, pp. 25—48. 1963

[2] S. P. Keeler and W. A. Backofen, "Determination of Forming Limits in
Automotive Stampings" Sheet Metal Indust. 42, pp. 683—691. 1965

[3] Z. Marciniak, and K. Kuczynski , "Limit Strains in the Process of
Stretch-Forming Sheet Metal" Int. J. Mech. Sci. 9, pp. 609-620. 1967

[4] William F. Hosford and Robert M. Caddell, “Metal forming: mechanics
and metallurgy" Englewood Cliffs, NJ. : Prentice Hall, c1993. ISBN-
0135885264

[5] A. Graf, WF. Hosford, “Calculations of Forming Limit Diagrams”,
Metall Trans A 21 (1): 87-94 Jan 1990,

[6] A. Graf, WF. Hosford, “The Effect of R-value on Calculated Forming
Limit Curves”, Forming limit diagrams: Concepts, Methods, and
applications: 153-163 TMS 1989 ISBN-0873390989

[7] F. Barlat “Crystallographic Texture, Anisotropic Yield Surfaces and
Forming Limits of sheet metals” Mater. Sci. Eng., 91, pp. 55-72. 1987

[8] A. Barata da Rocha, F. Barlat, J. M. Jalinier “Prediction of the Forming
Limit Diagrams of Anisotropic Sheets in Linear and Non-Linear
Loading" Mater Sci Eng 68 (2): 151-164 1985

[9] A. Graf, w. HOSFORD “Effect of Changing Strain Paths on Forming
Limit Diagrams of Al 2008-T4 Metall Trans A 24 (11): 2503-2512 Nov
1993

[10] A. Graf, w. HOSFORD “Calculations of Forming Limit diagrams for
Changing strain Paths" Metall Tans A 24 (11): 2497-2501 Nov 1993

[11] A. Graf, W. HOSFORD “The Inﬂuence of Strain-Path Changes on
Forming Limit Diagrams of A1 6111 T-4" Int J Mech Sci 36 (10): 897-
910 Oct 1994

[12] MC Butuc, A. Barata da Rocha, J.J. Gracio and J. Ferrereira Duarte, “A
more general model for forming limit diagrams prediction", J Mater
Process Tech 125: 213-218 Sp. Iss. SI Sep 9 2002.

[13] MG. Butuc, J.J. Gracio and A. Barata da Rocha, “A theoretical study

on forming limit diagrams prediction", J Mater Process Tech 142 (3):
714-724 Dec 10 2003.

68

 

V'-
IL"?

 

[14] F. Barlat “Forming Limit Diagrams — Predictions Based on Some
Microstructural Aspects of Materials" Forming limit diagrams: Concepts,
Methods, and applications: 275-301 TMS 1989 ISBN-0873390989

[15] MM. Moshksar, S. Mansorzadeh “Determination of the Forming
Limit Diagram for Al 3105 Sheet", J Mater Process Tech 141 (1): 138-
142 Oct 1 2003

[16] Z. Marciniak, K. Kuczynski and T. Pokora "Influence of Plastic
Properties of a Material on Forming Limit Diagram for Sheet-Metal in
Tension " Int. J. Mech. Sci. 15 (10), pp. 780-789. 1973

[17] K. S. Chan “Marciniak-Kuczynski Approach To Calculating forming
Limit Diagrams" Forming limit diagrams: Concepts, Methods, and
applications: 73-110 TMS 1989 ISBN-0873390989

[18] F. Barlat, Y. Maeda, K. Chung, M. Yanagawa, J. C. Brem, Y. Hayashida,
D. J. Lege, K. Matsui, S. J. Murtha, S. Hattori, R. C. Becker and S. Makosey
“Yield Function Development for Aluminum Alloy Sheets" J. Mech.
Phys. Solids 45, pp. 1727-1763 1997.

[19] F. Barlat, J. C. Brem, J.W. Yoon, K. Chung, R.E. Dick, D. J. Lege, F.
Pourboghrat, S.-H. Choi, E. Chu “Plane Stress Yield Function for
Aluminum Alloy Sheets—part 1: Theory" Int. J. of Plast. 19, pp. 1297-
1319. 2003.

69

PL f—ITiﬂi

h.-.

 

 

 

   

llllllllllllllllll