CHARACTERIZATION AND MULTISCALE CRYSTAL PLASTICITY MODELING OF
MULTIPHASE ADVANCED HIGH STRENGTH STEEL
By
Bassam Abdullah Mohammed

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering – Doctor of Philosophy
2017

ABSTRACT
CHARACTERIZATION AND MULTISCALE CRYSTAL PLASTICITY MODELING OF
MULTIPHASE ADVANCED HIGH STRENGTH STEEL
By
Bassam Abdullah Mohammed
The superior combination of strength and ductility of advanced high strength steel (AHSS) alloys
makes them particularly attractive for forming complex parts needed in the global automotive
industries. In recent years, the production of AHSS sheet has been of interest for lightweight
automotive applications without compromising either performance or cost standards. The federal
government has been making unceasing demands from automotive companies to produce
automobiles with fewer emissions in order to reduce pollution, and minimize petroleum
consumption. Consequently, developing lightweight vehicles has become increasingly important
for the automotive industry in recent years. New challenges in sheet metal forming processes
emerged with the application of the new generation of AHSS, resulting in additional research
efforts in experimental, theoretical, and constitutive modeling. Accordingly, the development and
application of new constitutive equations and failure criteria for AHSS sheets also increased. It is
particularly challenging to understand how the constituent phases in AHSS contribute to their
higher strength, and resistance to localized necking and fracture. Also, volume fraction, grain
morphology and distribution, and micromechanical properties of constituent phases contribute
significantly in the selection of the material to be used for automotive applications. To overcome
these challenges with the new generation of AHSS, the application of computer-aided material
design or integrated computational materials engineering (ICME) with multiscale modeling of
material processing is highly recommended.
ii

In this study, a computationally efficient, rate independent, crystal plasticity finite element (CPFE)
model was used to simulate a multiphase AHSS sheet undergoing a large plastic deformation. The
CPFE model was developed and implemented into ABAQUS as a user material subroutine,
VUMAT, to capture the mechanical properties of different phases of the steel sheet based on their
individual plastic deformation and slip systems. The macroscopic behavior of the polycrystalline
aggregate was then predicted based on the volume-averaged response of the representative phases,
and their volume fractions in the steel sheet.
The capability of the multiphase CPFE model was investigated based on experimental results
obtained from a quenched and partitioned Q&P980 steel sheet. The results from the CPFE model
were evaluated based on microscale uniaxial tension tests, and macroscale formability and
stamping tests. In all of the CPFE simulations, a random texture distribution was assumed for each
grain. However, to introduce the effect of inhomogeneity, two different methods were used to
model the initial volume fraction distribution of each phase. Inhomogeneity of the sheet metal
related to inhomogeneous distribution of phases and using a groove region in the specimen based
on the Marciniak-Kuczynski (MK) theory was also considered for the calculation of the forming
limit curve (FLC) using the CPFE model (MK-CPFE). The numerical study of the effect of
mechanical properties of constituent phases on the FLC suggests that it is possible to design
forming limit curves that would lead to textures that could enhance the resistance to localization.
Furthermore, the accurate prediction of strains at the onset of necking of a stamped part suggests
that a microstructure-sensitive model could be highly desirable for the stamping simulation of the
multiphase, third-generation advanced high strength steel (3GAHSS).

iii

This dissertation is dedicated to my daughters

iv

ACKNOWLEDGEMENTS

First, I want to express my sincere gratitude to my advisor, Dr. Farhang Pourboghrat for years of
support, leadership, assistance, and guidance during my PhD study. I would also like to thank the
members of my PhD committee; Dr. Thomas Bieler, Dr. Ronald Averill, and Dr. Thomas Pence for
their constructive criticism and helpful suggestions throughout my PhD research.
Thanks are due to the research members in our group, including Dr. Taejoon Park, Dr. Hyunki Kim,
Dr. Aboozar Mapar and Rasoul Esmaeilpour for their valuable suggestions during this work. My
special thanks are reserved for Dr. Fadi Abu-Farha and Dr. Jun Hu from Clemson University for
the successful collaborated work throughout the last three years.
Additionally, I would like to thank the Higher Committee for Education Development in Iraq
(HCED) for supporting my study. Furthermore, many thanks to the U.S. Department of Energy for
giving me the opportunity to participate in the “ICME Approach to Development of Lightweight
3GAHSS Vehicle Assembly” project.

Last, but not least, I would like to express sincere gratitude to my family for their persistence,
encouragement and support over the rough years. My warmest thanks to my parents for their
patience and everything they have done for me to become what I am today. They believed in me
and made me believe in myself during some painful times. Also, I would like to especially thank
my wife, without whose support and sacrifices I would not be here.

v

TABLE OF CONTENTS

LIST OF TABLES ....................................................................................................................... viii
LIST OF FIGURES ....................................................................................................................... ix
CHAPTER 1 INTRODUCTION ................................................................................................. 1
1.1 Background of AHSSs .........................................................................................................3
1.2 Why AHSSs .........................................................................................................................5
1.3 Generations of AHSSs .......................................................................................................10
1.4 Compositions and Microstructures of AHSSs ...................................................................12
1.5 Quenching and Partitioning (Q&P) Process Development to Replace Hot Stamping of
000000High-Strength Automotive Steel ........................................................................................15
1.6 Objective/Sub–Objectives of the Dissertation ...................................................................17
1.7 Overview of Dissertation ...................................................................................................18
CHAPTER 2 LITERATURE REVIEW..................................................................................... 21
2.1 Constitutive Models ...........................................................................................................21
2.1.1 Phenomenological models ........................................................................................22
2.1.2 Crystal Plasticity Models ..........................................................................................29
2.1.2.1 Crystal Plasticity Finite Element Methods (CPFEMs) ..................................32
2.1.2.2 CPFEM Numerical Integrations: ...................................................................33
2.2 Computational Modeling of AHSSs ..................................................................................36
2.2.1 Phenomenological Modeling of AHSSs ...................................................................36
2.2.1.1 Phenomenological Macroscale Modeling of AHSSs .....................................36
2.2.1.2 Phenomenological Microscale Modeling of AHSSs......................................40
2.2.2 Crystal Plasticity Modeling of AHSSs ......................................................................43
2.2.2.1 Macroscale Crystal Plasticity Modeling ........................................................43
2.2.2.2 Microscale Crystal Plasticity Modeling of AHSSs ........................................45
2.3 Failure and Forming Limits ...............................................................................................48
CHAPTER 3 THE CRYSTAL PLASTICITY MODEL ............................................................ 56
3.1 Crystal Plasticity Constitutive Model ................................................................................57
3.1.1 The Crystal Plasticity Yield Function .......................................................................57
3.1.2 Hardening ..................................................................................................................60
3.1.3 Homogenization ........................................................................................................61
3.1.4 Crystal Orientation ....................................................................................................62
3.2 Numerical Formulations ....................................................................................................64
3.3 Stress Integration Algorithm for Combined Constrains Single Crystal Plasticity Model .68
CHAPTER 4 MATERIAL CHARACTERIZATION ................................................................ 71
4.1 Q&P980 Microstructure.....................................................................................................71
4.2 Uniaxial Tension ................................................................................................................76
4.2.1 Experimental Setups .................................................................................................76
vi

4.2.2 Mechanical Properties of Q&P980 under Uniaxial Tension .....................................77
4.3 High Energy X-Ray Diffraction (HEXRD) Uniaxial Test .................................................81
4.3.1 Calibration of the Q&P980 Constituents Phases Based on CPFE Model ................81
CHAPTER 5 EXPERIMENTAL AND CRYSTAL PLASTICITY MODELING OF THE
HEMESPHERICAL PUNCH TESTS .......................................................................................... 87
5.1 Experimental Procedure of the Hemispherical Punch Test (Nakajima Test) .....................87
5.2 Numerical Simulations of the Nakajima Hemispherical Punch Stretching Tests ..............90
5.3 Assumptions for Implementation of Multiphase CPFE Model..........................................92
5.4 CPFE Model Verification and Discussions ........................................................................96
5.4.1 Punch Force-Displacement Analysis ........................................................................96
5.4.2 Strains Distribution Analysis ..................................................................................101
5.4.3 Identification of Necking Strains for FLC .............................................................. 110
CHAPTER 6 FORMING LIMIT CURVE PREDICTION BASED ON MK-CPFE MODEL 117
6.1 FLC Determination Procedures Based on Finite Element Simulation ............................ 118
6.2 The Finite Element Procedure and the Failure Criteria ................................................... 118
6.3 Results and Discussion ....................................................................................................125
6.3.1 Aspects Ratios Calibration with Analytical MK Model .........................................125
6.3.2 FLC Prediction Based on MK-CPFE ......................................................................131
CHAPTER 7 APPLICATION OF THE MULTIPHASE CPFE MODEL FOR SHEET METAL
STAMPING ................................................................................................................................ 142
7.1 T-Shape Stamping Experiment Setup ..............................................................................142
7.2 Assumptions for T-Shape Implementation .......................................................................144
7.3 FLC Validation and T-shape CPFE Simulation ...............................................................145
CHAPTER 8 THE ROLE OF CRYSTAL PLASTICITY MODELING IN CONTROLLING
THE DEFORMATION BEHAVIOR OF MULTIPHASE STEEL ............................................. 155
8.1 Representative Volume Element (RVE) Simulation (Hexahedral Element) ....................155
8.1.1 Effect of Phases Volume Fraction on the Stress-Strain Curve ................................159
8.1.2 Statistics of Stress and Strain Distributions ............................................................162
8.2 Microstructure Model RVE Simulation of Different Grains Sizes of Constituents
000000Phases ...............................................................................................................................168
8.2.1 RVE Construction and Meshing .............................................................................168
8.2.2 RVE Simulation Results..........................................................................................171
CHAPTER 9 CONCLUSIONS AND FUTURE WORKS ...................................................... 182
9.1 Conclusions ......................................................................................................................182
9.2 Recommendations for Future Work .................................................................................184
BIBLIOGRAPHY ........................................................................................................................186

vii

LIST OF TABLES

Table 2.1:

Mechanical parameters needed for the identification of several yield criteria
(Banabic, 2010a). ...................................................................................................... 28

Table 2.2:

Various yield criteria implemented in some commercial codes (Banabic, 2010a). .. 28

Table 4.1:

Chemical composition of Q&P980 sheet used in this study. .................................... 73

Table 4.2:

Mechanical properties of Q&P980. .......................................................................... 80

Table 4.3:

Slip systems for BCC crystals; mÎą and nÎą denote the slip direction and the slip
plane normal of the slip system Îą, respectively in the initial configuration. ............ 84

Table 4.4:

Slip systems for FCC crystals; mÎą and nÎą denote the slip direction and the slip
plane normal of the slip system Îą, respectively in the initial configuration. ............ 85

Table 4.5:

Material parameters for the CPFE model composed of four phases in Q&P980 steel
(unit: MPa). ............................................................................................................... 85

Table 8.1:

Material parameters for the Swift equation hardening model composed of four
phases in Q&P980 steel. ......................................................................................... 158

Table 8.2:

Material parameters for the CPFE model composed of four phases in Q&P980 steel
calibrated based on one orientation per gauss point (unit: MPa). ........................... 158

viii

LIST OF FIGURES

Figure 1.1:

Vehicle fuel economy improvement potential for various technologies (Stodolsky
et al., 1995). .............................................................................................................. 3

Figure 1.2:

Specific strength of various lightweight materials (Coates, 2012). .......................... 8

Figure 1.3:

The 2017 Kia Sportage body structure is 51% AHSSs (Courtesy of Kia Motors). .. 9

Figure 1.4:

Material average greenhouse gases emissions from primary production, in kg
CO2 e/kg material (Keeler et al., 2017). .................................................................... 9

Figure 1.5:

Overview of tensile strength and total elongation combinations for various classes
of conventional and advanced high strength sheet steel (AHSS) grades used in
automotive applications (Demeri, 2013; Matlock and Speer, 2009). ..................... 12

Figure 1.6:

Schematic show the composite phases of AHSSs (Keeler et al., 2017). (a) DP;
(b) TRIP; (c) CP; (d) MS; (e) FB; (f) TWIP. .......................................................... 14

Figure 1.7:

Schematic of the Q&P process for production of austenite-containing
microstructures starting with 100% austenite: Ci, CÎł and Cm respectively
represent the carbon contents of the initial alloy, austenite, and martensite
(Matlock et al., 2012). ............................................................................................. 16

Figure 2.1:

Twelve Slip system definition for FCC materials. .................................................. 31

Figure 2.2:

Geometry of slip in crystals (Hull and Bacon, 2001). ............................................. 31

Figure 2.3:

Schematics of forming limit diagram. ..................................................................... 55

Figure 3.1:

Diagram showing how rotation through the Euler angles φ1 ,, φ2 , in order 1, 2, 3
as shown describes the rotation between the specimen and crystal axes
(Engler and Randle, 2010). ..................................................................................... 64

Figure 4.1:

A representative micrograph of the Q&P980 microstructure showing the three
phases: ferrite (F), martensite (M) and retained austenite (RA). .............................. 73

Figure 4.2:

EBSD map and the associated inverse pole figures for both (a) BCC and (b) FCC
phases in the as-received steel. ............................................................................... 74

Figure 4.3:

Micrographs taken from cross sections of the Q&P980 sheet (a) perpendicular to
the transverse direction and (b) perpendicular to the rolling direction
(c) Microstructure of a long piece showing the size of the martensite
band (dark regions) in the microstructure. .............................................................. 75
ix

Figure 4.4:

Experimental setup for uniaxial tension testing. ..................................................... 77

Figure 4.5:

Uniaxial tensile stress-strain curves for the Q&P980 sheet steel with loading axis
aligned parallel, diagonal (DD) and normal (TD) to the rolling direction (RD). ... 79

Figure 4.6:

Optical microstructure examination of the adjacent areas reveals that strain bands
is due to microstructure inhomogeneity in multiphase AHSSs. ............................. 80

Figure 4.7:

Comparison of the experimental strain-stress curve with the CPFE model with
60 grain orientations assumed for each constituent phase and the Q&P980. ......... 86

Figure 4.8:

Comparison of the experimental strain-stress curve with the CPFE model in terms
of the number of grain orientations assumed per integration point. ....................... 86

Figure 5.1:

Bulge testing setups with DIC system. .................................................................... 88

Figure 5.2:

(a) Dimensional characteristics of the specimens used in the Nakajima tests; (b)
Principles and dimensions of the Nakajima test with a hemispherical punch;
(c) Nakajima punch. ................................................................................................ 89

Figure 5.3:

Formed samples after failure in the Nakajima test. ................................................. 89

Figure 5.4:

Schematic diagram for finite element model of a Nakajima test. ........................... 91

Figure 5.5:

The random distribution of phases orientations in the sheet metal elements (a)
180 mm (99310 elements), (b) 20 mm (36526 elements). ...................................... 95

Figure 5.6:

Comparison of punch force-displacement relationship computed with various
mesh sizes for 20 mm blank case. ........................................................................... 98

Figure 5.7:

A comparison of the experimental and numerically simulated punch force versus
punch displacement relationship during the Nakajima budging test with different
blanks sizes (a) 20 mm, (b) 60 mm, (c) 100 mm, (d) 180 mm. .............................. 99

Figure 5.8:
999999999
999999999
999999999
999999999
999999999
999999999
999999999
999999999

Comparison between the numerically simulated and experimental results for
rolling direction strain (Epsilon Y) of a 20 mm wide sheet at a punch depth of 20
mm. (a) CPFE with 140 grains with fixed distribution (b) CPFE with 140 grains
with random distribution (c) CPFE with 60 grains with fixed distribution (d) CPFE
with 60 grains with random distribution (e) von Mises (f) Hill48 (g) DIC. (Note: in
figure 5.8(g), the averaged minor and major strains histories at the fracture area are
calculated by DIC over a 2 mm diameter surface patch which is identified by the
red dot in each DIC image; colored lines are used to calculate the strain histories
along 60 mm section length at multiple locations) ............................................... 104

x

Figure 5.9:

Comparison between the numerically simulated and experimental results for the
rolling (Epsilon Y) and transvers (Epsilon X) strains of a 20 mm wide sheet at a
punch depth of 22 mm (a) CPFE with 140 grains with fixed distribution (b)
CPFEM with 140 grains with random distribution (c) CPFE with 60 grains with
fixed distribution (d) CPFE with 60 grains with random distribution (e) von Mises
(f) Hill48 (g) DIC. ................................................................................................. 105

Figure 5.10: Comparison between the numerically simulated and experimental results for the
rolling (Epsilon Y) and transvers (Epsilon X) strains of a 20 mm wide sheet at a
punch depth of 23.5 mm. (a) CPFEM with 140 grains with fixed distribution (b)
CPFE with 140 grains with random distribution (c) CPFE with 60 grains with
fixed distribution (d) CPFE with 60 grains with random distribution (e) von
Mises (f) Hill48 (g) DIC. ...................................................................................... 106
Figure 5.11: Comparison between the simulated and experimental results with blank width of
60 mm for rolling direction strain (Epsilon Y) at a punch depth of 24 mm (~1 mm
from experimental fracture point). (a) CPFE with 140 grains and fixed distribution
(b) CPFE with 140 grains and random distribution (c) von Mises (d) Hill48
(e) DIC. ................................................................................................................. 107
Figure 5.12: Comparison between the simulated and experimental results of the rolling (Epsilon
Y) and transvers (Epsilon X) strains of a 100 mm wide sheet at a punch depth of ~23
mm. (a) CPFE with 140 grains with fixed distribution (b) CPFE with 140 grains
with random distribution (c) von Mises (d) Hill48 (e) DIC. ................................ 108
Figure 5.13: Comparison between the simulated and experimental results for the rolling (Epsilon
Y) and transvers (Epsilon X) strains of a 180 mm wide sheet at a punch depth of ~30
mm (a) CPFE with 140 grains with fixed distribution (b) CPFE with 140 grains
with random distribution (c) von Mises (d) Hill48 (e) DIC. ................................ 109
Figure 5.14: Determination of instable necking strains using the linear best fitting (LBF)
method. ..................................................................................................................114
Figure 5.15: Identification of the necking area, 12 elements on the sheet’ surface of (a) 20 mm,
(b) 180 mm blank width. ........................................................................................114
Figure 5.16: Strain-based forming limit curve measured using the hemispherical dome
stretching Nakajima test (a) random phases distribution (b) fixed phases
distribution. ...........................................................................................................115
Figure 5.17: Comparison between experimental and simulated strain paths in balanced biaxial
888888888888case. .......................................................................................................................116

xi

Figure 6.1:

(a) Schematic of the unreformed finite element mesh, representing full and one
quarter of sheet blank. The sheet blank is subjected to biaxial straining by applying
proportional displacements, d1 and d2 and symmetric boundary conditions;
(b) Schematic of a sheet blank characteristics and dimensions; (c) A comparison
between the engineering and true strains histories on the prediction of FLC for
ρ=0.5. .................................................................................................................... 123

Figure 6.2:

Comparison between the predicted FLCs based on analytical MK model and
MK-FEM for different groove aspect ratios with fixed bulk aspect ratio (ARb =3/5):
(a) f0 =0.992; (b) f0 =0.994; (c) f0 =0.996. ............................................................. 128

Figure 6.3:

Bulk’s aspect ratio (ARb) effects on the of predicted MK-FEM forming limit curve
for fixed groove aspect ratio (ARg=4) and imperfection factor (f0 =0.992). .......... 130

Figure 6.4:

Comparison of the predicted MK-FEM and analytical MK model strain paths for
an imperfection factor f0 =0.992............................................................................ 130

Figure 6.5:

Comparison of the predicted FLCs by MK-CPFE (with f0 =0.992, ARg=4 and
ARb= 3/5) with experimental data: (a) 20 grains; (b) 40 grains; (c) 10 grains;
(d) Comparison of the fitted curves based on the number of grains. .................... 136

Figure 6.6:

Comparison of the predicted FLCs by MK-CPFE with different imperfection
factors. ................................................................................................................... 138

Figure 6.7:

(a) Initial random orientation; (b) Simulated textures of 2000 orientations at onset
of necking represented by (111) stereographic pole figure of an imperfection
factor f0 =0.992 and balanced biaxial tension ( ρ=1); (c) Equivalent plastic strain
contribution of the constituents phases at the onset of necking. ........................... 139

Figure 6.8:

FLC of Q&P980 predicted by MK-CPFE for the right- hand and LBF-CPFE of
the left-hand side. .................................................................................................. 141

Figure 7.1:

(a) The sheet metal stamping press tool components setup; (b) Dimensional
characteristics of the specimens used in the T-shape experiments. ...................... 143

Figure 7.2:

Schematic diagram for Finite element model of a T-shape stamping
simulation.............................................................................................................. 145

Figure 7.3:

CPFE simulation and experimental results of the T-shape stamping sample-1. (a)
Equivalent plastic strain distribution based on CPFE model; (b) Experimental final
shape results; (c) Comparison of punch force versus punch stroke relationship
during the stamping process.................................................................................. 148

Figure 7.4: The minor and major true strain distributions and FLC predicted by CPFE model
88888888888of T-shape sample-1 at the end of stamping process. ............................................ 150
xii

Figure 7.5:

CPFE simulation and experimental results of the T-shape stamping sample-2. (a)
Equivalent plastic strain distribution based on CPFE model; (b) Experimental
final shape results; (c) Comparison of punch force versus punch stroke
relationship during the stamping process. ............................................................. 151

Figure 7.6:

The minor and major true strain distributions and FLC predicted by CPFE
model of T-shape sample-2 at punch depth (a) 12.5 mm; (b) 13.5 mm and
(c) 14.0 mm. .......................................................................................................... 153

Figure 8.1:

(a) A representative volume element (RVE) of multiphase steel (Ferrite: 45%,
Martensite: 45%, New Martensite: 5%, Austenite: 5%); (b) Finite element meshes
and prescribed displacement boundary conditions for uniaxial tension
(U: translation, UR: rotation). ............................................................................... 157

Figure 8.2:

Macroscopic stress-strain responses of RVEs with varying the constituents phases
volume fraction. (a) CPFEM; (b) Isotropic FE. .................................................... 161

Figure 8.3:

Mises stress distribution of the RVE: (a) CPFE; (b) Isotropic FE. ....................... 164

Figure 8.4:

Max principle stress distribution of the RVE: (a) CPFE; (b) Isotropic FE. .......... 165

Figure 8.5:

Mises equivalent plastic strain distribution of the RVE: (a) CPFE; (b) Isotropic
FE. ......................................................................................................................... 166

Figure 8.6:

Stress triaxiality distribution of the RVE: (a) CPFE; (b) Isotropic FE. ................. 167

Figure 8.7:

Fitted log-normal grain size distribution used in DREAM.3D to generate equiaxed
grain structure. ...................................................................................................... 169

Figure 8.8:

A four-phase representative volume element of Q&P980 steel: (a) Synthetic
polycrystalline aggregate generated in DREAM.3D consisting of 3D voxels, and
generated from our integrated tool set; (b) full phase representation (Ferrite: 45%,
Martensite: 45%, New Martensite: 5%, Austenite: 5%); (b) Abaqus 3D RVE
undeformed meshes of the polycrystalline aggregate used in simulating uniaxial
tension. .................................................................................................................. 170

Figure 8.9: Simulation results of a four-phase representative volume element of Q&P980 steel
8……………at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution
8……………of maximum principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent
8……………plastic strain: (e) CPFE, (f) Isotropic FE. Distribution of stress triaxiality: (g)
8……………CPFE , (h) Isotropic FE. ....................................................................................... 172

xiii

Figure 8.10: Simulation results of a ferrite phase in the representative volume element of
88888888888Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic
88888888888FE. Distribution of maximum principle stress: (c) CPFE, (d) Isotropic FE.
88888888888Distribution of equivalent plastic strain: (e) CPFE, (f) Isotropic FE. Distribution
88888888888of stress triaxiality: (g) CPFE, (h) Isotropic FE. ................................................... 174
Figure 8.11: Simulation results of a martensite phase in the representative volume element of
Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic
FE. Distribution of maximum principle stress: (c) CPFE, (d) Isotropic FE.
Distribution of equivalent plastic strain: (e) CPFE, (f) Isotropic FE. Distribution
of stress triaxiality: (g) CPFE, (h) Isotropic FE. ................................................... 176
Figure 8.12: Simulation results of a new martensite phase in the representative volume element
of Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b)
Isotropic FE. Distribution of maximum principle stress: (c) CPFE, (d) Isotropic FE.
Distribution of equivalent plastic strain: (e) CPFE, (f) Isotropic FE. Distribution of
stress triaxiality: (g) CPFE, (h) Isotropic FE. ....................................................... 178
Figure 8.13: Simulation results of the austenite phase in the representative volume element of
Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic
FE. Distribution of maximum principle stress: (c) CPFE, (d) Isotropic FE.
Distribution of equivalent plastic strain: (e) CPFE, (f) Isotropic FE. Distribution
of stress triaxiality: (g) CPFE, (h) Isotropic FE. ................................................... 180

xiv

CHAPTER 1
INTRODUCTION
Numerous new functional materials have been developed in recent years. Nevertheless, steel
remains the most prosperous and cost-effective of all materials, and more than a billion tons are
used annually to enhance people’s quality of life and to continue to drive future progress.
Today’s vehicle programs must balance performance, safety, fuel efficiency, affordability, and
environmental concerns and still produce designs that are appealing to customers (Keeler and
Kimchi, 2014). Currently, the substitution of lightweight, high-strength materials for mild steel in
vehicle applications consists mostly of low-density materials (aluminum alloys, magnesium alloys,
fiber-reinforced composites, metal matrix composites) and high-density materials (high-strength
steels). Low-density functional materials still have a low effect on the overall dominant role of
steels, which have been established for more than a century in automotive manufacturing.
According to the U.S. Department of Energy’s (DOE’s) 2010 material technology report, with the
advent and development of new grades of high-strength steels, steel will still be the majority
material (at least 40%) in automotive components until 2035.
The estimated number of vehicles in the world exceeded 1.2 billion in 2014. According to
production statistics provided by the International Organization of Motor Vehicle Manufacturers
(OICA), the production of vehicles in 2013 exceeded 87.5 million vehicles, which was a 4%
increase over 2012, and about 95 million vehicles were produced in 2016, a 4.5% increase over
the production in 2015. Over the next five years, the production of vehicles will exceed 105 million
annually. With the increase in the number of vehicles produced globally, the demand for petroleum
all over the world will continue to increase. To achieve fuel-efficiency requirements, automakers
1

are developing new strategies and advanced technologies to improve engines, drivetrains,
transmissions, aerodynamics, tire rolling resistance, and vehicle weight. As shown in Figure 1.1,
reducing the weights of vehicles is the most effective means for improving fuel economy and
reducing energy consumption (Stodolsky et al., 1995). Reducing the weight of vehicles lowers the
inertial forces that the engine has to overcome to decelerate and stop them, and it also reduces the
power required to move and accelerate vehicles.
In an era of vehicle weight reduction, steel was not perceived as a lightweight or a high-tech
material because of its high density. This perception has changed with the introduction of new
grades of high-specific-strength steels, which offering thinner gauges to reduce vehicle weight.
The steel industry continues to develop new steel grades and innovative design concepts to meet
challenges related to vehicle performance, affordability, safety, fuel efficiency, and the
environment (Demeri, 2013).
The increased requirement by the federal government in the past few years on light weighting of
vehicle structures and improving the crashworthiness of automobiles have prompted the
application of advanced high strength steels (AHSSs). These AHSSs are characterized by their
multiphase structures. They have excellent mechanical, forming and even energy absorbing
properties.

2

Figure 1.1: Vehicle fuel economy improvement potential for various technologies (Stodolsky et al.,
1995).

1.1 Background of AHSSs
Steel has evolved over the years from mild steel in the early 1900s to high- strength, low- alloy
steels (HSLA) in the 1970s and to the first generation of AHSSs in the 1990s. In the 1970s, to
achieve the requirement to reduce the weight of vehicles as a result of the first and second world
oil crises, steel makers developed numerous grades of high-strength steels (HSSs) for the
automotive industry. Unfortunately, these HSSs, especially the HSLA and solid solution
strengthened (SSS) steels, turned out to be largely unsuccessful in the market due to their limited
improvement of strength and their impermissible decrease of ductility compared with conventional,
low-carbon steels (Horvath, 2010).
3

In 1994, Audi introduced to the market its A8 vehicle, the first model based on aluminum space
frame technology, and the dominant position of steels as materials for vehicle bodies began to
erode. The weight of the A8’s all-aluminum, body-in-white (BIW) vehicle, where BIW means that
the vehicle's sheet metal components have been welded together and they do not include moving
parts, e.g., the motor, trim, chassis parts, and paint, was considerably less than that of competing,
conventional cars with steel bodies. Consequently, in the same year, the world steel industry
reacted to the lightweight challenge by launching an international consortium of sheet steel
producers that comprised 35 companies from 18 countries. This consortium initiated the so-called
Ultra-Light Steel Auto Body (ULSAB) to explore opportunities for reducing the weight of cars by
using automotive components made of higher strength steels. The engineering efforts resulted in
the development of a steel car body that was 90% HSS. This car body had 25% less mass at 14%
less cost than the benchmarked sedans. In addition, the torsion and bending stiffness were
improved by 80% and 52%, respectively (Fonstein, 2015). Existing challenges coming from
carmakers and competitive materials resulted in intensive trials and implementation of earlier dualphase (DP) steels. By 1995, DP steels had been commercialized in Japan, the USA, and Europe
(Horvath, 2010).
Over the years since the ULSAB was introduced, the successes of AHSSs have motivated steel
companies to continue research on new types and grades of AHSS in order to bring these new
steels into the production process. The simultaneous development of new processes and equipment
to produce and form AHSSs has accelerated its use (Keeler et al., 2017).
More projects followed to demonstrate and communicate the capability of steel to meet demands
for increased safety and fuel efficiency through reducing the weights of various vehicle structures.
4

For example, the Ultra Light Steel Auto Closures (ULSAC) program mass produced efficient
hoods, doors, and deck lids from AHSSs. In conjunction with these new materials, the Ultra Light
Steel Auto Body Advanced Vehicle Concepts (ULSAB-AVC) program further refined fabrication
methods and achieved even greater weight reduction with AHSSs. In 2008, World Auto Steel, a
branch of the World Steel Association, began a program called Future Steel Vehicle (FSV), in
which steel members accelerated the development of new grades of AHSSs, further enhancing the
strength and ductility levels of the steels. The FSV program took automotive steel applications to
Giga-Pascal strength AHSSs and added the dimension of designing for reduced emissions over the
lifecycles of vehicles. Through engineering optimization, FSV achieved a 39% reduction in the
weight of the body structure (Keeler et al., 2017).
In 2013, The U.S. Department of Energy (DOE), in collaboration with multiple academic and
industrial partners, launched a project to develop third-generation AHSSs (3GAHSSs) for future
lightweight vehicles (Hector, 2013). At the same time, other 3GAHSS programs were
independently conducting related work, such as the projects related to nanosteel (Branagan, 2014).

1.2 Why AHSSs
Multiple factors drive the selection of materials for automotive applications, including safety, fuel
efficiency, environmental concerns, manufacturing concerns, durability, and quality. In the highlycompetitive automotive industry, cost also is extremely important in selecting materials. Cost is a
major driver for automakers, and reducing cost is the number top priority for the success of any
business. The cost of steels is always the most competitive in the market, and this is due mainly to
the tremendous magnitude of crude steel production worldwide. According to World Steel
Association, the global steel production reached 1.67 billion metric tons in 2014, far exceeding the
5

production of aluminum, which was in second place with 49 million metric tons of production,
followed by magnesium with 973,000 metric tons and carbon-fiber composites with 51,000 metric
tons. In the automotive industry, the most important lightweight material competing with steel is
aluminum alloys. Many studies have demonstrated the cost advantage of steel over aluminum
alloys, and the results have shown that the combined raw cost and conversion cost for aluminum
is three times that of steel, and the assembly cost of aluminum is 20 to 30 times that of steel. In
addition, the body structure of aluminum costs 60 to 80% more than steel (Demeri, 2013). Also, it
was estimated that, compared with steel, the unit cost of magnesium is 10 times higher, and the
cost of carbon-fiber composites is 14 times higher than steel (Fine and Roth, 2010).
Figure 1.2 shows that, with the introduction of AHSSs, the specific strength of those steels became
comparable to those of aluminum and magnesium alloys. Only the carbon fiber-reinforced
polymers have better specific strength values than steel, but their high cost and complicated
manufacturing preclude them from being considered for high-volume vehicle production (Coates,
2012, Demeri, 2013).
Although the densities of aluminum and magnesium alloys, at 2.7 and 1.8 g/cm3, respectively, are
lower than that of steel (~7.8 g/cm3), their yield strength and ductility combinations are much
lower than those of other grades of steel. Therefore, the weight reduction of AHSSs components
is achieved via increasing the strength of the material and decreasing the thicknesses of
components without compromising any safety and functional standards. In addition, due to the
high strength of AHSSs, exceeding 700 MPa, they are used extensively to manufacture strategic
structures for vehicles, especially the safety cage, to protect passengers from fatal injuries during
impact. For example, the structure of the 2017 Kia Sportage, Figure 1.3, was significantly
6

improved due to the extensive use of AHSSs. With 51% of the Sportage’s BIW consisting of AHSS,
versus the outgoing 2016 model with only 18%, torsional rigidity improved by 39%. As well, the
Sportage 2017 earned the highest safety designation possible (Top Safety Pick Plus (TSP+)) from
the Insurance Institute for Highway Safety (IIHS) when equipped with optional front crash
prevention systems. The rating reflects top scores in each of five crashworthiness tests as well as
the integration of driver assistance technologies to aid in crash prevention (Keeler et al., 2017).
In addition, AHSSs belong to the steel family, which consists of a wide range of grades, from highstrength up to 2.0 GPa for martensitic steels to high-ductility up to 80% for annealed austenitic
stainless steels. This provides more options to the designers of car bodies and also resists galvanic
corrosion when joining and assembling different steel components (Keeler et al., 2017). In addition,
AHSSs have their advantages all through the life cycle of a product, from the manufacturing of
the steel, through the mill producing the steel sheet to punching, bending, welding and assembling
the steel parts on the vehicle.
AHSSs are more environmentally-friendly than the low-density materials because the production
of these materials creates an offsetting emissions problem, i.e., their production is greenhouse gas
(GHG)-intensive and costly to the environment.
Figure 1.4 shows that the production of these low-density materials can result in 7 to 20 times
more emissions than the production of steel (World Auto Steel, 2011).
Another advantage is that steel is 100% recyclable. More than 92% of automotive steels can be
recycled without compromising their properties (Steel Recycle Institute, 2013).

7

Also, in forming sheet metal, in contrast to many aluminum and magnesium alloys, AHSSs panels
typically are cold-stamped to make automotive components, which further decreases the overall
cost of using AHSSs in automotive applications (Keeler et al., 2017).

Figure 1.2: Specific strength of various lightweight materials (Coates, 2012).

8

Figure 1.3: The 2017 Kia Sportage body structure is 51% AHSSs (Courtesy of Kia Motors).

Figure 1.4: Material average greenhouse gases emissions from primary production, in kg
CO2 e/kg material (Keeler et al., 2017).
9

1.3 Generations of AHSSs
AHSS steels are sometimes referred to by “generation,” and those steels in the “first generation”
include dual phase (DP), transformation induced plasticity (TRIP), complex-phase (CP), and
martensitic (MART) steels. The need to produce materials with significantly higher strengths has
led to the development of a “second generation” of AHSSs that consists of austenitic steels with
high manganese contents, which are closely related to conventional austenitic stainless steel, such
as twinning-induced plasticity (TWIP) steel and lightweight steel with induced plasticity (L-IP)
steels (Matlock and Speer, 2009).
Figure 1.5 shows the wide variety of deformation behaviors exhibited by first and second
generation AHSSs. Generally, the mechanical properties of the first generation follow the trend
formed by the conventional steels and high-strength steels (HSSs), which is higher strength and
lower ductility or vice versa. Therefore, formability problems are the main issue with first
generation AHSSs. The second generation of AHSSs possesses both high ductility and strength,
which are attractive in specific applications. These generations of AHSSs have used in car bodies
for such components as pillars, bumper beams, door rims, and floor cross-members (Keeler et al.,
2017).
For high production automotive applications, the second-generation AHSSs are attractive due to
their excellent formability. Second-generation AHSSs basically are austenitic steels. Austenite is
stabilized in the microstructure by significant additions of alloys, such as manganese (Mn) and
sometimes nickel (Ni), which are required to produce an austenitic microstructure. Despite their
advantages, second-generation AHSSs are in limited use because they have two negative aspects,
i.e., the high alloy content increases their cost, and they have a tendency for delayed cracking (i.e.,
10

the occurrence of fractures after the part is formed and stored) and poor weldability due to the high
alloy content (Nanda et al., 2016).
The third-generation of AHSSs (3GAHSSs) is the newest group of advanced steels, these steels
are still undergoing research and development (Demeri, 2013). The proposed development of
third-generation of AHSSs involves combinations of the properties of materials that bridge the gap
between the first and second generations of AHSSs (Figure 1.5). These future grades of steel are
intended to improve the ductility of the first generation of AHSSs and the affordability of the
second generation so that lightweight cars can be manufactured that provide improved fuel
economy and safety (Matlock and Speer, 2009).
The 3GAHSSs concept was pioneered by Speer et al. (2003) with the idea of utilizing a lowalloyed austenite-martensite dominant microstructure in 3GAHSSs rather than the conventional
ferrite-martensite microstructure that was dominant in HSSs and first-generation AHSSs. The basis
for this approach was that austenite would serve as the ductile constituent and martensite would
provide the desired strength. The possible strategies for producing 3GAHSSs were summarized
by De Moor et al. (2010), and they included 1) enhancing the properties of DP or TRIP steels, 2)
developing ultrafine bainitic microstructures, 3) the use of new fabrication approaches, such as
quenched-and-partitioned (Q&P) process, and 4) decreasing the Mn content in TWIP steels or
increasing the Mn content in TRIP steels. The promising 3GAHSS include, modified TRIP,
medium Mn (Fonstein, 2015; Hu et al., 2017) carbide free bainite (CFB) steels (Sugimoto et al.,
2002) and quenching and partitioning (Q&P) steels (Speer et al., 2003) are catalogs for the new
generation of AHSSs in the near future.

11

Figure 1.5: Overview of tensile strength and total elongation combinations for various classes of
conventional and advanced high strength sheet steel (AHSS) grades used in automotive
applications (Demeri, 2013; Matlock and Speer, 2009).

1.4 Compositions and Microstructures of AHSSs
Steel is an alloy that consists mostly of iron and carbon content between 0.0002% and 2.1% by
weight. Carbon is the most common alloying material for iron, but various other alloying elements,
such as manganese, chromium, nickel, molybdenum, and vanadium, are used in varying small
amounts. Some of the elements are residual, but others are added to impart useful properties, such
as strength, ductility, hardness, toughness, wear resistance, machinability, and weldability (Demeri,
2013).
12

The metallurgy associated with various grades of AHSSs and their processing are somewhat novel
compared to conventional steels. Their remarkable mechanical properties are related to the unique
processing used to produce them and to their structures. All AHSSs are produced by controlling
the chemistry and cooling rate from the austenite or austenite phase and the ferrite phase, either on
the runout table of the hot mill (for hot-rolled products) or in the cooling section of the continuous
annealing furnace (continuously annealed or hot-dip coated products). Research has provided
chemical and processing combinations that have created many additional grades and improved the
properties within each type of AHSSs.
The microstructure of AHSSs usually consists of four common phases, i.e., the austenite,
martensite, ferrite, and bainite phases. DP steels are a group of steels with a duplex microstructure
that consists of a soft ferrite matrix and a hard martensitic second phase in the form of islands
(Figure 1.6 (a)). TRIP steel contains 0.1 to 0.4 % C and other alloying elements, such as silicon,
aluminum, titanium, nickel, and vanadium. The microstructure of TRIP steels consists of austenite
embedded in a primary matrix of ferrite. In addition to a minimum of five volume percent of
retained austenite, hard phases, such as martensite and bainite, are present in varying amounts
(Figure 1.6(b)). The metastable retained austenite in TRIP steels transforms progressively to
martensite during plastic deformation. This combination of phases gives TRIP steels the high
formability of austenite during the initial stages of the stamping process and the high strength of
martensite at the end of the forming process.
Complex-Phase (CP) steels contain < 0.15% C and small amounts of martensite, retained austenite,
and pearlite in a ferrite-bainite matrix (Figure 1.6(c)). The martensitic steels (MSs) are
characterized by a martensitic matrix that contains small amounts of ferrite and/or bainite
13

(Figure 1.6(d)). Ferritic-bainitic (FB) steels have a microstructure of fine ferrite and bainite (Figure
1.6(e)). TWIP steels are austenitic steels with a high Mn content, i.e., 22 to 30%, that causes the
steel to be fully austenitic at room temperature (Figure 1.6(f)) (Keeler et al., 2017).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1.6: Schematic show the composite phases of AHSSs (Keeler et al., 2017). (a) DP; (b) TRIP;
(c) CP; (d) MS; (e) FB; (f) TWIP.

14

1.5 Quenching and Partitioning (Q&P) Process Development to Replace Hot Stamping
of High-Strength Automotive Steel
A novel steel processing concept called quenching and partitioning (Q&P) is being developed to
meet the demanding 3GAHSS cost and performance targets. The purpose of Q&P steel in the
context of automotive structures is to obtain a new type of ultrahigh-strength steel with good
ductility to improve fuel economy while promoting passenger safety. With a final microstructure
of ferrite, martensite, and retained austenite, Q&P steel exhibits an excellent combination of
strength and ductility. In general, the Q&P steels exhibit an ultrahigh strength in the range of 1000–
1400 MPa and a uniform elongation in the range of 10–20% (Wang and Speer, 2013). Also, steel
made via Q&P would save energy by avoiding high-temperature stamping. Replacement of hot
stamping would result in productivity gains for the steel industry by eliminating the time required
for heating the part before forming. By removing the need for a reheating furnace and dies,
automakers could reduce their energy consumption, as well as the footprint of their process lines
and their capital costs (U.S. Department of Energy, 2015).
In 2003, Speer et al. proposed an approach designated as the Q&P process to exploit novel
martensitic steels containing retained austenite (Q&P steel), based on the fact that carbon can
diffuse from supersaturated martensite into neighboring untransformed austenite and stabilize it to
room temperature. Since it first was proposed in 2003, Q&P steel has gained interest for its
potential to enhance steel’s strength and ductility properties with compositions similar to those of
transformation-induced plasticity (TRIP) steel, and it has been proposed as a third-generation
automotive steel.
Figure 1.7 shows the Q&P process, which consists of a two-step thermal treatment in which the
steel is quenched to a predetermined quench temperature (QT) in the martensite start temperature
15

to the martensite finish temperature range (Ms-Mf range) to produce a partially-martensitic,
partially-austenitic microstructure. The second step, the so-called partitioning step, aims at carbon
enrichment of the austenite by partial depleting the carbon in the martensite and transporting it to
the austenite. Thus, carbon stabilized austenite is retained in the microstructure after final
quenching to room temperature. Partitioning can be done at a higher temperature than the QT, i.e.,
the so-called two-step Q&P, or by holding at the quench temperature, i.e., the so-called one-step
(Matlock et al., 2012).

Figure 1.7: Schematic of the Q&P process for production of austenite-containing microstructures
starting with 100% austenite: Ci, CÎł and Cm respectively represent the carbon contents of the initial
alloy, austenite, and martensite (Matlock et al., 2012).

16

1.6 Objective/Sub–Objectives of the Dissertation
The aim of this dissertation was to provide an improved understanding of the modeling of
anisotropy in the elastic-plastic deformation behavior of multiphase steel from micro to macro
length scales by using a combination of microstructural and mechanical characterization and finite
element modeling approaches. Therefore, the development of a multiphase, computationallyefficient CPFE model is the prime objective and the unique contribution of the dissertation. This
prime objective includes the following specific sub–objectives:

1. To develop a model to study the mechanisms of deformation in the AHSS that idealizes the
constituent phases of the steel based on crystal plasticity constitutive relationships based on
crystallographic slip, which accounts for the initial anisotropy and for its gradual change due to
texture evolution during forming.
2. To integrate and implement the multiphase CPFE model with the commercial FE package,
ABAQUS, through a user materials subroutine, VUMAT. This will allow the user to conduct more
efficient CPFE simulations of multiphase materials at dramatically reduced computational cost
and with various selections for the properties of multiple phases. For this purpose, the following
tasks for the phases were established within the materials subroutine, i.e., (a) volume fraction; (b)
crystal structure (FCC or BCC); (c) volume fraction distribution (random, fixed, specified from
data file, and from EBSD); and (d) crystal orientations (random, specified from data file, and from
EBSD).
3. To examine the predictive capabilities of the multiphase CPFE model to capture the
inhomogeneous plastic deformation in complex microstructures that are commonly observed in
multiphase Q&P980 steel and to assess the relative influence of the distribution of the various
17

phases on the evolution of microstructural features for different loading conditions.
4. To focus on the prediction of the forming limit curve (FLC) of multiphase AHSSs using
microstructure-based CPFE simulations. The influence of inhomogeneities on the FLC of AHSS
sheet metal was studied in two categories in the framework of this research work, i.e., geometrical
and material inhomogeneities. The forming limits can be predicted by triggering strain localization
due to heterogeneous material in the hemispherical punch test (Nakajima test) or geometrical and
material heterogeneous in respect of Marciniak–Kuczynski (MK) model.
5. To investigate the capability of the multiscale and multiphase microstructure-based model to
capture the complex deformation behavior with a case study of stamping component-level finite
element forming.
6. To study the effects of the crystallographic orientations of the constituent phases on the
deformation and failure behaviors of multiphase steel during uniaxial tension.

1.7 Overview of Dissertation
The dissertation begins with an overall background of steels to justify the importance of pursuing
research of the new generations of AHSSs. The Introduction provides an overview of the literature
related to the historical development of AHSSs to justify the timing of this research. It also
provides an overview and information about the complex microstructure of AHSSs which defines
the dissertation principles. Also, the research objectives are presented in the introduction chapter.
Chapter 2 gives an overview of the available literature related to modelling the AHSSs. The chapter
begins with a brief background of the phenomenological and crystal plasticity models, and this is
18

followed by an overview of the extensive applications of these models in the recent years to model
the deformation behavior of AHSSs. The failure by forming limit curve also is described.
Chapter 3 introduces the numerical framework regarding the crystal plasticity finite element
method, and the homogenization method to simulate multiphase steels is presented. Also, the
formulations and the procedures of implementing the crystal plasticity constitutive model into the
user material subroutine (VUMAT) of ABAQUS software are outlined.
Chapter 4 presents the characterization of the materials and the calibration of parameters. The
results of the experimental measurements of microstructure, the results of the uniaxial tensile test,
and the results of high-energy X-ray diffraction (HEXRD) are presented and summarized. Also,
an approach is discussed for determining the parameters of the CPFE model along with the
HEXRD measurements.
Chapter 5 describes the experimental and numerical simulation procedures of the Nakajima
hemispherical punch test. The chapter focuses on the applicability of the multiphase CPFE model
to study the FLC and the effect of the inhomogeneity of materials on the predicted forming limits.
Chapter 6 describes the development of a methodology for simulating the FLC using the
multiphase CPFE model in conjunction with the Marciniak–Kuczynski (MK) approach. The
results of the forming limit curves simulation of a multiphase steel under various conditions are
presented and discussed.
Chapter 7 presents two case studies of component-level (T-shape) CPFE simulation and
experimental validation. Forming strains are validated with the FLCs predicted in Chapters 5
and 6.
19

Chapter 8 presents a series of numerical simulations of the representative volume element showing
the predictive capability of the proposed model.
Chapter 9 summarizes the results, provides conclusions, and makes some suggestions for future
research.

20

CHAPTER 2
LITERATURE REVIEW
The complex microstructures of AHSSs place a considerable demand on using advanced numerical
simulation analysis in the product design stage instead of the conventional experimental trial-anderror loops. Therefore, more attention has been paid to conducting accurate finite element (FE)
simulations of the sheet metal forming operations, one of the most extensively used production
processes in the automobile manufacturing industry. Currently, the accuracy of simulation results
is gaining credibility compared to experimental tests. Numerical simulation eliminates the need
for many experiments, and it reduces the time framework and manufacturing costs. In addition, it
provides very important information that is not readily available from experiments, such as local
stress fields and the initiation of damage.
The sheet metal forming simulations rely on an accurate description of the yielding behavior. In
general, sheet metals undergo severe plastic deformation during manufacturing processes, such as
cold rolling and stamping. This produces a preferential orientation to the grains and thus a variation
of mechanical properties at different orientations is to be expected. Due to the important impact of
deformation on the distribution of stresses and strains, the shapes of the formed parts are affected
by the anisotropic behavior of the material.

2.1 Constitutive Models
Over the last three decades, the application of FE simulations has undergone an increasing trend
in the various steps of the sheet forming processes, from design to testing. Consequently, further
improvements are required in terms of efficiency and accuracy. The accuracy and efficiency of
simulations have been improved dramatically due to the advent of fast computational resources
21

and parallel calculation techniques. However, due to the widespread application of AHSSs,
aluminum alloys, and complicated forming processes, improvements are still needed in realistic
material constitutive models, which are indispensable parts of any accurate simulation (Safaei,
2013).
The optimization of sheet metal forming processes requires in-depth knowledge of material
constitutive models and methods, so they can be implemented in user-friendly numerical tools. If
the discussion is restricted to polycrystalline materials, two different approaches can be used to
develop constitutive models. In the first approach, which is referred to as the phenomenological
approach, the average behavior of the material properties directly determines the global behavior
of the material. In the second approach, called crystal plasticity, the polycrystalline behavior is
described based on the behavior of each individual crystal. Without a doubt, both approaches have
numerous advantages.
2.1.1 Phenomenological models
Sheet metals exhibit either isotropic or anisotropic yielding behavior. An isotropic yield surface is
assigned to a material with identical mechanical properties at different orientations. Various
isotropic yield functions are available, such those provided by Tresca (1864), von Mises (1913),
Bishop and Hill (1951), Hershey (1954), Hosford (1972), Bassani (1977), Barlat and Richmond
(1987) and Budianski (1984). However, sheet metals generally exhibit a significant anisotropy of
mechanical properties. This is because they undergo severe plastic deformations during
manufacturing processes, such as cold rolling, which produces preferential orientations of the
microstructure. Consequently, the material obtains a direction-dependent mechanical behavior.
Material anisotropy highly affects the distribution of stresses and strains and, consequently, the
22

shapes of the final parts, their thicknesses, and possible instabilities, such as wrinkling for a deep
drawn part. With more anisotropic effects being observed in experiments, new anisotropic
phenomenological yield criteria were introduced. Starting from Hill’s quadratic anisotropy model,
various yield functions have been proposed to describe the initial anisotropy of metallic sheets.
In 1948, Hill added different parameters to the von Mises quadratic yield function and proposed
the well-known Hill48 model (Hill, 1948) as Eq. (2.1)
F ( yy   zz ) 2  G ( zz   xx ) 2  H ( xx   yy ) 2  2 L yz2  2 M  zx2  2 N xy2   2

(2.1)

In Eq. (2.1), the Hill48 quadratic yield function is represented by the six components of the stress
tensor (  ij ), which are defined based on the material coordinate system, i.e., x, y and z are the
rolling, transverse and thickness directions of the sheet. The term 

represents the effective

plastic stress. The terms F, G, H, L, M, and N are anisotropic coefficients, which can be determined
using the r-values or Lankford coefficients (Lankford et al., 1950) obtained from the uniaxial
tension experiments performed with tensile samples cut along three different directions (r0, r45,
and r90 along 0o, 45o, and 90o from the rolling direction, respectively).

F

r0
r0
(r  r )(2r45  1)
1
, G
, H
, N  0 90
r90 (1  r0 )
(1  r0 )
(1  r0 )
2r90 (1  r0 )

(2.2)

Usually, the material parameters for through-thickness shears are assumed to be same as the
isotropic von Mises yield function, i.e., L = M = 1.5 since results of those shear tests, which are
extremely difficult to perform, are not available.

23

The Hill48 model is still being used extensively in the industry, mainly due to its user-friendly
function and capability of representing the basic anisotropy. In addition, the Hill48 model
has been implemented in numerous finite element codes. However, this yield criterion cannot
represent materials that have severe anisotropy. Therefore, Hill made several improvements based
on the Hill48 model (Hill, 1993, 1990, 1979), as did Lin and Ding (1996) and Leacock (2006).
More sophisticated yield functions have been proposed to describe the initial anisotropy of metallic
materials. Balart and his coauthors proposed a series of ‘Yld’ criteria via substituting the stress
components in Hosford’s model by different stress transformation deviators. Barlat and Lian (1989)
proposed an anisotropic yield function (Yld89) based on a linear transformation of stress tensors.
A generalized version of the Yld89 yield function, i.e., the Yld91 yield function, considers threedimensional (3D) problems (Barlat et al., 1991). Barlat et al. (1997) proposed the Yld94 yield
function that improved the prediction of pure shear yield stress without affecting other in-plane
yield stresses. In order to improve the performance of the Yld94 function, Barlat et al. (1997a)
proposed another yield function, i.e., Yld96, for both three-dimensional and plane stress problems.
To overcome the relative convexity of the Yld96 function and to obtain a better prediction of inplane variation of yield stresses and Lankford coefficients, Barlat et al. (2003) proposed a new
two-dimensional (2D) plane stress yield function, i.e., Yld2000-2d. This model gained
considerable popularity mainly because of its accurate predictions of yield stresses and Lankford
coefficients at rolling, diagonal, and transverse directions but also because of its balanced biaxial
stress state. In 2005, Barlat and coworkers (Barlat et al., 2005) generalized the plane stress
Yld2000-2d yield function to consider six-component stress states for three-dimensional (3D)
problems yield function, and the generalized function was identified as Yld-2004.

24

Banabic and coauthors used the classical method of developing yield functions for adding weight
coefficients and extended Hershey’s formulation to a series of Banabic-Balan-Comsa (BBC) yield
criteria. Banabic et al. (2005) improved the plane stress BBC2000 yield function (Banabic et al.,
2000) to BBC2005, which takes an additional experimental input for parameter identification.
Comsa and Banabic (2008) proposed the BBC2008 yield function that shares some features with
Barlat’s Yld2004 function in that both require the same experimental data for parameter
identification and can be used in 3D and plane stress cases. However, the BBC2008 function is
based on more parameters and consequently requires more experimental inputs.
Cazacu and Barlat (2004) proposed a new yield function based on the linear transformation
approach to consider the effect of strength differential observed as yielding asymmetry in pressureinsensitive materials, such as hexagonal close-packed (HCP) metals, including titanium, textured
magnesium, and magnesium alloys.
In 1950, Hill proposed a general framework of plane stress yield functions based on a polynomial
expression (Hill, 1950). Polynomial functions were very attractive for many researchers due to
convenient calculation of derivatives that is needed for implementation into the finite element. The
basic idea of polynomial functions is to use a fourth order polynomial function to represent
complex anisotropic yielding (Gotoh, 1977). Moreover, some fourth-, sixth-, and even eighth-order
polynomial formulations were developed with 9, 16, and 25 model parameters, respectively (Soare
et al., 2008; Soare and Barlat, 2010). These models were developed originally for plane stress
problems, but they can be extended conveniently for 3D conditions.
The main advantage of using phenomenological models in the finite element method is the short
simulation times. However, this approach neglects the effect of texture and its evolution during the
25

deformation process. The fact that the same yield surface is used at every material point in the
specimen makes the predictions from this approach questionable in forming operations, which
involve complex deformation modes. This is because the yield surfaces in different regions of the
specimen are expected to evolve differently due to differences in texture evolution caused by the
difference in local deformation histories.
Phenomenological models can account for the effect of the initial plastic anisotropy, but they
neglect the effect of the evolution of anisotropy during plastic deformation, which has a prominent
role in large strain sheet metal forming operations. They still have some difficulties in accounting
for yield surface evolution, which requires the incorporation of complex equations to account for
the variation of anisotropy in the finite element simulations. In addition, the calibration of
phenomenological models is very arduous because many experiments are required, and they are
only valid for monotonic loading. For complex metal forming, assuming monotonic loading for
predicting limit strains is not adequate. Table 2.1 shows the mechanical parameters needed for the
identification of several yield criteria. On the basis of this list, we can estimate the amount of
experimental tests and the costs required for the identification of various yield criteria. The main
question of interest is whether a biaxial yield stress and biaxial anisotropy coefficient must be
determined since this requires a special apparatus, i.e., an apparatus for either cross tensile tests,
hydraulic bulge tests, or disk compression. In addition, Table 2.1 refers to plane-stress models. The
3D notations that are used in the table signify that the model is extendable to spatial stress states.
The yield criteria that belong to the Hershey family use an exponent chosen in accordance with
the crystallographic structure of the material (Banabic, 2010a).

26

Usually, the main criterion for selecting the yield function for implementation in FE code is the
precision of the prediction of anisotropic behavior. However, commercial software package for FE
simulation are not always equipped with the state-of-the-art material constitutive laws. Table 2.2
presents the main commercial FE software programs and the anisotropic yield criteria
implemented in them. The Yld89 and Yld2000-2d models have been used by various users in the
material subroutines of both ABAQUS and LS-DYNA commercial software programs (Banabic,
2010a).

27

Table 2.1: Mechanical parameters needed for the identification of several yield criteria (Banabic,
2010a).

Table 2.2: Various yield criteria implemented in some commercial codes (Banabic, 2010a).

28

2.1.2 Crystal Plasticity Models
The elastic–plastic deformation response of crystalline aggregates depends on the direction of
loading. This phenomenon is due to the anisotropy of the elastic tensor and to the orientation
dependence of the activation of the crystallographic deformation mechanisms, which take place
by dislocation slip on crystal planes. A consequence of crystalline anisotropy is that the associated
mechanical phenomena, such as changes in shape, crystallographic texture, strength, strain
hardening, and deformation-induced surface roughening and damage also are orientationdependent (Roters et al., 2010). In this context, physics-based crystal plasticity theories have
shown remarkable success in predicting both the anisotropic mechanical response of crystalline
materials and the evolution of the underlying texture in plastic deformation.
G. I. Taylor’s Bakerian Lectures at the Royal Society in 1923 formally initiated the systematic
study of the finite deformation of single cubic crystals. Taylor and Elam (1923) primarily
concentrated on the orientation of the crystal lattice. The directions of deformation in a single
crystal were determined by X-ray analysis and geometrical measurements. However, the concept
of crystal plasticity is based mainly on the pioneering work of Taylor (1938, 1934a, 1934b). He
observed that slip in FCC aluminum crystals occurred on octahedral {111} planes, parallel to the
<110> edges of each plane. An octahedron has four pairs of parallel planes, so FCC crystals have
four slip planes with three slip directions on each, meaning that a total of 12 slip systems are
available in each crystal (Figure 2.1). Then, Taylor resolved the shear stress on these planes parallel
to each direction and observed that the slip direction with the highest resolved shear stress
accommodates the deformation.

29

For slip to occur, a characteristic shear stress is required. When a mechanical load is applied on a
single crystal, deformation occurs by activation of possible slip system/systems in the crystal.
Which slip systems are activated is determined by the orientation of the applied stress and the
resultant resolved shear stress acting on the slip planes along the slip directions. For example, in
Figure 2.2 a crystal is deformed under tension by an applied force F along the axis of the cylinder.
Assuming the cross sectional area to be A, the stress developed parallel to the applied force is
σ = F/A. This force has a component in the slip direction FcosΝ with Ν being the angle between F
and the slip direction. This force component acts over the slip surface with an area A /cosϕ, where
angle ϕ is the rotation from F to the slip plane normal. Therefore, the resolved shear stress acting
on the slip plane along the slip direction can be represented as



F
cos  cos 
A

(2.3)

τ is the shear stress resolved on to a slip system. Furthermore, if a critical force Fc is required to
initiate slip, the corresponding shear stress is denoted τc or critical resolved shear stress (CRSS).
The translation cosλcosϕ from σ to τ is known as the Schmid factor.
Taylor proposed the shear stress that is required to move two parallel plane edge dislocations
separated by l past each other to be
  b / l

(2.4)

where Îą is a constant related to the strength of obstacles in the matrix, Îź is the shear modulus, and
b is the Burgers vector. The average dislocation spacing l is proportional to the inverse root of the
dislocation density.
30

Figure 2.1: Twelve Slip system definition for FCC materials.

Figure 2.2: Geometry of slip in crystals (Hull and Bacon, 2001).
31

2.1.2.1 Crystal Plasticity Finite Element Methods (CPFEMs)
Early methods to describe plastic anisotropy with simple boundary conditions, such as those
formulated by Taylor (1938) or Bishop and Hill (1951a, 1951b), were developed by simplifying
the assumptions of strain or stress homogeneity. They were not designed to consider explicitly the
mechanical interactions among crystals in a polycrystal or to respond to complex internal and
external loading boundary conditions. To overcome these issues, methods in the form of FE
approximations have attracted a lot of attention in the field of crystal plasticity. These crystal
plasticity finite-element (CPFE) models are based on the variational solution of equilibrium
between forces and the compatibility of displacements using a weak form of the principal of virtual
work within a finite-volume element. The CPFEMs use the extensive information gained from
experimental and theoretical studies of the deformation of single crystals and dislocation dynamics
to guide the development of more advanced continuum field theories. One main advantage of
CPFE models is their ability to solve crystal mechanical problems based on complicated internal
and/or external boundary conditions. This aspect is not a mere computational advantage; rather, it
is an inherent part of the physics of crystal mechanics, since it enables one to understand and use
the boundary conditions that are imposed by inter- and intra-grain micromechanical interactions
(Roters et al., 2010).
The first CPFE simulations were performed by Peirce et al. (1982). Due to computational
restrictions, they used a simplified setup of two symmetrical slip systems in order to study the
tensile behavior of a single crystal. Later, Harren and Asaro (1989) and Harren et al. (1988)
extended these simulations to a polycrystalline arrangement using a 2D representative volume
element (RVE) with two or three slip systems. In 1991, Becker was the first to perform simulations
on the basis of the 12 slip systems of an FCC crystal. Using a 3D model for the crystallographic
32

texture, Becker simulated the channel–die deformation of a columnar polycrystal aggregate
(Becker, 1991) and that of a single crystal (Becker et al., 1991). Since that time, ever-increasing
numbers of CPFE simulations have been performed, enabled by the increase in computational
power. In the field of direct or one-to-one crystal plasticity models, numerous grain-scale and
subgrain-scale problems have been investigated using meshes with subgrain resolutions and, in
part, complex 2D and 3D grain arrangements. Since its first use in finite element method, the
CPFEM has been expanded and enhanced into a whole family of constitutive and numerical
formulations, and they have been applied to a broad variety of crystal mechanical problems. In this
context, Roters et al. (2010) presented a list of more than 200 major contributions associated with
different applications of the CPFEM in the last two decades.
2.1.2.2 CPFEM Numerical Integrations:
Crystal plasticity numerical modeling remains a controversial subject, especially with the
development of CPFE codes. These codes generally are used to simulate the mechanical behavior
of metallic components and structures or to predict the overall behavior of crystal aggregates (by
using some homogenization techniques). Therefore, developing numerical integration algorithms
with improved accuracy, speed, and efficiency has become more desirable.
Two different modeling frameworks have been proposed in the literature for the evolution of slip
rates, i.e., rate-dependent and rate-independent approaches. Limiting the discussion to the rateindependent approaches, there have been two long-standing problems, i.e., 1) how to integrate the
set of highly non-linear constitutive equations and 2) how to resolve the ambiguity arising from
the non-uniqueness in the activity of slip systems (Ben Bettaieb et al., 2012).

33

Non-linearity has two distinct sources, i.e., 1) the behavior of the material, which is related to the
expression of the hardening law (when hardening is considered) and 2) the geometrical
configuration, which is related to the evolution of the rotation of the crystal lattice. Several
numerical algorithms have been developed and tested in order to solve this set of non-linear
equations. Typically, these numerical algorithms reach the solution incrementally and can be
classified into two main families, i.e., 1) explicit algorithms based on a forward Euler scheme and
2) implicit algorithms based on a backward Euler scheme. A third smaller, intermediate family that
exists between implicit and explicit algorithms is referred to as semi-implicit algorithms.
The development of explicit integration algorithms is straightforward, but it generally requires
very small times or loading increments to avoid numerical instabilities (Amirkhizi and NematNasser, 2007; Anand and Kothari, 1996; Knockaert et al., 2002, 2000; Maniatty et al., 1992; Peirce
et al., 1982; Rossiter et al., 2010; Zamiri et al., 2007; Zamiri and Pourboghrat, 2010). The critical
shear stresses and the rotation are assumed to be constant over the time increment, and their values
are chosen to be equal to their values at the beginning of the time step (Ben Bettaieb et al., 2012).
The implementation of implicit schemes is more complicated than the implementation of explicit
schemes. In these implicit algorithms, the critical shear stresses and the rotation of the crystal
lattice are unknown, producing a set of nonlinear equations. This set of equations must be solved
by traditional, iterative methods, such as the Newton–Raphson procedure or the fixed-point
method. In semi-implicit integration algorithms, the crystal lattice’s rotation is evaluated at the
beginning of the time increment. However, the computation of the critical shear stresses is based
on an implicit scheme. So, this algorithm class is implicit in the hardening rates and explicit with
regards to the increment of lattice rotation (McGinty and McDowell, 2006; Watanabe et al., 2010).
34

The second long-standing difficulty in the rate-independent numerical approach concerns the
method used to solve the well-known ambiguity of the non-uniqueness of the set of active slip
systems and the corresponding slip rates. The literature provides several contributions that
attempted to solve this difficulty.
Hamelin et al. (2011) developed a first-order selection criterion based on the minimization of the
change of the internal work with respect to the von Mises strain. Driver et al. (1984), Fortunier
and Driver (1987), and Skalli et al. (1985) simulated the behavior of single crystals of aluminum
using a method in which a combination of the active set of slip systems was selected based on the
fact that some combinations rotate towards orientations that require less dissipation of internal
work than others. Franciosi and Zaoui (1991) used a minimizing, hardening criterion in order to
select the set of active slip systems. However, several authors have used a rate-dependent
formulation based on power-type creep laws without differentiation of the slip systems into active
and inactive sets via a loading function (Asaro and Needleman, 1985; Mathur and Dawson, 1998;
Peirce et al., 1983, 1982). Anand and Kothari (1996) and Knockaert et al. (2000) used the singular
value decomposition technique when inverting the Jacobian matrix of the active slip systems.
Miehe and SchrĂśder (2001), SchrĂśder and Miehe (1997) suggested methods of general inversion
of the Jacobian of the active yield criterion functions in a reduced space. Gambin (1992), Guan et
al. (2006), and Montheillet et al.(1985) proposed a rate-independent approach based on the concept
of crystals with smooth yield surfaces. In this approach, only one yield function is used to calculate
the spin of the crystal and the incremental shear strains on the active slip systems. Therefore, the
ambiguity of the slip system in crystal plasticity is avoided. Zamiri and Pourboghrat (2010)
developed an effective methodology to define a smooth yield function based on the logarithms of
Schmid’s law for all slip systems and incorporating a regulation parameter.
35

2.2 Computational Modeling of AHSSs
This section presents the numerical approaches used to describe the mechanical behavior of AHSSs
at different length scales. From an engineering perspective, the high interest in simulating the
macroscale flow response of AHSSs is fueled by the requirement to adequately predict the
material’s behavior when exposed to sheet forming and crash loading. From a micromechanical
perspective, the interest is based on the substantial challenges associated with modeling highmechanical-contrast composite materials, such as the elastoplastic deformation of phases,
mechanics of interfaces, size effects, and damage.
2.2.1 Phenomenological Modeling of AHSSs
For many years, researchers have worked on modeling the deformation behavior of AHSSs using
different phenomenological models. They have been used different methods to capture their
macroscale and microscale behaviors.
2.2.1.1 Phenomenological Macroscale Modeling of AHSSs
Component-scale or forming FE simulations of AHSSs usually omit the crystalline nature of AHSS
to decrease the computation time while trying to predict certain metal-forming and crash-relevant
properties, such as strain hardening and forming limits. Many researchers have used this approach
to model the deformation behavior of AHSSs.
Banu et al. (2006) applied the quadratic Hill48 formulation in combination with the Swift law for
isotropic hardening to simulate the stamping process of a DP steel rail-shaped part. Tarigopula et
al. (2008) also used an isotropic yield surface, but with Voce hardening to study the localization in
DP steels under dynamic loading using digital image correlation and FE analysis. Tarigopula et al.,
(2008a) presented an approach of fitting parameters for yield surface approximation and
36

constitutive laws to measured data. The constitutive model includes isotropic and kinematic
hardening. Also Tarigopula et al. (2008c) studied fracture in a shear test and compared the findings
with experimental results obtained by digital image correlation (DIC) technique at a resolution
above the grain size.
Luo and Wierzbicki (2010) applied a developed Modified Mohr–Coulomb (MMC) ductile fracture
criterion to analyze the failure behavior of a DP780 steel sheet during stretch-bending operations.
The plasticity and ductile fracture of the sheet was were fully characterized fully by a the Hill48
model and an MMC fracture model. Hu et al. (2010) conducted an experimental study of cyclic
loading of DP and TRIP steels and a model that was derived to predict failure.
Instead of using a shear test, Hoon et al. (2011) reported how a displacement-controlled, drawbending test was developed, conducted, and simulated DP steel using a constitutive law that
accounts for thermomechanical effects, which some have claimed are important for the prediction
of shear failure. However, fracture was not modeled explicitly in this approach. Chung et al. (2011)
introduced a modified damage model to study the hole expansion formability of TWIP940 and
TRIP590 as AHSS grades and 340R as a high-strength grade. The damage model, along with the
anisotropic yield function Hill48 incorporated into the ABAQUS/Explicit FEM code, performed
well to predict hole expansion ratios and the sensitivity of the surface condition. Wu-rong et al.
(2011) compared the capabilities of three different yielding models (Hill48, YLD89 and BBC2005)
for describing certain averaged directional properties for four types of AHSSs.
Toros et al. (2012) used the FE analysis method to investigate and characterize the formability and
springback of TRIP800 steel based on three different material models (Hill48, YLD89, and
YLD2000-2D). Firat (2012) proposed an FE technique based on plastic instability and used Swift’s
37

diffuse necking and Hill’s localized necking concepts to predict the FLC of DP600 steel. To
compare simulated predictions with experimental results, Lim et al. (2012) used isotropic yielding
and the Hollomon hardening law to simulate the elastic springback for multiple types of AHSSs
subjected to sheet metal–forming operations.
Lou and Huh (2013) predicted the onset of ductile fracture by a newly-proposed ductile fracture
criterion of DP980 steel in various stress states, from shear to plane strain tension, where most
ductile fractures occur in sheet metal-forming processes. Kim et al. (2013) evaluated the sheet
formability of DP590 and TRIP590 steels by the limiting dome height (LDH) test and FE
simulations using three constitutive models (von Mises, Hill48, and Yld2000-2d). Either the Swift
or the Voce type hardening laws were used to fit the flow curves beyond the occurrence of
instability. Gruben et al. (2013) used a modified Mohr–Coulomb and the extended Cockcroft–
Latham fracture criteria in explicit FE simulations to study ductile crack propagation in a DP steel
sheet. Huh et al. (2013) studied the effect of strain rate on yielding anisotropy for two AHSS sheets
of TRIP590 and DP780 in addition to a conventional steel sheet of SPCC for comparison. They
compared the results based on three different yield descriptions. Qin et al. (2013) investigated the
effects of strain rate on the tensile properties of DP700 and DP500 commercial steels.
Chung et al. (2014) studied the formability of DP980 using FE simulation of the hemispherical
punch tests as well as the circular cup drawing test with isotropic hardening of the Hill48 yield
function and strain rate sensitivity. Lian et al. (2014) used the Lemaitre damage model that
accounts for the Lode angle effect with isotropic hardening to study the ductile fracture of DP600.
Panich et al. (2014) used different yield criteria (von Mises, Hill48 and Yld2000-2d) with the Swift
and Voce hardening models to model the plastic behavior of TRIP780 steel. Requena et al. (2014)
38

studied the evolution of damage in a DP980 steel. The experimental results are correlated with
those obtained by finite element analysis using a Gurson–Tvergaard–Needleman (GTN)
framework with a Chu and Needleman formulation to introduce the effect of the nucleation of
cavities. Habibi et al. (2015) investigated the formability of the automotive TWIP steel sheets
experimentally and numerically using various ductile/shear fracture criteria. Weiss et al. (2015)
studied the effect of volume fraction and hardness of martensite on the Bauschinger effect of five
different grades of DP steels.
Majidi et al. (2016) studied the springback behavior of DP780 and DP980 after different forming
conditions using experimental results and FE simulation. Panich et al. (2016) predicted both the
initiation of ductile damage and fracture of dual-phase steel grade DP780.
Anderson et al. (2017) implemented a hybrid experimental-numerical procedure to determine the
failure surface of a DP780 steel sheet as a function of the effective plastic strain, triaxiality, and
Lode parameter using butterfly specimens with in situ digital image correlation strain measurement
and supporting FE calculations. Erice et al. (2017) used the non-associated quadratic plasticity
model with Swift-Voce hardening and Johnson–Cook type of strain rate and temperature
dependency to model the large deformation response of three AHSSs (DP980, CP980, and
CP1180). Gruben et al. (2017) used an experimental-numerical approach to determine the strain
localization and ductile fracture of high-strength DP and martensitic steel sheet materials.
Charoensuk et al. (2017) predicted the ductile failure for DP780 and DP1000 by using a combined
approach of experiments and numerical simulations. FE simulations of the tensile tests were
conducted and the stress triaxialities, equivalent plastic strains, and Lode angles were evaluated
for critical areas.
39

2.2.1.2 Phenomenological Microscale Modeling of AHSSs
The development of a computational framework for predicting the mechanical performance
(overall stress–strain response) of AHSSs through modeling the AHSS’s their microstructure using
the virtually- generated, representative volume element (RVE) was the dominant approach for the
majority of researchers. Of the various AHSSs, DP steels, introduced in the mid-1970s, have
experienced the most rapid rate of utilization fastest growth in the automotive industry. Numerous
computational methods have been utilized to model deformation behavior and damage nucleation
of DP steels because of the simple two-phase microstructure that develops in these steels.
In the early works of Karlsson and SundstrĂśm (1974), the 2D-RVE was conducted to study the
deformation of two phase steel using FEM. They indicated that models that are developed to
investigate plastic deformation behavior in two phase metals use idealized assumptions that are
restricted to their models. Also, they concluded that very inhomogeneous strain distributions exist
in the two-phase structure and that the inhomogeneity is strongly dependent on the hardness ratio
of the phases and that develop at the very early stages of straining.
Tomota and Tamura, 1982 studied the deformation fields in different DP steels. They reported that
the degree inhomogeneity of the plastic deformation is influenced significantly by three factors
i.e., 1) the volume fraction of the martensite phase, 2) the yield stress ratio of the ferrite–martensite
phases, and 3) the shape of the martensite phase. The results of their FE simulation by a 2D-RVE
indicated that martensite was more plastically deformed and that the concentration of strain in the
ferrite matrix was reduced when the martensite grains formed a chain-like networked structure.

40

The generic building block–type approaches presented by Karlsson and Sundström (1974) and
Tomota and Tamura (1982) usually are derived from somewhat unrealistic assumptions regarding
the nonlinear phase and constitutive response of the interface, so they typically fail to reproduce
the underlying microscale mechanisms in metallic composite materials loaded into the plastic
regime. However, they showed that FE analysis could be used to examine local micromechanical
phenomena, but the determination of global properties was the main goal of their studies. In these
early works, the constitutive models were simple, e.g., an isotropic response was assumed for both
elastic and plastic regimes.
Improved isotropic models for multiphase metals have been developed over the last three decades,
and they provide better results. The composition of the steel, volume fraction, size, morphology,
distribution of phases, grain size, and boundary interface are some of the main parameters that
have been studied of AHSSs.
In most isotropic models, the flow curve of individual phases for ferrite and martensite is quantified
using the dislocation-based strain hardening model (Rodriguez and Gutierrez, 2003). This model
has been used extensively because of the simplicity of predicting the flow curves of phases based
on parameters that are dependent on the chemical composition of the material and the crystal size.
Many researchers have used real microstructures RVEs, which usually are determined from
scanning electron microscopy (SEM) images, to simulate the plastic flow and damage behavior of
AHSSs, especially for two-phase steels (Alaie et al., 2015; Alharbi et al., 2015; Asgari et al., 2009;
Han et al., 2014; Kadkhodapour et al., 2014, 2011; Latypov et al., 2016; Madej et al., 2014; MarviMashhadi et al., 2012; Paul, 2013a; Ramazani et al., 2014, 2013a, 2013b, 2012; Sirinakorn et al.,
2014; Sodjit and Uthaisangsuk, 2012; Vajragupta et al., 2012).
41

The mechanical properties of the constituent phases were determined by Sun et al. (2009a, 2009b)
and Choi et al. (2010, 2009) by using synchrotron-based in-situ high-energy X-ray diffraction
(HEXRD) experiments. They predicted the macroscopic stress–strain response of DP and TRIP
steels based on simulating real 2D-RVE obtained from SEM. Those authors used isotropic
elastoplastic constitutive formulations without considering the evolution of damage or the effect
of the crystallographic orientation on microstructure-level heterogeneity. Based on these
simulations, the main conclusion was that the softening and failure of steels is governed by plastic
strain localization that results from the incompatible deformation between the harder martensite
phase and the softer ferrite matrix rather than from void nucleation, growth, and coalescence.
A major disadvantage of using real microstructures for the determination of global properties is
the limited size of the material sample and its properties. Also, the FE simulations are still done in
2D to decrease computation time and, more importantly, because 3D data of real AHSSs
microstructures are difficult to obtain (Calcagnotto et al., 2010). For these reasons, the use of
idealized RVE (statistical models) has been promoted rather than the use of real microstructures
RVE for the determination of the global mechanical properties, such as plastic flow behavior and
influence of individual phase properties on the average material response (Abid et al., 2017, 2015;
Amirmaleki et al., 2016; Delannay et al., 2007; Govik et al., 2014; Liedl et al., 2002; Paul, 2013b;
Paul and Kumar, 2012; Uthaisangsuk et al., 2011, 2009). However, all of the aforementioned
studies have used the classical (local) plasticity theory or local plasticity–damage theory for
conducting micromechanical simulations.

42

2.2.2 Crystal Plasticity Modeling of AHSSs
In general, crystal plasticity finite element models can also be classified as macroscale and
microscale, depending upon the scale of volume that they represent at a material point.
2.2.2.1 Macroscale Crystal Plasticity Modeling
Understanding AHSS sheet metal responses undergoing many deformation modes is important for
industrial metal forming processes, particularly for further improvement of product design as well
as for optimizing the forming methods. Currently, the stamping simulation of complex automotive
parts with standard material models available in commercial finite element analysis codes lack the
sophistication necessary to accurately predict the forming and failure behavior of multiphase
AHSSs. These simple empirical constitutive models do not properly account for the complex
microstructure of AHSS sheets, and the evolution of their texture during the deformation process
(Kraska et al., 2009). Alternatively, macroscale metal forming simulations can be performed using
microscale material models, e.g. single crystal plasticity theory, within a multiscale analysis
framework. This approach is more appropriate, since material anisotropy in sheet metal forming
is, in fact, strongly related to the crystallographic orientation distribution of material grains
(Tjahjanto et al., 2015). Therefore, applying CPFEM for AHSSs is highly recommended for the
reason that it bridges the gap between the polycrystalline texture and macroscopic mechanical
properties, and it permits for a more profound consideration of metal anisotropy in the stamping
process simulation (Kraska et al., 2009). Macroscale use of such promising multiphase steels in
structural applications strongly depends on the ability to simulate the rather complex mechanical
response on the continuum scale. Employment of a direct finite element analysis on the scale of
individual grains is an accurate means but computationally far too expensive and impractical to be
used in component engineering. Therefore, suitable methods to coarse-grain (or homogenize) the
43

plastic response from the microstructural scale to a composite response valid at a material point of
the continuum are indispensable (Tjahjanto et al., 2007).
Macroscale or so-called Large scale CPFEM typically refers to the modeling of more than one
crystal orientation assigned to one integration point. Usually, the number of grains assigned to one
integration point are large enough to calculate the average behavior of the polycrystalline materials
using a homogenization scheme, as well as the material anisotropy through texture evolution
during the plastic deformation (Roters et al., 2010).
Macroscale CPFEM has been widely used to simulate the forming of sheet metals and extruded
tubes, especially with face centered cubic (FCC) materials such as aluminum and copper alloys.
These metals are relatively simple to implement since they have a single-phase microstructure and
their initial texture can be measured and implemented into CPFEM, and the texture evolution can
be tracked using a single crystal constitutive model (Zamiri and Pourboghrat, 2010; Li et al., 2009;
Zamiri et al., 2007; Guan et al., 2006; Zhao et al., 2004; Grujicic and Batchu, 2002). Moreover,
many researchers have been applying CPFEM to model body centered cubic crystals (BCC) such
as single phase steels (Kim et al., 2012; Tikhovskiy et al., 2008; Raabe et al., 2005, Xie and
Nakamachi, 2002).
Macroscale application of CPFEM to sheet metal forming processes with more than one phase are
scarce in the literature. In a recent study by Tjahjanto et al. (2015), the authors investigated the
analysis of deep drawing of DP800 with only BCC structure (ferrite and martensite) using
previously developed homogenization scheme called relaxed grain cluster (RGC) and material
simulation platform DAMASK (Roters et al., 2012).

44

2.2.2.2 Microscale Crystal Plasticity Modeling of AHSSs
CPFEM microscale simulations were performed with the goal of developing an in-depth
understanding of the microstructures of AHSSs. Hence, their fidelity improves as more of the
underlying physical mechanisms are incorporated. Examples of the effects that matter at this scale
are textures, morphology, constitutive strain-hardening details, damage models and void formation.
Modeling studies of the behavior of AHSSs at the grain level are presented in this subsection. In
contrast to Section 2.2.1.2 when the consideration is given to the phase contrast between martensite
and ferrite, additional constitutive effects are provided by using crystal plasticity formulations to
correctly model the anisotropic crystalline deformation modes of AHSSs.
In the past few years, industrial demand for the use of AHSS to manufacture lightweight parts has
placed new demands on the development of more sophisticated CPFE models for metals with more
than one phase. Many of the newly-developed models are based on microscale CPFEM of DP steel.
In microscale CPFEM simulations, the size of the FE mesh is equal to or smaller than the grain
size. With such sub-grain resolutions, the heterogeneity of deformation within a single crystal can
be taken into account.
Chen et al. (2014) used a combination of micropillar compression tests and CPFE with 3D-RVE
and random orientation to determine the flow strength and strain rate partitioning in uniaxial
tension in two commercial DP sheet steels. The work by Ghassemi-Armaki et al. (2013) was
similar to the work of Chen et al. (2014) but it placed more focus on modeling the constitutive
response of fully-martensitic sheet steel M190. J. H. Kim et al. (2012) proposed a model to predict
the cyclic behavior of DP steels. They incorporated a simplified dislocation density model into the
crystal plasticity model. Verma et al. (2014) used a CPFE model to analyze the two-stage
45

deformation behavior of DP780 steel based on 3D-RVE. They predicted the cross-effect observed
during the two-stage deformation fairly accurately using the CP model.
The failure of DP steels based on microscale has been the dominant approach among some
researchers. Choi et al. (2013) used a real 2D microstructure of DP980 to study the effect of the
orientation of crystals on the deformation and failure behaviors of steel during uniaxial tension.
They found that the incorporation of crystallographic orientations into DP steel was important,
even though the mechanical contrast between ferrite and martensite clearly has a more dominant
role in strain partitioning. That effect was proved by Alharbi et al. (2015) when they compared
their isotropic simulation results of DP980 with those of Choi. Their results showed an increase of
14% for the maximum stress in martensite after using a CPFE model, while an increase of only 8%
was recorded for the maximum stress in the ferrite phase. Therefore, the stress values computed in
the microstructure of the two-phase steel investigated in their work are likely to be inaccurate given
the isotropic assumption made in their model. Also, they mentioned that, apart from the effect of
crystallographic orientation and given the lack of knowledge about the geometry of the
microstructure in the third dimension, the assumption of 2D plane stress conditions is another
source of uncertainty for the computed stress values even though this approach has been used in
most similar studies in the literature.
A similar work of Choi, coupled CPFE simulations to experiments, Woo et al. (2012) focused on
simple geometrical features to investigate the importance of ferrite crystallographic orientation of
DP980 steel. Tasan et al. (2014a, 2014b) investigated strain localization and damage in DP steels
by performing in-situ experiments and crystal plasticity simulations based on 2D real
microstructure RVE. The strain distribution was measured by microscopic-digital image
46

correlation techniques (ÂľDIC). The results showed that the characteristics of the DP microstructure
strongly affect the strain localization process. It was emphasized that larger ferrite grains deform
plastically earlier than smaller grains.
In addition to above studies, there are a few publications in which the small scale CPFEM is used
to model more than two phases. In one of these studies, Srivastava et al. (2015) investigate the
plastic deformation of Q&P steel experimentally using CPFEM. The stress-strain curve for ferrite
and martensite phases were experimentally obtained by performing micropillar compression tests
at different crystal orientations. Due to the small sizes of retained austenite grains, the authors
calculated phase properties and phase transformation behavior of this phase based on the
macroscopic stress-strain behavior of the steel, and the previously measured properties of the
ferrite and martensite phases using CPFEM. Although feasible, an approach based on
microstructure simulation is highly computationally expensive for simulation of the large
deformation of engineered components. Moreover, one must be careful with dimensional
measurement inaccuracies and surface defects incurred during the milling process when using the
micropillar compression results. In addition, since the deformation is free from neighboring grains
during micropillar compression test, the dislocations will disappear when they glide to the free
surfaces of the micropillar and the pileup mechanisms of dislocations may not occur in such tests
(Hu et al., 2016a). Therefore, the flow and work hardening behaviors of a single phase obtained
from such a compression test can be different from that of a grain surrounded by other phases and
grains in a polycrystalline material under tension test such as HEXRD.

47

2.3 Failure and Forming Limits
The sheet metal formability is a measure of its capability to deform plastically during a forming
process to produce a part with definite requirements on mechanics, dimension and appearance,
being mainly limited by the occurrence of flow localization or instability.
Forming limit is actually a very broad concept, since it refers to plastic deformation to a given
shape without defects. Nevertheless, the primary defect concern in sheet metal forming is tearing
which means a sheet has lost its load carrying capacity. Typically, tearing is preceded by localized
necking or so-called thickness/local necking, which is a manifestation of critical strain localization
and a confirmed stage of component strength reduction. Thus, onset of localized necking is usually
considered the forming limit in sheet metal stamping (Keeler and Kimchi, 2014). Forming limits
from multiple deformations modes are usually connected to identify the forming limit curve (FLC)
which is plotted on a forming limit diagram (FLD) in Figure 2.3.
FLC experiments are actually based on various strip and cup forming methods with different punch
and/or specimen geometry (Banabic, 2010a). The out-of-plane hemispherical dome stretching test
or so-called Nakajima test (Nakajima et al., 1971) is widely used to obtain FLC for sheet metals.
The simplicity of the tools, the simple shape of the specimens and the possibility of covering the
entire domain of the FLDs are the main advantages of Nakajima test. Specimens with different
widths are usually employed in this test to generate proportional strains such as uniaxial tension,
plane strain, and biaxial stretch. Hence, the strain distribution in the FLC is relatively uniform.
Moreover, the punch force-displacement relationship obtained from collective Nakajima tests
provide important information about the material behavior that is relevant to the stamping process.

48

Failure in sheet metal forming processes are often evaluated based on the onset of localized
necking, which is characterized by the FLC in the principal strain space. FLC separates the region
of uniform sheet deformation from the region of slightly greater deformation, where the sheet will
likely develop a local deformation instability or neck. The concept of FLC was pioneered by Keeler
(1965) for the tension-tension domain of the right hand side of the FLC and extended by Goodwin
(1968) to the tension-compression domain. The FLC concept has proved to be quite useful for
representing conditions at the onset of sheet necking, and now is a standard tool for characterizing
the material forming behavior. It has since been widely used as a standard tool to study the
formability of metal sheets.
Although the concept of FLC is simple to understand and apply, it has limited application. The
actual response and failure of a material during plastic deformation depends on a large number of
coupled effects, making experimental determination of failure non trivial and time consuming
(Serenelli et al., 2011). In addition, scatter in the experimental data for a given sheet is often so
large that researchers sometimes question the accuracy and precision with which the FLC has been
determined (Janssens et al., 2001). Experimental determination of the onset of diffuse necking
becomes more important when it comes to stamping AHSS sheets, since they are relatively new to
stamping engineers. Typically, there are limited number of FLC available for the control of
stamping failures in the press shop applications, as well as to designers to validate formability
predictions with FE analysis.
In view of the continuous effort to reduce development costs, considerable effort has been made
to develop reliable theoretical predictive models, based on the continuum theory of plasticity and
different instability criteria for the assessment of a sheet metal's formability. Early theoretical
49

models for the calculation of FLCs assuming homogeneous sheet metals were proposed by Swift
(1952) for describing diffuse necking and Hill (1952) for localized necking. Marciniak and
Kuczyński (1967) proposed the MK model based on the premise that sheet metals are nonhomogeneous from both geometrical and microstructural point of view. According to MK model,
an infinite sheet metal contains a local imperfection region, which is a band having an initial
thickness less than that of the sheet. These authors showed that the presence of such an
imperfection could lead to unstable deformation in the thinner region and subsequent localized
necking and failure. The MK model has been one of the most widely applied model for FLC
prediction, and well established for predicting sheet necking among researchers (Banabic, 2010b).
In the last five decades, authors have been implementing various yield criteria and hardening rules
in the MK model. A common observation is that the FLC predictions are very sensitive and depend
crucially on the constitutive equation being used (Barlat and Lian, 1989b). The quadratic Hill48
criterion (Hill, 1948) and the von Mises plasticity theory in conjunction with the MK model are
frequently used, however, they overpredict the limit strains for biaxial paths. Sowerby and Duncan
(1971) used Hill48 yield criterion and Swift type isotropic hardening law (Swift, 1952) to calculate
the FLC. They showed that the right-hand side of FLC increases with an increase in stain hardening
exponent (n-value), and strongly depends on the normal anisotropic parameter (r-value).
Several non-quadratic yield functions have been developed and subsequently implemented in the
computational models of the limit strains prediction. The right-hand side of the FLC using MK
analysis with the non-quadratic Hosford yield criterion (Hosford, 1972) was predicted by Graf and
Hosford (1993, 1994). The computed FLC was found to increase by increasing the n-value, and
was unaffected by the r-value. The predicted FLC, although still higher, was in better agreement
50

with the experimental FLC compared to the results obtained by using the Hill48 yield criterion.
Barlat et al. (1997) calculated the FLCs of one aluminum and two steel sheet materials using the
Yld94 yield function with reasonable prediction of FLC. Cao et al. (2000) predicted the FLC using
MK model together with K-B yield function developed by Karafillis and Boyce (1993) and
analyzed the effect of changing the strain-paths on the FLC. An algorithm proposed by Ávila and
Vieira (2003) was used to predict the right side of the FLC using the methodology proposed by
Marciniak and Kuczynski using five different yield criteria. They confirmed that the type of yield
criterion used in the analysis strongly influences the FLC. More recently, much valuable
contributions have been made to predict FLC curves using MK analysis with more accurate
phenomenological models (Stoughton and Zhu, 2004, Soare and Banabic, 2009 , Eyckens et al.,
2009, Ahmadi et al., 2009, Zadpoor et al., 2009, Eyckens et al., 2011, Panich et al., 2013, Zhang
et al., 2014, Eyckens et al., 2015).
In the MK model, instability is viewed as a process in which the strain state in a region of local
weakness evolves to one of plane strain as deformation proceeds uniformly in the non-defected
region (bulk). While this characteristic agrees with that observed for many sheet failures, the
analytical MK approach is based on a simple model. The MK defect in analytical model is infinite
in length, with no end effects taken into account; such defects do not exist in real sheet materials.
The analytical MK analysis, therefore, implicitly imposes nonphysical boundary conditions by
way of the uniform minor strains throughout the sheet sample. These two simplifications allow for
the derivation of closed form one-dimensional solutions (allow variations in field quantities in one
direction only) that MK sought, but overstate the impact of a real material defect (Wagoner et al.,
1989). Finite element modeling allows relaxation of these two restricting assumptions by solving
for nonuniform internal stresses and strains under prescribed displacement boundary conditions.
51

Therefore, researches have used phenomenological models and finite element method for the
numerical determination of limit strains within the framework of the MK model (MK-FEM).
Narasimhan and Wagoner (1991) used the FEM approach to predict the FLC of the sheet specimen
with finite defect in the middle of the sheet far from the applied boundary condition. They studied
numerous material aspects and their effects on the predicted FLC, as well as validated their results
with analytical MK model. However, with an imperfection factor less than 0.98, the strain paths
in the groove region showed an earlier deviation from the applied strain path of the bulk and
eventually tended to plane strain. This is in contrast with the MK model where the strain path of
the groove should be proportional at the beginning and the deviation could be achieved later. In
the MK analysis based finite element method (MK-FEM), in addition to the imperfection factor,
the aspect ratios of the initial groove and bulk have a significant influence on the rate of strain
localization and thus on the predicted FLCs. These effects were not taken into account by some
researchers in their MK-FEM model, and the imperfection factor was the only parameter used for
predicting FLC using different phenomenological models (Evangelista et al., 2002; Boogaard and
Hu, 2003; Zhang et al., 2010). Recently, Srivastava et al. (2016) predicted the FLC of two DP
steels based on the standard rate-dependent, isotropic plasticity model (von Mises), and MK-FEM.
The flow behavior of DP steel used for FLC prediction was calculated with 3D finite element RVE
simulation coupled with a microstructure based crystal plasticity model. In their study, the
imperfection factor and the aspect ratio of the groove were only determined from the measurements
of the surface roughness of the sheets using 3D and DIC. It was found that the prediction of FLC
via their method was in good agreement with experimental results with only limited strains range.
Beyond phenomenological constitutive applications, it is well known and widely recognized that
the crystallographic texture evolution strongly affects the FLC and the macroscopic anisotropy of
52

polycrystalline sheet metals. Research has shown that the localization of plastic flow is influenced
by deformation-induced textures and anisotropy (Asaro and Needleman, 1985). Thus, more
realistic constitutive models, often involving crystal plasticity, are preferable. These
microstructure-based models provide a framework to better understand the relation between flow
localization and the materials microstructure, and subsequently on forming limits. Issues such as
material anisotropy, shape of the yield-surface and crystal orientation can be directly addressed
within a microstructure-based model.
The earlier analytical polycrystalline model based Taylor homogenization scheme and analytical
MK model was proposed by Zhou and Neale (1995) to predict FLCs for annealed FCC metal sheet.
Many recent studies have employed crystal plasticity models for predicting FLCs for sheet metals
based on a single-phase assumption. It is well documented that factors such as the initial
imperfection, texture, hardening rate, temperature, and the choice of the constitutive model affect
the prediction of the FLC (Wu et al. (1997, 1998, 2004, 2005); Inal et al., 2005; Yoshida et al.,
2007; Signorelli et al., 2009; John Neil and Agnew, 2009;Yang et al., 2010; Serenelli et al., 2011;
Wang et al., 2011; Yoshida and Kuroda, 2012; Chiba et al., 2013; Kim et al., 2013; Mohammadi
et al., 2014; LĂŠvesque et al., 2015; Cyr et al., 2017).
The analytical prediction of FLC based on crystal plasticity and MK analysis with more than one
phase are scarce in the literature. In a recent studies by Schwindt et al. (2015a, 2015b) the authors
investigated formability prediction of DP780 with only BCC structure (ferrite and martensite). In
their study, due to the lack of experimental flow curves for ferrite and martensite, four simulated
uniaxial flow curves of the phases were calibrated with the experimental flow curve for DP780.
Also, because of the high computational cost of the Visco-Plastic Self-Consistent (VPSC) crystal53

plasticity scheme adopted in their work, only limited number of orientations were assumed.
However, in this study, it was noticed that even with the assumption of thousands of orientations
for multiple phases in steels, the FLCs predicted by MK model of multiphase steels show, to some
extent, sensitivity to the phase distribution.
The typical FLCs for cold rolled, Q&P980, DP780 and DP980 steels have been reported by Wang
and Speer (2013). The authors revealed that the formability of Q&P980 is superior to that of DP980
steel, and reaches similar levels as that of DP 780. Recently, the formability of three dual phase
DP980 steels and one Q&P980 steel was evaluated by Zhang et al. (2016) using hemispherical
punch tests and DIC technique. Inconsistent and scattered results were observed when they
compared the formability of the four AHSS sheets based on the commonly used FLC method.
According to their observations of strain distribution from hemispherical tests and metallurgical
analysis of the Q&P980 steel, the inconsistency and scatter in the formability was attributed to the
microstructural features of the sheet metal. They observed that strain valleys occur in the regions
with higher martensite volume fraction, and conversely strain peaks occur in regions with a lower
volume fraction of martensite. The authors proposed an analytical method based on mathematical
expectation and variance to evaluate the formability of sheet metals and correlate them to the
global strain distribution observed in DIC. Also, they concluded that traditional formability
evaluation methods might not be applicable to the newly developed complex materials with
multiphase.

54

Figure 2.3: Schematics of forming limit diagram.

55

CHAPTER 3
THE CRYSTAL PLASTICITY MODEL
To model the mechanical behavior of multiphase steel under complex loading, the computationally
efficient crystal plasticity model introduced by Zamiri and Pourboghrat (2010) was used. This
model uses a yield function to describe the elastic-plastic dislocation slip of a single crystal, and
was originally developed to model the plastic deformation of single phase FCC materials, such as
aluminum and copper alloys. To study the mechanical behavior of multiphase steels, the effect of
the body centered cubic (BCC) slip systems for ferrite, martensite and the new martensite which
is transformed from pre-existing austenite were added. The hardening evolution of each phase,
which represents the slip resistance of slip systems were calculated from measured lattice strain
distributions obtained from high energy X-ray diffraction (HEXRD). These curves were
subsequently fitted to the extended Voce hardening equation with four parameters in Eq. (3.11).
In the proposed model, for computational efficiency and due to the low volume fraction of the
newly formed martensite in Q&P980 steel, it was assumed that the phase transformation of
austenite can be ignored in the crystal plasticity model. Therefore, the newly formed martensite
was treated as an additional phase independent of ferrite, pre-existing martensite and austenite
phases, and its shear stress-shear strain hardening curve was obtained from the in-situ high energy
X-ray diffraction (HEXRD).
The mechanical response of materials using macroscale 2D CPFEM has been widely used to
simulate plane strain or plane stress deformation states (Zamiri and Pourboghrat, 2010; Li et al.,
2009; Zamiri et al., 2007; Guan et al., 2006). However, 2D models are not capable of predicting
strain patterns that are out of plane. Furthermore, it has been shown by Simha et al. (2008) that the
56

state of stress during localized deformation (necking/shear banding) is triaxial, which cannot be
modeled using a 2D approach. Also, Rossiter et al. (2013) showed that 2D CPFEM is inadequate
for accurate modeling of the state of strain and stress as well as many other phenomena in sheet
metal forming such as bending and surface roughness. Therefore, this study utilized threedimensional (3D) CPFEM to predict localized deformation after necking, which results in out-ofplane deformation in sheet metal.

3.1 Crystal Plasticity Constitutive Model
The formulation of the crystal plasticity model is described below to show the fundamental
relationship between slip systems, shear rates, and hardening parameters as a function of the rate
of plastic deformation.
3.1.1 The Crystal Plasticity Yield Function
A velocity gradient in the plastic deformation with respect to the material coordinate system can
be decomposed into a rate of deformation D P and a spin tensor W P as:
LP  DP  WP

(3.1)

An elastoplastic problem is usually defined as a constrained optimization problem aimed at finding
the optimum stress tensor and internal variables for a given strain increment. In such a problem,
the objective function is defined based on the principle of maximum dissipation. This function has
terms describing the incremental release of the elastic strain energy, and dissipation due to
incremental plastic work, and the constraint is the yield function (Simo and Hughes, 1998):
trial
1
trial
1

 Min [l (x)  (σ j 1  Σ) : C : (σ j 1  Σ)  (q  q j ) : E : (q  q j )]

(σ j 1 , q j 1 )  

Subjected to : f ( Σ, q)  0)





57

(3.2)

where ÎŁ is the design variable (stress tensor) to be found, q is a vector containing internal
variables such as strain hardening and kinematic hardening parameters that need to be found, C is
the material stiffness matrix, f (ÎŁ, q) is the yield function, and E is the so-called matrix of
generalized hardening moduli. One of the solutions to the above problem is the following equation
for the plastic rate of deformation:

DP  

f (σ, q)
σ

(3.3)

here,  is a Lagrange multiplier, or as it is called in plasticity theory the consistency parameter.
In a crystal plasticity problem, deformation is defined by yield functions for each slip system in a
crystal. Assuming the Schmid law is valid for plastic deformation of a single crystal, then for any
slip system a yield function can be defined as:

f (σ, q) 

| σ : P |

 y

1

(3.4)

The constraints of problem (3.2) can be combined and replaced by an equivalent single constraint
defined by Zamiri and Pourboghrat (2010) as:
   | σ : P |   
1 
N

f (σ, q)  ln  exp  
1  


   1
 m   y
 



where  is the slip system,

N

(3.5)

is the number of slip systems,  y is the critical resolved shear

stress (CRSS),  is a parameter to determine the closeness to the condition, | σ : P |   y , and m
is a parameter for a flexible yield function shape. Also, Pι is a matrix of the Schmid tensor and I
describes the orientation of a slip system, defined as:
58

1
1
P  [(I )  (I )T ]  [m  n  (m  n )T ]
2
2

(3.6)

where n  is a unit vector representing the normal to the slip plane, and m  is a unit vector
denoting the slip direction. The plastic deformation gradient can be expressed as:

N

D P     P

(3.7)

 1

where   are the accumulated slip rates. The plastic spin tensor, which represents the material
rotation due to slip, can be expressed as:

N

W P     ω

(3.8)

 1

where ω matrix is the anti-symmetric part of I , defined as:

1
1
ω  [(I )  (I )T ]  [m  n  (m  n )T ]
2
2

(3.9)

Using Eq. (3.3), (3.5), and (3.7), it can be shown that during the plastic deformation of a single
crystal, the slip rate on any slip system can be expressed by:

   σ : P




exp

1


 y
 m   y


  σ : P
N


m  exp  

1


 m   y
 1


sgn(σ : P )

  

59

(3.10)

where  represents the index of summation for slip systems, m and  are material parameters
that control the shape of the single-crystal yield surface, which have been shown to have a direct
relationship with the stacking fault energy (SFE) of the material. In the CPFE model, the material
parameter  =80 was found to be the best value to represent the anisotropy of the steel sheet in
terms of the shape of the yield function, and m was assumed to be equal to one.
3.1.2 Hardening
The generalized Schmid factor can be found for each slip system  by finding the scalar product
of the normalized stress tensor (using the Frobenius norm) and the Schmid matrix. In order to
define the resistance against shear stress, an extended Voce type strain hardening rule (Tome et al.,
1984) was adopted


 0   
 ()   0  1  1   1  exp     
 1 









(3.11)

where  is the total accumulated shear strain on a slip system  over the entire deformation,
and   is the instantaneous flow stress on the slip system as a function of  . The initial
deformation is characterized by  0 , and  0 is the hardening rate, while the asymptotic
hardening is characterized by  1 and 1 . The extended Voce law represents the exponential
hardening which consists of the initial critical resolved shear stress,  0 , and the work hardening
rate,  0 . The excess stress  1 can be obtained before the critical resolved shear stress reaches
saturation where the sum of  0 and  1 is the saturation stress at which the work hardening rate
approaches zero. The 1

represents a constant work hardening rate for the linear hardening part

of extended Voce-law. The accumulated total shear strain on all slip systems  , can be defined
as:
60

 N
   
  1


 dt


(3.12)

3.1.3 Homogenization
In crystal plasticity, the term “homogenization” refers to the transition between the microscale and
the macroscale. In the present study, the Taylor type homogenization scheme (Taylor, 1938) was
used for the calculation of the homogenized stress in case multiple grains were considered for each
integration point in the finite element simulations. According to the Taylor model, which is an
upper bound model, each grain of the polycrystal aggregate undergoes the same deformation and
deformation rate. Therefore, the deformation gradient in Eq. (3.7) in all grains are the same and
equal to the overall deformation gradient.

By assuming that grains have the same volume V c for simplicity, the homogenized stress, σ ,
can be calculated from the crystal stress σc :

σ 

Vc
V

σ

c



c

1
M

σ

c

(3.13)

c

where V is the total volume of grains, and M is the number of grains (Zamiri et al., 2007). For
multiphase materials, Eq. (3.13) can be rewritten in the form:

 1
σ   fi 
i 1
 Mi

σ
c

i
c





(3.14)

Where fi is the volume fraction of each phase. Eq. (3.14) is then used in the explicit algorithm
to calculate the macroscopic stress tensor from the crystal stress tensors for all phases.

61

3.1.4 Crystal Orientation
For the material texture representation, the Euler angles refer to three rotations that, when
performed in the correct sequence, transform the specimen coordinate system onto the crystal
coordinate system-in other words, specify the orientation matrix Q. There are several different
conventions for expressing the Euler angles. The most commonly used are those formulated by
Bunge, as shown in Figure 3.1 (Bunge, 1982, Engler and Randle, 2010). The rotations are as
follows: φ1 about the normal direction ND, transforming the transverse direction TD into TD՛ and
the rolling direction RD into RD՛.  about the axis RD՛ (in its new orientation). φ2 about ND″ (in
its new orientation).
where 1 ,  , 2 are the Euler angles (Bunge definition). The effect of the operation sequence of
these three rotations can be seen in Figure 3.1. Analytically, the three rotations are expressed as:

 cos 1
Q1    sin 1

 0

sin 1
cos 1
0

0
0


1 

0
0 
1

Q  0 cos  sin  


0  sin  cos  

 cos  2
Q2    sin  2

 0

62

sin  2
cos  2
0

(3.15)

(3.16)

0

0


1 

(3.17)

By multiplication of these three matrices in order, an expression is obtained, which links the
rotation matrix Q to define the orientation of a crystal by the three Euler angles:

 cos 1 cos 2  sin 1 sin 2 cos  sin 1 sin  2  cos 1 sin  2 cos  sin  2 sin  
Q (1 ,  ,  2 )  Q2 .Q .Q1    cos 1 cos  2  sin 1 sin  2 cos   sin 1 sin  2  cos 1 cos  2 cos  cos  2 sin  



sin 1 sin 
 cos 1 sin 
cos  
(3.18)
In the CPFE model, the grain orientation matrix Q can be updated with the following equation:

Qt0 t  Qt0  Q  [ I  Ωt0 t t ]Qt0

(3.19)

where t0 t is the increment of lattice spin rate tensor in the finite element coordinate system can
be defined as:

N

  W     ω

(3.20)

 1

where W is the total spin tensor. More details on updating crystal orientation can be found in
Guan et al. (2006) and Zamiri et al. (2007).

63

Figure 3.1: Diagram showing how rotation through the Euler angles φ1 ,, φ2 , in order 1, 2, 3 as
shown describes the rotation between the specimen and crystal axes (Engler and Randle, 2010).

3.2 Numerical Formulations
For the numerical formulation, the predictor–corrector scheme (Simo and Hughes 1998) was
employed for the derivation of the equations and numerical implementation of the rateindependent crystal plasticity model.

Assuming an elastic stress–strain relationship, for a given discrete strain increment,

ε n   ε n1  ε n  , the trial Cauchy stress is calculated by assuming purely elastic deformation,
preserving the state variables from the previous n-th discretized time step ( ε n  εen and
 kn10   n ),

σ kn10  σ n  σ nk 10  σ n  Ce  ε n

64

(3.21)

Here, the superscript, k, denotes the iteration number ( k  0 : trial step), the subscripts, n and n  1 ,
denote the discrete time step number. The (n+1)-th step is purely elastic if the following yield
condition is satisfied:
N
   σ kn10 : P
 




ln  exp
 1    0

k 0
 
  1
 m   y ( n )
 


1

(3.22)

Otherwise, the (n+1)-th step is elasto-plastic. Then, the accumulated total shear strain, n1 , and
the Cauchy stress, σ n 1 , should be iteratively updated until the consistency condition in Eq. (3.23)
is satisfied.

nk1  Tol

(3.23)

where,



k
n 1

N
   σ kn 1 : P
 




 ln  exp
 1  

k
 
  1
 m   y ( n )
 


1

(3.24)

Here, Tol is the tolerance, which is numerically an infinitesimal value chosen to guarantee the
solution accuracy.
The Lagrangian multiplier, nk1 , the accumulated total shear strain,  kn 1 , and Cauchy stress,
σ kn 1 , are iteratively updated until the consistency condition in Eq. (3.23) is satisfied. The increment

of the Lagrangian multiplier is calculated by

65

 ( )kn 1  nk11  nk1  

 kn 1

(3.25)

  


  n 1
k

  
The denominator, 
 , can be linearized and decomposed by the chain rule:
  n 1
k

  
    σ 
     

 
 
 
 

  n 1  σ  n 1    n 1    n 1    n 1
k

k

k

k

k

(3.26)

Here,

 N sign(σ : P )
   σ : P




  P
exp

1


k
k



m


 
y
y
  1
  
 f (σ ) 









  σ : P
N

 σ  n 1  σ  n 1 



m
exp

1


  y

m



1

 


k







n 1

(3.27)

and,

 
  

   
   n 1   y
k

k

   y
 
n 1  

k



 n 1

(3.28)

where

k

 σ : P
   σ : P
 



 
exp

1

k
k
 2


m



y

  
 f (σ ) 
  y
  
 

    

 

   σ : P
N
 
  y  n 1   y  n 1 
 1  
 m exp  


m




1
y

 
 

n 1

66

(3.29)

Also, from Eq. (3.11),

  y

 

k

 

  1
n 1 

k


  0      0 
  0    

1

exp







exp
     



1

  1



1
 1

1
 n 1


(3.30)

The following relationship is derived under the associated flow rule:

d   d 

f (σ)
 y

(3.31)

Therefore,

k

k
 f (σ ) 
      
  
 k

  sign   n 1 


  
 



  n 1    n 1    n 1
y  n 1


(3.32)

 σ 
 e f (σ) 

  C 

σ n 1
  n 1


(3.33)

k

k

and,

k

k

which is derived from the following linear elastic constitutive law:
f (σ ) 

d σ  C e  d ε e  C e  ( d ε  d ε p )  C e   dε  d 

σ 


(3.34)

Then, the accumulated slip rate is updated by

  
k

  ( )n 1
  n 1
k



k 1
n 1

 

k
n 1

  ()

k
n 1

 

67

k
n 1

(3.35)

and the plastic slip rate is updated by

  kn 11    nk 1   (  ) kn 1    kn 1   ( ) kn 1

f (σ )
 

(3.36)

where

k

 sign(σ : P )
   σ : P
 



 
exp

1

  y
 

m


y
f (σ ) 

 






  σ : P
N



m exp  
 1 




 m   y
 1

n 1


(3.37)

Finally, from Eq. (3.33) the Cauchy stress can be updated by:

 σ 
k

  ( )n 1
  n 1
k

σ

k 1
n 1

σ

k
n 1

 σ

k
n 1

σ

k
n 1

(3.38)

The above formulations were used in the implementation of polycrystal CPFE model into the
commercial Abaqus/Explicit finite element code with aid of the user-defined material subroutine
VUMAT, which allows the user to define the mechanical behavior of a material and to interface
with any externally defined programs.

3.3 Stress Integration Algorithm for Combined Constrains Single Crystal Plasticity
Model
Summarizes for the algorithm that is used to find the admissible stress tensor for stress integration.
Here, the n subscript indicates a finite element increment and the k superscript refers to a plasticity
return mapping iteration as shown in the previous section. At each integration point the finite
element solver calls the plasticity algorithm. This scheme then calculates the admissible stress
68

tensor for the integration point and updates the solver. The details of this algorithm are given below.
1. The σˆ n , n , ˆ, Qn , qn variables are known at a finite element increment n for a
crystal in the element coordinate system. σ̂ n is the stress in the element coordinate
system, n is the accumulated shear strain, ˆn increment of strain in the element
coordinate system, Qn is the crystal orientation matrix and qn are the internal
variables. The stress and strain tensors are rotated into the crystal coordinate system.
σ n  Qnσˆ nQnT

 n  Qn ˆn QnT

2. Initialize variables.

 nk10   n  Ce :  n
qnk10  qn

 kn10   n
nk10  0

3. Evaluate the yielding condition.



k
n 1

N
   σ kn 1 : P
 


 ln  exp    k  1  
 
  1
 m   y ( n )
 


1

69

If

 kn1  0

go to step 9.

4. Compute the plastic consistency parameter.

 ( )kn 1  

 kn 1
k

k
k
 
    σ 

 
   
 σ n 1   n 1   y

   y
 
n 1  

k

    k
 

n 1    n 1

5. Update the stress, accumulated shear strain, and the shear strain on the slip system.

 σ 
k

  ( ) n 1
  n 1
k

σ

k 1
n 1

σ

k
n 1

 σ

k
n 1

σ

k
n 1

  
k

  ( )n 1
  n 1
k



k 1
n 1

 

k
n 1

  ()

k
n 1

 

k
n 1

  kn 11    nk 1   (  ) kn 1    kn 1   ( ) kn 1

f (σ )
 

6. If  kn 11  Tol go to step 8, otherwise go to step 7.
7. Update the consistency parameter.
nk1  nk11   ( ) kn 1
k  k 1

then go to 3.
8. Update the internal variables such as the crystal orientation using Eq. (3.19).
9. Rotate the stress tensor back to the element coordinate system.

70

CHAPTER 4
MATERIAL CHARACTERIZATION
In this Chapter, the target multiphase Q&P980 steel was investigated experimentally on its
mechanical properties based on different aspects, including microstructure and uniaxial tension
tests. Also, an approach to determine the CPFE model parameters along with as-received HEXRD
measurements is presented.

4.1 Q&P980 Microstructure
A commercial Q&P980 steel sheet (~1 mm thickness) from Baosteel was investigated in the asreceived condition. The chemical compositions of the steel sheet are detailed in Table 4.1. The
Q&P980 microstructure consists mainly of ferrite, martensite and retained austenite (Figure 4.1).
The martensitic phase, formed after quenching process, has a body-centered tetragonal (BCT)
structure with a unit cell that is a distortion of the body centered cubic (BCC) ferrite phase, a
consequence of the anisotropic carbon distribution. The reduction of carbon in the martensite
matrix during partitioning process alters the lattice parameter of martensite to make the BCT
crystal structure closer to the cubic ferrite (BCC), the equilibrium phase at the partitioning
temperature. These changes render the martensitic phase softer after partitioning relative to the
martensite after quenching.
Experimental microstructural observations of the Q&P980 steel sheet under consideration revealed
that this material consists of 44±2% by volume of the martensite phase, ~8–10% by volume of the
retained austenite phase, and the rest is ferrite.

71

The crystallographic texture of the steel was analyzed based on electron backscatter diffraction
(EBSD) measurements. Figure 4.2(a) and (b) show representative EBSD map for the steel phases
indexed as BCC for both martensite and ferrite, and FCC for retained austenite. The orientation
scatter of the phases in the inverse pole figures show that the Q&P980 sheet steel does not possess
a pronounced texture and almost appears to have a random crystallographic orientation.
Optical micrographs taken from the middle of the sheet for cross sections perpendicular to the
transverse direction (TD) and perpendicular to the rolling direction (RD) are shown in
Figure 4.3(a) and (b) respectively. The spatial distribution of martensite (dark regions) is not
uniform, and the microstructure exhibits martensite banding located in the mid-thickness of the
sheet and parallel to the rolling plane Figure 4.3(c). Such martensite morphology and
inhomogeneous distribution have a significant influence on the accumulation of strains in soft
phases such as ferrite (light regions) and subsequently on damage.

72

Table 4.1: Chemical composition of Q&P980 sheet used in this study.
Alloying Element

Mn

P

Al

Si

C

% Weight

1.8

0.02

0.05

1.5

0.2

Figure 4.1: A representative micrograph of the Q&P980 microstructure showing the three phases:
ferrite (F), martensite (M) and retained austenite (RA).

73

Figure 4.2: EBSD map and the associated inverse pole figures for both (a) BCC and (b) FCC
phases in the as-received steel.
74

Figure 4.3: Micrographs taken from cross sections of the Q&P980 sheet (a) perpendicular to the
transverse direction and (b) perpendicular to the rolling direction (c) Microstructure of a long piece
showing the size of the martensite band (dark regions) in the microstructure.

75

4.2 Uniaxial Tension
The objective of this section is to evaluate the macroscopic mechanical properties of Q&P980 steel
during macroscopic uniaxial tension required for phenomenological models (von Mises and
Hill48).
4.2.1 Experimental Setups
To determine the macroscopic flow behavior, specimens with a 50 mm gauge length and a flat
dog-bone geometry were deformed in uniaxial tension. The tension specimens were prepared by
water-jet cutting according to ASTM E8 (ASTM Int., 2009). The uniaxial tensile tests were
performed with the electro-mechanical INSTRON 5985 universal load frame, fitted with a 250 kN
load cell as shown in Figure 4.4. A 3D system DIC technique ARAMIS 5M was used to monitor
the deformation of the test specimens. For strain measurements using the DIC technique, one
surface of the specimen was painted with white paint and then black paint was sprinkled on it to
produce a random speckled pattern. The universal load frame and the DIC systems were triggered
simultaneously at the beginning of each test, thus all strains images were tagged by the
corresponding force during testing. All tests were performed at a constant crosshead speed of 10
mm/min, yielding an average strain rate of ~0.002 s-1, and the DIC cameras were set to capture
images at a fixed rate of 2 fps.

76

Figure 4.4: Experimental setup for uniaxial tension testing.

4.2.2 Mechanical Properties of Q&P980 under Uniaxial Tension
The mechanical properties of the Q&P980 were obtained from specimens along the rolling
direction (RD), the diagonal direction (DD) and the transverse direction (TD). Each
setup/configuration was tested a total of five times. The tensile test was conducted at a constant
speed for both the yield strength and r-values determination. The r-value, or the Lankford
coefficient (Lankford et al., 1950), reveals the directionality of plastic deformation of sheet metals,
which is induced by statistical orientation tendency of grains produced in rolling. It can be
identified as the ratio of the true plastic width strain to the true plastic thickness strain.

r

77

w
t

(4.1)

The uniaxial tensile stress–strain curves for the Q&P980 at room temperature with different
alignment to the rolling direction are shown in Figure 4.5. The measured mechanical properties of
Q&P980 (yield strength (YS), ultimate tensile strength (UTS) and elasticity constants) by the
uniaxial tensile tests are summarized in Table 4.2.
Here, the uniaxial tension data along the rolling direction was considered as the reference state for
characterization of the Swift type isotropic hardening law (Swift, 1952):

  K ( 0   )h

(4.2)

where  and  are the effective stress and equivalent plastic strain, respectively, while K,  0
and h are the hardening coefficients. The hardening coefficients were determined using the least
square method based on the experimental true stress-plastic strain curve along the RD. The
following material parameters have been obtained in this manner: K = 1413.64 MPa,

 0 = 2.1192×10-3 and h = 0.099.
Figure 4.5 also shows the Q&P980 sheet’s inherent microstructure inhomogeneity as bands
appearing on its surface detected by DIC along diagonal and transverse directions at 4% strain. It
can be seen that the rolling direction develops a preferred crystallographic orientation in the
microstructure that results in a higher strength. In this regard, mechanical properties tend to be
lower in the transverse direction which is perpendicular to the martensite band direction
(Figure 4.3).
From optical microstructure examination shown in Figure 4.6, in a strain valley, there are more
martensite grains, which are strong but not ductile, while in a strain peak, there are more ferrite
78

grains, which are comparatively more ductile. Such a microstructure inhomogeneity was primarily
caused in sheet metal rolling processes. It has also been reported by (Zhang et al., 2016) and
observed in many other sheet metals, especially those multiphase AHSSs.

Figure 4.5: Uniaxial tensile stress-strain curves for the Q&P980 sheet steel with loading axis
aligned parallel, diagonal (DD) and normal (TD) to the rolling direction (RD).

79

Table 4.2: Mechanical properties of Q&P980.
Material Elasticity Constants
C11 (GPa)

C12 (GPa)

C44 (GPa)

269.23

115.38

76.92

Direction

YS (MPa)

UTS (MPa)

r-value

RD (0o)

814.5

1019.3

0.8721

DD (45o)

778.4

1001.6

0.9465

772.3

1005.8

0.9357

o

TD (90 )

Figure 4.6: Optical microstructure examination of the adjacent areas reveals that strain bands is
due to microstructure inhomogeneity in multiphase AHSSs.

80

4.3 High Energy X-Ray Diffraction (HEXRD) Uniaxial Test
In-situ volumetric diffraction measurements of a sample in a tensile test with HEXRD have been
recently used to provide experimental data of the constituent phases of AHSSs as inputs to various
constitutive models (Cong et al., 2009; Hu et al., 2016b; Jia et al., 2009, 2008; Sun et al., 2009a).
The in-situ HEXRD test produces diffraction patterns that are indicative of changes in the steel
microstructure as a specimen is incrementally strained, such as in tensile testing (Jia et al., 2008).
In-situ HEXRD test under uniaxial tensile loading along rolling direction was performed with the
Q&P980 steel at Advanced Photon Source (APS) - Argonne National Laboratory by Pacific
Northwestern National Laboratory researchers (PNNL) to determine the mechanical properties of
the constituent phases as well as to obtain the austenite volume fraction evolution during the
deformation process. Through the detailed data analysis of the diffraction pattern during tensile
loading, the lattice strains of various crystal planes were first calculated as a function of the
macroscopic strain. The properties of the constituent phases in the Q&P980 were obtained using
the procedure described by Hu et al. (2016), which can be referred to for more details on the
calibration and determination of phase properties of the multiphase Q&P steels with in-situ
HEXRD. The as-received flow curves of the constituents phases are used in the next section to
calibrate with CPFE model.
4.3.1 Calibration of the Q&P980 Constituents Phases Based on CPFE Model
The as-received data from the in-situ HEXRD were used for the calibration of the CPFE
parameters of the hardening model. For the sake of computational efficiency, the CPFE model was
run with the reduced integration gauss point of just one at each element. The hardening parameters
for each phase were calibrated using a specimen with 50 mm gauge length and 1 mm thickness.
81

The specimen was meshed using the eight node C3D8R brick elements from the
ABAQUS/Explicit code. Separate calibrations were conducted based on the number of random
orientations (grains) assumed per integration point for each element.
The slip systems used for the deformation of BCC crystals (ferrite, old martensite, and new
martensite) were assumed to be (12 {110} <111> and 12 {112}<111> slip systems) (Table 4.3),
while the following 12 slip systems ({111} <110> slip systems) were considered for the FCC
crystal (retained austenite) (Table 4.4).
In order to accurately predict individual phase parameters for the steel, the fit involved using stressstrain data for each phase, as well as the global stress strain curve for the Q&P980. The commercial
optimization software LS-OPT (Stander et al., 2010) was used to calibrate the material parameters
using CPFE model against the tension tests. The LS-OPT software was integrated with ABAQUS
software using a Python script written specifically for the calibrated geometry to extract the data
and find the parameters based on the number of grains per integration point. The iterative
calibration procedure to find the hardening parameters for each phase was carried out until the
simulation and the as-received curves obtained from the in-situ HEXRD matched with each other.
After roughly fitting the hardening parameters for each phase to its corresponding HEXRD stressstrain curve, another method was applied to further optimize those parameters. In this method, we
tried to find the most accurate parameters for the hardening equation by only fitting two phases
and finding the parameters for the remaining phases such that it provides the best fit to the
macroscopic stress-stain curve of the Q&P980 steel. In this case, the fitted parameters for austenite
and the new martensite were kept similar to the fitted values in the first attempt since together they
constituted only 10% of the steel.
82

The final hardening parameters which were fitted to each phase using the above calibration method
are listed in Table 4.5. Figure 4.7 shows the comparison between the stress-strain curve for
individual phases and the best fitted curves with 60 grain orientations assumed per individual phase.
Using the microstructure information, the overall flow curve of Q&P980 steel sheet was modeled
with CPFE model based on 45% ferrite, 45% pre-existing martensite, 5% retained austenite, and
5% transformed martensite. Figure 4.7 also shows the comparison between predicted and
experimental macroscopic stress-strain curve for Q&P980 steel.
The effect of changing the number of grain orientations assumed per integration point on computed
stress-strain curves for the case of uniaxial tensile test is shown in Figure 4.8. It can be seen that
changing the number of grain orientations assumed per integration point does not have a significant
impact. Therefore, the same calibrated parameters obtained with 60 grain orientations were also
used in all CPFEM simulations with different number of grain orientations. It was also shown by
others that increasing or reducing the number of grain orientations will have little impact on the
stress-strain curve (Zamiri et al., 2007; Guan et al., 2006).

83

Table 4.3: Slip systems for BCC crystals; mÎą and nÎą denote the slip direction and the slip plane
normal of the slip system Îą, respectively in the initial configuration.



m



n

1

-1

1

1

1

1

0

2

1

-1

1

1

1

0

3

1

1

1

-1

1

0

4

1

1

-1

-1

1

0

5

1

-1

-1

1

0

1

6

1

1

-1

1

0

1

7

1

1

1

-1

0

1

8

-1

1

-1

-1

0

1

9

-1

-1

1

0

1

1

10

1

-1

1

0

1

1

11

1

1

1

0

-1

1

12

-1

1

1

0

-1

1

84

m

n

13

-1

-1

1

1

1

2

14

-1

1

1

1

-1

2

15

1

1

1

-1

-1

2

16

1

-1

1

-1

1

2

17

1

1

1

-1

2

-1

18

-1

1

1

1

2

-1

19

-1

1

-1

1

2

1

20

1

1

-1

-1

2

1

21

-1

1

1

2

1

1

22

1

-1

1

2

1

-1

23

1

1

1

2

-1

-1

24

1

1

-1

2

-1

1

Table 4.4: Slip systems for FCC crystals; mÎą and nÎą denote the slip direction and the slip plane
normal of the slip system Îą, respectively in the initial configuration.



m

n

1

1

-1

0

1

1

1

2

1

0

-1

1

1

1

3

0

1

-1

1

1

1

4

1

0

1

-1

1

1

5

1

1

0

-1

1

1

6

0

1

-1

-1

1

1

7

1

0

-1

1

-1

1

8

0

1

1

1

-1

1

9

1

1

0

1

-1

1

10

1

-1

0

1

1

-1

11

1

0

1

1

1

-1

12

0

1

1

1

1

-1

Table 4.5: Material parameters for the CPFE model composed of four phases in Q&P980 steel
(unit: MPa).
Phase

𝝉𝜶𝟎

𝝉𝜶𝟏

𝜽𝜶𝟎

𝜽𝜶𝟏

Ferrite

348

130

1100

117

Martensite

487

55

1994

185

Austenite

600

2

1600

550

New Martensite

496

39

2200

1788

85

Figure 4.7: Comparison of the experimental strain-stress curve with the CPFE model with 60 grain
orientations assumed for each constituent phase and the Q&P980.

Figure 4.8: Comparison of the experimental strain-stress curve with the CPFE model in terms of
the number of grain orientations assumed per integration point.
86

CHAPTER 5
EXPERIMENTAL AND CRYSTAL PLASTICITY MODELING OF THE
HEMESPHERICAL PUNCH TESTS
The experimental and numerical simulation of Nakajima hemispherical punch tests for Q&P980
steel are presented in this chapter. For CPFE model validation purposes, the calibrated parameters
from uniaxial tension test were adopted for the simulation of Nakajima tests with different sheet
widths. In addition, the comparison between the multiphase CPFE model and two
phenomenological models (isotropic von Mises and anisotropic Hill48 yield criteria) was
conducted for validation of the proposed model. Also, the numerical results were further evaluated
by comparing the predicted against experimental results.

5.1 Experimental Procedure of the Hemispherical Punch Test (Nakajima Test)
The Nakajima bulging test was conducted using Interlaken Servo-Press 150 formability test system
with the roof DIC system Figure 5.1. As shown in Figure 5.2(a), different blanks with varying
width (w) ranging from 20 mm up to 180 mm, and a 1 mm gauge thickness were stretched by the
hydraulic displacement of the Nakajima punch at a constant tool velocity of 1 mm/s. A constant
blank-holder force was applied to the sheet metal via an annular cylindrical die securely tightened
under the action of a 600 kN force. A schematic of the specimens’ dimensions and test setup are
shown in Figure 5.2(a) and (b), respectively. A double layer of Teflon sheet (0.1 mm sheet thickness)
were placed between the blank and the hemispherical punch Figure 5.2(c) to prevent premature
fracture of the specimen and to limit the effects of friction, particularly at the pole (International
Standard ISO 12004-2, 2008).

87

A 3D DIC system similar to the one used in the uniaxial tensile test was employed to measure
strain histories on the blank’s surface. A load cell built into the machine was used to record the
reaction force on the Nakajima punch. The measured strains on the surface of the sheet metal were
synchronized with the measured punch stroke and punch force for comparisons with numerical
predictions. Also, to study the formability of the Q&P980 based on FLC, the evolution of major
and minor strains in a 2 mm diameter area containing the failure location was tracked. The tracking
of major and minor strains was performed for every Nakajima specimen to represent different
deformation modes. Figure 5.3 depicts dimensions of the Nakajima test specimens, and the
deformed samples after failure.

Figure 5.1: Bulge testing setups with DIC system.

88

Figure 5.2: (a) Dimensional characteristics of the specimens used in the Nakajima tests; (b)
Principles and dimensions of the Nakajima test with a hemispherical punch; (c) Nakajima punch.

Figure 5.3: Formed samples after failure in the Nakajima test.

89

5.2 Numerical Simulations of the Nakajima Hemispherical Punch Stretching Tests
Numerical simulations of the Nakajima hemispherical punch stretching tests were performed with
the calibrated CPFE model for the 1 mm thickness Q&P980 steel to evaluate the effect of strain
conditions on the forming processes. Figure 5.4 shows the finite element model of the die, punch,
blank and blank holder used in the Nakajima test. In contrast with many FEM studies, which only
assumed one quarter of the full blank with symmetric boundary conditions, a full blank was used
in this study to predict the strain localization location. However, in order to avoid unnecessary
calculations and enhance the computational speed, a fully constrained boundary condition along
the draw bead edge of the sheet was assumed (elements outside the draw bead were removed from
the model). This type of boundary condition constrains all combinations of displacements and
rotations of a node.
Analytical rigid body surfaces were utilized for the die, blank holder and punch. The Coulomb
friction model was used to describe the contact interaction between the rigid body tooling and the
deformable blanks. Because the Teflon sheets were applied to the steel surfaces of the tools during
the experiment, a friction coefficient of Îź=0.04 was assumed for the frictional parts in the
simulation (Oberg et al., 2004). Furthermore, because the strain-rate independent material model
was applied in finite element simulations, and in order to enhance the calculation speed, massscaling technique was used by setting the minimum step increment time 10-7 seconds. Also, a
smooth step amplitude was assumed for the punch to minimize the dynamic effects for bulging
process in the finite element simulation (ABAQUS User’s Manual, 2000).

90

Figure 5.4: Schematic diagram for finite element model of a Nakajima test.

91

5.3 Assumptions for Implementation of Multiphase CPFE Model
The polycrystal CPFE model was implemented into the commercial Abaqus/Explicit finite element
code with aid of the user-defined material subroutine VUMAT. The model was then calibrated
based on the weight fraction of each phase in the 1 mm thick Q&P980 sheets used in the Nakajima
test with different widths.
The experimental microstructural observations of the Q&P980 steel sheet under consideration
revealed that this material consists of 44±2% by volume of the martensite phase, ~8–10% by
volume of the retained austenite phase, and the rest is ferrite. Using this information, the overall
microstructure of Q&P980 steel sheet was modeled with CPFE model based on 45% ferrite, 45%
pre-existing martensite, 5% retained austenite, and 5% transformed martensite. In the absence of
an austenite to martensite phase transformation model, it was assumed that a 5% new (transformed)
martensite exists in the undeformed sheet from the beginning of deformation. However, in the real
material, the volume faction of the new martensite evolves with increasing plastic strain. The
assumption of having 5% new martensite available from the beginning of the deformation will not
affect the material behavior in a significant way, because the newly transformed martensite will
have low volume fraction compared to many TRIP steels. This minor effect was also reported by
other researchers (Srivastava et al., 2015).
For the complex microstructure of the Q&P980 steel sheet and based on the EBSD data
in Figure 4.2, a random crystallographic orientation was assumed for grains in an aggregate. The
number of grain orientaions assumed per aggregate were 140 or less in order to avoid large
computational time, while still representing the actual behaviour of the steel sheet. Furthermore,
using very large numbers of grain orientations as input for CPFE simulations is often neither
92

practical nor scientifically rewarding, since it provides more texture information than is usually
necessary for predicting plastic anisotropy of formed parts and the evolution of texture (Raabe and
Roters, 2004).
Mapping the grains and texture distributions observed in the sheet metal to the finite element
analysis (FEA) integration points requires a method which is physically meaningful. At the same
time, the texture information must be reduced to an optimal level in which the CPFEM accurately
models the complex forming condition at reasonable computational costs. Therefore, in this study,
two methods were adopted to represent the initial texture of the Q&P980 steel. In the first method,
we assumed that all FE integration points have fixed numbers of randomly oriented grains
proportional to their phase volume fraction in the material. For example, when using 140 grain
orientations to represent the Q&P980 sheet at an integration point, those 140 were distributed as
63 ferrite (45%), 63 old martensite (45%), 7 austenite (5%) and 7 new martensite (5%) orientations.
In the second method, a random phase distribution for the grains was assumed at each integration
point, but the total distribution for each phase in the whole sheet was kept at 45% ferrite, 45% old
martensite, 5% austenite, and 5% new martensite. In the second method, not all integrations points
have equal volume fractions of phases. Figure 5.5 shows examples of phase distributions in 180
mm and 20 mm Nakajima test blank sizes with 140 grains per integration point. In Figure 5.5(a),
for the 180 mm circular blank, 99,310 elements with one integration point per element were
assumed. Therefore, the entire model of the sheet had a total of 13,903,400 grain orientations to
represent the blank. Of those orientations, there were 6,256,530 ferrite and 6,256,530 old
martensite grain orientations, representing 45% of each phase. Also, 695,170 grain orientations
were assumed for austenite and the new martensite representing a volume fraction of 5% for each
93

of them. As the distribution curve in Figure 5.5(a) shows, there are approximately 6,723 elements
with 63 grain orientations (out of 140 total) representing the ferrite phase, and almost similar
number for the martensite phase. Those same phases, in extreme cases, are represented by as little
as 37 grain orientations in 2 elements, and by as many as 88 grain orientations in 2 other elements.
In contrast, the other distribution curve in Figure 5.5(a) shows that there is no representation of the
austenite or the new martensite phases in about 80 elements. Also, there are about 15,096 elements
with 8 grain orientations representing austenite and the new martensite phases. Using such
distribution schemes were necessary to approximately represent the material microstructure and
inhomogeneity in terms of different phase distributions. A similar scheme was adopted for all the
Nakajima test blank widths in this study. Figure 5.5(b) shows an example of different phase
distributions assumed for the 20 mm wide blank.

94

(a)

(b)
Figure 5.5: The random distribution of phases orientations in the sheet metal elements (a) 180 mm
(99310 elements), (b) 20 mm (36526 elements).

95

5.4 CPFE Model Verification and Discussions
5.4.1 Punch Force-Displacement Analysis
Using the calibrated material parameters, FE simulations were performed until the strain
localization (necking) occurred. In order to optimize the element size in the Nakajima
hemispherical punch stretching test, the effect of the element size on the punch force-displacement
was investigated. As can be seen from Figure 5.6, the simulation results did not change
significantly for element sizes smaller than or equal to 0.5 mm. Therefore, to reduce the
computational time, the blanks were meshed with C3D8R brick elements and the element size was
set as 0.5 mm in all cases.
Figure 5.7(a) – (d) show the comparison between punch force-displacement curves obtained from
multiphase CPFE model, Hill48 and von-Mises models with experimental curves for four different
specimen widths. Using different specimen widths is valuable since as the width of the specimen
increases, the deformation changes from the uniaxial tension for the thinnest specimen to plane
strain for the wider specimen, and finally to the biaxial stretch for the circular specimen.
In the narrow blank with a 20 mm width in Figure 5.7(a), the predicted force-displacement curves
based on the crystal plasticity model with different number of grain orientations and phase
distributions perfectly match the experimental curve, including the predicting of the localized
failure point at the punch depth of ~25 mm. In the case of 20 grains per integration point, however,
the punch force-displacement curve drops a little earlier compared with 60 and 140 grains, which
perfectly match each other. The von-Mises and Hill48 yield criteria show a little higher force than
the CPFE model, and an earlier localization of around ~24 mm in the case of Hill48, and at ~22
mm in the case of von Mises. The CPFE simulation results are in better agreement with the
96

experimental data partly because the deformation mode is very close to the uniaxial tension, which
was used along the rolling direction as the reference state to calibrate the crystal plasticity material
parameters.
In the case of intermediate blank width of 60 mm, Figure 5.7(b), the deformation mode is in
between the uniaxial tension and plane strain tension. As can be seen, the simulated punch forcedisplacement curves for all CPFE models coincide with each other and show an excellent fit to the
experimental results up to about 23 mm punch displacement, but they overestimate the localized
necking associated with a drop in the curve. The curves for the von Mises and Hill48 models are
somewhat higher than the experimental curve, but the localized necking showed almost a good
match at 25 mm punch displacement with higher punch force.
The deformation mode for 100 mm blank width, Figure 5.7(c), is almost similar to the plane strain
case. The CPFE results in this case again show very good match with the experimental forcedisplacement curve up to a displacement of about 24 mm, while the curves for the von Mises and
Hill48 show somewhat higher forces, as well as an earlier necking. All the material models in this
case failed to correctly predict the necking point of the specimen in comparison with the
experimental result.
As for the 180 mm circular blank, Figure 5.7(d), where the deformation mode is close to the
balanced biaxial tension, the von Mises and Hill48 models show almost similar results for the
punch force-displacement, but predict a delayed localization point. The CPFE model’s prediction
of the punch force-displacement curve with different number of grain orientations and phase
distributions match the experimental curve very well up to a punch displacement of about 26 mm.

97

Figure 5.6: Comparison of punch force-displacement relationship computed with various mesh
sizes for 20 mm blank case.

98

(a)

(b)
Figure 5.7: A comparison of the experimental and numerically simulated punch force versus punch
displacement relationship during the Nakajima budging test with different blanks sizes (a) 20 mm,
(b) 60 mm, (c) 100 mm, (d) 180 mm.

99

Figure 5.7 (cont’d)

(c)

(d)
100

5.4.2 Strains Distribution Analysis
Further verification of the validity of the proposed multiphase CPFE model was performed through
strain analysis, in which the strain distribution from the CPFE was compared with experimental
strains measured by a DIC system in an area that included the fracture location. The red dots shown
in DIC images in Figure 5.8-Figure 5.13 are the fracture location.
Figure 5.8-Figure 5.13 show the comparison of predicted and measured strain distributions and
strain evolution along the rolling and transverse directions of the narrow blank (20 mm) at punch
depths of 20 mm, 22 mm, and 23.5 mm. It can be seen that compared with DIC results, the CPFE
model with different grain orientations and phase distributions can predict the evolution of strains
and their distributions, with increasing applied punch displacement (Figure 5.8). In contrast, the
von-Mises and Hill48 results show a different strain distribution.
In Figure 5.9 at a punch depth of 22 mm, for the 20 mm wide blank, it can be seen that the CPFE
strain distribution prediction is remarkably close to those measured from DIC. The DIC and the
CPFE images show almost similar symmetry in the strain pattern, homogeneity, and smoothness
especially with 140 grains per integration point. Interestingly, the multiphase CPFE model was
also able to show remarkable results with respect to the localized necking strain values, location
and direction, as seen in Figure 5.10, at a punch displacement of ~23.5 mm (1 mm away from
fracture).
Also, for the 60 mm blank width, the CPFE model was able to accurately predict the strain values
and distributions at a punch displacement of 24 mm (Figure 5.11). From this figure, the CPFE
model with random phase distribution was more accurate than that with fixed phase distribution
in predicting the strain distributions and the necking location. The CPFE model also shows good
101

prediction in terms of strain values and distributions for the 100 mm blank width at a punch
displacement of ~23 mm (Figure 5.12). Also, the CPFE simulation based on the random phase
distribution successfully captures the significant strain localization which was observed in DIC,
Figure 5.12(e).
It is interesting to note that although the CPFE model predicts a higher force-displacement at the
necking point, Figure 5.7(d), compared with experimental and phenomenological models, it
predicts a more realistic strain distribution compared with Hill48 model which had a closer fit in
terms of punch force displacement at the necking point Figure 5.13. Moreover, the CPFE model
was remarkably able to predict the localized necking location and direction at the correct punch
displacement of 30 mm, even though Figure 5.7(d) shows that the experimental force-displacement
curve for the 180 mm blank is lower than the CPFE model’s prediction.
Overall, the strain analysis based on the multiphase CPFE model shows a very good agreement
with the experimental results (obtained from DIC) compared with phenomenological models. The
CPFE model was able to accurately predict strain values and localized necking locations and
directions compared with experimental DIC results.
As mentioned in chapter 4, both the CPFE model and the phenomenological models were
calibrated based only on the uniaxial tension data. Although the strain range in the uniaxial tension
data used for the calibration were relatively small compared to the strain range in the Nakajima
bulging test, the simulated strain distributions based on the CPFE model agreed very well with
experiments in the Nakajima bulging test as shown in Figure 5.8-Figure 5.13. Also, in all cases,
the simulated punch force-displacement based on the CPFE model matched well with experiments,
while the phenomenological models over predicted the punch force in the 60 mm and 100 mm
102

wide blanks as shown in Figure 5.7(b) and (c), respectively. In general, quadratic yield functions
have higher yield stress in the plane strain tension mode compared to the non-quadratic yield
functions. Therefore, quadratic yield functions, such as von Mises and Hill48 functions, which
were considered for the phenomenological models in this study, lead to the over-prediction of the
punch force for the 60 mm and 100 mm wide specimens. As shown in Figure 5.16 of the next
section, the deformation modes for the 60 mm and 100 mm wide specimens during the Nakajima
bulging tests are near the plane strain tension mode, in which quadratic yield functions over predict
the stress level.
It can be noted that strain localizations observed from the DIC measurements (Figure 5.8-Figure
5.13) occurred at earlier punch strokes compared with the onset of necking from the punch force
versus punch stroke curves in Figure 5.7. As shown in Figure 5.8-Figure 5.13, the onset of the
strain localization observed from DIC measurements were at punch strokes of 23.5 mm, 24 mm,
23 mm, and 30 mm, for specimens with 20 mm, 60 mm, 100 mm, and 180 mm widths, respectively.
The evolution of the strain localization in the FE simulation is affected by the applied material
model. However, the CPFE model in this work was not intended to account for the material
behavior beyond the strain localization. In other words, strain rate sensitivity and damage behavior
were not considered in this study. Also, the extended Voce type hardening rule shown in Eq. (3.11)
converges to the linear hardening rule as the exponent part vanishes with the increasing
accumulated shear strain,  . Therefore, the material models could not accurately predict the
necking of the punch force versus punch stroke curves.

103

(a)

(e)

(b)

(c)

(f)

(d)

(g)

Figure 5.8: Comparison between the numerically simulated and experimental results for rolling
direction strain (Epsilon Y) of a 20 mm wide sheet at a punch depth of 20 mm. (a) CPFE with 140
grains with fixed distribution (b) CPFE with 140 grains with random distribution (c) CPFE with
60 grains with fixed distribution (d) CPFE with 60 grains with random distribution (e) von Mises
(f) Hill48 (g) DIC. (Note: in figure 5.8(g), the averaged minor and major strains histories at the
fracture area are calculated by DIC over a 2 mm diameter surface patch which is identified by the
red dot in each DIC image; colored lines are used to calculate the strain histories along 60 mm
section length at multiple locations).

104

(a)

(b)

(e)

(f)

(c)

(d)

(g)

Figure 5.9: Comparison between the numerically simulated and experimental results for the rolling
(Epsilon Y) and transvers (Epsilon X) strains of a 20 mm wide sheet at a punch depth of 22 mm (a)
CPFE with 140 grains with fixed distribution (b) CPFEM with 140 grains with random distribution
(c) CPFE with 60 grains with fixed distribution (d) CPFE with 60 grains with random distribution
(e) von Mises (f) Hill48 (g) DIC.
105

(a)

(b)

(e)

(f)

(c)

(d)

(g)

Figure 5.10: Comparison between the numerically simulated and experimental results for the
rolling (Epsilon Y) and transvers (Epsilon X) strains of a 20 mm wide sheet at a punch depth of 23.5
mm. (a) CPFEM with 140 grains with fixed distribution (b) CPFE with 140 grains with random
distribution (c) CPFE with 60 grains with fixed distribution (d) CPFE with 60 grains with random
distribution (e) von Mises (f) Hill48 (g) DIC.
106

(a)

(b)

(d)

(c)

(e)

Figure 5.11: Comparison between the simulated and experimental results with blank width of
60 mm for rolling direction strain (Epsilon Y) at a punch depth of 24 mm (~1 mm from experimental
fracture point). (a) CPFE with 140 grains and fixed distribution (b) CPFE with 140 grains and
random distribution (c) von Mises (d) Hill48 (e) DIC.

107

(a)

(b)

(d)

(c)

(e)

Figure 5.12: Comparison between the simulated and experimental results of the rolling (Epsilon Y)
and transvers (Epsilon X) strains of a 100 mm wide sheet at a punch depth of ~23 mm. (a) CPFE
with 140 grains with fixed distribution (b) CPFE with 140 grains with random distribution (c) von
Mises (d) Hill48 (e) DIC.

108

(a)

(b)

(d)

(c)

(e)

Figure 5.13: Comparison between the simulated and experimental results for the rolling (Epsilon
Y) and transvers (Epsilon X) strains of a 180 mm wide sheet at a punch depth of ~30 mm (a) CPFE
with 140 grains with fixed distribution (b) CPFE with 140 grains with random distribution (c)
von Mises (d) Hill48 (e) DIC.
109

5.4.3 Identification of Necking Strains for FLC
The Nakajima test provides the data necessary to construct a FLC for a sheet metal. As explained
in Section 5.1, during the Nakajima test the strain evolution in a small surface patch is traced by
DIC for each deformation path (uniaxial tension, plane strain tension, balanced biaxial, etc.) from
the beginning until the failure occurs. The averaged minor and major strains are calculated by DIC
over a 2 mm diameter surface patch (total area ~3.14 mm2) which includes the fracture point. The
red dot in each DIC image in Figure 5.8-Figure 5.13 show this area.
The FLC for Q&P980 multiphase steel was constructed using DIC, as well as with CPFE
simulation results obtained from a similar size area and approximately in the same position of
fracture in the DIC. After plotting the history of major and minor strains, as well as the thinning
rate, as a function of the punch stroke (Figure 5.14), the time dependent linear best fitting (LBF)
method was employed to determine the major and minor strains at the onset of instability.
The LBF technique was first proposed by Volk and Hora (2011) and later modified by Hotz et al.
(2013), for calculating the thinning rate of the material in the vicinity of the necking point. When
the thinning rate at a local area on the surface of the sheet metal grows faster than the thinning rate
in the surrounding area, it signals the onset of necking. In this study, the onset of necking was
evaluated with LBF method using two best fit straight lines; one is calculated as a straight fitting
line through the area of stable deformation (see Stable Line Fit, Figure 5.14), and the second is a
gliding straight line (see Unstable Line Fit, Figure 5.14) fitted through the thinning rate and
calculated for each point in time based on the change in the punch position (Hotz et al., 2013). The
intersection of these two lines identifies the punch stroke at which the unstable necking starts and
the corresponding major and minor strains are taken to construct the FLC.
110

Figure 5.7(a)–5.8(d) show that compared with experimental curves, the necking strains predicted
by von Mises and Hill48 models for various modes of deformation are inconsistent. That is, for 20
mm and 100 mm wide specimens the necking point are predicted to occur sooner, while for the
180 mm wide specimen it is predicted to occur later. Figure 5.8-Figure 5.13 also show that in
general the surface strain distributions predicted by these two phenomenological material models
do not accurately match DIC results. Based on these observations, it was decided to not apply the
LBF method to find strains at the onset of necking for these two phenomenological material models.
In order to determine the onset of necking with multiphase CPFE model, 12 elements (3x4
arrangement) measuring 0.5 mm per side (3 mm2 total area) were identified on the surface of the
specimen where strain localization occurred (Figure 5.15). The averaged logarithmic major (𝜀1)
and minor (𝜀2) strains from these 12 elements were plotted as a function of punch stroke. The
thickness direction strain (thinning strain) at the necking area was calculated based on 𝜀3= -(𝜀1+𝜀2).
Since our calculations are based on rate independent (time in simulation does not match the time
in experiments), all strains and thinning rate in the simulation were determined based on the
corresponding punch displacement instead of time. The value of the thinning rate was calculated
at each point using:

d  3  3i 1   3i
 i 1
dZ
Z  Zi

(5.1)

Here, the subscripts, i and i+1, denote the discrete time step number, and Z is the punch position.

111

In Figure 5.14, the unstable line fit goes through one of the three points on the thinning rate curve,
just before the specimen fails, with the largest slope. Note that the same method was applied to
both experimental and numerical simulation results to determine the onset of strain localization in
this study.
The experimental and numerical FLC calculated using the LBF method are plotted in Figure 5.16.
Forming limit strains numerically obtained by random and fixed phases distribution matched very
well with experiments for the 20 mm wide blank, and reasonably well matched for the 180 mm
wide blank, as shown in Figure 5.16(a) and (b). However, the forming limit strains obtained
numerically for other specimen widths were higher compared with experimental FLC. Compared
with the random phase distribution, the forming limit strains predicted with the fixed phase
distribution were slightly higher.
Overall, the inhomogeneity due to phase distributions in the Q&P980 sheet plays an important role
in controlling the localized strain values for the FLC. As mentioned in Section 5.4.1, for the
20 mm wide specimen the uniform deformation is dominated from the beginning, and the
deformation mode stays very close to the uniaxial tension. Given that the uniaxial tension test data
along the rolling direction was used to calibrate the material parameters for each phase, the
inhomogeneity did not have a considerable impact on the predicted necking strain value in FLC
for the 20 mm wide sheet as a function of assumed phase distributions.
For the balanced biaxial case, the formability of the sheet is less affected by inherent defects such
as metal inhomogeneity. Hence, the difference between necking strains based on random and fixed
phase distributions for the 180 mm wide specimen was not as large as in other deformation modes.
Also, it was observed that the simulated strain path for the 180 mm wide sheet with random phase
112

distributions is very close and in good agreement with the experimental strain path (Figure 5.17).
This implies that major and minor strains obtained by the random phase distribution accurately
represents the actual state of deformation, which is almost monotonously proportional to the
experimental strain path. In contrast, it can be seen from Figure 5.17 that the strain path from the
simulation with the fixed phase distribution is slightly different than the experimental deformation
path.
It can be clearly seen that FLC localized necking points are highly affected by inhomogeneity in
blank sizes ranging between 20 mm and 180 mm, as shown in Figure 5.16(a) and (b). The lower
necking strain values predicted by the random phase distribution are attributed to early localization
of strains in the elements with larger volume fraction of softer ferrite or austenite phases, compared
with the stronger martensite.
In summary, the inhomogeneity of the steel sheet at the micro level has significant impact on the
macro level performance, and the correct prediction of stress and strain due to the inhomogeneity
can lead to better estimation of FLC strains. The inhomogeneity effect due to the assumption of
random phase distribution showed better results in comparison with fixed phase distribution.
However, the assumption of random phase distribution was not sufficient to bring the FLC points
in deformation modes between the uniaxial and balanced biaxial deformation closer to the
experimental FLC. By considering more of the material inhomogeneity related to the phase
distribution, it will be possible to bring those strains much closer to the experimental FLC, without
adversely affecting the uniaxial and balanced biaxial necking strains.

113

Figure 5.14: Determination of instable necking strains using the linear best fitting (LBF) method.

(a)

(b)

Figure 5.15: Identification of the necking area, 12 elements on the sheet’ surface of (a) 20 mm, (b)
180 mm blank width.

114

(a)

(b)
Figure 5.16: Strain-based forming limit curve measured using the hemispherical dome stretching
Nakajima test (a) random phases distribution (b) fixed phases distribution.
115

Figure 5.17: Comparison between experimental and simulated strain paths in balanced biaxial
case.

116

CHAPTER 6
FORMING LIMIT CURVE PREDICTION BASED ON MK-CPFE MODEL
In this chapter, a multiscale computational framework based on MK theory (Marciniak and
Kuczyński, 1967) and crystal plasticity finite element CPFE model (i.e., MK-CPFE) is developed
to predict FLC for multiphase steel sheets.
Generally, the most reliable method to numerically generate FLC for sheet metals is to simulate
the Nakajima hemispherical punch test with similar boundary conditions and specimen sizes used
in the experiments. However, given the complex geometry of the Nakajima hemispherical punch
test and specimen shapes, and highly intensive computational time, especially with CPFEM, it was
opted to use the simpler rectangular geometry of the specimen used in the MK analysis for FLC
prediction.
It was found that the analytical prediction of FLC using MK model in conjunction with crystal
plasticity and Taylor homogenization scheme does not accurately predict the forming limits. The
main shortcoming of this analytical approach was identified to be its assumption that all crystal
orientations exist and evolve at only one material point, and therefore crystal orientation
distribution of the specimen cannot be accounted for. In contrast, by using FEA to predict forming
limits, multiple integration points with different orientation distributions can be assumed to more
accurately represent microstructural properties of the sheet metal. Also, by using FEA to predict
FLC using MK model, it is possible to impose the more physical boundary conditions, such as
displacement or velocity, rather than strains, which are imposed in analytical models.

117

6.1 FLC Determination Procedures Based on Finite Element Simulation
Similar to MK theory, the right-hand side of FLC is determined with CPFE using a rectangular
specimen with a groove in the mid-section, as shown in Figure 6.1(a). To determine the appropriate
geometry of the rectangular specimen, however, a number of finite element simulations were first
performed with von Mises yield function (i.e., MK-FEM) to study the sensitivity of predicted FLC
to geometrical parameters such as bulk and groove thicknesses ( t0b and t 0g ), the transition region
between bulk and groove sections, and other geometrical aspect ratios ( ARb and ARg ). Using
the optimum specimen geometry and the calculated imperfection factor f 0 (where t0g  f 0 t0b ),
finite element simulations with CPFE model were performed to predict the right-hand side of FLC
for the multiphase steel sheet.

6.2 The Finite Element Procedure and the Failure Criteria
In the MK model, an imperfection, or a groove, is assumed to exist initially in the sheet metal as
illustrated in Figure 6.1(a). In real materials, however, numerous imperfections are present at
arbitrary angles, and therefore, imperfections aligned at different angles should be considered for
obtaining more accurate results. However, due to the inhomogeneous phase distributions in the
Q&P980 sheet, the localized strain and fracture were mostly observed to take place along the
rolling direction in the experiments. Therefore, for the sake of simplicity and computational cost
reduction, only an imperfection parallel to the rolling direction (RD), which corresponds to the x1
axis in Figure 6.1(a), was considered for the calculation of the right-side of FLC based on the MKFEM model.
In this study, a sheet specimen with the bulk thickness of t0b , width (dimension along x1) of 2W,
and a total length (dimension along x2) of 2L was considered, as shown in Figure 6.1(a). The
Cartesian tensor notation was used, and the origin of the coordinate system was taken to be at the
118

center of the specimen (i.e. sheet specimen was oriented so that the x1, x2, and x3-axes were
parallel to the RD, TD, and ND of the sheet material, respectively). In order to reduce the
computational cost, symmetric boundary conditions along x1 and x2 axes were assumed for
regions 0  x1  W and 0  x2  L , so that only one quarter of the sheet specimen was needed for
computational modeling, as shown in Figure 6.1(a).
The commercial software MATLAB was used to generate different dimensions for groove and
bulk regions. The dimensional values were then used in the mesh generation, creation of a finite
element model, and subsequent simulations with the commercial software ABAQUS. Solid
elements were used for FEA simulations as illustrated in Figure 6.1(a), with fine symmetric mesh
(0.05 mm × 0.05 mm elements) in the groove region. A gradual increase in the mesh size in the
bulk region, along x2 direction from 0.1 mm to 0.4 mm, was assumed in order to increase
computational efficiency. To have a smooth transition of strains between the bulk and groove, a
transition region with one element (0.1 mm thickness) was adopted, as shown in Figure 6.1(b).
The generated strains in the transition region were neglected in the analysis since they did not
represent either the bulk or the groove region.
In all the MK-FEM simulations, the nominal thickness of the specimen in the bulk region was
assumed to be t0b =1.0 mm, which corresponds to the measured average thickness of the Q&P980
sheets used in all the experiments. Also, the specimen width (W) was fixed to 3 mm length. The
other two in-plane dimensions, Lg and Lb were iteratively adjusted until the best values for the
bulk and groove aspect ratios ( ARb and ARg ) were found, such that the MK-FEM predicted FLC
matched the FLC calculated by analytical MK model. The selected parameters were chosen such
that a direct comparison with the analytical MK model could be made. Therefore, the following
three-groove aspect ratios ( ARg ) were considered, 4, 6, and 7.5, in the MK-FEM simulations.
119

In addition, three aspect ratios for the bulk region were assumed for calibration. To determine the
thickness of the groove region, it was gradually reduced until strain localization occurred inside of
the imperfection area. The initial thickness of the groove region is given by:
t0g  f 0 t0b

(6.1)

where f 0 is the inhomogeneity factor, t 0g and t0b are initial thicknesses of the groove and bulk,
respectively.
The sheet specimen was deformed under in-plane biaxial straining by prescribing the
displacements d1 and d2. Therefore, the state of biaxial stretching was imposed as following:



1
2

 0    1

(6.2)

where  is the strain ratio corresponding to the deformation path, 1 and  2 are strains caused
by imposing displacements d1 and d2 on planes A1 and A2, respectively, as shown in Figure 6.1(b).
It is important to mention that by imposing uniform displacements d1 and d2 in the finite element
analysis, only proportional engineering strains will be developed, and not the proportional true
strains needed to construct FLC. To get around this problem, non-uniform displacements d1 and d2
had to be generated analytically first, and then imposed in the MK-FEA simulations to generate
proportional true strains throughout the deformation. Figure 6.1(c) visually shows the difference
in proportionality between engineering and true quantities generated by assuming that   0.5 in
engineering quantities (blue line). From the non-proportional engineering strain history, true strain
history was converted and illustrated in red line, which can clearly show that the proportionality
in true strains were maintained.
120

During initial stages of the biaxial stretching of the sheet simulations, the sheet deformed
uniformly in both the bulk and groove regions. However, as plastic deformations increased, a strain
gradient developed and strains inside the groove region grew much faster than in the bulk region,
resulting in a localized neck. To predict the limit strains from the groove region corresponding to
the formation of a neck, the following failure criterion was used. A FLC was then constructed by
repeating the simulations for several different strain paths  (i.e., 0    1 ), and plotting the
principal strains in the bulk region once the assumed failure criterion was reached. In this study,
the failure criterion was based on the onset of necking, identified as:
g
  22
b
  22

 10

(6.3)

where,



g
22

1

Ng

 22b 

1
Nb

Ng

  
n 1

Nb

g n
22

  
n 1

b n
22

(6.4)

(6.5)

The sum of  22 in equations (6.4) and (6.5) is taken over all elements within groove and bulk
regions, respectively. Ng and Nb are the total number of elements in groove and bulk regions,
respectively.
The value of  22b

at the onset of necking in Eq. (6.5) will be referred to as the ‘major strain’

and the corresponding value of 11b

(analogous to  22b

‘minor strain’ for construction of FLCs.
121

in Eq. (6.5)) will be referred to as the

To predict FLCs of the steel sheet based on phenomenological von-Mises model, finite element
simulations of in-plane biaxial straining of the steel sheet specimen characterized by varying aspect
ratios of the groove and the bulk were carried out. These simulations used the calibrated hardening
parameters based on the Swift type model shown in Eq. (4.2).

122

Figure 6.1: (a) Schematic of the unreformed finite element mesh, representing full and one quarter
of sheet blank. The sheet blank is subjected to biaxial straining by applying proportional
displacements, d1 and d2 and symmetric boundary conditions; (b) Schematic of a sheet blank
characteristics and dimensions; (c) A comparison between the engineering and true strains histories
on the prediction of FLC for ρ=0.5.
123

Figure 6.1 (cont’d)

(c)

124

6.3 Results and Discussion
6.3.1 Aspects Ratios Calibration with Analytical MK Model
Forming limits of the steel sheet under biaxial stretching were calculated using the framework
described in Section 6.2. The calculated forming limits based on von Mises and MK-FEM were
compared with the calculated data based on the analytical MK model. The predicted FLCs based
on changing the aspect ratio of the groove, ARg , and imperfection factor, f 0 , for a fixed bulk
aspect ratio of ARb  3 5 , are shown in Figure 6.2(a)-(c). From the comparison of the relative
positions of the different curves in Figure 6.2, it is apparent that the size of the groove as
characterized by the aspect ratio of the groove ARg , for a given f 0 , influences the level of the
FLC. The conclusion is that a defect (assumed groove thickness) can cause a strain localization
that is controlled by both f 0 and its aspect ratio ARg ( W Lg ) .
Compared with the analytical MK model, the MK-FEM model overestimates the magnitude of
FLC0 (plane strain) as the assumed groove aspect ratio becomes larger than 4. For example, as
shown in Figure 6.2(c), for f 0 =0.996 and ARg =7.5, the FLC0 computed by FEM is on average
13% larger than those computed by the analytical MK model. This is an expected result as the
assumed analytical MK defect is more severe for a given sheet thickness than the defect assumed
in the MK-FEM model, and hence the plastic deformation localizes faster. Therefore, the ARg
for the MK-FEM model must be adjusted to match the specified FLC0 for a given material.
The FLC calculated by the MK-FEM model with the bulk aspect ratio of ARb 1.0 , predicts
forming strains that are slightly larger than those predicted by the analytical MK model for strain
ratios ρ less than 0.5, as shown in Figure 6.3. Less sensitivity was noticed when the values of ARb
reduced to 3/5 and 3/7. The effect of the aspect ratios on FLC predictions can be summarized as
follows; increasing the value of ARg increases the values of the limit strains, whereas the values
125

of the limit strains decrease when the values of ARb increases. However, the values of the limit
strains are less sensitive to the changes in ARb compared with the changes in ARg values.
The shape of analytical MK and MK-FEM FLCs are nearly identical for strain paths in the range
of 0 ≤ ρ ≤ 0.5, and FLC0, except for strain paths in the range of 0.6 ≤ ρ ≤ 1, as shown in Figure 6.2.
The comparison between various strain paths in groove and bulk regions is also shown in
Figure 6.4 for five different cases. The respective strain paths for the groove and bulk calculated
from the analytical MK model generally match the strains paths predicted by MK-FEM. However,
small deviation can be observed in the cases of ρ > 0.6.
As mentioned, the deformation behavior in the groove and the adjacent bulk with MK-FEM and
analytical MK model is different because of the nonuniformity introduced by the presence of finite
dimensions of groove and bulk with MK-FEM. Therefore, all limit strains from MK-FEM are not
expected to have similar values to those predicted by analytical MK model due to the uniqueness
of each model. Hence, the calibration method with analytical MK model was proposed here to find
the best parameters that give approximate limit strains to that calculated by MK-FEM, and
subsequently avoid choosing arbitrary dimensions that may affect the FLC prediction by adopting
different constitutive models.
As expected, the FLCs calculated by MK models (analytical MK and MK-FEM), using von Mises
and Swift type model for Q&P980 showed low plane strain (FLC0) value, and high balanced
biaxial strain limits, which is inconsistent with what was observed from experimental results
(Wang and Speer, 2013; Zhang et al., 2016). It should be mentioned that FLC predicted by MK
models (analytical MK and MK-FEM) is the combined result of the assumed failure model and
the consitututive model that describes material properties and behaviour under plastic deformation.
126

In the next section, it will be shown that by adopting a more sophisticated material model (crystal
plasticity), but the same failure model described in Section 6.2, more accurate results that better
fit experimental measurments can be obtained.

127

(a)

(b)
Figure 6.2: Comparison between the predicted FLCs based on analytical MK model and MK-FEM
for different groove aspect ratios with fixed bulk aspect ratio (ARb =3/5): (a) f0 =0.992; (b)
f0 =0.994; (c) f0 =0.996.
128

Figure 6.2 (cont’d)

(c)

129

Figure 6.3: Bulk’s aspect ratio (ARb) effects on the of predicted MK-FEM forming limit curve for
fixed groove aspect ratio (ARg=4) and imperfection factor (f0 =0.992).

Figure 6.4: Comparison of the predicted MK-FEM and analytical MK model strain paths for an
imperfection factor f0 =0.992.
130

6.3.2 FLC Prediction Based on MK-CPFE
It is well known that a slight change in the shape of the yield surface for a sheet metal can have a
rather large effect on predicted FLC (Lian et al., 1989). In the multiphase MK-CPFE forming limit
curve analysis carried out here, beside the mechanical properties of various phases, the shape of
the ‘yield surface’ for the sheet metal is specified by the assumed initial texture, as well as its
evolution.
The MK analysis predicts the FLC based on the growth of an initial imperfection factor f 0 .
However, the strength of the imperfection cannot be directly measured by physical experiments.
Thus, it is possible to "predict" the experimentally observed FLC0 by a suitable choice of the f 0 .
In other words, since f 0 has no physical meaning, it is simply a fitting parameter that allows a
model to be calibrated by matching the experimentally measured FLC0. Here, the lowest initial
imperfection factor ( f 0 0.992) and the initial aspects ratios of ARg = 4 and ARb = 3/5 were
chosen for the FLC prediction based on MK-CPFE. These values showed the best fit for the lowest
strain limit in plane strain condition (FLC0).
The initial texture (represented by 912,000 grains) with 20 grains per integration point was
employed in the FLC simulations, as shown in Figure 6.5(a). The orientations and volume fractions
of phases were distributed randomly in all integration points. Figure 6.5(a) shows two different
strategies to model FLC based on random distribution of phases. FLC1 was predicted by using the
same randomly generated volume fraction and orientation distributions in all integration points for
all strain paths (0 ≤ ρ ≤ 1). On the other hand, FLC2 was predicted by using different randomly
generated volume fraction and orientation distributions in all integration points for all strain paths
(0 ≤ ρ ≤ 1). In the second method, the material subroutine (multiphase CPFE-VUMAT) randomly
generated a new volume fraction of phases and orientation distributions for each simulation run.
131

In all the simulations, however, the volume fraction of phases were kept at 45% ferrite, 45% old
martensite, 5% austenite, and 5% new martensite in the bulk and the groove sections of the
specimen. It should be noted that for each strain path, ρ, the MK-CPFE simulation was repeated
five times, and the predicted limit strains using this method is plotted in Figure 6.5(a). The FLC2
curve was then obtained by applying the best fit to these data. Using such distribution scheme was
necessary to approximately represent the material microstructure and inhomogeneity in terms of
different phase distributions with well-known MK approach, and the results are in better agreement
with the scattered experimental limit strains of the Q&P980.
To study the sensitivity of the total number of grain orientations assumed per integration point on
FLC2 curve, additional simulations with 40 grains (represented by 1,824,000 grains) and 10 grains
(represented by 456,000 grains) were also performed and compared with predicted FLC2 using 20
grains. The fitted FLC2 curves for 40 grains and 10 grains are shown in Figure 6.5(b) and Figure
6.5(c), respectively. These figures also show that compared with the case of using 20 grain
orientations, scattering in the data for each strain path, ρ, significantly increases for the case of 10
grains, but it also significantly decreases when 40 grain orientations is used. Increasing the number
of grain orientations per integration point will reduce the overall material inhomogeneity and
creates less scatter. Figure 6.5(d) shows that FLC2 curves predicted by 20 and 40 grains are similar,
while that for 10 grains predicts a lower and more conservative FLC2 curve.
Although the incorporation of crystallographic texture and its evolution increases the modeling
effort and the computational time for each simulation, it is necessary to consider the phase
distribution features in order to improve the prediction of the plastic anisotropy. On the other hand,
a realistic distribution of phases has to be considered as well. This can be seen when the simulation
132

results in Figure 6.5, indicate that the distribution of constituent phases affect the anisotropy and
subsequently the forming strains of the multiphase AHSS. Such approaches can be used to
facilitate a more detailed material modeling.
Figure 6.6 shows the effect of initial imperfection factor f 0 on the predicted FLC with the same
procedure and material parameters adopted as before for f 0 =0.992 and 20 grains per integration
point. As expected, the larger the initial imperfection f 0 , the larger the critical strain for sheet
necking. As f 0 increases from 0.992 to 0.994, the FLC0 increases by ~5% with no significant
change in the FLC shape. However, FLC0 computed with f 0 =0.996 is on average 15% larger than
that predicted with f 0 =0.992. Also, the shape of FLC predicted by f 0 =0.996 is different than
those predicted by f 0 =0.992 and f 0 =0.994. Such change in the shape of FLC can also be observed
when using a non-quadratic yield function as f 0 approaches the value of one (Zhang et al., 2014).
The influence of texture evolution on FLC prediction has already been shown to be important on
the necking in single phase metals (Wu et al., 1998). In addition to initial texture effect, for
multiphase steels, complexity arises due to a difference in the flow characteristics of constituent
phases. For example, Figure 6.7(b) shows the evolution of 2000 randomly oriented grains per
phase at the onset of necking for multiphase Q&P980 under balanced biaxial tension in the bulk
and groove region with f 0 = 0.992. As expected, ferrite has more texture intensity and tendency
than other phases to accommodate the balanced biaxial deformation mode in the bulk. One may
attribute this to the fact that ferrite, beside the high initial volume fraction of 45%, is the softest
phase with respect to other three phases. Therefore, the strain heterogeneity in elements that
contain high volume fraction of ferrite is much higher than that in other elements in CPFEM
simulation. Austenite with 5% volume fraction shows a low intensity, and even a little lower
intensity can be observed for the harder new martensite when it is compared with the initial texture,
as shown in Figure 6.7(a).
133

The equivalent plastic strain distribution at the onset of necking in each individual phase inside
and outside of the groove is plotted in Figure 6.7(c). The ferrite phase shows a little higher plastic
strain compared with other phases in the groove and bulk regions. Based on experimental
observations, the flow strength of ferrite is lower than that of other phases; therefore, it is expected
that plastic deformation initiate in the softer ferrite phase, compared with other harder phases
which might even be deforming elastically. From crystal plasticity perspective, however, the
plastic deformation could also simultaneously occur in a harder phase, as long as the harder phase
is more aligned with the loading direction. For that reason, the plastic strain in the ferrite phase is
only slightly larger than that in martensite, which is a harder phase.
As pointed out before, in the FLC analysis carried out based on MK-CPFE model, the shape of the
yield surface for the steel sheet is determined based on the initial texture evolution, and the
mechanical properties of constituent phases. Consequently, the individual phase property has one
of the most important impacts on the FLC prediction. It can be noticed that ferrite may have an
effect on the predicted limit strain and subsequently on FLC more than other phases. This agrees
well with DIC measurements presented in chapters 4 and 5 when localization in many cases tend
to concentrate where ferrite phases are.
It is well known that for linear strain paths, ρ < 0, the strain localization occurs at an angle inclined
to the prescribed major strain rate direction (Hutchinson and Neale, 1978). Hence, the assumption
of normal groove to major strain rate direction will overestimate the strain limits for the left-hand
side of FLC. Therefore, the left-hand side of FLC was constructed using the limit strain obtained
from the time dependent linear best fitting (LBF) method from chapter 5, as shown in Figure 6.8.
LBF is employed to determine the major and minor strains at the onset of necking using the
134

simulated hemispherical punch test of a narrow specimen of 20 mm wide. That method showed
good results for determining the limit strains for uniaxial and balanced biaxial cases, however,
overestimated the experimental results for in-between deformation paths. That overestimation was
related to the inhomogeneity effect, which was solved in the current study by incorporating the
MK approach to construct the right-hand side of FLC.
The simulation methods carried out in the current study successfully predicted the FLC of the
Q&P980 steels. In particular, it is interesting to note that the predicted FLC is in agreement with
the experimental measurements, as shown in Figure 6.8.

135

(a)

(b)
Figure 6.5: Comparison of the predicted FLCs by MK-CPFE (with f0 = 0.992, ARg=4 and
ARb= 3/5) with experimental data: (a) 20 grains; (b) 40 grains; (c) 10 grains; (d) Comparison of
the fitted curves based on the number of grains.
136

Figure 6.5 (cont’d)

(c)

(d)
137

Figure 6.6: Comparison of the predicted FLCs by MK-CPFE with different imperfection factors.

138

Ferrite
Martensite
Austenite

New Martensite

(a)

Bulk

Groove

(b)
Figure 6.7: (a) Initial random orientation; (b) Simulated textures of 2000 orientations at onset of
necking represented by (111) stereographic pole figure of an imperfection factor f0 =0.992 and
balanced biaxial tension ( ρ=1); (c) Equivalent plastic strain contribution of the constituents phases
at the onset of necking.

139

Figure 6.7 (cont’d)

(c)

140

Figure 6.8: FLC of Q&P980 predicted by MK-CPFE for the right-hand and LBF-CPFE of the lefthand side.

141

CHAPTER 7
APPLICATION OF THE MULTIPHASE CPFE MODEL FOR SHEET METAL
STAMPING
The T-shape stamping was used to evaluate the formability of the Q&P980 steel sheet under
complex forming conditions. The maximum punch stroke without tearing or splitting the sheet was
taken as a relative formability index for different test setups.
The T-shape was used in the present study to validate proposed methods for predicting the FLC
based on multiphase CPFE model. This test was specically designed to closely replicate actual
forming conditions experienced by the sheet metal during the continuous stamping production.

7.1 T-Shape Stamping Experiment Setup
The experiments of T-shape stamping were carried out by using Interlaken Servo-Press system.
Figure 7.1(a) shows the tooling and the experimental setup. T-shape tests of 1.0 mm sheet thickness
of Q&P980 steel were performed for two different blanks dimensions as shown in Figure 7.1(b).
The formability of Q&P980 steel was investigated under a constant binder load of 20.0 Ton, and a
punch speed of 0.1 mm/s. The stamping process took place in two steps: first a constant binder
load was applied on the blank, and then the T-shape was formed under the actions of die and punch.

142

(a)

(b)
Figure 7.1: (a) The sheet metal stamping press tool components setup; (b) Dimensional
characteristics of the specimens used in the T-shape experiments.

143

7.2 Assumptions for T-Shape Implementation
Figure 7.2 shows the finite element model of the die, punch, blank and blank holder used in the Tshape simulation. Analytical rigid body surfaces were utilized for the die, blank holder and the
punch. Blanks were meshed with 0.333 mm×1 mm elements (three layers through the thickness)
and C3D8R element type from the ABAQUS/ Explicit element library. It is well known that
computational times increase almost quadratically by increasing element number (Tekkaya, 2000).
Also, reducing the size of the blank element will require additional mesh refinement for the Tshape part, and subsequently increase the computational time required by the microstructure-based
CPFE model. Hence, it is important to mention here that multiphase CPFE model was less sensitive
to changes made in the element size compared with phenomenological models under large
deformation applications. Therefore, additional mesh refinement below 1 mm (surface element
size) and 0.333 mm (thickness element size) was not considered in this study.
The Coulomb friction model was used to describe the contact interaction between the rigid body
tooling and the deformable blank. A friction coefficient of 0.12 was assumed in the CPFE
simulation, as it showed the best fit to the experimental punch force-stroke relationship.
Furthermore, since a strain-rate independent material model was assumed in the finite element
simulations, and to enhance the calculation speed, mass-scaling was used by setting the minimum
step increment time to 10-7 seconds. Stamping simulation of a T-shape carried out using CPFE
model with 120 orientations per integration point. Also, the random volume fraction distribution
and orientations of the phases were assumed for the sheet metal.

144

Figure 7.2: Schematic diagram for Finite element model of a T-shape stamping simulation.

7.3 FLC Validation and T-shape CPFE Simulation
Figure 7.3(a) and (b) illustrate the contour plot of the equivalent plastic strain and the
corresponding experimental results, respectively, when the punch force reaches the peak value in
the T-shape stamping of sample-1, as shown in Figure 7.3(c). Visual inspection of experimentally
stamped T-shape, Figure 7.3(b), showed that the part did not fail. It can be seen from simulation
results, Figure 7.3(a), that plastic strains are highest in the corner areas of the T-shape part where
there is severe bending and stretching deformations. However, failure of the T-shape part could
not be ascertained by visual inspection of the simulation results alone. Therefore, major and minor
strains at the final punch depth, Figure 7.3(c), were calculated based on plane stress conditions and
plotted on the forming limit diagram, as shown in Figure 7.4. Strains with critical values from
multiple regions (elements) with different colors (upper-right corner of Figure 7.4) for the
145

simulated T-shape were projected onto the forming limit diagram as shown in Figure 7.4.
According to the simulation results, and compared with the FLC predicted with CPFE model, some
of the major and minor strains plotted on the forming limit diagram indicate that the stamped Tshape might be experiencing necking. From the experimental perspective, cracks and wrinkles are
the type of failures that can be predicted from simulation results using fracture forming limit curve
(FFLC) instead of forming limit curve (FLC) (Chung et al., 2014). Hence, the onset of necking
cannot be identified just from visual inspection, without experimentally measured forming strains.
In addition, the critical strains that crossed the FLC in Figure 7.4 may have been due to the slightly
larger punch force predicted at the end of the simulation result, as shown in Figure 7.3(c). However,
overall, the predicted punch force-stroke curve based on the CPFE model and 0.12-friction
coefficient is in good agreement with the experimental curve, Figure 7.3(c).
The simulation results with the contour of the equivalent plastic strain and the experimental final
shape of the T-shape (sample-2) with holes, are shown in Figure 7.5(a) and (b), respectively. Here,
the simulation results from Figure 7.5(a) at ~17.5 mm punch stroke cannot be used for formability
assessment because the experimental punch force-stroke curve showed earlier fracture at
approximately 14.3 mm punch depth, Figure 7.5(c). Therefore, the strains were collected for three
punch depths prior to the maximum load occurring at the 14.3 mm punch depth. For the safe Tshape stamping at 12.5 mm punch depth, Figure 7.5 (c), it was noticed that true strains everywhere
in the sheet were less than the critical limit strains as defined by FLC in Figure 7.6(a). Also, it was
observed that the maximum plastic deformation generated around the holes with no contact
surfaces, gradually extended to adjacent regions around the holes as the stamping process
progressed to 13.5 mm punch depth, as shown in Figure 7.6 (b), when the sheet started to
experience the onset of necking around the holes. At approximately 14 mm punch depth, the CPFE
146

punch force-stroke curve overestimates the experimental curve. The calculated strains plotted on
the forming limit diagram in Figure 7.6(c) at 14 mm punch depth, indicate that the sheet metal is
close to fracture since many of the strain points are above the FLC with additional critical strains
also occurring at locations away from the holes.
As can be seen from these figures, the CPFE model was able to predict the safe strains, as well as
necking strains and failure locations in the sheet at various conditions. However, it can be seen
that fracture observed from the experimental test for sample-2, Figure 7.5(c), occurred at an earlier
punch stroke of 14.3 mm compared with the fracture from the punch force versus punch stroke
curves predicted by the CPFE model. The CPFE model in this work was not intended to account
for the material behavior beyond the strain localization. In other words, strain rate sensitivity and
damage behavior with element deletion were not considered in this study. Also, the extended Voce
type hardening rule shown in Eq. (3.11) converges to the linear hardening rule as the exponent
part vanishes with the increasing accumulated shear strain. Therefore, the material model could
not accurately predict the necking of the punch force versus punch stroke curves of sample-2.

147

(a)

(b)
Figure 7.3: CPFE simulation and experimental results of the T-shape stamping sample-1. (a)
Equivalent plastic strain distribution based on CPFE model; (b) Experimental final shape results;
(c) Comparison of punch force versus punch stroke relationship during the stamping process.

148

Figure 7.3 (cont’d)

(c)

149

Figure 7.4: The minor and major true strain distributions and FLC predicted by CPFE model of
T-shape sample-1 at the end of stamping process.

150

(a)

(b)

Figure 7.5: CPFE simulation and experimental results of the T-shape stamping sample-2.
(a) Equivalent plastic strain distribution based on CPFE model; (b) Experimental final shape
results; (c) Comparison of punch force versus punch stroke relationship during the stamping
process.

151

Figure 7.5 (cont’d)

(c)

152

(a)

(b)
Figure 7.6: The minor and major true strain distributions and FLC predicted by CPFE model of
T-shape sample-2 at punch depth (a) 12.5 mm; (b) 13.5 mm and (c) 14.0 mm.
153

Figure 7.6 (cont’d)

(c)

154

CHAPTER 8
THE ROLE OF CRYSTAL PLASTICITY MODELING IN CONTROLLING THE
DEFORMATION BEHAVIOR OF MULTIPHASE STEEL
Experimental observations using digital image correlation (DIC) based on scanning electron
microscopy (SEM) topography images indicated that both the morphology and the spatial
distribution of the phases in AHSSs can significantly affect the heterogeneity of the strain
distribution within the microstructure and, subsequently, the initiation of damage during uniaxial
tension (Ghadbeigi et al., 2010). Therefore, in this chapter, virtual representative volume elements
(RVEs) that depict the heterogeneous microstructure of multiphase Q&P980 steel were used to
perform a parametric study of the effects of the multiphase microstructure. The effect of the
crystallographic orientation on the deformation of Q&P980 steel was investigated using the
CPFEM. For comparison with the CPFEM results, an isotropic, elastoplastic von Mises FEM also
was used to simulate the heterogeneity of the strain–stress partitioning of Q&P980 steel without
considering the crystallographic orientation of the constituent phases.

8.1 Representative Volume Element (RVE) Simulation (Hexahedral Element)
The FE simulations were performed using RVE discretized with 20×20×20, eight-node continuum
elements of the type C3D8R. Each element represents a single grain. Therefore, a total of 8000
randomly-oriented grains are simulated. Based on the volume fraction of the phases in the RVE,
3600 grain orientations (out of a total of 8000) were assumed for the ferrite phase, and a similar
number was assumed for the martensite phase orientations, representing 45% of each phase. Also,
400 grain orientations were assumed for austenite and the new martensite, representing a volume
fraction of 5% for each of them, as shown in Figure 8.1(a). Within each grain (element), the
155

properties are assumed to be homogeneous and the RVE is initially stress free. The multiphase
steel microstructure information (volume fraction, phase type, crystal orientation) were embedded
within material subroutine VUMAT and defined at each integration point.
The boundary conditions considered for the uniaxial tension simulation were applied to the four
planes that comprised the 3D-RVE mesh, as shown in Figure 8.1(b). A displacement-controlled
loading was applied on the left-hand side of the RVE to a prescribed constant nominal strain rate
parallel to the x3 axis, while maintaining zero resultant forces on planes perpendicular to the x2 and
x3 axes. This allows the computation of the stress and strain distribution in the microstructure and
the prediction of the uniaxial tensile stress-strain curves for the multiphase steel microstructures.
For a comparison with the results of CPFEM, isotropic von Mises constitutive formulations also
were implemented into another VUMAT user subroutine. The stress–strain relationship of the
constituent phases were fitted for isotropic von Mises by using a Swift type (Swift, 1952)
hardening law, Eq. (4.2). Hardening parameters, which were fitted to each phase, are listed in
Table 8.1. In addition, the hardening parameters of the CPFE model, which were calibrated
assuming one grain per integration point, are listed in Table 8.2.

156

(a)

(b)
Figure 8.1: (a) A representative volume element (RVE) of multiphase steel (Ferrite: 45%,
Martensite: 45%, New Martensite: 5%, Austenite: 5%); (b) Finite element meshes and prescribed
displacement boundary conditions for uniaxial tension (U: translation, UR: rotation).

157

Table 8.1: Material parameters for the Swift equation hardening model composed of four phases
in Q&P980 steel.
Phase

𝑲 (MPa)

𝜺𝟎

𝒉

Ferrite

1504

0.01

0.2092

Martensite

1532

0.0001

0.0911

Austenite

2174

0.031

0.3325

New Martensite

4794

0.02

0.4416

Table 8.2: Material parameters for the CPFE model composed of four phases in Q&P980 steel
calibrated based on one orientation per gauss point (unit: MPa).
Phase

𝝉𝜶𝟎

𝝉𝜶𝟏

𝜽𝜶𝟎

𝜽𝜶𝟏

Ferrite

278

130

450

80

Martensite

440

60

1994

185

Austenite

472

1

800

450

New Martensite

440

40

2200

1300

158

8.1.1 Effect of Phases Volume Fraction on the Stress-Strain Curve
The roles of the constituent phases in controlling the overall response of the multiphase steel are
of particular interest. Therefore, we created several RVEs to probe the effects of the phases on the
uniaxial stress–strain response of the microstructure using CPFEM, and we compared the results
with isotropic FEM.
The simulations of the RVE, as detailed in the previous sections, were conducted to predict the
macroscopic flow curves for model materials in which the retained austenite was removed or
replaced with either ferrite, martensite, new martensite, or an equal volume fraction of martensite
and ferrite.
Figure 8.2(a) shows that the levels of the stress-strain curves predicted by CPFEM based on
changing the volume fraction of the constituents phase were, to some extent, different from those
calculated by the isotropic FEM shown in Figure 8.2(b). In Figure 8.2, the flow curve for the
material that contained no austenite and 10% new martensite has a higher initial yield stress and
exhibits high strength following macroscopic yielding. In contrast, the flow curves for the material
that contain 5% or 10% of austenite show that the retained austenite decreased the yield stress at
applied strains, thereafter increasing the flow stress by continuous work-hardening. This
improvement in strain hardening would be expected to delay the onset of localization and possibly
result in better ductility.
Figure 8.2(a) and (b) show how the variations in the volume fractions of the martensite and ferrite
phases influence the macroscopic flow behavior. The difference in the flow behavior due to a 10%
variation in the volume fraction of the martensite phase decreases with progressive tensile
deformation.
159

Collectively, the fitted flow curves based on CPFE and isotropic FE did not identify the effect of
the constituent phases on the total behavior of the material under plastic deformation. In the next
sections, we show how the initial crystallographic orientation and phase type of the constituent
phases significantly affected the micromechanical deformation and failure behaviors of multiphase
steel during simple uniaxial tension.

160

(a)

(b)
Figure 8.2: Macroscopic stress-strain responses of RVEs with varying the constituents phases
volume fraction. (a) CPFEM; (b) Isotropic FE.
161

8.1.2 Statistics of Stress and Strain Distributions
The internal stress and strain distributions in the RVEs with 45% ferrite, 45% martensite, 5% new
martensite, and 5% austenite were examined in more detail. A particular focus was on the statistics
of stress and strain measures, which are known to be related to the nucleation and growth of cracks
and voids and, subsequently, to damage. Statistical figures of the Mises stress, maximum principal
stress (which is related to nucleation of cracks), equivalent plastic strain (which is related to
nucleation of voids in a ductile matrix), and the stress triaxiality (the ratio of the first and second
stress invariants which is related to void growth in ductile fracture) for an RVE are presented in
this section.
Figure 8.3 shows the distribution of Mises stress within the constituent phases in Q&P980 steel
during uniaxial tension. Figure 8.3(a) shows that the range of Mises stress distribution, based on
CPFE simulation, in the constituent phases was broader than that determined by the isotropic FE,
Figure 8.3(b), which was relatively narrow to limited stresses values during uniaxial tension. The
ranges of Mises stress in the martensite phases and austenite were broader than that in the ferrite
phase during uniaxial tension. The relatively higher hardening rates combined with crystal
orientation of these phases seem to have contributed to the broader range of Mises stress
distribution with CPFE simulation. In contrast, the range of Mises stress distribution in the
martensite phase was similar to that in the ferrite phase, and the range of the new martensite phase
was similar to that of the austenite phase during uniaxial tension.
Figure 8.4 shows the distribution of the maximum principle stress within the constituent phases in
Q&P980 steel during uniaxial tension. Figure 8.4(b) shows that the range of the distribution of the
maximum principle stress in the constituent phases determined by the isotropic FE results was
162

relatively narrow compared with that determined by the CPFEM, Figure 8.4(a). Figure 8.4(a) and
(b) also indicate that the martensite phase was subjected to a much higher value of maximum
principle stress than the ferrite phase during tension.
The distribution of the equivalent plastic strain within the constituent phases in Q&P980 steel
during uniaxial tension is shown in Figure 8.5. The range of equivalent plastic strain in the ferrite
phase covered the entire range of the other phases in CPFE simulation, Figure 8.5(a). This result
indicates that all phases undergo a non-negligible plastic deformation to satisfy the strain
compatibility during uniaxial tension. However, it should be noted that the range of the equivalent
plastic strain distribution was relatively narrow in isotropic FE simulation, Figure 8.5(b), compared
with that of the equivalent plastic strain distribution, as shown in Figure 8.5(a).
Figure 8.6 shows the distribution of the triaxiality of the stress within the constituent phases in
Q&P980 steel during uniaxial tension. For uniaxial tensile loading, the triaxiality should
be 0.333 for a homogeneous material, but a wide range of triaxiality values can be observed for a
multiphase microstructure. The local triaxiality in the ferrite phase depends on the positions of the
adjacent phases. The deformation in the ferrite phase is constrained locally by the neighboring
hard phases, resulting in the buildup of triaxiality.
Figure 8.6(a) shows that, based on CPFE simulation, the range of the distribution of stress
triaxiality in the constituent phases was broader than that determined by the isotropic FE, Figure
8.6(b). The ferrite phase was much broader than the other phases in the CPFE simulation, Figure
8.6 (a). Figure 8.6 indicates that the stress triaxiality in ferrite can have higher values than in other
phases, therefore the probability of void nucleation is highest in the ferrite phase.

163

(a)

(b)
Figure 8.3: Mises stress distribution of the RVE: (a) CPFE; (b) Isotropic FE.

164

(a)

(b)
Figure 8.4: Max principle stress distribution of the RVE: (a) CPFE; (b) Isotropic FE.

165

(a)

(b)
Figure 8.5: Mises equivalent plastic strain distribution of the RVE: (a) CPFE; (b) Isotropic FE.

166

(a)

(b)
Figure 8.6: Stress triaxiality distribution of the RVE: (a) CPFE; (b) Isotropic FE.

167

8.2 Microstructure Model RVE Simulation of Different Grains Sizes of Constituents
Phases
The microstructure model was generated to analyze the effect of crystallographic orientation on
the micromechanical deformation and failure behaviors of multiphase steel during uniaxial tension
with distinguishable microstructural characteristics. The generation of the microstructure of the
Q&P980 was conducted using DREAM.3D software (Groeber and Jackson, 2014), which is
capable of creating synthetic microstructures using experimentally-determined microstructure
morphological and crystallographic statistics.
8.2.1 RVE Construction and Meshing
The DREAM.3D software package was used in conjugate with a developed MATLAB code to
build the finite element model of the multiphase metal microstructure. The MATLAB code that is
developed in this study utilizes the inbuilt microstructure using DREAM.3D software from EBSD
data into an ABAQUS input file. The code supplies the information of the generated microstructure
based on the desired parameters of the material subroutine (e.g., crystallographic orientation,
crystal structures, slip system, and hardening parameters).
The grain sizes of the distribution of the three phases computed from the EBSD image using
MTEX software were used to generate the synthetic microstructures of the three-dimensional
grain-size distribution using log-normal fitting (Figure 8.7). The newly-formed martensite was
assumed to have existed from the beginning, so the grains sizes of the new martensite phase were
assumed to be similar to the sizes of the retained austenite grains. Also, random orientations of all
grains were assumed in the Dream.3D RVE.

168

A cubical, synthetic microstructure was built in 3D using DREAM.3D, as shown Figure 8.8(a), for
four constituents phases, Figure 8.8(b), that consisted of 2226 grains each with a specific initial
lattice orientation. The 2226 grains were meshed with 1,728,000 brick elements with reduced
integration (C3D8R), as shown in Figure 8.8(c). For uniaxial tension, a similar boundary condition,
as described in Section 8.1, was applied to the RVE constructed using Dream.3D.

Figure 8.7: Fitted log-normal grain size distribution used in DREAM.3D to generate equiaxed
grain structure.

169

(a)

(b)

(c)

Figure 8.8: A four-phase representative volume element of Q&P980 steel: (a) Synthetic
polycrystalline aggregate generated in DREAM.3D consisting of 3D voxels, and generated from
our integrated tool set; (b) full phase representation (Ferrite: 45%, Martensite: 45%, New
Martensite: 5%, Austenite: 5%); (b) Abaqus 3D RVE undeformed meshes of the polycrystalline
aggregate used in simulating uniaxial tension.
170

8.2.2 RVE Simulation Results
Figure 8.9 compares the CPFE and isotropic-FE simulation results of the RVE multiphase steel
during uniaxial tension. Figure 8.9(a) and (c) clearly show the effect of the grain orientation on
Mises and maximum principle stresses compared to the isotropic FE results in Figure 8.9(b) and
(d). When comparing Figure 8.9(a) and (b) to Figure 8.8(b), it is clear that the martensite phases
were subjected to a much higher stress than the ferrite and austenite phases during uniaxial tension.
However, note that in Figure 8.9(a), some ferrite grains (red circles) have high stress values with
the same level of the martensite phase, while some martensite grains (orange circles) show less
stress, with values similar to the majority of ferrite phase stress. That effect is related to the grain
and its neighbor’s orientations, which can have a significant effect on the stress and strain level of
that individual grain. When Figure 8.9(a) is compared to Figure 8.9(b), it is clear that this effect
was neglected in the isotropic FE simulations with a homogenous distribution within the individual
phases. Figure 8.9(a) and (c) show that the deformed microstructure exhibited higher
heterogeneous distribution of stresses than were shown in Figure 8.9(b) and (d). This effect of
heterogeneity also is evident in Figure 8.9 (e), (f), (g), and (h) for the equivalent plastic strain and
triaxiality distribution.
In order to explain the individual behaviors of the phases and the crystallographic effect in the
multiphase steel, individual regions of constituent phases were separated and shown in four figures.
Figure 8.10, 8.11, 8.12 and 8.13 show the distribution of Mises stress, maximum principle stress,
equivalent plastic strain and traiaxiality within the constituent phases during uniaxial tension based
on the CPFE and isotropic FE simulations. From these figures, it can be clearly show that isotropic
FE results which are based on phases contrast only cannot show the significant contributions of
the crystallographic orientation as in CPFE simulations.
171

Martensite
Ferrite

(a)

(b)

(c)

(d)

Figure 8.9: Simulation results of a four-phase representative volume element of Q&P980 steel at
12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution of maximum
principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent plastic strain: (e) CPFE, (f)
Isotropic FE. Distribution of stress triaxiality: (g) CPFE, (h) Isotropic FE.
172

Figure 8.9 (cont’d)

(e)

(f)

(g)

(h)

173

(a)

(b)

(c)

(d)

Figure 8.10: Simulation results of a ferrite phase in the representative volume element of Q&P980
steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution of
maximum principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent plastic strain: (e)
CPFE, (f) Isotropic FE. Distribution of stress triaxiality: (g) CPFE, (h) Isotropic FE.
174

Figure 8.10 (cont’d)

(e)

(f)

(g)

(h)

175

(a)

(b)

(c)

(d)

Figure 8.11: Simulation results of a martensite phase in the representative volume element of
Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution
of maximum principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent plastic strain:
(e) CPFE, (f) Isotropic FE. Distribution of stress triaxiality: (g) CPFE, (h) Isotropic FE.

176

Figure 8.11 (cont’d)

(e)

(f)

(g)

(h)

177

(a)

(b)

(c)

(d)

Figure 8.12: Simulation results of a new martensite phase in the representative volume element of
Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution
of maximum principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent plastic strain:
(e) CPFE, (f) Isotropic FE. Distribution of stress triaxiality: (g) CPFE, (h) Isotropic FE.

178

Figure 8.12 (cont’d)

(e)

(f)

(g)

(h)

179

(a)

(b)

(c)

(d)

Figure 8.13: Simulation results of the austenite phase in the representative volume element of
Q&P980 steel at 12% strain. Distribution of Mises stress: (a) CPFE, (b) Isotropic FE. Distribution
of maximum principle stress: (c) CPFE, (d) Isotropic FE. Distribution of equivalent plastic strain:
(e) CPFE, (f) Isotropic FE. Distribution of stress triaxiality: (g) CPFE, (h) Isotropic FE.

180

Figure 8.13 (cont’d)

(e)

(f)

(g)

(h)

181

CHAPTER 9
CONCLUSIONS AND FUTURE WORKS
9.1 Conclusions
The microstructures of AHSSs have become increasingly complex, incorporating typically more
than one microstructural feature to adjust the material’s properties according to specific needs.
This is especially true when more than one phase is present, and multiple phases usually are
introduced to improve strength and ductility. To simulate the behavior of such complex materials,
the underlying effects, such as the interactions of the different phases, hardening of the slip systems,
and the distribution of the phases, must be considered to produce useful results. In this context, the
multiphase CPFE model, which is believed to be the most efficient and fastest crystal plasticity
model currently available, has proven to be a powerful tool for determining the local and global
mechanical response of multiphase AHSSs.
The experimental results of Q&P980 indicate that, for these new AHSSs, especially those with
complex microstructures, such as Q&P980, a comprehensive experimental investigation is
required to determine the effects of the features of the microstructure. According to all of the
observations of the strain distribution in uniaxial tensile tests, the formability tests, and the
metallurgical analysis of the Q&P980 steel, it was found that the inconsistency of formability was
attributable to the microstructural features of the sheet metal. Forming limit strains based on
experiments associated with the onset of localized necking within a local area for a specified
deformation path cannot accurately represent the formability involved during a forming process.
A multiscale and multiphase CPFE model that accounted for the grain rotation and the evolution
of anisotropy was applied to large deformations in the Nakajima hemispherical punch test. Two
182

different methods were considered to represent the inhomogeneous distribution of phases in steel,
and it was found that the random phase distribution had more accurate results than the fixed phase
distribution. The results of the CPFE simulation of various deformation modes, including uniaxial
tension, plane strain tension, and balanced biaxial tension, showed good fits with the experimental
data in terms of punch force-displacement in the Nakajima hemispherical punch test. The FLC
predicted based on the LBF method showed good results for determining the limit strains for the
uniaxial and balanced biaxial cases, but the method provided results that were slightly higher than
the experimental results for the in-between deformation paths. That overestimation was related to
the inhomogeneity effect, which was solved by incorporating the Marciniak-Kuczynski (MK)
approach and finite element methodology to construct the FLC based on the CPFE simulation.
The numerical study on the effect of the mechanical properties of the constituent phases on the
FLC suggested the possibility of designing forming processes that would lead to textures that
enhance the resistance to localization. In addition, the accurate prediction of strains at the onset of
necking of the T-shape stamped part offers a promising microstructure-sensitive model that will
be highly desirable for the stamping simulation of the third-generation, advanced, high-strength
steels.
Crystal plasticity, in particular, has been used extensively as a highly-capable framework for
studying the effects of the microstructure on plastic deformation in metals, but most attempts to
model AHSSs still rely on phenomenological models and simple 2D-RVE simulations. This is not
sufficient for predicting the complex behaviors of AHSSs because it ignores the crystallographic
texture and the phase type that influences the mechanisms of deformation of the constituent phases
in AHSSs.
183

9.2 Recommendations for Future Work
The work presented in this study provides an excellent tool for computer-aided engineering and
material characterization. Unlike phenomenologically-based material models, this tool could be
used to develop a better understanding of the behavior of multiphase steel based on a sophisticated,
multiphase microstructure model, rather than simply fitting data for modeling. Several
recommendations for further improvement of the multiphase CPFE model and future
developments in this area are provided below:
1- The forming limit curve predicted by the hemispherical punch test can be improved to yield
better results by increasing the inhomogeneity of the material. That can be achieved by using
suitable distribution functions for the distribution of the phases in the sheet metal that can be
implemented into the material subroutine VUMAT.
2- When the steel has a high-volume fraction of austenite, the multiphase CPFE model can be
modified to account for the transformation of the martensitic phase that is induced by the plastic
deformation of the retained austenite. That can be achieved by considering the invariant shear
deformation of the lattice and the relationships between the orientations of the parent austenite and
the transformed martensite. For the computational efficiency of the CPFE model and because of
the small sizes of austenite grains, another method can be used in which the isotropic FE phase
transformation for the parent austenite is used in conjunction with the CPFE deformation for ferrite
and martensite.
3- Another important future research topic is to consider the prediction of damage using the
multiphase CPFE model. Additional work can be conducted to account for isotropic damage in the
proposed CPFE model. However, for more sophisticated cases, damage can be accounted for in
184

the CPFE model by using several micro-scale damage parameters of the constituent’s phases,
which evolve on the slip systems with deformation or time. These micro-scale parameters can be
a function of temperature and microstructural parameters. Then, these modified CPFE models can
be used to obtain macro-scale damage parameters for the material.
4- A 3D-RVE developed using Dream.3D of multiphase microstructure can be used to study the
effect of grain size and crystallographic orientation to provide more information about the effects
of the features of the microstructure. Also, it can be used to calculate the anisotropy coefficients
of phenomenological yield functions (such as the Yld2000 or Yld2004 yield functions) using the
CPFE modeling.
5- The CPFE model also could be extended to non-cubic crystalline materials (HCP structure),
such as magnesium, titanium, or zinc alloys.

185

BIBLIOGRAPHY

186

BIBLIOGRAPHY

ABAQUS User’s Manual, 2000. version 6.14, Hibbit, Karlsson & Sorenson, Providence, RI.
Abid, N.H., Abu Al-Rub, R.K., Palazotto, A.N., 2017. Micromechanical finite element analysis of
the effects of martensite morphology on the overall mechanical behavior of dual phase steel.
Int. J. Solids Struct. 104–105, 8–24.
Abid, N.H., Abu Al-Rub, R.K., Palazotto, A.N., 2015. Computational modeling of the effect of
equiaxed heterogeneous microstructures on strength and ductility of dual phase steels.
Comput. Mater. Sci. 103, 20–37.
Alaie, A., Ziaei Rad, S., Kadkhodapour, J., Jafari, M., Asadi Asadabad, M., Schmauder, S., 2015.
Effect of microstructure pattern on the strain localization in DP600 steels analyzed using
combined in-situ experimental test and numerical simulation. Mater. Sci. Eng. A 638, 251–
261.
Alharbi, K., Ghadbeigi, H., Efthymiadis, P., Zanganeh, M., Celotto, S., Dashwood, R., Pinna, C.,
2015. Damage in dual phase steel DP1000 investigated using digital image correlation and
microstructure simulation. Model. Simul. Mater. Sci. Eng. 23, 85005.
Amirkhizi, A. V., Nemat-Nasser, S., 2007. A framework for numerical integration of crystal elastoplastic constitutive equations compatible with explicit finite element codes. Int. J. Plast. 23,
1918–1937.
Amirmaleki, M., Samei, J., Green, D.E., van Riemsdijk, I., Stewart, L., 2016. 3D micromechanical
modeling of dual phase steels using the representative volume element method. Mech. Mater.
101, 27–39.
Anand, L., Kothari, M., 1996. A computational procedure for rate-independent crystal plasticity. J.
Mech. Phys. Solids 44, 525–558.
Anderson, D., Butcher, C., Pathak, N., Worswick, M.J., 2017. Failure parameter identification and
validation for a dual-phase 780 steel sheet. Int. J. Solids Struct. 124, 89–107.
Asaro, R.J., Needleman, A., 1985. Overview no. 42 Texture development and strain hardening in
rate dependent polycrystals. Acta Metall. 33, 923–953.
Asgari, S.A., Hodgson, P.D., Yang, C., Rolfe, B.F., 2009. Modeling of Advanced High Strength
Steels with the realistic microstructure-strength relationships. Comput. Mater. Sci. 45, 860–
866.
ASTM Int., 2009. Standard Test Methods for Tension Testing of Metallic Materials 1. Astm 1–27.
187

Banabic, D., 2010a. Sheet Metal Forming Processes: Constitutive Modelling and Numerical
Simulation. Springer Berlin Heidelberg.
Banabic, D., 2010b. A review on recent developments of Marciniak-Kuczynski model. Comput.
Methods Mater. Sci. 1–24.
Banabic, D., Aretz, H., Comsa, D.S., Paraianu, L., 2005. An improved analytical description of
orthotropy in metallic sheets. Int. J. Plast. 21, 493–512.
Banabic, D., Balan, T., Comsa, D.., 2000. A new yield criterion for orthotropic sheet metals under
plane-stress conditions, in: 7th Cold Metal Forming Conference, Cluj Napoca, Romania. pp.
217–224.
Banu, M., Takamura, M., Hama, T., Naidim, O., Teodosiu, C., Makinouchi, A., 2006. Simulation
of springback and wrinkling in stamping of a dual phase steel rail-shaped part. J. Mater.
Process. Technol. 173, 178–184.
Barlat, F., Aretz, H., Yoon, J.W., Karabin, M.E., Brem, J.C., Dick, R.E., 2005. Linear
transfomation-based anisotropic yield functions. Int. J. Plast. 21, 1009–1039.
Barlat, F., Brem, J.C., Lege, D.J., Chung, K., 1997a. Advanced Methods in Materials Processing
Defects. Stud. Appl. Mech. 45, 265–272.
Barlat, F., Brem, J.C., Yoon, J.W., Chung, K., Dick, R.E., Lege, D.J., Pourboghrat, F., Choi, S.H.,
Chu, E., 2003. Plane stress yield function for aluminum alloy sheets—part 1: theory. Int. J.
Plast. 19, 1297–1319.
Barlat, F., Lege, D.J., Brem, J.C., 1991. A six-component yield function for anisotropic materials.
Int. J. Plast. 7, 693–712.
Barlat, F., Lian, K., 1989a. Plastic behavior and stretchability of sheet metals. Part I: A yield
function for orthotropic sheets under plane stress conditions. Int. J. Plast. 5, 51–66.
Barlat, F., Lian, K., 1989b. Plastic behavior and stretchability of sheet metals. Part I: A yield
function for orthotropic sheets under plane stress conditions. Int. J. Plast. 5, 51–66.
Barlat, F., Maeda, Y., Chung, K., Yanagawa, M., Brem, J.C., Hayashida, Y., Lege, D.J., Matsui, K.,
Murtha, S.J., Hattori, S., Becker, R.C., Makosey, S., 1997b. Yield function development for
aluminum alloy sheets. J. Mech. Phys. Solids 45, 1727–1763.
Barlat, F., Richmond, O., 1987. Prediction of tricomponent plane stress yield surfaces and
associated flow and failure behavior of strongly textured f.c.c. polycrystalline sheets. Mater.
Sci. Eng. 95, 15–29.
Bassani, J.L., 1977. Yield characterization of metals with transversely isotropic plastic properties.
Int. J. Mech. Sci. 19, 651–660.
188

Becker, R., 1991. Analysis of texture evolution in channel die compression-I. Effects of grain
interaction. Acta Metall. Mater. 39, 1211–1230.
Becker, R., Butler, J.F., Hu, H., Lalli, L.A., 1991. Analysis of an aluminum single crystal with
unstable initial orientation (001) [110] in channel die compression. Metall. Trans. A 22, 45–
58.
Ben Bettaieb, M., DĂŠbordes, O., Dogui, A., DuchĂŞne, L., Keller, C., 2012. On the numerical
integration of rate independent single crystal behavior at large strain. Int. J. Plast. 32–33, 184–
217.
Bishop, J.F.W., Hill, R., 1951a. A theory of the plastic distortion of a polycrystalline aggregate
under combined stresses. Philos. Mag. Ser. 7 42, 414–427.
Bishop, J.F.W., Hill, R., 1951b. A theory of the plastic distortion of a polycrystalline aggregate
under combined stresses. Philos. Mag. Ser. 7 42, 414–427.
Bishop, J.F.W., Hill, R., 1951c. CXXVIII. A theoretical derivation of the plastic properties of a
polycrystalline face-centred metal. London, Edinburgh, Dublin Philos. Mag. J. Sci. 42, 1298–
1307.
Boogaard, A.H. Van Den, Hu, J., n.d. Prediction of sheet necking with shell finite element models.
Branagan, D.J., 2014. NanoSteel 3rd Generation AHSS : Auto Evaluation and Technology
Expansion. Gt. Des. Steel.
Budianski, B., 1984. Anisotropic Plasticity of Plane-Isotropic Sheets. Mech. Mater. Behav. 15–29.
Bunge, H. J., 1982. Texture analysis in materials science. Butterworths Publ., London.
Calcagnotto, M., Ponge, D., Demir, E., Raabe, D., 2010. Orientation gradients and geometrically
necessary dislocations in ultrafine grained dual-phase steels studied by 2D and 3D EBSD.
Mater. Sci. Eng. A 527, 2738–2746.
Cao, J., Yao, H., Karafillis, A., Boyce, M.C., 2000. Prediction of localized thinning in sheet metal
using a general anisotropic yield criterion. Int. J. Plast. 16, 1105–1129.
Cazacu, O., Barlat, F., 2004. A criterion for description of anisotropy and yield differential effects
in pressure-insensitive metals. Int. J. Plast. 20, 2027–2045.
Charoensuk, K., Panich, S., Uthaisangsuk, V., 2017. Damage initiation and fracture loci for
advanced high strength steel sheets taking into account anisotropic behaviour. J. Mater.
Process. Tech. 248, 218–235.

189

Chen, P., Ghassemi-Armaki, H., Kumar, S., Bower, A., Bhat, S., Sadagopan, S., 2014. Microscalecalibrated modeling of the deformation response of dual-phase steels. Acta Mater. 65, 133–
149.
Chiba, R., Takeuchi, H., Kuroda, M., Hakoyama, T., Kuwabara, T., 2013. Theoretical and
experimental study of forming-limit strain of half-hard AA1100 aluminium alloy sheet.
Comput. Mater. Sci. 77, 61–71.
Choi, K.S., Liu, W.N., Sun, X., Khaleel, M.A., 2009. Microstructure-based constitutive modeling
of TRIP steel: Prediction of ductility and failure modes under different loading conditions.
Acta Mater. 57, 2592–2604.
Choi, K.S., Soulami, A., Liu, W.N., Sun, X., Khaleel, M.A., 2010. Influence of various material
design parameters on deformation behaviors of TRIP steels. Comput. Mater. Sci. 50, 720–
730.
Choi, S.H., Kim, E.Y., Woo, W., Han, S.H., Kwak, J.H., 2013. The effect of crystallographic
orientationonthe micromechanical deformation and failure behaviors of DP980 steel during
uniaxial tension. Int. J. Plast. 45, 85–102.
Chung, K., Lee, C., Kim, H., 2014. Forming limit criterion for ductile anisotropic sheets as a
material property and its deformation path insensitivity, Part II: Boundary value problems.
Int. J. Plast. 58, 35–65.
Chung, K., Ma, N., Park, T., Kim, D., Yoo, D., Kim, C., 2011. A modified damage model for
advanced high strength steel sheets. Int. J. Plast. 27, 1485–1511.
Coates, G., 2012. Importance of Materials and Manufacturing Emissions for Future Vehicle
Considerations.
Comsa, D.S., Banabic, D., 2008. Plane-stress yield criterion for highly-anisotropic sheet metals,
in: NUMISHEET 2008, Switzerland.
Cong, Z.H., Jia, N., Sun, X., Ren, Y., Almer, J., Wang, Y.D., 2009. Stress and strain partitioning of
ferrite and martensite during deformation. Metall. Mater. Trans. A Phys. Metall. Mater. Sci.
40, 1383–1387.
De Moor, E., Gibbs, P.J., Speer, J.G., Matlock, D.K., Schroth, J.G., 2010. Strategies for thirdgeneration advanced high-strength steel development.
Delannay, L., Doghri, I., Pierard, O., 2007. Prediction of tension-compression cycles in multiphase
steel using a modified incremental mean-field model. Int. J. Solids Struct. 44, 7291–7306.
Demeri, M.Y., 2013. Advanced high-strength Steels. ASM International.

190

Cyr, E., Mohammadi, M., Brahme, A., Mishra, R.K., Inal, K., 2017. Modeling the formability of
aluminum alloys at elevated temperatures using a new thermo-elasto-viscoplastic crystal
plasticity framework. Int. J. Mech. Sci. 128–129, 312–325.
Engler, O., Randle, V., 2010. Descriptors of Orientation. Introd. to Texture Anal. Macrotexture,
Microtexture, Orientat. Mapp. 15–50.
Erice, B., Roth, C.C., Mohr, D., 2017. Stress-state and strain-rate dependent ductile fracture of
dual and complex phase steel. Mech. Mater. 0, 1–22.
Eyckens, P., Mulder, H., Gawad, J., Vegter, H., Roose, D., van den Boogaard, T., Van Bael, A., Van
Houtte, P., 2015. The prediction of differential hardening behaviour of steels by multi-scale
crystal plasticity modelling. Int. J. Plast. 73, 1–23.
Eyckens, P., Van Bael, A., Van Houtte, P., 2009. Marciniak-Kuczynski type modelling of the effect
of Through-Thickness Shear on the forming limits of sheet metal. Int. J. Plast. 25, 2249–2268.
Fine, C., Roth, R., 2010. Lightweight Materials for Transport: Developing a Vehicle Technology
Roadmap for the Use of Lightweight Materials.
Fonstein, N., 2015. Advanced High Strength Sheet Steels, Advanced High Strength Sheet Steels.
Franciosi, P., Zaoui, A., 1991. Crystal hardening and the issue of uniqueness. Int. J. Plast. 7, 295–
311.
Gambin, W., 1992. Refined analysis of elastic-plastic crystals. Int. J. Solids Struct. 29, 2013–2021.
Ghadbeigi, H., Pinna, C., Celotto, S., Yates, J.R., 2010. Local plastic strain evolution in a high
strength dual-phase steel. Mater. Sci. Eng. A 527, 5026–5032.
Ghassemi-Armaki, H., Chen, P., Bhat, S., Sadagopan, S., Kumar, S., Bower, A., 2013. Microscalecalibrated modeling of the deformation response of low-carbon martensite. Acta Mater. 61,
3640–3652.
Goodwin, G.M., 1968. Application of strain analysis to sheet metal forming in the press shop. SAE
Pap. 680093, 767–774.
Gotoh, M., 1977. A theory of plastic anisotropy based on yield function of fourth order (plane
stress state)-II. Int. J. Mech. Sci. 19, 513–520.
Govik, A., Rentmeester, R., Nilsson, L., 2014. A study of the unloading behaviour of dual phase
steel. Mater. Sci. Eng. A 602, 119–126.
Graf, A., Hosford, W., 1993. Effect of Changing Strain Paths on Forming Limit Diagrams of Al
2008-T4. Metall. Trans. a-Physical Metall. Mater. Sci. 24, 2503–2512.
191

Graf Alejandro, H.W., 1994. The influence of strain-path changes on formint limit diagrams of Al
6111-T 4. Int. J. Mech. Sci. 36, 897–910.
Groeber, M. a, Jackson, M. a, 2014. DREAM.3D: A Digital Representation Environment for the
Analysis of Microstructure in 3D. Integr. Mater. Manuf. Innov. 3, 5.
Gruben, G., Hopperstad, O.S., Børvik, T., 2013. Simulation of ductile crack propagation in dualphase steel. Int. J. Fract. 180, 1–22.
Gruben, G., Morin, D., Langseth, M., Hopperstad, O.S., 2017. European Journal of Mechanics A
/ Solids Strain localization and ductile fracture in advanced high-strength steel sheets. Eur. J.
Mech. / A Solids 61, 315–329.
Grujicic, M., Batchu, S., 2002. Crystal plasticity analysis of earing in deep-drawn OFHC copper
cups. J. Mater. Sci. 37, 753–764.
Guan, Y., Pourboghrat, F., Barlat, F., 2006. Finite element modeling of tube hydroforming of
polycrystalline aluminum alloy extrusions. Int. J. Plast. 22, 2366–2393.
Habibi, N., Zarei-Hanzaki, A., Abedi, H.R., 2015. An investigation into the fracture mechanisms
of twinning-induced-plasticity steel sheets under various strain paths. J. Mater. Process.
Technol. 224, 102–116.
Hamelin, C.J., Diak, B.J., Pilkey, A.K., 2011. Multiscale modelling of the induced plastic
anisotropy in bcc metals. Int. J. Plast. 27, 1185–1202.
Han, Q., Asgari, A., Hodgson, P.D., Stanford, N., 2014. Strain partitioning in dual-phase steels
containing tempered martensite. Mater. Sci. Eng. A 611, 90–99.
Harren, S. V., Asaro, R.J., 1989. Nonuniform deformations in polycrystals and aspects of the
validity of the Taylor model. J. Mech. Phys. Solids 37, 191–232.
Harren, S. V., Dève, H.E., Asaro, R.J., 1988. Shear band formation in plane strain compression.
Acta Metall. 36, 2435–2480.
Hector, L.G., 2013. The Next Generation of Advanced High Strength Steels – Computation ,
Product Design and Performance.
Hershey, A., 1954. The Plasticity of an Isotropic Aggregate of Anisotropic Face Centered Cubic
Crystals. J. Appl. Mech. Trans ASME 21, 241.
Hill, R., 1993. A user-friendly theory of orthotropic plasticity in sheet metals. Int. J. Mech. Sci. 35,
19–25.
Hill, R., 1990. Constitutive modelling of orthotropic plasticity in sheet metals. J. Mech. Phys.
Solids 38, 405–417.
192

Hill, R., 1979. Theoretical Plasticity of Textured Aggregates. Math. Proc. Cambridge Philos. Soc.
85, 179–191.
Hill, R., 1952. On discontinuous plastic states, with special reference to localized necking in thin
sheets. J. Mech. Phys. Solids 1, 19–30.
Hill, R., 1950. The mathematical theory of plasticity. Clarendon Press, Oxford.
Hill, R., 1948. A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proc. R. Soc.
London A Math. Phys. Eng. Sci. 193, 281–297.
Hoon, J., Hyun, J., Piao, K., Wagoner, R.H., 2011. The shear fracture of dual-phase steel. Int. J.
Plast. 27, 1658–1676.
Horvath, c. . d., 2010. Advanced Steels for Lightweight Automotive Structures. Materials,
Design and Manufacturing for Lightweight Vehicles.
Hosford, W.F., 1972. A Generalized Isotropic Yield Criterion. J. Appl. Mech.
Hotz, W., Merklein, M., Kuppert, A., Friebe, H., Klein, M., 2013. Time dependent FLC
determination – Comparison of different algorithms to detect the onset of unstable necking
before fracture 549, 397–404.
Hu, B., Luo, H., Yang, F., Dong, H., 2017. Recent progress in medium-Mn steels made with new
designing strategies, a review. J. Mater. Sci. Technol. 6–13.
Hu, X., Choi, K.S., Sun, X., Ren, Y., Wang, Y., 2016a. Determining Individual Phase Flow
Properties in a Quench and Partitioning Steel with In Situ High-Energy X-Ray Diffraction
and Multiphase Elasto-Plastic Self-Consistent Method. Metall. Mater. Trans. A.
Hu, X., Choi, K.S., Sun, X., Ren, Y., Wang, Y., 2016b. Determining Individual Phase Flow
Properties in a Quench and Partitioning Steel with In Situ High-Energy X-Ray Diffraction
and Multiphase Elasto-Plastic Self-Consistent Method. Metall. Mater. Trans. A Phys. Metall.
Mater. Sci. 1–17.
Hu, Z.G., Zhu, P., Meng, J., 2010. Fatigue properties of transformation-induced plasticity and dualphase steels for auto-body lightweight: Experiment, modeling and application. Mater. Des.
31, 2884–2890.
Huh, J., Huh, H., Lee, C.S., 2013. Effect of strain rate on plastic anisotropy of advanced high
strength steel sheets. Int. J. Plast. 44, 23–46.
Hull, D., Bacon, D.J., 2011. Introduction to Dislocations, Fifth Edit. ed. Butterworth-Heinemann.
Hutchinson, J. W., Neale, K.W., 1978. Sheet Necking II: time independent behaviour, in: Koistinen
K. P., Wang N. -M. (Eds.). Mech. Sheet Met. Forming, Plenum, New York 127–153.
193

Inal, K., Neale, K.W., Aboutajeddine, a., 2005. Forming limit comparisons for FCC and BCC
sheets. Int. J. Plast. 21, 1255–1266.
International Standard ISO 12004-2, 2008. Metallic Materials-sheet and Strip Determination of
Forming Limit Curves. Part 2: Determination of Forming Limit Curves in the Laboratory.
International Organization for Standardization, Geneva, Switzerland.
Janssens, K., Lambert, F., Vanrostenberghe, S., Vermeulen, M., 2001. Statistical evaluation of the
uncertainty of experimentally characterised forming limits of sheet steel. J. Mater. Process.
Technol. 112, 174–184.
Jia, N., Cong, Z.H., Sun, X., Cheng, S., Nie, Z.H., Ren, Y., Liaw, P.K., Wang, Y.D., 2009. An in
situ high-energy X-ray diffraction study of micromechanical behavior of multiple phases in
advanced high-strength steels. Acta Mater. 57, 3965–3977.
Jia, N., Lin Peng, R., Wang, Y.D., Johansson, S., Liaw, P.K., 2008. Micromechanical behavior and
texture evolution of duplex stainless steel studied by neutron diffraction and self-consistent
modeling. Acta Mater. 56, 782–793.
John Neil, C., Agnew, S.R., 2009. Crystal plasticity-based forming limit prediction for non-cubic
metals: Application to Mg alloy AZ31B. Int. J. Plast. 25, 379–398.
Kadkhodapour, J., Anbarlooie, B., Hosseini-Toudeshky, H., Schmauder, S., 2014. Simulation of
shear failure in dual phase steels using localization criteria and experimental observation.
Comput. Mater. Sci. 94, 106–113.
Kadkhodapour, J., Butz, A., Ziaei Rad, S., 2011. Mechanisms of void formation during tensile
testing in a commercial, dual-phase steel. Acta Mater. 59, 2575–2588.
Karafillis, A.P., Boyce, M.C., 1993. A general anisotropic yield criterion using bounds and a
transformation weighting tensor. J. Mech. Phys. Solids 41, 1859–1886.
Karlsson, B., SundstrĂśm, B.O., 1974. Inhomogeneity in plastic deformation of two-phase steels.
Mater. Sci. Eng. 16, 161–168.
Keeler, S., Kimchi, M., 2014. Advanced High-Strength Steels Application Guidelines Version 5.0.
World AutoSteel.org 511.
Keeler, S., Kimchi, M., J Mooney, P., 2017. Advanced High-Strength Steels Apllication Guidelines
Version 6.0 314.
Keeler, S.P., 1965. Determination of Forming Limits in Automotive Stampings. SAE Tech. Pap.
14B, 656–669. doi:10.4271/650535
Kim, E.-Y., Choi, S.-H., Shin, E.-J., Yoon, J., 2012. Simulation of earing behaviors in bake
hardening steel exhibiting a strong off-γ-fiber component. Int. J. Solids Struct. 49, 3573–3581.
194

Kim, J.H., Kim, D., Barlat, F., Lee, M.-G., 2012. Crystal plasticity approach for predicting the
Bauschinger effect in dual-phase steels. Mater. Sci. Eng. A 539, 259–270.
Kim, J.H., Lee, M.G., Kim, D., Barlat, F., 2013. Numerical procedures for predicting localization
in sheet metals using crystal plasticity. Comput. Mater. Sci. 72, 107–115.
Kim, S., Lee, J., Barlat, F., Lee, M.G., 2013. Formability prediction of advanced high strength
steels using constitutive models characterized by uniaxial and biaxial experiments. J. Mater.
Process. Technol. 213, 1929–1942.
Knockaert, R., Chastel, Y., Massoni, E., 2002. Forming limits prediction using rate-independent
polycrystalline plasticity. Int. J. Plast. 18, 231–247.
Knockaert, R., Chastel, Y., Massoni, E., 2000. Rate-independent crystalline and polycrystalline
plasticity, application to FCC materials. Int. J. Plast. 16, 179–198.
Kraska, M., Doig, M., Tikhomirov, D., Raabe, D., Roters, F., 2009. Virtual material testing for
stamping simulations based on polycrystal plasticity. Comput. Mater. Sci. 46, 383–392.
Lankford, W.T., Snyder, S.C., Bausher, J.A., 1950. New criteria for predicting the press
performance of deep drawing sheets. Trans. Am. Soc. Met. 42, 1197–1232.
Latypov, M.I., Shin, S., De Cooman, B.C., Kim, H.S., 2016. Micromechanical finite element
analysis of strain partitioning in multiphase medium manganese TWIP+TRIP steel. Acta
Mater. 108, 219–228.
Leacock, A.G., 2006. A mathematical description of orthotropy in sheet metals. J. Mech. Phys.
Solids 54, 425–444.
LĂŠvesque, J., Mohammadi, M., Mishra, R.K., Inal, K., 2015. An extended Taylor model to simulate
localized deformation phenomena in magnesium alloys. Int. J. Plast. 78.
Li, D.Y., Peng, Y.H., Zhang, S.R., Liu, S.R., 2009. Numerical simulation of sheet metal stamping
by using ODF data. Int. J. Mech. Sci. 51, 41–51.
Lian, J., Barlat, F., Baudelet, B., 1989. Plastic behaviour and stretchability of sheet metals. Part II:
Effect of yield surface shape on sheet forming limit. Int. J. Plast. 5, 131–147.
Lian, J., Feng, Y., MĂźnstermann, S., 2014. A Modified Lemaitre Damage Model
Phenomenologically Accounting for the Lode Angle Effect on Ductile Fracture. Procedia
Mater. Sci. 3, 1841–1847.
Liedl, U., Traint, S., Werner, E.A., 2002. An unexpected feature of the stress-strain diagram of
dual-phase steel. Comput. Mater. Sci. 25, 122–128.

195

Lim, H., Lee, M.G., Sung, J.H., Kim, J.H., Wagoner, R.H., 2012. Time-dependent springback of
advanced high strength steels. Int. J. Plast. 29, 42–59.
Lin, S.B., Ding, J.L., 1996. A modified form of Hill’s orientation-dependent yield criterion for
orthotropic sheet metals. J. Mech. Phys. Solids 44, 1739–1764.
Lou, Y., Huh, H., 2013. Prediction of ductile fracture for advanced high strength steel with a new
criterion: Experiments and simulation. J. Mater. Process. Technol. 213, 1284–1302.
Luo, M., Wierzbicki, T., 2010. Numerical failure analysis of a stretch-bending test on dual-phase
steel sheets using a phenomenological fracture model. Int. J. Solids Struct. 47, 3084–3102.
Madej, L., Wang, J., Perzynski, K., Hodgson, P.D., 2014. Numerical modeling of dual phase
microstructure behavior under deformation conditions on the basis of digital material
representation. Comput. Mater. Sci. 95, 651–662.
Majidi, O., Barlat, F., Lee, M.G., 2016. Effect of slide motion on springback in 2-D draw bending
for AHSS. Int. J. Mater. Form. 9, 313–326.
Maniatty, A.M., Dawson, P.R., Lee, Y. ‐S, 1992. A time integration algorithm for elasto‐
viscoplastic cubic crystals applied to modelling polycrystalline deformation. Int. J. Numer.
Methods Eng. 35, 1565–1588.
Marciniak, Z., Kuczyński, K., 1967. Limit strains in the processes of stretch-forming sheet metal.
Int. J. Mech. Sci. 9, 609–620.
Marvi-Mashhadi, M., Mazinani, M., Rezaee-Bazzaz, A., 2012. FEM modeling of the flow curves
and failure modes of dual phase steels with different martensite volume fractions using actual
microstructure as the representative volume. Comput. Mater. Sci. 65, 197–202.
Mathur, E.B., Dawson, P.R., 1998. Elastoplastic finite element analyses of metal deformations
using polycrystal constitutive models. Comput. Methods Appl. Mech. Eng. 165, 23–41.
Matlock, D.K., Speer, J.G., 2009. Third generation of AHSS: microstructure design concepts,
Microstructure and Texture in Steels.
Matlock, D.K., Speer, J.G., De Moor, E., Gibbs, P.J., 2012. Jestech Recent Developments in
Advanced High Strength Sheet Steels for Automotive Applications: an Overview. Jestech 15,
1–12.
McGinty, R.D., McDowell, D.L., 2006. A semi-implicit integration scheme for rate independent
finite crystal plasticity. Int. J. Plast. 22, 996–1025.
Miehe, C., SchrĂśder, J., 2001. Comparative study of stress update algorithms for rate-independent
and rate-dependent crystal plasticity. Int. J. Numer. Methods Eng. 50, 273–298.
196

Mohammadi, M., Brahme, A.P., Mishra, R.K., Inal, K., 2014. Effects of post-necking hardening
behavior and equivalent stress-strain curves on the accuracy of M-K based forming limit
diagrams. Comput. Mater. Sci. 85, 316–323.
Montheillet, F., Gilormini, P., Jonas, J.J., 1985. Relation between axial stresses and texture
development during torsion testing: A simplified theory. Acta Metall. 33, 705–717.
Nakajima, K., Kikuma, T., Hasuka, K., 1971. Study on the formability of steel sheets. Yawata Tech.
Rep. No. 284, 678–680.
Nanda, T., Singh, V., Singh, V., 2016. Third generation of advanced high-strength steels :
Processing routes and properties. Mater. Des. Appl. 0, 1–30.
Narasimhan, K., Wagoner, R.H., 1991. Finite element modeling simulation of in-plane forming
limit diagrams of sheets containing finite defects. Metall. Trans. A 22, 2655–2665.
Oberg, E., Jones, F.D., Horton, H.L., Ryffel, H.H., 2004. Machinery’s Handbook (27th Edition) &
Guide to Machinery’s Handbook, Industrial Press. Industrial Press.
Panich, S., Suranuntchai, S., Jirathearanat, S., Uthaisangsuk, V., 2016. A hybrid method for
prediction of damage initiation and fracture and its application to forming limit analysis of
advanced high strength steel sheet. Eng. Fract. Mech. 166, 97–127.
Panich, S., Uthaisangsuk, V., Suranuntchai, S., Jirathearanat, S., 2014. Investigation of anisotropic
plastic deformation of advanced high strength steel. Mater. Sci. Eng. A 592, 207–220.
Paul, S.K., 2013a. Real microstructure based micromechanical model to simulate microstructural
level deformation behavior and failure initiation in DP 590 steel. Mater. Des. 44, 397–406.
Paul, S.K., 2013b. Effect of material inhomogeneity on the cyclic plastic deformation behavior at
the microstructural level: micromechanics-based modeling of dual-phase steel. Model. Simul.
Mater. Sci. Eng. 21, 55001.
Paul, S.K., Kumar, A., 2012. Micromechanics based modeling to predict flow behavior and plastic
strain localization of dual phase steels. Comput. Mater. Sci. 63, 66–74.
Peirce, D., Asaro, R.J., Needleman, A., 1983. Material rate dependence and localized deformation
in crystalline solids. Acta Metall. 31, 1951–1976.
Peirce, D., Asaro, R.J., Needleman, A., 1982. An analysis of nonuniform and localized deformation
in ductile single crystals. Acta Metall. 30, 1087–1119.
Qin, J., Chen, R., Wen, X., Lin, Y., Liang, M., Lu, F., 2013. Mechanical behaviour of dual-phase
high-strength steel under high strain rate tensile loading. Mater. Sci. Eng. A 586, 62–70.

197

Raabe, D., Roters, F., 2004. Using texture components in crystal plasticity finite element
simulations. Int. J. Plast. 20, 339–361.
Raabe, D., Wang, Y., Roters, F., 2005. Crystal plasticity simulation study on the influence of texture
on earing in steel. Comput. Mater. Sci. 34, 221–234.
Ramazani, A., Ebrahimi, Z., Prahl, U., 2014. Study the effect of martensite banding on the failure
initiation in dual-phase steel. Comput. Mater. Sci. 87, 241–247.
Ramazani, A., Mukherjee, K., Prahl, U., Bleck, W., 2012. Transformation-induced, geometrically
necessary, dislocation-based flow curve modeling of dual-phase steels: Effect of grain size.
Metall. Mater. Trans. A Phys. Metall. Mater. Sci. 43, 3850–3869.
Ramazani, A., Mukherjee, K., Quade, H., Prahl, U., Bleck, W., 2013a. Correlation between 2D and
3D flow curve modelling of DP steels using a microstructure-based RVE approach. Mater.
Sci. Eng. A 560, 129–139.
Ramazani, A., Schwedt, A., Aretz, A., Prahl, U., Bleck, W., 2013b. Characterization and modelling
of failure initiation in DP steel. Comput. Mater. Sci. 75, 35–44.
Requena, G., Maire, E., Leguen, C., Thuillier, S., 2014. Separation of nucleation and growth of
voids during tensile deformation of a dual phase steel using synchrotron microtomography.
Mater. Sci. Eng. A 589, 242–251.
Rodriguez, R.M., Gutierrez, I., 2003. A unified formulation to predict the tensile curves of steels
with different microstructures. Mater. Sci. Forum 426–432, 4525–4530.
Rossiter, J., Brahme, A., Inal, K., Mishra, R., 2013. Numerical analyses of surface roughness
during bending of FCC single crystals and polycrystals. Int. J. Plast. 46, 82–93.
Rossiter, J., Brahme, a., Simha, M.H., Inal, K., Mishra, R., 2010. A new crystal plasticity scheme
for explicit time integration codes to simulate deformation in 3D microstructures: Effects of
strain path, strain rate and thermal softening on localized deformation in the aluminum alloy
5754 during simple shear. Int. J. Plast. 26, 1702–1725.
Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D.D., Bieler, T.R., Raabe, D., 2010. Overview
of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity
finite-element modeling: Theory, experiments, applications. Acta Mater. 58, 1152–1211.
Roters, F., Eisenlohr, P., Kords, C., Tjahjanto, D.D., Diehl, M., Raabe, D., 2012. DAMASK: The
dĂźsseldorf advanced material simulation kit for studying crystal plasticity using an fe based
or a spectral numerical solver. Procedia IUTAM 3, 3–10.
Safaei, M., 2013. Constitutive modelling of anisotropic sheet metals based on a non-associated
flow rule. PhD disseration, Ghent University.
198

Saleh, A.A., Pereloma, E. V., Clausen, B., Brown, D.W., TomĂŠ, C.N., Gazder, A.A., 2014. Selfconsistent modelling of lattice strains during the in-situ tensile loading of twinning induced
plasticity steel. Mater. Sci. Eng. A 589, 66–75.
SchrĂśder, J., Miehe, C., 1997. Aspects of computational rate-independent crystal plasticity.
Comput. Mater. Sci. 9, 168–176.
Schwindt, C., Schlosser, F., Bertinetti, M.A., Stout, M., Signorelli, J.W., 2015. Experimental and
Visco-Plastic Self-Consistent evaluation of forming limit diagrams for anisotropic sheet
metals: An efficient and robust implementation of the M-K model. Int. J. Plast. 73, 62–99.
Serenelli, M.J., Bertinetti, M.A., Signorelli, J.W., 2011. Study of limit strains for FCC and BCC
sheet metal using polycrystal plasticity. Int. J. Solids Struct. 48, 1109–1119.
Signorelli, J.W., Bertinetti, M.A., Turner, P.A., 2009. Predictions of forming limit diagrams using
a rate-dependent polycrystal self-consistent plasticity model. Int. J. Plast. 25, 1–25.
Simha, C.H.M., Inal, K., Worswick, M.J., 2008. Orientation and Path Dependence of Formability
in the Stress- and the Extended Stress-Based Forming Limit Curves. J. Eng. Mater. Technol.
130, 41009.
Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Interdisciplinary Applied
0000Mathematics. Springer-Verlag, Berlin.
Sirinakorn, T., Wongwises, S., Uthaisangsuk, V., 2014. A study of local deformation and damage
of dual phase steel. Mater. Des. 64, 729–742.
Soare, S., Banabic, D., 2009. A discussion upon the sensitivity of the MK model to input data. Int.
J. Mater. Form. 2, 503–506.
Soare, S., Barlat, F., 2010. Convex polynomial yield functions. J. Mech. Phys. Solids 58, 1804–
1818.
Soare, S., Yoon, J.W., Cazacu, O., 2008. On the use of homogeneous polynomials to develop
anisotropic yield functions with applications to sheet forming. Int. J. Plast. 24, 915–944.
Sodjit, S., Uthaisangsuk, V., 2012. Microstructure based prediction of strain hardening behavior of
dual phase steels. Mater. Des. 41, 370–379.
Sowerby, R., Duncan, J.L., 1971. Failure in sheet metal in biaxial tension. Int. J. Mech. Sci. 13,
217–229.
Speer, J., Matlock, D.K., De Cooman, B.C., Schroth, J.G., 2003. Carbon partitioning into austenite
after martensite transformation. Acta Mater. 51, 2611–2622.

199

Srivastava, A., Bower, A.F., Hector, L.G., Carsley, J.E., Zhang, L., Abu-Farha, F., 2016. A
multiscale approach to modeling formability of dual-phase steels. Model. Simul. Mater. Sci.
Eng. 24, 25011.
Srivastava, A., Ghassemi-Armaki, H., Sung, H., Chen, P., Kumar, S., Bower, A.F., 2015.
Micromechanics of plastic deformation and phase transformation in a three-phase TRIPassisted advanced high strength steel: Experiments and modeling. J. Mech. Phys. Solids 78,
46–69.
Stander, N., Roux, W., Goel, T., Eggleston, T., Craig, K., 2010. LS-OPT ® User’s Manual - A
Design Optimization and Probabilistic Analysis Tool.
SteelRecycleInstitute, 2013. Steel: North America’s Most Recycled Material. Steel North Am.
Most Recycl. Mater.
Stodolsky, F., Vyas, A., Cuenca, R., 1995. Lightweight materials in the light-duty passenger vehicle
market: their market penetration potential and impacts.
Sugimoto, K., Kanda, A., Kikuchi, R., Hashimoto, S., Kashima, T., Ikeda, S., 2002. Ductility and
Formability of Newly Developed High Strength Low Alloy TRIP-aided Sheet Steels with
Annealed Martensite Matrix. ISIJ Int. 42, 910–915.
Sun, X., Choi, K.S., Liu, W.N., Khaleel, M. a., 2009a. Predicting failure modes and ductility of
dual phase steels using plastic strain localization. Int. J. Plast. 25, 1888–1909.
Sun, X., Choi, K.S., Soulami, A., Liu, W.N., Khaleel, M.A., 2009b. On key factors influencing
ductile fractures of dual phase (DP) steels. Mater. Sci. Eng. A 526, 140–149.
Swift, H.W., 1952. Plastic instability under plane stress. J. Mech. Phys. Solids 1, 1–18.
Tarigopula, V., Hopperstad, O.S., Langseth, M., Clausen, A.H., 2008a. Elastic-plastic behaviour
of dual-phase, high-strength steel under strain-path changes. Eur. J. Mech. A/Solids 27, 764–
782.
Tarigopula, V., Hopperstad, O.S., Langseth, M., Clausen, A.H., Hild, F., 2008b. A study of
localisation in dual-phase high-strength steels under dynamic loading using digital image
correlation and FE analysis. Int. J. Solids Struct. 45, 601–619.
Tarigopula, V., Hopperstad, O.S., Langseth, M., Clausen, A.H., Hild, F., Lademo, O.G., Eriksson,
M., 2008c. A study of large plastic deformations in dual phase steel using digital image
correlation and FE analysis. Exp. Mech. 48, 181–196.
Tasan, C.C., Diehl, M., Yan, D., Zambaldi, C., Shanthraj, P., Roters, F., Raabe, D., 2014a.
Integrated experimental-simulation analysis of stress and strain partitioning in multiphase
alloys. Acta Mater. 81, 386–400.
200

Tasan, C.C., Hoefnagels, J.P.M., Diehl, M., Yan, D., Roters, F., Raabe, D., 2014b. Strain
localization and damage in dual phase steels investigated by coupled in-situ deformation
experiments and crystal plasticity simulations. Int. J. Plast. 63, 198–210.
Taylor, G.I., 1938. Plastic strain in metals. Twenty-eighth May Lect. to Inst. Met.
Taylor, G.I., 1934a. The Mechanism of Plastic Deformation of Crystals. Part I. Theoretical. Proc.
R. Soc. A Math. Phys. Eng. Sci. 145, 362–387.
Taylor, G.I., 1934b. The mechanism of plastic deformation of crystals. Part II—comparison with
observations. Proc. R. Soc. A Math. Phys. Eng. Sci. 145, 388–404.
Taylor, G.I., Elam, C.F., 1923. Bakerian Lecture. The Distortion of an Aluminium Crystal during
a Tensile Test. Proc. R. Soc. London. Ser. A 102, 643 LP-667.
Tekkaya, A.E., 2000. State-of-the-art of simulation of sheet metal forming. J. Mater. Process.
Technol. 103, 14–22.
Tikhovskiy, I., Raabe, D., Roters, F., 2008. Simulation of earing of a 17% Cr stainless steel
considering texture gradients. Mater. Sci. Eng. A 488, 482–490.
Tjahjanto, D., 2007. Iso-work-rate weighted-taylor homogenization scheme for multiphase steels,
assisted by transformation-induced plasticity effect. Steel, 1–7.
Tjahjanto, D.D., Eisenlohr, P., Roters, F., 2015. Multiscale deep drawing analysis of dual-phase
steels using grain cluster-based RGC scheme. Model. Simul. Mater. Sci. Eng. 23, 45005.
Tome, C., Canova, G.R., Kocks, U.F., Christodoulou, N., Jonas, J.J., 1984. The relation between
macroscopic and microscopic strain hardening in F.C.C. polycrystals. Acta Metall. 32, 1637–
1653.
TOMOTA, Y., TAMURA, I., 1982. Mechanical Behavior of Steels Consisting of Two Ductile
Phases. Trans. Iron Steel Inst. Japan 22, 665–677.
Toros, S., Polat, A., Ozturk, F., 2012. Formability and springback characterization of TRIP800
advanced high strength steel. Mater. Des. 41, 298–305.
Tresca, H., 1864. Mémoire sur l’écoulement des corps solides soumis à de fortes pressions. C.R.
Acad. Sci. Paris 59, 754.
U.S. Department of Energy, 2015. Quenching and Partitioning Process Development to Replace
Hot Stamping of High-Strength Automotive Steel.
Uthaisangsuk, V., Prahl, U., Bleck, W., 2011. Modelling of damage and failure in multiphase high
strength DP and TRIP steels. Eng. Fract. Mech. 78, 469–486.
201

Uthaisangsuk, V., Prahl, U., Bleck, W., 2009. Characterisation of formability behaviour of
multiphase steels by micromechanical modelling. Int. J. Fract. 157, 55–69.
Vajragupta, N., Uthaisangsuk, V., Schmaling, B., MĂźnstermann, S., Hartmaier, A., Bleck, W., 2012.
A micromechanical damage simulation of dual phase steels using XFEM. Comput. Mater. Sci.
54, 271–279.
Verma, R.K., Biswas, P., Kuwabara, T., Chung, K., 2014. Two stage deformation modeling for DP
780 steel sheet using crystal plasticity. Mater. Sci. Eng. A 604, 98–102.
Volk, W., Hora, P., 2011. New algorithm for a robust user-independent evaluation of beginning
instability for the experimental FLC determination. Int. J. Mater. Form. 4, 339–346.
von Mises, R., 1913. Mechanics of solid bodies in the plastically-deformable state. Math.-phys.
Klasse 4, 1–10.
Wagoner, R., Chan, K., Keeler, S., 1989. Forming Limit Diagrams: Concepts, Methods and
Applications.
Wang, H., Wu, P.D., Boyle, K.P., Neale, K.W., 2011. On crystal plasticity formability analysis for
magnesium alloy sheets. Int. J. Solids Struct. 48, 1000–1010.
Wang, L., Speer, J.G., 2013. Quenching and Partitioning Steel Heat Treatment. Metallogr.
Microstruct. Anal. 2, 268–281.
Watanabe, I., Setoyama, D., Iwata, N., Nakanishi, K., 2010. Characterization of yielding behavior
of polycrystalline metals with single crystal plasticity based on representative characteristic
length. Int. J. Plast. 26, 570–585.
Weiss, M., Kupke, A., Manach, P.Y., Galdos, L., Hodgson, P.D., 2015. On the Bauschinger effect
in dual phase steel at high levels of strain. Mater. Sci. Eng. A 643, 127–136.
Woo, W., Em, V.T.T., Kim, E.-Y., Han, S.H.H., Han, Y.S.S., Choi, S.-H., 2012. Stress–strain
relationship between ferrite and martensite in a dual-phase steel studied by in situ neutron
diffraction and crystal plasticity theories. Acta Mater. 60, 6972–6981.
World Auto Steel, 2011. Future Steel Vehicle.
Wu-rong, W., Chang-wei, H., Zhong-hua, Z., Xi-cheng, W., 2011. The limit drawing ratio and
formability prediction of advanced high strength dual-phase steels. Mater. Des. 32, 3320–
3327.
Wu, P.D., Graf, A., MacEwen, S.R., Lloyd, D.J., Jain, M., Neale, K.W., 2005. On forming limit
stress diagram analysis. Int. J. Solids Struct. 42, 2225–2241.

202

Wu, P.D., MacEwen, S.R., Lloyd, D.J., Neale, K.W., 2004. Effect of cube texture on sheet metal
formability. Mater. Sci. Eng. A 364, 182–187.
Wu, P.D., Neale, K.W., Giessen, E., Jain, M., MacEwen, S.R., Makinde, a., 1998. Crystal
plasticity forming limit diagram analysis of rolled aluminum sheets. Metall. Mater. Trans. A
29, 527–535.
Wu, P.D., Neale, K.W., Van Der Giessen, E., 1997. On Crystal Plasticity FLD Analysis. Proc. R.
Soc. London A 453, 1831–1848.
Xie, C.L., Nakamachi, E., 2002. Investigations of the formability of BCC steel sheets by using
crystalline plasticity finite element analysis. Mater. Des. 23, 59–68.
Yang, M., Dong, X., Zhou, R., Cao, J., 2010. Crystal plasticity-based forming limit prediction for
FCC materials under non-proportional strain-path. Mater. Sci. Eng. A 527, 6607–6613.
Yoshida, K., Ishizaka, T., Kuroda, M., Ikawa, S., 2007. The effects of texture on formability of
aluminum alloy sheets. Acta Mater. 55, 4499–4506.
Yoshida, K., Kuroda, M., 2012. Comparison of bifurcation and imperfection analyses of localized
necking in rate-independent polycrystalline sheets. Int. J. Solids Struct. 49, 2073–2084.
Zamiri, A., Pourboghrat, F., Barlat, F., 2007. An effective computational algorithm for rateindependent crystal plasticity based on a single crystal yield surface with an application to
tube hydroforming. Int. J. Plast. 23, 1126–1147.
Zamiri, A.R., Pourboghrat, F., 2010. A novel yield function for single crystals based on combined
constraints optimization. Int. J. Plast. 26, 731–746.
Zhang, C., Leotoing, L., Zhao, G., Guines, D., Ragneau, E., 2010. A methodology for evaluating
sheet formability combining the tensile test with the M-K model. Mater. Sci. Eng. A 528,
480–485.
Zhang, F., Chen, J., Chen, J., 2014. Effect of through-thickness normal stress on forming limits
under Yld2003 yield criterion and M-K model. Int. J. Mech. Sci. 89, 92–100.
Zhang, L., Lin, J., Min, J., Ye, Y., Kang, L., 2016. Formability Evaluation of Sheet Metals Based
on Global Strain Distribution. J. Mater. Eng. Perform. 25, 1–11.
Zhao, Z., Mao, W., Roters, F., Raabe, D., 2004. A texture optimization study for minimum earing
in aluminium by use of a texture component crystal plasticity finite element method. Acta
Mater. 52, 1003–1012.
Zhou, Y., Neale, K.W., 1995. Predictions of forming limit diagrams using a rate-sensitive crystal
plasticity model. Int. J. Mech. Sci. 37, 1–20.
203