Linköping Studies in Science and Technology

Thesis No. 1529

Modelling of Failure
in High Strength Steel Sheets

Oscar Björklund

LIU–TEK–LIC–2012:14

Department of Management and Engineering, Division of Solid Mechanics

Linköping University,
SE–581 83, Linköping, Sweden
http://www.solid.iei.liu.se/

Linköping, April 2012

Cover:
Results from finite element simulation of a forming limit test.

Printed by:
LiU-Tryck, Linköping, Sweden, 2012
ISBN 978–91–7519–895–5
ISSN 0280–7971

Distributed by:
Linköping University
Department of Management and Engineering
SE–581 83, Linköping, Sweden

c 2012 Oscar Björklund

This document was prepared with LATEX, April 30, 2012

No part of this publication may be reproduced, stored in a retrieval system, or be

transmitted, in any form or by any means, electronic, mechanical, photocopying,
recording, or otherwise, without prior permission of the author.
Preface

The work presented in this thesis has been carried out at the Division of Solid
Mechanics, Linköping University with financial support from the VINOVA PFF
project ”Fail” and the SFS ProViking project ”SuperLight Steel Structures”. The
industrial partners Dynamore Nordic, Outokumpu Stainless, Saab Automobile,
Scania, SSAB, SvereaIVF and Volvo Car Corporation are also gratefully acknowl-
edged for their support.
I would like to thank my supervisor, Professor Larsgunnar Nilsson for all his sup-
port and guidance during course of this work. I also greatly appreciate all my
colleagues at the Division of Solid Mechanics for their help, in particular Rikard
Larsson, Lic. Eng., for all his support and guidance concerning constitutive mod-
elling.
For their assistance during the mechanical testing throughout this project I would
like to thank Andreas Lundstedt at Outokumpu Stainless and Bo Skoog, Ulf
Bengtsson and Sören Hoff at Linköping University.
I am also grateful to all my friends and my family for their support over the course
of all these years. I could not have done it without you.

Oscar Björklund
Linköping, April 2012

”Results! Why, man, I have gotten a lot of results. I know several thousand things
that won’t work”
Thomas Edison

iii
Abstract

In this theses the high strength steel Docol 600DP and the ultra high strength steel
Docol 1200M are studied. Constitutive laws and failure models are calibrated and
verified by the use of experiments and numerical simulations. For the constitu-
tive equations, an eight parameter high exponent yield surface has been adopted,
representing the anisotropic behaviour, and a mixed isotropic-kinematic hardening
has been used to capture non-linear strain paths.
For ductile sheet metals three different failure phenomena have been observed: (i)
ductile fracture, (ii) shear fracture, and (iii) instability with localised necking. The
models for describing the different failure types have been chosen with an attempt
to use just a few tests in addition to these used for the constitutive model. In this
work the ductile and shear fracture have been prescribed by models presented by
Cockroft-Latham and Bressan-Williams, respectively. The instability phenomenon
is described by the constitutive law and the finite element models. The results
obtained are in general in good agreement with test results.
The thesis is divided into two main parts. The background, theoretical framework,
mechanical experiments and finite element models are presented in the first part.
In the second part, two papers are appended.

v
List of Papers

In this thesis, the following papers have been included:

I. R. Larsson, O. Björklund, L. Nilsson, K. Simonsson (2011). A study of
high strength steels undergoing non-linear strain paths - Experiments and
modelling, Journal of Materials Processing Technology, Volume 211, pp. 122-
132.
II. O. Björklund, R. Larsson, L. Nilsson (2012). Failure of high strength steel
sheets - Experiments and modelling, Submitted.

Own contribution

In the first paper, Rikard Larsson and I jointly performed the modelling and ex-
perimental work. However, Rikard Larsson was responsible for writing the paper.
In the second paper, the modelling and experimental work were once again per-
formed by Rikard Larsson and myself, but I was responsible for evaluating the
failure phenomena and for writing the paper.

vii
Contents

Preface iii

Abstract v

List of Papers vii

Contents ix

Part I Theory and Background

1 Introduction 3

2 Steels 5
2.1 Crystal Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Dual Phase Steel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Martensitic Steel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Deformation and Fracture 9

4 Constitutive Modelling 11
4.1 Isotropic Hardening . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.2 Kinematic Hardening . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.3 Mixed Hardening . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.4 Effective Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5 Fracture Modelling 19
5.1 Ductile Fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2 Shear Fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.3 Phenomenological Fracture Models . . . . . . . . . . . . . . . . . . 23

6 Modelling Instability 29
6.1 Instability in Plane Strain . . . . . . . . . . . . . . . . . . . . . . . 29
6.2 Analytical Instability Models . . . . . . . . . . . . . . . . . . . . . 31
6.3 The Marciniak and Kuczynski Model . . . . . . . . . . . . . . . . . 33
6.4 Finite Element Model . . . . . . . . . . . . . . . . . . . . . . . . . . 35

ix
 6.5 Evaluation of Instability . . . . . . . . . . . . . . . . . . . . . . . . 36

7 Mechanical Experiments 39
7.1 Pre-Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
7.2 Tensile Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.3 Simple Shear Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.4 Plane Strain Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.5 Balanced Biaxial Test . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.6 Nakajima Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

8 Finite Element Modelling 45

8.1 Pre-Straining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
8.2 Tensile Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
8.3 Shear Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
8.4 Plane Strain Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
8.5 Nakajima Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

9 Conclusions and Discussion 51

10 Review of Appended Papers 53

Bibliography 55

Part II Appended Papers

Paper I: A study of high strength steels undergoing non-linear strain
paths - Experiments and modelling . . . . . . . . . . . . . . . . . . 63
Paper II: Failure of high strength steel sheets - Experiments and modelling 77

x
 Part I

Theory and Background

Introduction
1

Due to the desire from industries to shorten product development lead time, the
use of Simulation Based Design (SBD) is rapidly increasing. The development of
finite element (FE) methods and the rapid growth of computational power have
made it possible for SBD to advance from research to industrial utilisation. In
addition to reduced development time and cost, the possibility of determining the
product’s properties at an early stages of the design process is of utmost impor-
tance e.g. to find out if a car is safe in a crash situation. Simultaneously with
the development of SBD, material suppliers have developed more advanced steels
with improved mechanical properties. However, these advanced steels often show
anisotropic behaviour and more complex hardening compared to traditional low-
carbon steels. Consequently, the demand for research on material models, including
failure models suitable for these steels, has increased.
The constitutive model that has been used throughout this work includes an eight
parameter high exponent yield surface with mixed isotropic-kinematic hardening.
The model has the ability to represent the anisotropic behaviour that is shown in
mechanical tests, even for a complex loading path. To further increase the ability
to predict the product’s potential and thereby enable a more optimised product,
the necessity of improving failure models has arisen. At present, progressive fail-
ure is not included in automotive FE simulations. In order to predict the likeli-
hood of failure, the predicted strain states are compared to experimental forming
limit curves (FLC). However, the experimental FLC is constructed for linear strain
paths, although its locus has been shown to depend on the strain path. Therefore,
the prediction of failure of a component, e.g. a car body part subjected to a crash
simulation, lacks information from its previous design history. Thus, an improved
phenomenological model for prediction of failure is desired. A phenomenological
failure model must have the ability to inherit the properties of a component from
its prior manufacturing process.
Macroscopic fractures have always been of great interest. As early as the beginning
of 16th century, Leonardo da Vinci explained fractures in terms of mechanical
variables. He established that the load an iron wire could carry depended on the
length of the wire, as a consequence of the amount of voids in the material. The
longer the wire, the more voids, which led to lower load-carrying capacity. Even
if material fracture has been studied for a long time, the underlying microscopical
fracture mechanisms are hard to translate to a phenomenological model. Since the

3
CHAPTER 1. INTRODUCTION

underlying mechanisms need to be represented on a length scale, which can be used

in an FE simulation, it would be computationally expensive, and hence too time
consuming, to represent the fracture on a micromechanical scale. Throughout this
work, the expressions failure and fracture are used extensively:
• Failure - loss of load-carrying capacity of a component or a structure.
• Fracture - material separation into two, or more, pieces.
The term failure incorporates the term fracture but may also be other structural
phenomena which do not include a material separation, e.g. material and geometri-
cal instabilities. In ductile sheet metals three different failure phenomena have been
observed: (i) ductile fracture, (ii) shear fracture and (iii) material instability with
localisation. The third term instability with localisation is also sometimes denoted
as purely-plastic failure. In this work the different phenomena have been modelled
separately and the focus has been placed on using failure models which can easily
be calibrated from simple mechanical tests. For modelling of the ductile and shear
fractures the Cockroft-Latham and Bressan-Williams models, respectively, have
been used. Each of their models contains only one material parameter and are
easy to calibrate. The instability phenomenon is described by the FE model and
the constitutive law and needs no further mechanical experimental calibration.
Even if the objective of this project is to identify improved material and failure
models suited for high strength sheet metals, the benefit in a longer perspective
may be the opportunity to produce lighter and safer structures. The environmental
profit of producing lighter products is obvious for the automotive industry, but also
other industrial applications can profit from reducing weight i.e. lighter containers
will result in an ability to increase the pay-load. The safety of occupants in cars
has been an increasingly important marketing issue for the automotive industry.
By using new materials to the very edge a more optimised product can be obtained
and thereby many fatalities may be prevented.

4
Steels
2
The classification of different steel grades is often made by the amount and type
of alloying substances. However, it is also common to talk about different groups
of steels, e.g. High-Strength Low Alloy (HSLA) steel, High-Strength Steel (HSS),
Ultra High-Strength Steel (UHSS), Advance High-Strength Steel (AHSS), Dual
Phase (DP) Steel, and Transformation-Induced Plasticity (TRIP) Steel. A specific
steel can be a part of more than one group as there is no unequivocal division.
However, in this work the materials studied are categorised into the two groups
HSS and UHSS. These groups are grouped according to the yield strength, and
steels with a yield strength between 210 to 550 MPa are classified as HSS and
steels with a yield strength over 550 MPa are classified as UHSS, see Opbroek
(2009).
The steels covered in this study are Docol 600DP and Docol 1200M. The name Do-
col is a SSAB trademark. Docol 600DP is an HSS with a dual face structure consist-
ing of about 75% ferrite and 25% martensite, where the two-phase microstructure
is produced by heat treatment. Docol 1200M is an UHSS with a fully martensitic
steel produced by water quenching from an elevated temperature in the austenitic
range, see Olsson et al. (2006). The nominal thicknesses of the steel sheets studied
were 1.48 mm and 1.46 mm for Docol 600DP and Docol 1200M, respectively. For
details on the chemical composition of the Docol 600DP and Docol 1200M, see
Table 1.

Table 1: Chemical composition of Docol 600DP and Docol 1200M.

C Si Mn P S Nb Al Fe
[%] [%] [%] [%] [%] [%] [%] [%]
Docol 600DP 0.10 0.20 1.5 0.01 0.002 - 0.04 balance
Docol 1200M 0.11 0.2 1.6 0.015 0.002 0.015 0.04 balance

5
CHAPTER 2. STEELS

2.1 Crystal Structure

The crystal structure of the atoms in a metal is composed of regular three dimen-
sional patterns called a crystal lattice. There are a number of different types of
crystal lattices and it is common to talk about unit cells of different types. For
metals, the most common types of unit cells are body-centered cubic (BCC), face-
centered cubic (FCC), hexagonal close-packed (HCP), and body-centered tetrago-
 $ 
nal (BCT), see Figure 1. For a more complete overview of different types of unit
cells see e.g. Askeland (1998). For pure iron there are three different types of
1 ,
phases called α, γ and δ. The α-iron,  known
also  as ferrite, is present at temper-
atures up to 912◦ C and has a BCC phase structure. Between 912◦ C and 1394◦ C
the pure iron is known as γ-iron, or austenite, which has an FCC phase structure.
Above 1394◦ C and up to the melting point at 1538◦ C the iron is known as δ-iron,
or delta ferrite, which also has a BCC phase structure. By adding carbon or other
alloying elements the range of the different phases will be altered, for more detailed
information see e.g. Krauss (2005). By controlled cooling of the metal, the phases
can be prevented from changing, thereby a phase which is present under high tem-
peratures can be obtained at room temperature. However, also  new
 '((phases such
as martensite can be obtained by rapidly cooling austenite.

(a) BCC (b) FCC (c) HCP (d) BCT

'(

Figure 1: Different types of unit cells common for metals, (a) body-centered cubic,
(b) face-centered cubic, (c) hexagonal close-packed and (d) body-centered tetragonal.

2.2 Dual Phase Steel

  )((
Dual phase steels mainly consist of ferrite and a minor part of martensitic second
phase particles. The second phase particles can be seen as channels and darker
areas between the ferritic islands in Figure 2(a). An increase of the amount of hard
martensitic particles generally increases strength and reduces ductility. The dual
phase is produced either by controlled cooling from the austenitic phase or from
the ferrite-austenite phase. The carbon (C) enables the formation of martensite at
practical cooling rates. By adding other alloying element such as manganese (Mn),

6 #
 2.3. MARTENSITIC STEEL

chromium (Cr), molybdenum (Mo), vanadium (V) and nickel (Ni) individually or
in combination, an increase in the hardenability of the steel can be obtained.

2.3 Martensitic Steel

Martensitic steels are produced by a quenching from a temperature in the austenitic
range. The steels are characterised by a martensitic matrix with small parts of fer-
rite or bainite. The martensite can be seen as small needle shaped particles in
Figure 2(b). The hardenability can be increased by adding carbon (C) or other al-
loying elements such as manganese (Mn), silicon (Si), chromium (Cr), molybdenum
(Mo), boron (B), vanadium (V) and nickel (Ni).

(a)

(b)

Figure 2: Microstructure for: (a) Docol 600DP and (b) Docol 1200M.

7
Deformation and Fracture
3
The deformation of steel can be divided into elastic and inelastic parts. The elastic
(reversible) deformation, which occurs at the atomic level, does not cause any per-
manent deformation of the material. Inelastic or plastic deformation, on the other
hand, causes permanent deformation. This deformation occurs within the crystal
structure by a procedure known as slip. The slip produces crystal defects, called
dislocations, and when the deformation continues a large number of these defects
occurs. The dislocations in the crystal lattice cause obstacles, which prevent the
mobility of subsequent dislocations and thereby contribute to a hardening mech-
anism. Elastic or plastic deformations do not destroy the atomic arrangement of
the material while fractures, on the other hand, cause discontinuities within the
material. These discontinuities cause stress concentrations which will increase the
’rate’ of fracturing. There are two main types of fracture: brittle and ductile.
Brittle fracture is the breaking of interatomic bonds without noticeable plastic de-
formation. These fractures occur when the local strain energy becomes larger than
the energy necessary to pull the atom layers apart. Brittle fracture occurs mainly
in high-strength metals with poor ductility and toughness. However, even metals
that have normal ductility may exhibit brittle fracture at low temperatures, in
thick sections or at high strain-rates. The surface of a brittle fracture is charac-
terised by its flat appearance and the fracture surface is most often perpendicular
to the applied load.
Ductile fracture, on the other hand, is caused by an instability which is the result
of very extensive plastic deformations occurring in the surroundings of crystalline
defects. The global deformation in a ductile fracture may be either large or small,
depending on the density of the defects. The ductile fracture is described as the
initiation, growth and coalescence of voids in the material, and finally the loaded
area has been reduced and the material fails. The surface of the ductile fracture is
characterised by the presence of shear lips, which results in the form of a cup and
a cone on the two fracture surfaces. It is also often possible to see the dimples that
are caused by the micro-voids. The material examined in this study is assumed
to fail only in a ductile manner, and consequently the brittle fracture type is not
considered.

9
CHAPTER 3. DEFORMATION AND FRACTURE

In ductile sheet metals or thin walled structures, failure can be caused by one
or a combination of the following; (i) ductile fracture, (ii) shear fracture or (iii)
instability with localisation, see Figure 3.

) * +   ! 

4(a)
# $    "  (b)

 

     

(c)

Figure 3: Failure types in sheet metals, (a) ductile fracture, (b) shear fracture, (c)
instability with localisation.

The ductile fracture arises since mostmaterials

    contain inclusions or particles, and
during deformation voids are likely to nucleate at these locations. Further plastic
deformation will cause them to grow until they link together to produce a ductile
fracture. Hydrostatic tension will increase the volume of the voids and thereby
favours this type of fracture.
Shear fracture can be caused either 
by extensive
!    slip on the activated slip planes,
see Dieter (1986), or as a result of void nucleation in the slip bands. Both of these
mechanisms are favoured by shear stresses. When voids nucleate in the slip bands,
the loaded area is reduced such that plastic flow localises there. Continued shear
increases the area of the voids until separation occurs. Furthermore, as stated by
Teirlinck et al. (1988) ”Voids which extend in shear need not increase in volume, so
shear fracture is less pressure-dependent than ductile fracture, though it remains
more pressure-dependent than purely-plastic failure”.
Instability or plastic failure occurs if no other mechanism intervenes. In the ordi-
nary tensile test of a sheet metal sample two types of instability can occur: diffuse
necking and localised necking. Both mechanisms occur when the strain hardening
can no longer compensate for the increased load. Diffuse necking in a tensile test
occurs as a reduction of the waist over a length which has about the same size
as the width of the waist. Localised necking takes place inside the diffuse zone
as a thickness reduction along a narrow band which has about the same size as
the thickness. In a more general loading case, localisation can occur without a
preceding diffuse necking. Since the strain localises inside the neck the material
fails there due to ductile or shear fracture.

10
Constitutive Modelling
4
In order to be able to formulate the hypothesis on which the constitutive models
are based, it is necessary to understand the basic physical deformation mechanisms
of the material as described in Chapter 3. Since applications in sheet metal forming
and vehicle collision involve large deformations and rotations of the material, it is
necessary to formulate the constitutive equations in a consistent manner. By using
a co-rotational material framework, see e.g. Belytschko et al. (2000), these large
deformations and rotations can be accounted for. The co-rotated Cauchy stress
tensor can the be expressed as

σ̂ = RT σR (1)

where R is the orthogonal transformation tensor and σ̂ denotes the co-rotated

Cauchy stress tensor, σ. From now on (ˆ·) denotes the co-rotated quantity of
(·). The co-rotated rate of deformation tensor is assumed to follow an additive
decomposition of the elastic and plastic parts

e p
D̂ = D̂ + D̂ (2)

where superscript e and p denotes elastic and plastic parts, respectively. In cases
of small elastic deformations, a hypo-elastic stress update can be assumed, i.e.

˙ = Ĉ : (D̂ − D̂ p )
σ̂ (3)

where Ĉ is the co-rotated fourth order elastic stiffness tensor. Henceforth, the
co-rotated superscript, (ˆ·), will be excluded and all subsequent relations will be
related to co-rotated configuration. Plastic deformation will not occur as long as
the stress state is considered as elastic. Whether or not the stress state is in the
elastic regime is determine by the yield function, f
y
f = σ̄ − σiso (4)

11
CHAPTER 4. CONSTITUTIVE MODELLING

y
where σ̄ is the effective stress and σiso is the current yield stress. The yield function
determines the elastic limit of the material and the hypersurface, when f = 0 is
denoted as the yield surface. As long as f < 0 no plastic deformations will occur
and when f = 0 and f˙ = 0 plastic flow will occur. The plastic part of the rate of
deformation tensor will be determined from the plastic potential, g, according to

∂g
D p = λ̇ (5)
∂σ

where λ̇ is the plastic multiplier and the gradient of g determines in which direction
the material will flow. In an associated, in contrast to a non-associated, flow
rule the flow is directed normal to the yield surface and in that case the plastic
potential and the yield function coincide, i.e. g = f . This, assumption has been
used throughout this work. Three different modes of loading can occur, (i) elastic
case then λ̇ = 0 and f < 0, (ii) plastic loading λ̇ ≥ 0 and f = 0 , and (iii)
neutral loading if the direction of loading is tangential to the yield surface λ̇ = 0
and f = 0. These three cases may be expressed using the Karush-Kuhn-Tucker
conditions (KKT-conditions).

λ̇ ≥ 0, f ≤ 0, λ̇f = 0 (6)

4.1 Isotropic Hardening

As discussed earlier, plastic deformation is caused by dislocations in the material
and the growth of these dislocations will prevent the mobility of new ones and
thereby increase the stress that is necessary to produce an increased plastic strain.
y
Most models for describing the current yield stress σiso are given as functions of
p
the equivalent plastic strain ε̄ defined as

Zt
p
ε̄ = ε̄˙p dt (7)
0

where t is the time. However, other quantities, e.g. the plastic strain-rate and
temperature, may also affect the hardening of the material. A large amount of
analytic functions for describing the hardening can be found in literature. It has
also been quite common to use the experimental hardening curves directly. The
tensile test is often used for describing hardening up to diffuse necking after which
some kind of extrapolation is needed. In this work, an extended Voce law, see Voce
(1948), was fitted to the hardening data from the tensile test up to diffuse necking.

12
 4.1. ISOTROPIC HARDENING
  

After diffused necking an extrapolation using the Hollomon law (see Hollomon
(1945)) i.e.
 a
power-law, was adopted by the use of inverse modelling of a shear
test. The analytic hardening
σ2 function can be expressed σas
2
"  

 2
 !
 X
 σ + p
QRi (1 − e−CRi ε̄ ) ε̄p ≤ εt
y p 0
σiso (ε̄ ) = (8)


i=1
σ1 σ1
A + B(ε̄p )C ε̄p > εt
"  
where σ0 , Q

Ri , CRi , A, B and C are material constants and εt is the transition
 
   !
 !the extended
point between
1
Voce and Hollomon  !hardening. The requirements of
a smooth curve, i.e. a C transition between the Voce and Hollomon expressions,
makes it possible to determine the constants A, B and C with only one new pa-
  
0
rameter denoted σ100 , which can be interpreted as the stress at 100% plastic strain.

2
X t
A + B(εt )C = σ0 + QRi (1 − e−CRi ε )
i=1
2
X (9)
t
CB(εt )C−1 = CRi QRi e−CRi ε
i=1
0
A + B = σ100

The isotropic hardening can be illustrated as an expansion of the yield surface in

all directions, see Figure 4(a).     

σ1 σ1 σ1

σ2 σ2 σ2

(a) (b) (c)

    
Figure 4: Different types of hardening (a) isotropic, (b) kinematic and (c) mixed.

13
CHAPTER 4. CONSTITUTIVE MODELLING
,  - %

4.2 Kinematic Hardening

Plastic deformation
Σ = σ −mayα also introduce anisotropic behaviour of the material, & e.g.
as manifested by the Bauschinger effect, which is the phenomenon of a lower yield
stress in the case of reverse loading. For modelling this effect it is useful to introduce
a kinematic hardening model. The kinematic hardening enables the yield surface
2 2  
to move in the stress space, see Figure
Σ 4(b), and thereby obtain a reduced yield
α̇ =
stress in the case of
α̇i =
reverse
CXi QXi − αi
loading, see
σ̄ Figure 5.
&
i=1 i=1

σ σRD
 

 

ε σT D   

Σ=σ−α *
Figure 5: Ilustraition of the Bauschinger effect in the case of reverse loading.
   ! 
The motion 2of the 2yield surface Σ is described by the backstress tensor, α. The
α̇ = α̇i =
related overstress
CXi QXi − αi
tensor,
-
Σ = σσ̄ − α, replaces the Cauchy stress tensor, σ, in the
i=1 i=1
effective stress function, σ̄ = σ̄(σ − α). An interpretation can be seen in Figure 6.

σT D

σ
Σ

α σRD

Figure 6: Interpretation of the overstress, Σ,%backstress, α, and Cauchy stress

tensor, σ.

In this work the evolution of

thebackstress has been described by a two-component
  
law, presented by Frederick and Armstrong (2007)

14

.
 4.3. MIXED HARDENING

2
X 2
X  
Σ
α̇ = α̇i = CXi QXi − αi (10)
i=1 i=1
σ̄

where QXi and CXi are material constants determined from pre-strained tensile
tests.

4.3 Mixed Hardening

Throughout this work a combination of isotropic and kinematic hardening has
been used, i.e. the yield surface can both expand and move, see Figure 4(c). In the
special case of uniaxial tension in the φ-direction the only non-zero components
in the overstress, backstress, and Cauchy stress tensors will be Σφ , αφ , and σφ ,
respectively. The evolution of the backstress can then be expressed as

2
X p
p
αφ (ε̄ ) = rφ QXi (1 − e−CXi ε̄ ) (11)
i=1

where rφ = Σyφ /σ̄ is constant and αφ (0) = 0 has been used. The yield stress, σφy ,
can now be expressed as a function of the equivalent plastic strain, ε̄p .

 2

 X  p p 

 σ + QRi (1 − e−CRi ε̄ ) + QXi (1 − e−CXi ε̄ ) ε̄p ≤ εt
 0

σφy (ε̄p ) = rφ i=1

2 (12)

 X

 A + B(ε̄ p C
) +
p
QXi (1 − e−CXi ε̄ ) p
ε̄ > ε t

i=1

As can be seen from Equation 12, the mixed hardening may, in the special case
of monotonic loading, be interpreted as a summation of isotropic and kinematic
hardening, see Figure 7. The total stress at 100% plastic strain, σ100 , can now be
0
expressed as the function of the isotropic stress at 100% plastic strain, σ100 , and
the saturation of the kinematic hardening.

2
X
0
σ100 = σ100 + QXi (13)
i=1

15
CHAPTER 4. CONSTITUTIVE MODELLING

1000
σ100

800 ′
σ100
Stress, σ [MPa]

600

400

200
σyy
t σiso
ε
α
0
0 0.2 0.4 0.6 0.8 1
Equivalent plastic strain, ε̄p [-]

Figure 7: Hardening decomposition in the case of monotonic loading.

In order to simplify the parameter identification procedure, a linear mixture be-

tween the isotropic and kinematic hardening has been assumed according to

QRi = (1 − βi )Qi , QXi = βi Qi , i = 1, 2

(14)
CRi = CXi = Ci , i = 1, 2

where βi determines the amount of kinematic hardening, and is restricted to a

value between 0 and 0.9 in order to avoid infeasible values of Ci .

4.4 Effective Stress

Many expressions for effective stress have been presented over the years. The
expression for the effective stress according to von Mises is well known and often
used. However, since the von Mises expression only shows isotropic behaviour, and
since the rolling process used for producing sheet metals usually leads to different
properties in different directions, an anisotropic effective stress function is needed.
The rolling direction (RD), transversal direction (TD) and normal direction (ND)
are often used to describe the axes of orthotropy. Anisotropy can occur both in
the yield stress and in the plastic flow. The anisotropy in plastic flow is described
by the Lankford parameter, R, see e.g. Hosford and Cadell (1993), or the plastic
strain ratios, k, defined as

16
 4.4. EFFECTIVE STRESS

dεpT dεpT −Rφ

Rφ = p , kφ = = (15)
dεN dεpL Rφ + 1

where the indices on the logarithmic plastic strain increment indicate the transver-
sal, T , normal N , and longitudinal, L, direction, respectively. The anisotropy in
yield stress is described by

σφy
rφ = (16)
σref

where σφy is the yield stress in the φ direction and σref is a reference yield stress.
 
Equations 15 and 16 are expressed for the tensile test case but similar equations
can be obtained for the balanced biaxial test

dεpT D σby
Rb = kb = , r b = (17)
dεpRD σref

where the subscript b denotes balance biaxial test. The interpretation of the plastic
strain ratios, k, and the yield stress ratios, r, are shown in Figure 8.

   σT D
 
     1
    
  1

kb
 (rb σref , rb σref )
k90
(0, r90 σref )


  

(r00 σref , 0) σRD

k00
1


Figure 8: Yield locus for a plane stress case with τRD = 0.

  
17


CHAPTER 4. CONSTITUTIVE MODELLING

Expressions for representing anisotropic behaviour have been presented by numer-

ous authors and a good overview is given in Banabic et al. (2010). In the present
study the model presented by Aretz (2005) has been used

 1/a
1
σ̄(Σ) = (|∆01 |a + |∆02 |a + |∆001 − ∆002 |a ) (18)
2

 s 2
∆01 A8 Σ∗11 + A1 Σ∗22 A2 Σ∗11 − A3 Σ∗22
= ± + A24 Σ12 Σ21 (19)
∆02 2 2

 s 2
∆001 Σ∗ + Σ∗22 A5 Σ∗11 − A6 Σ∗22
= 11 ± + A27 Σ12 Σ21 (20)
∆002 2 2

where A1 , ..., A8 and a are material constants. The model was originally derived for
a plane stress case, where the 11-direction and 22-direction correspond to the RD
and TD, respectively. However, a regularisation that enables a C 0 continuous thick-
ness across the element edges has been carried out. Thus the through-thickness
normal stress has been included in the effective stress measures, Σ∗11 = Σ11 − Σ33
and Σ∗22 = Σ22 − Σ33 , where the 33-direction corresponds to the ND of the sheet.
The identification of the parameters has been implemented both by direct methods,
using the strain and stress ratios, and inverse modelling of the shear test.

18
Fracture Modelling
5
Fracture can generally be divided into ductile and brittle fracture, as mentioned in
Chapter 3. The ductile fracture is categorised by a considerable amount of plastic
deformation prior to material separation. The brittle fracture, on the other hand,
occurs at small or no plastic deformation. However, the steels examined in this
study are considered as ductile, therefore only this type of fracture will be consid-
ered. From a micromechanical point of view the ductile fracture is characterised
by nucleation, growth and coalescence of voids in the material until finally the
load-bearing area has been reduced and material separation occurs. The reduction
in load-bearing area due to voids will lead to material softening. In damage me-
chanics this effect is described as damage accumulation affecting the constitutive
behaviour and the models presented by Gurson (1977) and Lemaitre (1985) are
two examples of such models. One simple way to consider this damage is proposed
by Lemaitre and Chaboche (1990), given as a relationship between the initial and
the damaged area in a certain direction, see Figure 9.

Figure 9: Initial area, S, and damaged area, SD , in the normal direction, n̄.

The damage variable, Dn , can then be interpreted as the relationship between the
damaged area, SD , and the initial area, S, respectively.

SD
Dn = (21)
S
In this definition of damage, ultimate fracture is expected when Dn reaches the
value of 1, i.e. when the entire surface is damaged and there is no material left to

19
CHAPTER 5. FRACTURE MODELLING

keep the parts together. The area that can carry the load on the material is the
area given by the difference between the initial and damaged area (S − SD ). If the
stress far away from the damaged region, σ∞ , is considered and the effective stress
working on the material in the damaged region is evaluated, this stress is given by

Sσ∞ σ∞
σ̄ = = (22)
S − SD 1 − Dn

Regarding the fracture models, these do not affect the constitutive law before frac-
ture occurs. Most fracture models consider fracture to occur when some state value
is reached. Numerous fracture models have been presented and an overview of some
of these is presented in Wierzbicki et al. (2005). Several different phenomena may
contribute to the failure process. In Teirlinck et al. (1988) four failure phenomena,
observed in uniaxial tension specimens, are described: (i) plastic failure, (ii) duc-
tile fracture, (ii) shear fracture and (iv) cleavage and brittle intergranular failure.
It should be noted that failure type (iv) cleavage and brittle intergranular failure
is considered as a brittle fracture, which is not included in this study. The first
three failure types presented are also identified for sheet metals in Hooputra et al.
(2004), were the term plastic failure is generalised to represent any sheet instabil-
ity. Plastic failure or sheet instability is often the primary mechanism leading to
failure, even if no material separation occurs at the point of instability. At insta-
bility the deformation is localised in a small region and further deformation will
rapidly lead to fracture. Failure due to sheet instability is examined separately in
Chapter 6. In this chapter the ductile and shear fractures, together with some of
the phenomenological models used for their representation in an FE environment,
are presented. First some important variables and quantities are specified which
are frequently used in the phenomenological models. As stated by Teirlinck et al.
(1988) the fracture is strongly connected to hydrostatic pressure, p, which can be
expressed in terms of the first invariant of the Cauchy stress, I1 , see Equation 23.

I1 = tr(σ) = σ1 + σ2 + σ3 = 3σm = −3p

1 1
J2 = s : s = tr(s2 ) (23)
2 2
1
J3 = det(s) = tr(s3 )
3

Another important stress measure is the deviatoric stress tensor, s,

I1
s=σ− I (24)
3

where I is the unit tensor. Also the second and third invariants of the deviatoric
stress tensor, J2 and J3 , see Equation 23, are frequently used. From the first

20
invariants of the Cauchy stress, I1 , the second and third invariants of the deviatoric
stress, J2 and J3 , the Haigh-Westergaard coordinates, ξ, ρ and θ can be defined

1
ξ = √ I1
3
p
ρ = 2J2 = σvM (25)
√
27 J3
cos(3θ) =
3 J23/2

where σvM is the equivalent von Mises stress. These coordinates can be explained
 
from Figure 10, where the ξ-axis is aligned with the hydrostatic axis (σ1 = σ2 = σ3 ),
ρ-axis is the radius of the von Mises cylinder and θ is the Lode angle, cf. also
Figure 11. The point B in the principal stress space can now be defined in the
Haigh-Westergaard coordinates as

    r  
σ1 1 cos(θ)
ξ
 σ2  = √  1  + 2
ρ  cos(θ − 2π/3)  (26)
3 3
σ3 1 cos(θ + 2π/3)

σ2

ξ
σ3
=
σ2
=
σ1
O
θ B

A
π ,   ρ

O σ1

σ3

Figure 10: Von Mises cylinder shown in the principal Cauchy stress space.
  

21

'
CHAPTER 5. FRACTURE MODELLING

σ3

A θ
σ1 σ2
B

Figure 11: Von Mises cylinder shown in the π-plane.

The importance of the Lode angle, θ, on fracture has been reported, e.g. by Bai
and Wierzbicki (2008). Also the importance of the stress triaxiality, η, which is the
relationship between hydrostatic pressure and the equivalent von Mises stress, see
Equation 27, has been demonstrated by several authors e.g. Oyane et al. (1980)
  
and Bao and Wierzbicki (2004).

p
η=− (27)
σvM

5.1 Ductile Fracture

Even if the fracture models do not affect the elasto-plastic or elasto-viscoplastic
constitutive laws and thus produce a material softening, it is useful to consider
them as a type of growing function. One useful representation is

Zεf
F (state variables)dε̄p ≤ C (28)
0

where the function, F , may depend on any state variable and is integrated on
the effective plastic strain, ε̄p . The state variables can be divided into observable
variables, e.g. temperature, T , total strain, ε, or internal, e.g. plastic strain, εp , or
stress state, σ, among others. The simplest form of Equation 28 is when F ≡ 1, in
which case the model predicts the equivalent plastic strain to fracture. However, a
fracture criterion given by constant equivalent plastic strain to fracture is contrary
to experimental observations.

22
 5.2. SHEAR FRACTURE

5.2 Shear Fracture

Shear fracture is a result of extensive slip on the activated slip planes. Consequently
this fracture is favoured by shear stresses. The shear fracture can also be found
under certain conditions when voids nucleate in the slip band and reduce the load-
bearing area, so that continuous plastic flow localises there. Since sheet metal
applications often are modelled by shell elements, which most often are restricted
to plane stress assumptions, it is necessary to construct a fracture model that can
predict fracture due to transversal shear.

5.3 Phenomenological Fracture Models

A short review of some of the phenomenological fracture models presented in lit-

erature is given below.

The Cockroft and Latham Criterion

Cockroft and Latham (1968) suggested a criterion based on an accumulated stress

and plastic strain. More precisely they argued that the plastic work must be an
important factor for the fracture. The total amount of plastic work done per unit
volume at the fracture point is

Zεf
p
W = σ̄dε̄p (29)
0

y
where σ̄ is the current equivalent stress (σ̄ = σiso ), εf is the fracture strain, and ε̄p
the equivalent plastic strain. However, the current stress σ̄, unlike the peak stress
σ1 , is not influenced by the shape of the necked region. A criterion based on the
total amount of plastic work will therefore predict a fracture independent of this
shape, which is contrary to experiments. Therefore, the total amount of plastic
work cannot provide a good criterion by itself, since necking is important according
to experimental observations. A more reasonable fracture criterion would be to take
the magnitude of the highest tensile stress into account. Therefore, Cockroft and
Latham proposed that fracture occurs in a ductile material when

Zεf 
< σ1 >  p
W = σ̄ dε̄ (30)
σ̄
0

23
CHAPTER 5. FRACTURE MODELLING

reaches a critical value, Wc , for a given temperature

and strain rate. The non-
dimensional stress concentration factor, <σσ̄1 > , represents the effect of the highest
tensile stress, < σ1 >= max(σ1 , 0). The reduced form

Zεf
W = < σ1 > dε̄p ≤ Wc (31)
0

is used for the fracture evaluation, and fracture is expected when the integral
reaches a critical value, Wc , which is determined from experiments. The Cockroft
and Latham model implies that fracture in a ductile material depends both on the
stresses and plastic strain states, i.e. neither stress nor strain alone can describe
ductile fracture. The benefit of using the largest principal stress, σ1 , is that it can
be expressed as a function of the hydrostatic pressure, p, the second invariant of
the deviator stress, J2 and the Lode angle, θ, as

r
4J2
σ1 = −p + cos θ (32)
3

The Wilkins Criterion

The model by Wilkins et al. (1980), also known as the Rc Dc model, states that
two factors increase the fracture threshold: the hydrostatic stress and the asym-
metric stress. The hydrostatic stress accounts for the growth of holes by spalling.
Interrupted tension tests have shown initiation and growth of voids which form
a fracture surface. The asymmetric stress accounts for the observation that the
elongation at fracture decreases as the shear load increases in fracture tests with
combined stress loads. The accumulated damage is expressed as

Zεf
D= ω1 ω2 dε̄p ≤ 1 (33)
0

and the weight factors for the hydrostatic stress and the asymmetric stress are
expressed as

 C1
1
ω1 = , ω2 = (2 − AD )C2 (34)
1 − C3 σ m

24
 5.3. PHENOMENOLOGICAL FRACTURE MODELS

where ω1 is the hydrostatic pressure weight, ω2 is the asymmetric stress weight, dε̄p
is the equivalent plastic strain increment, σm is the hydrostatic pressure and C1 ,
C2 , and C3 are material constants. The parameter AD included in the asymmetric
stress weight, ω2 is defined as

 
s2 s2
AD = min , (35)
s3 s1

where s1 , s2 and s3 are the principal stress deviators. The parameter AD ranges
from 0 to 1, and when AD = 1 the stress field is symmetric (and asymmetric when
AD = 0). Fracture is expected when the damage parameter, D, in Equation 33
reaches 1.

The Oyane Criterion

Oyane et al. (1980) suggested a fracture model to consider the effect of stress
triaxiality on ductile fracture. The model is derived from the plasticity theory for
porous media as

Zεf
(C1 + η) dε̄p ≤ C2 (36)
0

where C1 and C2 are material constants and η is the stress triaxiality.

The Johnson and Cook Criterion

The fracture model presented by Johnson and Cook (1985) is a purely phenomeno-
logical model, which is based on a similar relation as the hardening model presented
by Johnson and Cook (1983). The model uses a damage parameter D, and when
this parameter reaches the value of 1, ultimate fracture is expected. The definition
of the damage parameter is

Zεf
1 p
D= dε̄ ≤ 1 (37)
Φ
0

where εf is the equivalent strain to fracture and dε̄p is the increment of equivalent
plastic strain. The expression for the function Φ is given as
  p 
−C3 η

ε̄˙
Φ = C1 + C2 e 1 + C4 ln (1 + C5 T ) (38)
ε̇0

25
CHAPTER 5. FRACTURE MODELLING

where C1 , ..., C5 are material constants, which can be determined from experiments,
η is the stress triaxiality, ε̄˙p is the equivalent plastic strain-rate, ε̇0 is a reference
strain-rate and T is the temperature. As noted from the Equations 37 and 38 the
model depends on strain, strain-rate, temperature and stress triaxiality.

The Bressan and Williams Criterion

Shell elements are often used for modelling sheet metal applications, but are fre-
quently restricted to plane stress conditions. Thus, they are not able to describe
transversal shear stresses. It is therefore necessary to consider a fracture criterion
which can be used to predict this phenomenon. Bressan and Williams (1983) sug-
gested a model for predicting instability through the thickness of a sheet metal.
They suggested that the plastic strain increment should be equal to zero at an in-
clined direction through the thickness of the sheet, i.e. dεpt = 0 in the et -direction
shown in Figure 12. By assuming that the directions of the principal  . * +
stress and
strains coincide, and by using the stress and strain tensor components according
to Equation 39, together with the rotation tensor to rotate the stress and strains
4incline
to the *+ , directions
$  (en , e2 , et ),
e3

et
e2

en

e1

 $ 
 plane
Figure 12: Inclined through   the
 thickness
%   of
the
sheet.
  The stress and strain
on the et -plane is determined from the values in the coordinate axis ei , where
i = 1...3, and the angle, θ.

dεpt = sin2 θ dεp1 + cos2 θ dεp3 = 0 &'

   
σ1 0 0 dεp1 0 0
σ =  0 σ2 0  , dεp =  0 dεp2 0  (39)
0 dε0p1 +0dεp3 0 0 dεp3
cos 2θ =
dεp1 − dεp3
&)
the strain increment in the et -direction and the shear stress on the inclined surface
can be found as
−dεp2
cos 2θ =
2dεp1 + dεp2
=−
β
2+β
&3
26

τnt = cosθsinθ(σ1 − σ3 ) = sin2θ

(σ1 − σ3 )
2
&4


τnt =
σ1 − σ3
2
1−(
β 2
2+β
) ≤ τc ':

:
 5.3. PHENOMENOLOGICAL FRACTURE MODELS

dεpt = sin2 θdεp1 + cos2 θdεp3 = 0 (40)

and

σtn = −(σ1 ) sin θ cos θ (41)

Assuming that fracture occurs when the shear stress, σtn , on the inclined surface
reach a critical value τc , and that the material is plastically incompressible, the
following is obtained

−dεp2 −β
cos 2θ = = (42)
2dεp1 + dεp2 2+β

and

2τc
sin 2θ = − (43)
σ1

Here β = dεp2 /dεp1 is the relationship between the in-plane principal plastic strains
and τc is a material constant. From Equations 42 and 43 the following inequality
can be obtained

s  2
σ1 −β
τ= 1− ≤ τc (44)
2 2+β

where σ1 /2 is the shear stress on the inclined surface, cf. Figure 12. As long as
the inequality is fulfilled no fracture is expected due to through-thickness shear.
However, in this work a normal stress through the thickness has been introduced to
obtain C 0 continuity over the element edges. Therefore, the criterion in Equation 44
is modified such that

s  2
σ1 − σN D −β
τ= 1− ≤ τc (45)
2 2+β

where σN D is the through-thickness normal stress, and thus (σ1 − σN D )/2 is the
shear stress on the inclined surface.

27
CHAPTER 5. FRACTURE MODELLING

The Han and Kim Criterion

Han and Kim (2003) suggested a criterion for ductile fracture, which combines
the model by Cockcroft-Latham with a maximum shear stress criterion and the
through thickness strain. Thus,
Zεf
< σ1 > dε̄p + C1 (σ1 − σ3 )/2 + C2 εt ≤ C3 (46)
0

where C1 , C2 , and C3 are material parameters, determined from experiments, and

εt is the thickness strain. Here it has been used so that the principal stresses follows
σ1 > σ2 > σ3 . Then the term (σ1 − σ3 )/2 is the maximum shear stress.

28
Modelling Instability
6
Even if there is no material separation during an instability, it is most often consid-
ered as a limit of the operational use of the product. At instability the deformation
is localised as a narrow band of the same width as the thickness of the sheet, see
e.g. Dieter (1986). Due to the localisation, further deformation will rapidly lead
to a material fracture. The prediction of instability limits is therefore important
in order to find the correct failure behaviour.

6.1 Instability in Plane Strain

The instability limit can, for the special case of plane strain, be derived analytically.
However, the elastic strains are neglected here and all strains are assumed to be
plastic. Consider a rod with a length, l, width, w, and thickness, t. The force
in the length direction can be expressed as F11 = σ11 wt, where the subscript 11
corresponds to the length direction. The instability arises when the force reaches
a maximum viz. dF11 = 0. In terms of stress and strains, the following is obtained

dσ11
= σ11 (47)
dεp11

where dll = dεp11 , dw

w
= dεp22 = 0 (plane strain), dtt = dεp33 , has been used and the
fact that the material is assumed to be plastically incompressible, i.e. dεp11 + dεp22 +
dεp33 = 0. By assuming a plane stress assumption, the stress and strain components
can be expressed as

   
1 0 0 1 0 0
ε= εp11  0 0 0  , σ = σ11  0 γ 0  (48)
0 0 −1 0 0 0

where γ is the stress ratio producing plane strain deformation, see Figure 13.

29
 ⎣ ⎦ 3 3 ⎣ ⎦
σ3 cos(θ + 2π/3)

CHAPTER 6. MODELLING INSTABILITY

n
σ11


γ
σ22

  
Figure 13: Stress path for a plane strain condition.

By using the principle of equivalent plastic work,

σ : εp = σ11 εp11 = σ̄ ε̄p (49)

where σ̄ and ε̄p are the effective stress and plastic strain, respectively. Since the
effective stress, σ̄, is a homogeneous function of degree one, i.e. σ̄(χσ) = |χ|σ̄(σ),
the stress in the 11-direction can be expressed as
%
p(γ)
z }| {
1
σ11 = σ̄ (50)
σ̄(σ11 = 1, σ22 = γ, σ12 = 0)

From Equation 49 the strain can be expressed as

1 p
ε11 = ε̄ (51)
p(γ)

y
From the yield function σ̄ = σiso (ε̄p ) is known. Therefore, the instability condition,
y
Equation 47, can be expressed in terms of the yield stress, σiso , where the only
p
unknown is the equivalent strain, ε̄ .

y
dσiso (ε̄p ) 1 y p
= σ (ε̄ ) (52)
dε̄ p p(γ) iso

The strains causing localisation can then be evaluated from Equation 51. The
instability at plane strain is therefore strongly affected by the hardening after
necking, and by increasing σ100 from Equation 13 the strain at instability can be
increased, see Figure 14.

30
 6.2. ANALYTICAL INSTABILITY MODELS

1000 y
σiso
y
dσiso
dε̄p p(γ)
900
Stress, σ [MPa]

800

700

600 σ100

500
0.15 0.2 0.25 0.3
Equivalent plastic strain, ε̄p [-]

Figure 14: Instability under plane strain.

6.2 Analytical Instability Models

The theory in Section 6.1 can be used to obtain similar expressions for more gen-
eral instability models. In Aretz (2004) a strategy is given for the computational
implementation of the models according to Hill (1952), Swift (1952) and Hora et al.
(1996). These models are presented in terms of principal stresses and strains, see
Table 2. Some general assumptions for all these models are
• Elastic strains are neglected, i.e. ε = εp and dε = dεp .
• Plane stress is assumed, i.e. σ11 , σ22 are the only non-zero principal stresses
(σ11 ≥ σ22 ).
• No shear strains are assumed, i.e. εp11 , εp22 and εp33 are the only non-zero
strains. εp11 and εp22 are the major and minor principal strains, respectively.
• The material is assumed to be plastically incompressible, i.e.
dεp11 + dεp22 + dεp33 = 0.
• Linear strain paths are assumed, i.e. εp11 /εp22 = dεp11 /dεp22 = const.
• Only isotropic hardening is considered and the material follows an associative
flow rule.
• The anisotropy axes coincide with the principal strain and stress axes.

31
CHAPTER 6. MODELLING INSTABILITY

Table 2: Instability models according to Hill (1952), Swift (1952) and Hora et al.
(1996) expressed in terms of principal stresses and strains.
dσ11 dεp
Hill p = σ11 (1 − β) β = 22
dε11 dεp11
dσ11
Swift = σ11
dεp11
∂σ11 ∂σ11 ∂β dεp
Hora et al. p + p = σ11 β = 22
∂ε11 ∂β ∂ε11 dεp11

The stress and strain-rate ratios are expressed by γ and β, respectively. How-
ever, the strain-rate ratio can be expressed in terms the stress ratio, γ, since an
associative flow rule is assumed.

dεp22 ∂ σ̄/∂σ22
β= = (53)
dεp11 ∂ σ̄/∂σ11 σ11 =1,σ22 =γ

Thus, the strain-rate ratio can be expressed as β = β(γ). From Equation 50 the
stress can be expressed in terms of the equivalent measurement and the positive
function p(γ). The equivalence of principle plastic work can now be used to express
the strain in terms of the equivalent strain measurement

1 p
σ̄ ε̄p = σ11 εp11 + σ22 εp22 = σ11 εp11 (1 + γβ(γ)) ⇒ εp11 = ε̄ (54)
q(γ)

where the function q(γ) = p(γ)(1+γβ(γ)) has been introduced. By prescribing the
y
yield stress, σiso = Y , as a function of the accumulated equivalent plastic strain,
p
ε̄ , the localisation criteria can be expressed as shown in Table 3.

Table 3: Instability models according to Hill (1952), Swift (1952) and Hora et al.
(1996) expressed in terms of yield stress functions.
Hill Y 0 (ε̄p )q(γ) = Y (ε̄p )(1 + β(γ))

Swift Y 0 (ε̄p )q(γ) = Y (ε̄p )

p(γ)q(γ)β(γ)
Hora et al. Y 0 (ε̄p )p(γ)q(γ) − Y (ε̄p ) = p(γ)Y (ε̄p )
β 0 (γ)ε̄p

In Table 3 new derivatives have been introduced defined as Y 0 (ε̄p ) = dY (ε̄p )/dε̄p ,
β 0 (γ) = dβ(γ)/dγ and p0 (γ) = dp(γ)/dγ. From these equations the only unknown

32
 6.3. THE MARCINIAK AND KUCZYNSKI MODEL

quantity is the accumulated equivalent plastic strain, ε̄p , if the stress path is as-
sumed to be prescribed and proportional. By prescribing different strain paths (i.e.
stress ratios) a FLC can be constructed.

6.3 The Marciniak and Kuczynski Model

The model described by Marciniak and Kuczynski (1967), often referred to  as
the 
MK model, considers an imperfection or defect in a sheet metal, most often realised
 reduction
as a thickness   over a narrow band, see Figure 15. The thickness difference
is described by an imperfection factor F = tB /tA .
σ11

tA
A

σ22 tB
B

  
Figure 15: The geometry for the MK-model.
σ11
The elastic deformations are neglected in thisBanalysis.
1
Due to equilibrium in the
region B the stress in the imperfection can be expressed
A1 as

B0
tA A σA 1
B
σ11 = σ11 = 11 A0 (55)
tB F αB
Compatibility in the 22-direction requires that the strains in this direction must
be identical in both the regions 1

αA

εA
22 = εB
22 ⇔ dεA
22 = dεB
22 σ22 (56)
  
A plane stress state is used and the material is considered as plastically incom-
pressible. Furthermore, a proportional deformation of the region A is specified
as

   
1 0 0 1 0 0
A 
σ A = σ11 0 γA 0  , ε A = ε A
11
 0 βA 0  (57)
0 0 0 0 0 −(1 + βA )

33


CHAPTER 6. MODELLING INSTABILITY

where γA and βA describe the deformation. The strain ratio, βA , can be expressed
as,
!"   

∂ σ̄/∂σ22
βA = (58)
∂ σ̄/∂σ11 σ11 =1,σ22 =γA ,σ12 =0

since that the material is assumed to follow an associated flow rule. The stress
B
σ11 can then be found from  Equation 55. Since the stress in the length direction
is slightly larger inside the band, the yield surface is first reached here. However,
dεB
11
βB condition, Equation
due to the compatibility 56, no plastic deformation can occur
before the yield surface is also reached
A
dε11 for region A. If further deformation occurs,
the state of stress in the band is moved on the yield surface towards point B1 ,   
see Figure 16(a). At βA a certain point the yield limit is also reached in the uniform
region, point A1 , and plastic deformation can occur. However, due to the restriction
that the strain in the 22-direction needs to be equal in both regions A and B,
A B
the equivalent dε 22 = dεstrain
plastic 22 will be different, i.e. ε̄pB > ε̄pA , see Figure 16(b). A
continuous deformation will increase the difference in equivalent plastic strain and
 
move the stress state in region B towards the state of plane strain, see Bf in Figure
16(a). When the state of plane strain is reached, instability is assumed to occur.

σ11
Bf

Af

B1
A1

B0 
A0
dεB
11
βB

 dεA
11

1 βA

γA

σ22 dεA B
22 = dε22

(a)  (b)  



Figure 16: MK model (a) stress state for region A and B and (b) relationship
between the strain paths. &

34


        

δ22

δ11 δ11

δ22 6.4. FINITE ELEMENT MODEL

6.4 Finite Element Model   

δ22
Detailed FE models can generally be used to capture the instability phenomena,
B A
e.g. Lademo et al. (2004a). To predict the instability, a square patch of elements,
see Figure 17, is used with has an inhomogeneity in its thickness distribution.
δ11
The patch is stretched in different δ11 linear strain paths and to
directions to obtain
produce the FLC in the forming limit diagram (FLD). The strain paths have been
prescribed such that
C
δ22



δ11 = w0 eεcosθ − 1 , δ22 = w0 eεsinθ − 1 (59)
  
where w0 is the width of the plate, θ = tanβ is the relationship between the strains,
and ε has been given as a smooth function of time.

δ22

δ11 δ11
w0

δ22
w0

  
Figure 17: Patch of elements for the instability prediction.

The focus of the study is on the relationship (between the local thickness strain
increment and the average thickness strain increment, on an area Ω.

∆ε33,i
ζi = (60)
∆εΩ33

The average thickness strain increment is

N
1X
∆εΩ
33 = ∆ε33,i (61)
N i=1

where N is the number of elements in the domain Ω. Instability is assumed to

occur when ζi reaches a critical value in any element i. However, since the strain-
rate shows an oscillating behaviour pattern, it is necessary to define the point of

35
CHAPTER 6. MODELLING INSTABILITY

instability over a period of time, viz. ζi needs to exceed a critical value ζc for a
number of time steps nt . Lademo et al. (2004a) describe the thickness variation
as a normal distributed random field with mean value µ and standard deviation,
usually denoted σ. However, according to Fyllingen et al. (2009) one drawback
with this method is that the variation depends on the number of nodes. Hence, a
refinement of the FE mesh will lead to a different random field. Consequently the
variation in the thickness is described here independently of the FE mesh

t(x, y) = µ(x, y) + Z(x, y) (62)

where both the mean value, µ, and the residual term, Z, may depend on the global
coordinates x and y. In this work the mean value has been set to a constant value,
i.e. µ(x, y) = µ, and for the residual term, Z, a Gaussian zero mean homogeneous
random field according to Shinozuka and Deodatis (1996) has been used.

6.5 Evaluation of Instability

The different methods presented above can be summarized in an FLD, see Fig-
ure 18. As can be seen all models predict about the same instability limit in
their domain of usage. However, all instability models used here, except the FE-
based model which can be used for any material model, use a plane stress as-
sumption. The material model presented in Chapter 4 is modified to also include
a through-thickness normal stress, σN D , to obtain a C 0 continuous element for-
mulation. Therefore, the FE-based model has been used throughout this study
for instability predictions. A square patch of finite elements is used with an in-
homogeneous thickness distribution in order to predict the instability. The patch
is stretched such that linear strain paths are obtained in different directions. In
Lademo et al. (2004a) and Lademo et al. (2004b) the limit strains causing local-
isation are considered as the total strains in the patch when an instability has
occurred. However, in this study the local strains in an element within the locali-
sation area are considered. In order to find the strain limit at localised necking in
the patch, two elements are considered: one inside the localisation zone and one
some distance away from it, see Figure 19(a) elements A and B, respectively. The
limit value is then obtained as the strains in the finite element within the locali-
sation zone at a stage when the strains in the distant element do not increase, see
Figure 19(b).

36
 0.5
Swift 1952
Hill 1952
Strain in 11-direction, ε11 [-] Hora et al. 1996
0.4 Marciniak and Kuczynski 1967
Lademo et al. 2004

0.3

0.2

0.1
       

0
δ22 −0.1 0 0.1 0.2 0.3
Strain in 22-direction, ε22 [-]

Figure 18: Different models to predict instability.

δ11 δ11

δ22 0.6
A (element inside localisation zone)
B (element far away from localisation zone)
   0.5 C (elements across the localisation zone)
Local strain, ε∗ [-]

0.4
δ22
B A
0.3

δ11 δ11 0.2

0.1

C
0
δ22 0 0.05 0.1 0.15 0.2
Global strain, ε = ln( δ+w
w0 ) [-]
0


(a)   (b)

Figure 19: Evaluation procedure for instability limits, (a) chosen elements and
(b) local vs. global strain.
δ22

δ11 δ11
w0

δ22
w0

  

(
Mechanical Experiments
7
A number of mechanical tests have been performed to calibrate both the constitu-
tive relations and the fracture models. In this study six different mechanical tests
have been conducted:
• Pre-deformation
• Tensile
• Shear
• Plane strain
• Bulge
• Nakajima
The pre-deformation specimens have been performed in an MTS hydraulic machine
with a 250 kN load cell. The tensile, shear and plane strain tests have been carried
out in an INSTRON 5582 machine with a 10 kN load cell.

7.1 Pre-Deformation
Pre-deformation of large specimens, according to Figure 20(a), was performed in
order to produce broken strain paths. Due to limitations of the test rig, the speci-
men length, width and thickness were 1100 mm, 125 mm and 1.5 mm, respectively.
Since all material tests in this work were performed under quasi-static conditions
the deformation was performed at a constant crosshead speed of 5 mm/min, which
resulted in a strain-rate of approximately 10−4 s−1 . The unloading was set to last
for 1 min. The design of the specimen was made in an attempt to maximize the
zone of homogenous strain. The homogeneity at the centre of the specimen was
confirmed both by FE simulation, see Figure 20(b), and by ocular measurements.
The mid section of the specimen was divided into 16 equal sections, in which the
thickness, width and length were measured both before and after completed load-
ing. The pre-deformations were performed both in the RD and in the TD, i.e.
ψ = 0◦ and ψ = 90◦ according to Figure 21. For each direction the materials have
been pre-deformed to two significant strain levels; one level chosen close to the
strain causing diffused necking and the other level chosen to approximately half of

39
CHAPTER 7. MECHANICAL EXPERIMENTS

the first one. The levels of pre-straining for each material and direction are dis-
played in Table 4. After the pre-straining operation tensile, shear and plane strain
specimens have been machined out from the centre part of the larger specimens,
see Figure 20(c).

(a)

(b)

(c)

Figure 20: a) Geometry of pre-strain specimen. b) Sketch including some cut out
specimens. c) Distribution of longitudinal plastic strain. The iso-curves represent
strain levels ranging from 9% to 11 % with 0.2% steps. Dimensions are given in
mm.

RD

Figure 21: Angles defining pre-strain directions, ψ, and directions of subsequent

testings, φ = ψ + θ.

40
 7.2. TENSILE TEST

Table 4: Approximative plastic strains after pre-straining.

Material Docol 600DP Docol 1200M
Direction RD TD RD TD
p
ε [%] 5 10 4.5 8 0.5 1 0.4 0.7

7.2 Tensile Test

The tensile tests have been used to represent the hardening up to diffuse necking
and also to represent the anisotropy of the yield function and the kinematic hard-
ening. Tensile test specimens according to Figure 22 have been performed in the
φ = 0◦ , 45◦ and 90◦ directions, where φ is the angle to the RD, see Figure 21, both
for the as-received and the pre-deformed material specimens. During the tensile
tests the load, elongation and width contraction were recorded. The elongation, εL ,
has been measured by an INSTRON 2620-601 extensometer with a gauge length
of L0 = 12.5 mm, and the width contraction, εW , has been measured with an
MTS 632.19B-20 extensometer over the entire width of the specimen. The defor-
mation was performed with a constant crosshead speed of 0.45 mm/min, which
resulted in a strain-rate of approximately 10−4 s−1 . The anisotropy in yield stress
and plastic flow were evaluated directly from the tensile tests of the as-received
specimens. Also the hardening behaviour up to diffused necking is obtained di-
rectly from the experimental data. Numerical simulations and inverse modelling
are used to identifying the kinematic hardening of the pre-strained materials.

Figure 22: Geometry of the tensile test with dimensions in mm.

7.3 Simple Shear Test

The geometry of the simple shear test specimen used in this study is shown in
Figure 23. The shear tests are performed on both as-received and on pre-deformed
specimens in the φ = 0◦ , 45◦ and 90◦ directions. The two attachments for the
shear test have been in the form of a pin at both ends to prevent the specimens
from rotational loading. In order to obtain a strain-rate of approximately 10−4 s−1 ,
the displacement rate of 0.03 mm/min has been used. The shear test results are
used in combination with inverse modelling to gain an accurate description of the

41
CHAPTER 7. MECHANICAL EXPERIMENTS

hardening behaviour for large strains, to find the yieldsurface exponent, and to
calibrate the Cockroft-Latham fracture model.

Figure 23: Geometry of the shear specimen with dimensions in mm.

7.4 Plane Strain Test

The geometry of the plane strain test (notched tensile test) is shown in Figure 24.
Similar designs have been utilised in other studies, e.g. Lademo et al. (2008).
The plane strain tests were performed in φ = 0◦ , 45◦ and 90◦ directions, with a
displacement rate of 0.06 mm/min. During the test the load, grip motion, width
reduction and midsection elongation were measured. The latter was measured by
an INSTRON 2620-601 extensometer width a gauge length of 23 mm, while the
width reduction was measured by an MTS 632.19B-20 extensometer. The plane
strain test is used in combination with inverse modelling in order to calibrate the
Bressan-William fracture model.

7.5 Balanced Biaxial Test

A balanced biaxial bulge test was performed on the virgin materials. The pressure
was applied with a punch made of silicon and a pattern of randomly placed dots
was sprayed on the sheet surface. The specimen was video recorded by two cameras
during the subsequent testing. The in-plane strains and the radius of the bulge
were then evaluated from the recording of the pattern motion using an ARAMIS
strain measurement system. Balanced biaxial tests have been used in order to
obtain the stress and strain ratios rb and kb , respectively, cf. Sigvant et al. (2009).

42
 7.6. NAKAJIMA TEST

Figure 24: Geometry of the plane strain specimen with dimensions in mm.

7.6 Nakajima Test

A number of Nakajima tests, see ISO (2008), have been conducted in order to
evaluate the failure behaviour. The tests were made on virgin material and the
geometries, see Figure 25, were chosen so that the first quadrant (ε1 ≥ 0, ε2 ≥ 0)
of the FLD was covered.

(a) (b)

Figure 25: Nakajima specimen geometry dimmensions in mm: (a) circular (b)
with different waists (L=130,140,150,160,165 and 170 mm).

The Nakajima tests were performed in an Interlaken ServoPress 150. The machine
has a punch diameter of 100 mm and a load cell of 550 kN. The clamping force
was limited to 700 kN and the punch motion was set to 1 mm/s. During the
test the force and displacement of the punch was recorded. The specimens were

43
CHAPTER 7. MECHANICAL EXPERIMENTS

treated with three layers of oil and plastic film to reduce friction. Before the
test a mesh with a 2 mm grid size was etched onto the specimens. The strains
- -
were then evaluated optically by an AutoGrid 4.1 Strain Analysis System, which
used information from four cameras, each of which recorded 30 images per second
during the test. The image just before the fracture was used to evaluate the limit
strains. Since the image just before fracture is used, the strains may be well past
the limit of localised necking. Therefore a method similar to the one described by
Bragard et al. (1972) has been used to evaluate the limit strains at the onset of
localised necking. This method is straightforward and the major and minor strain
distributions on a few lines across the instability are considered. Afterwards the
strains inside the localised region are excluded and a polynomial fit is adopted to
the remaining strains. The maximum strain from this polynomial is chosen as the
major strain limit, see Figure 26. For the minor strain a linear fit is performed at
the location of the maximum.

Excluded Points Major strain

Minor strain
0.4 Polynomial Fit
Major Strain at
0.3 Localised Necking
Strain, ε [-]

x
0.2

0.1 Minor Strain at

Localised Necking
0

−0.1
−15 −10 −5 0 5 10 15
Position along, x [mm]

(a) (b)
  $      %                
                          
Figure 26: Illustration of the Bragard et al. (1972) method; (a) polynomial fit
to strains and (b) lines of elements across the localisation zone used for strain
evaluation.

%

44
 Finite Element Modelling
8
    

Finite element analyses of the tensile, shear, plane strain and Nakajima experiments
have been performed using LS-DYNA,
3 ! seeHallquist (2009). All FE meshes have
been produced in TrueGrid, see Rainsberger (2006). A non-local treatment of the
5
fracture parameters is needed in order to reduce the mesh dependency. However,
<    
the C 0 continuity of the thickness across the element edges as obtained by the
use of an extended shell element formulation contributes to a regularisation of
the thickness strain, and thus reduces the mesh dependency. The introductionεxz = 0 of
0
the normal through-thickness stress, σN D , enables a C continuity in the element
formulation, cf. Borrvall and Nilsson (2003). However, the inclusion of the normal
    
stress also introduces transversal shear stresses and strains, see Figure 27. The
transversal shear stresses are not entering into the elasto-plastic constitutive model,
3 !butare
 updated elastically.      
5 5
<     <    
εxz = 0 εxz = 0

(a) (b)

Figure 27: Element    (a) ordinary plane

 formulation   
 element
stress (b)&a continuous
thickness based element.
5
<    
⎧ prevent ⎫ ⎡the localisation and will thus contribute ⎤⎧ ⎫
The transversal shear stresses will
⎪
⎪ σ ⎪
⎪ λ + 2μ λ λ 0 0 0 ⎪
⎪ εxx ⎪
⎪
⎪
to a stiffening effect. An improved
⎪
⎪
xx ⎪
result
⎪ ⎢
⎪ can be obtained by using an effective ⎥⎪⎪
⎪
⎪
⎪
⎪
⎪
ε
⎪ xz 
= 0 ⎪
⎪ ⎢ ⎥ ⎪
⎪ ⎪
⎪
⎪
stress, which takes all stresses ⎪into ⎪
σyyaccount ⎢ or
λ by reducing
λ + 2μ the
λ transversal
0 0 0 ⎥ ⎪ εyy
shear ⎪
⎪
⎪
⎪
⎪ ⎢
⎪ ⎥⎪⎪
⎪
⎪
⎪
⎪
⎪
stresses by a correction factor according ⎪ ⎢ ⎥ ⎪ ⎪
⎨ σzz ⎬ ⎢to the
λ theory λ of λReissner-Mindlin,
+ 2μ 0 0 0 ⎥see
⎨ ε ⎬
Hughes (2000). The latter method ⎪ has
⎪
= ⎢
been
⎢ adopted here, i.e. the
⎥
transversal
⎥ ⎪
zz
⎪
%1
⎪ σxy ⎪ ⎪ ⎢ 0 0 0 μ 0 0 ⎥⎪ ⎪ ε ⎪
shear stresses have been reduce⎪ ⎪
⎪by a shear
⎪ ⎪ ⎢
⎪ ⎢correction factor, κ, according to ⎥ ⎪
⎪ ⎥ ⎪ xy
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪ σyz ⎪
⎪
⎪ ⎢ 0 0 0 0 κμ 0 ⎥⎪
⎪
⎪ εyz
⎪
⎪
⎪
⎪
⎪ ⎪
⎪ ⎣ ⎦⎪
⎪ ⎪
⎪
⎪ ⎪ ⎪ ⎪
⎩ σzx ⎭ 0 0 0 0 0 κμ ⎩ εzx ⎭
     &
E E
σyz = κ εyz , σzx = κ εzx (63)
2(1 + ν) 2(1 + ν)
⎧
⎪
⎪
⎫
⎪
⎪
⎡
μ= E ⎤⎧ ⎫Eν
λ = (1+ν)(1−2ν)
⎪
⎪ ⎪ %2
⎪ εxx ⎪
σxx λ + 2μ λ λ 0 0 2(1+ν)
0
⎪
⎪ ⎪
⎪ ⎢ ⎥⎪ ⎪ ⎪ 45
⎪
⎪ ⎪
⎪ ⎢ ⎥⎪⎪ ⎪
⎪
⎪
⎪ σyy ⎪
⎪ ⎢ λ λ + 2μ λ 0 0 0 ⎥⎪⎪ εyy ⎪
⎪
⎪
⎪ ⎪
⎪ ⎢ ⎥⎪⎪ ⎪
⎪
⎪
⎪ ⎪
⎪ ⎢ ⎥⎪⎪ ⎪
⎪
⎨ σ ⎬ ⎢ λ λ λ + 2μ 0 0 0 ⎥⎨ εzz ⎬
⎪
zz
⎪
=⎢
⎢
⎥
⎥⎪ %1
εxy ⎪
⎪
⎪
⎪
⎪
⎪
σxy ⎪
⎪
⎪
⎪
⎪
⎢
⎢
⎢
0 0 0 μ 0 0 ⎥⎪
⎥⎪
⎥⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
&
⎪
⎪ σyz ⎪
⎪ ⎢ 0 0 0 0 κμ 0 ⎥⎪⎪ εyz ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪ ⎣ ⎦⎪⎪
⎪
⎪
⎪
⎪
⎪
⎩ σzx ⎪
⎭ ⎪
⎩ ⎪
0 0 0 0 0 κμ εzx ⎭

μ= E
2(1+ν) λ= Eν
(1+ν)(1−2ν)
%2
CHAPTER 8. FINITE ELEMENT MODELLING

where E and ν are the Young’s modulus and Poisson’s ratio, respectively. However,
κ will also affect the hardening after diffused necking. Thus, the hardening after
diffused necking has been modified compared to what was presented in Larsson
et al. (2011) in order to maintain a good agreement between simulations and test
results. The parameters of the hardening after necking and the shear correction
factor, κ, have been obtained from inverse modelling of the tensile, shear and plane
strain tests.

8.1 Pre-Straining
FE analyses of the tensile and shear experiments have been performed both for
virgin and pre-strained specimens. The pre-straining operation, which produced
the non-proportional strain paths, was made by stretching one single " #
element to $
the strain levels shown in Table 4. The back-stress tensor, α, equivalent plastic
strain, ε̄p , and fracture parameter, W , are subsequently mapped on the models
3  
of the small  used in the subsequent analyses, see Figure 28. Since the
specimens
strain field from stretching one single element is mapped on to all elements of the
model of each small specimen, no strain inhomogeneities will be present in the
subsequent simulations.

'#  "  7    "

( p
α, ε̄ , W (
F
ψ
F ⇒ F
φ
F

y - y -
x x

Figure 28: One pre-strained element from which

   the
backstress, α, equivalent plastic
strain, ε̄p , and fracture parameter, W , are mapped to the models for the subsequent
analyses.

8.2 Tensile Test

The FE model of the tensile specimen is shown in Figure 29. It consists of about
5500 shell elements. The characteristic element size is approximately 0.5 mm. A
good agreement between the stress-strain relationship simulation and tensile test
results can be observed up to diffused necking, as reported in Larsson et al. (2011).
By introducing a shear correction factor κ = 0.05 in the shell formulation and mod-
ifying the hardening curve after diffuse necking, improved agreement was obtained
between the simulation and test results, even for the crosshead displacement after
diffuse necking, see Figure 30.

46

&
 8.3. SHEAR TEST

Figure 29: FE model of the tensile test specimen.

15 30

25

10 20

Force, F [kN]
Force, F [kN]

15

5 φ = 0◦ exp. 10 φ = 0◦ exp.
φ = 0◦ sim. φ = 0◦ sim.
φ = 45◦ exp. φ = 45◦ exp.
φ = 45◦ sim. 5 φ = 45◦ sim.
φ = 90◦ exp. φ = 90◦ exp.
φ = 90◦ sim. φ = 90◦ sim.
0 0
0 5 10 15 0 1 2 3 4
Displacement, δ [mm] Displacement, δ [mm]

(a) (b)

Figure 30: Results from tensile test of virgin material in different material direc-
tion (a) Docol 600DP (b) Docol 1200M.

8.3 Shear Test

The FE model of the shear specimen is shown in Figure 31. It consists of about
12500 shell elements. The characteristic element size in the shear zone is approxi-
mately 0.06 mm, see Figure 31(b). A good agreement between simulation results
and the force-displacement relationships of the shear experiment was found after
the correction of the shear factor and modification of the hardening curve, see
Figure 32.

(a) (b)

Figure 31: (a) FE model of the shear specimen. (b) Details of the mesh in the
shear zone.

47
CHAPTER 8. FINITE ELEMENT MODELLING

3 4

3.5
2.5
3
2
Force, F [kN]

Force, F [kN]
2.5

1.5 2

1.5
1 φ = 0◦ exp. φ = 0◦ exp.
φ = 0◦ sim. 1 φ = 0◦ sim.
φ = 45◦ exp. φ = 45◦ exp.
0.5 φ = 45◦ sim. φ = 45◦ sim.
φ = 90◦ exp. 0.5 φ = 90◦ exp.
φ = 90◦ sim. φ = 90◦ sim.
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5
Displacement, δ [mm] Displacement, δ [mm]

(a) (b)

Figure 32: Results from shear test of virgin material in different material direction
(a) Docol 600DP (b) Docol 1200M.

8.4 Plane Strain Test

The FE model of the plane strain test specimen can be seen in Figure 33. The
model consists of about 7900 shell elements. The element size in the centre of
the specimen is approximately 0.2 mm. The C 0 continuous element formulation
used through this work enabled an improved correspondence between FE simu-
lation and experimental results compared to the conventional shell element for-
mulation, see Figure 34. A good agreement between simulation results and the
force-displacement relations for the plane strain tests was found, see Figure 35.

Figure 33: FE model of the plane strain specimen.

48
 8.4. PLANE STRAIN TEST

25

20
Force, F [kN]

15

10

5
Experiment
Plane Stress
With Normal Stress
0
0 0.5 1 1.5 2 2.5 3
Displacement, δ [mm]

Figure 34: Different element formulation used for simulations of the plane strain
test.

25 50

20 40
Force, F [kN]

Force, F [kN]

15 30

10 20
φ = 0◦ exp. φ = 0◦ exp.
φ = 0◦ sim. φ = 0◦ sim.
5 φ = 45◦ exp. 10 φ = 45◦ exp.
φ = 45◦ sim. φ = 45◦ sim.
φ = 90◦ exp. φ = 90◦ exp.
φ = 90◦ sim. φ = 90◦ sim.
0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5
Displacement, δ [mm] Displacement, δ [mm]

(a) (b)

Figure 35: Results from plane strain test of virgin material in different material
direction (a) Docol 600DP (b) Docol 1200M.

49
CHAPTER 8. FINITE ELEMENT MODELLING

8.5 Nakajima Test

The characteristic element size in the FE models of the Nakajima specimens was
approximately 0.75 mm. Since the Nakajima tests are used for prediction of the
failure behaviour, it was necessary to introduce a small disturbance of the thickness
in order to capture the instability correctly. The thickness variation has here been
given using the same method as for the patch discussed in Chapter 6. The Nakajima
tests were used in order to evaluate the failure behaviour and a good agreement
was obtained. However it has been seen that, the friction chosen between the
specimen and punch surface in the simulation model is of a great importance in
the prediction of the instability. The friction used in all simulation models has
therefore been chosen such that the instability is predicted as close as possible to
the specimen with a waist of 60 mm, see Figure 36.

80 2000
Experiments
Sim. µ=0.01

Fracture parameter, W [MPa]

70 Sim. µ=0.05
Sim. µ=0.10
60 1500
Punch force, F [kN]

50

40 1000
F
30

20 500
W
10

0 0
0 5 10 15 20 25 30
Punch displacement, δ [mm]

Figure 36: Different friction for prediction of instability of the Nakajima test of
Docol 600DP in RD with a waist of 60 mm.

50
Conclusions and Discussion
9

One major objective of the proposed fracture models was that they should be
easy to calibrate. Since both the models investigated in this thesis have only one
material parameter, just one basic mechanical experiment and the corresponding
FE simulation are sufficient for calibration. However, in several loading situations,
it is hard to identify whether the fracture is initiated prior to the instability or
after. In the plane strain test a localisation is observed before fracture, both in
the FE simulations and the experiments. The fracture parameters calibrated from
these experiments will then depend on how well the deformation is captured in
the localised zone. It may be argued that the plane strain test should not be used
for this calibration, since the instability arises before the fracture. The fracture
model calibrated from the plane strain test can, since the instability arises before
the fracture, be considered as an upper limit for the fracture within the FLD.
Even if one single experiment and one FE simulation are sufficient to calibrate
each of the fracture models used in this study, one fracture parameter value was
evaluated from specimens in all material directions of the virgin material. From
the evaluation of these results it can be seen that there is a large variation in
the parameters obtained from the different material directions and also within the
different tests in the same direction. The behaviour of the material in different
material directions cannot be predicted by the fracture models chosen for this
study, since the fracture models adopted are isotropic, although the constitutive
model is anisotropic. In order to predict the different fracture behaviour observed
in different directions, an anisotropic fracture model is needed. The variation of
the fracture appearance within the same material direction can be explained by the
fact that the fracture is strongly influenced by inclusions and impurities within the
material, and since these are located differently in different specimens the fractures
will differ. It may thus be preferable to specify the fracture parameters with lower
and upper limits rather than as a specific value.
Most of the Nakajima experiments and simulations show an instability before frac-
ture. Consequently the prediction of this instability is very important. However,
in this study the prediction of the instability is made by using the constitutive
model and FE models with significantly dense meshes. In order to obtain a good
prediction it is therefore necessary to carefully calibrate the material and fracture
models before subsequent use. Since the shape of the yield surface, the harden-
ing law and also the element formulation all affect the instability. Therefore it

51
CHAPTER 9. CONCLUSIONS AND DISCUSSION

is important to choose suitable mechanical experiments and element formulation

for the calibration. The validation of the failure behaviour has been performed
in experiments and FE simulations of the Nakajima tests. In the evaluation of
the quality of the fracture and instability predictions, their representation is very
important. One direct method is to compare the force-displacement curves from
simulations and experiments directly. However, since the FLD is a commonly-used
tool for failure prediction, it is also useful to represent failure predictions as FLC
in an FLD. The curves representing the locus of different failures, i.e. ductile frac-
ture, shear fracture and instability, are easy to interpret. Their representation of
the strains from the experiments causing localisation is, however, not that trivial.
There are many approaches on how to evaluate the forming limit. In experiments
in this study, the deformation was monitored during the deformation, and thus the
strains were captured close to the instance of fracture. However, the strains may
then be well past the onset of localisation and the limit of instability needs to be
evaluated using a different method. In some methods a few grid points outside the
fracture zone are evaluated and the strains there are considered to be the limit.
In this study the strains of failure are evaluated with a polynomial fit across the
instability zone.
The results obtained from the failure prediction are generally in a good agreement
with experiments. However, a large variation of fracture parameters is observed
and an interval in its values is to be preferred, since some variation is always
obtained in mechanical experiments. The variation in the mechanical experiments
may be caused by a number of different factors such as variation in setup, variation
in material and measurement error.

52
Review of Appended Papers
10

Paper I

A study of high strength steels undergoing non-linear strain paths -

Experiments and modelling

This paper presents an evaluation of the constitutive behaviour, including plastic

anisotropy and mixed isotropic-kinematic hardening of two high strength steels,
Docol 600DP and Docol 1200M, during strain path changes. A series of tensile
and shear tests was performed on both virgin and pre-strained materials. The ini-
tial anisotropy and work hardening parameters were obtained from tensile tests,
shear tests and a bulge test of the virgin material, whereas the kinematic harden-
ing parameters were identified by comparing numerical predictions to experimental
results related to the pre-strained materials. Numerical predictions using the ob-
tained parameters agree well with the experimental results, both in the case of
proportional, and under non-proportional strain paths.

Paper II

Failure of high strength steel sheet - Experiments and modelling

Failure in thin sheet metal structures of ductile material is usually caused by one
of, or a combination of, the following: ductile fracture, shear fracture or localised
instability. In this paper the failure of the high strength steel Docol 600DP and the
ultra high strength steel Docol 1200M is explored. The constitutive model used
in this study includes plastic anisotropy and mixed isotropic-kinematic hardening.
For modelling of the ductile and shear fracture the models presented by Cockroft-
Latham and Bressan-Williams have been used. The instability phenomenon is
described by the constitutive law and the finite element (FE) models. For calibra-
tion of the failure models and validation of the results, an extensive experimental
series has been conducted including shear tests, plane strain tests and Nakajima
tests. The geometries of the Nakajima tests have been chosen so that the first

53
CHAPTER 10. REVIEW OF APPENDED PAPERS

quadrant of the forming limit diagram (FLD) were covered. The results are pre-
sented both in an FLD and using prediction of force-displacement response of the
Nakajima test employing element erosion during the FE simulations. The classical
approach for failure prediction is to compare the FE simulations with experimental
determined forming limit curves (FLC). Here phenomenological models have been
used which show promising results. The benefit of using phenomenological models
is the ability to follow failure behaviour for more complex loading paths.

54
Bibliography

Aretz, H., 2004. Numerical restrictions of the modified maximum force crite-
rion for prediction of forming limits in sheet metal forming. Modelling and
Simulation in Materials Science and Engineering 12 (4), 677–692.
Aretz, H., 2005. A non-quadratic plane stress yield function for orthotropic sheet
metals. Journal of Materials Processing Technology 168 (1), 1 – 9.
Askeland, D. R., 1998. The Science and Engineering of Materials. Nelson Thornes
Ltd, Sheffield.
Bai, Y., Wierzbicki, T., 2008. A new model of metal plasticity and fracture with
pressure and Lode dependence. International Journal of Plasticity 24 (6),
1071 – 1096.
Banabic, D., Barlat, F., Cazacu, O., Kuwabara, T., 2010. Advances in anisotropy
and formability. International Journal of Material Forming 3, 165–189.
Bao, Y., Wierzbicki, T., 2004. On fracture locus in the equivalent strain and
stress triaxiality space. International Journal of Mechanical Sciences 46 (1),
81 – 98.
Belytschko, T., Liu, W. K., Moran, B., 2000. Nonlinear Finite Elements for
Continua and Structures. Wiley, Chichester.
Borrvall, T., Nilsson, L., 2003. Revision of the implementation of material 36 in
LS-DYNA. Engineering Reaserch Nordic AB, Linköping, private communi-
cation.
Bragard, A., Baret, J.-C., Bonnarnes, H., 1972. A simplified technique to deter-
mine the FLD at the onset of necking. Report no. 33, Rapport Centre de
Recherche de la Mètallurgie, Liège.
Bressan, J., Williams, J., 1983. The use of a shear instability criterion to predict
local necking in sheet metal deformation. International Journal of Mechanical
Sciences 25 (3), 155–168.
Cockroft, M. G., Latham, D. J., 1968. Ductility and the workability of metals.
Journal of the Institute of Metals 96, 33–39.
Dieter, G. E., 1986. Mechanical Metallurgy. McGraw-Hill, New York.

55
CHAPTER 10. REVIEW OF APPENDED PAPERS

Frederick, C., Armstrong, P., 2007. A mathematical representation of the multi-

axial Bauschinger effect. Materials at High Temperatures 24 (1), 1–26.
Fyllingen, O., Hopperstad, O., Lademo, O.-G., Langseth, M., 2009. Estimation
of forming limit diagrams by the use of the finite element method and Monte
Carlo simulation. Computers and Structures 87 (1-2), 128 – 139.
Gurson, A., 1977. Continuum theory of ductile rupture by void nucleation and
growth: Part I - yield criteria and flow rules for porous ductile media. Journal
of Engineering Materials and Technology 99 (1), 2 – 15.
Hallquist, J., 2009. LS-DYNA Theory Manual. Livermore Software Technology
Corporation, Livermore.
Han, H. N., Kim, K.-H., 2003. A ductile fracture criterion in sheet metal forming
process. Journal of Materials Processing Technology 142 (1), 231 – 238.
Hill, R., 1952. On discontinuous plastic states, with special reference to localized
necking in thin sheets. Journal of the Mechanics and Physics of Solids 1 (1),
19 – 30.
Hollomon, J. H., 1945. Tensile deformation. Transaction of the American Insti-
tute of Mining, Metallurgical and Petroleum Engineers 162, 268–290.
Hooputra, H., Gese, H., Dell, H., Werner, H., 2004. A comprehensive failure
model for crashworthiness simulation of aluminium extrusions. International
Journal of Crashworthiness 9 (5), 449 – 63.
Hora, P., Tong, L., Reissner, J., 1996. A prediction method for ductile sheet metal
failure in FE-simulation. Proceedings of the 3’rd International Conference
Numisheet’96, Dearborn, Michigan, 252–256.
Hosford, W., Cadell, R., 1993. Metal Forming Mechanics and Metallurgy. Prentis-
Hall, New York.
Hughes, T. J. R., 2000. The Finite Element Method - Linear Static and Dynamic
Finite Element Analysis. Dover Publications, Mineola, New York.
ISO, 2008. Metallic materials -sheet and strip - determination of forming limit
curves - Part 2: Determination of forming limit curves in the laboratory.
12004-2:2008, Austrian Standards Institute, Wien.
Johnson, G. R., Cook, W. H., 1983. Fracture characteristics of three metals sub-
jected to various strains, strain rates, temperatures and pressures. Seventh
International Symposium on Ballistics, Hague, 1–7.
Johnson, G. R., Cook, W. H., 1985. Fracture characteristics of three metals
subjected to various strains, strain rates, temperatures and pressures. Engi-
neering Fracture Mechanics 21 (1), 31 – 48.
Krauss, G., 2005. Steels Processing, Structure, and Preformance. ASM Interna-
tional Materials Park, Ohio.

56
Lademo, O.-G., Berstad, T., Hopperstad, O. S., Pedersen, K. O., 2004a. A nu-
merical tool for formability analysis of aluminium alloys. Part I: Theory.
Steel Grips 2, 427–432.
Lademo, O.-G., Pedersen, K., Berstad, T., Furu, T., Hopperstad, O., 2008. An
experimental and numerical study on the formability of textured AlZnMg
alloys. European Journal of Mechanics - A/Solids 27 (2), 116 – 140.
Lademo, O.-G., Pedersen, K. O., Berstad, T., Hopperstad, O. S., 2004b. A numer-
ical tool for formability analysis of aluminium alloys. Part II: Experimental
validation. Steel Grips 2, 433–437.
Larsson, R., Björklund, O., Nilsson, L., Simonsson, K., 2011. A study of high
strength steels undergoing non-linear strain paths - experiments and mod-
elling. Journal of Materials Processing Technology 211, 122–132.
Lemaitre, J., 1985. A continuous damage mechanics model for ductile fracture.
Journal of Engineering Materials and Technology 107, 83–89.
Lemaitre, J., Chaboche, J.-L., 1990. Mechanics of Solid Materials. Cambridge
University Press, Cambridge.
Marciniak, Z., Kuczynski, K., 1967. Limit strains in the processes of stretch-
forming sheet metal. International Journal of Mechanical Sciences 9 (9),
609–620.
Olsson, K., Gladh, M., Hedin, J.-E., Larsson, J., 2006. Microalloyed high-strength
steels. Advanced Materials and Processes 164 (8), 44–46.
Opbroek, E. G., 2009. Advanced high strength steel application guidelines. Tech.
Rep. WorldAutoSteel.
Oyane, M., Sato, T., Okimoto, K., Shima, S., 1980. Criteria for ductile fracture
and their applications. Journal of Mechanical Working Technology 4 (1), 65
– 81.
Rainsberger, R., 2006. TrueGrid User’s Manual. XYZ Scientific Applications,
Inc., Livermore.
Shinozuka, M., Deodatis, G., 1996. Simulation of multi-dimensional Gaussian
stochastic fields by spectral representation. Applied Mechanics Reviews 49,
29–53.
Sigvant, M., Mattiasson, K., Vegter, H., Thilderkvist, P., 2009. A viscous pressure
bulge test for the determination of a plastic hardening curve and equibiaxial
material data. International Journal of Material Forming 2 (4), 235–242.
Swift, H. W., 1952. Plastic instability under plane strain. Journal of the Mechan-
ics and Physics of Solids, 1, 1–18.
Teirlinck, D., Zok, F., Embury, J., Ashby, M., 1988. Fracture mechanism maps
in stress space. Acta Metallurgica 36 (5), 1213 – 1228.

57
CHAPTER 10. REVIEW OF APPENDED PAPERS

Voce, E., 1948. The relationship between stress and strain for homogenus defor-
mation. Journal of the Institute of Metals 74, 537–562.
Wierzbicki, T., Bao, Y., Lee, Y.-W., Bai, Y., 2005. Calibration and evaluation of
seven fracture models. International Journal of Mechanical Sciences 47 (4-5),
719 – 743.
Wilkins, M. L., Streit, R. D., Reaugh, J. E., 1980. Cumulative-strain-damage
model of ductile fracture: Simulation and prediction of engenering fracture
tests. Tech. Rep. Lawrence Livermore National Laboratory, Livermore.

58