Multiscale crystal plasticity modeling of multiphase advanced high strength steel

Accepted Manuscript

Multiscale Crystal Plasticity Modeling of Multiphase Advanced High Strength Steel Bassam Mohammed , Taejoon Park , Farhang Pourboghrat , Jun Hu , Rasoul Esmaeilpour , Fadi Abu-Farha PII: DOI: Reference:

S0020-7683(17)30211-1 10.1016/j.ijsolstr.2017.05.007 SAS 9569

To appear in:

International Journal of Solids and Structures

Received date: Revised date: Accepted date:

30 November 2016 18 April 2017 6 May 2017

Please cite this article as: Bassam Mohammed , Taejoon Park , Farhang Pourboghrat , Jun Hu , Rasoul Esmaeilpour , Fadi Abu-Farha , Multiscale Crystal Plasticity Modeling of Multiphase Advanced High Strength Steel, International Journal of Solids and Structures (2017), doi: 10.1016/j.ijsolstr.2017.05.007

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

ACCEPTED MANUSCRIPT

Highlights  New multiphase crystal plasticity model embedded in Abaqus explicit for advanced high strength steel. Properties of all phases of steel are considered in the model with two proposed types of phase volume fraction distribution in the steel sheet.



The model is computationally efficient and could be utilized for formability and stamping simulations of multiphase steel.



The model was validated with independent experimental tests.

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT

Multiscale Crystal Plasticity Modeling of Multiphase Advanced High Strength Steel Bassam Mohammed a, b, Taejoon Parkb, Farhang Pourboghratb,c *, Jun Hud, Rasoul Esmaeilpourb, Fadi Abu-Farhad

Department of Mechanical Engineering, Michigan State University, 428 S Shaw Ln, East Lansing, MI 48824, United States b

CR IP T

a

Department of Integrated Systems Engineering, The Ohio State University, 210 Baker Systems, 1971 Neil Avenue, Columbus, OH 43210, United States

c

Department of Mechanical and Aerospace Engineering, The Ohio State University, 201 W 19th Ave, Columbus, OH 43210, United States

Department of Automotive Engineering, Clemson University, Greenville, SC 29607, United States

AN US

d

April 17, 2017

M

* Corresponding author. Tel.: +1 (614) 292 3124; E-mail: [email protected]

Abstract

ED

Designing multiphase metals based on their constituent phase micromechanical properties and using these metals for manufacturing lightweight automotive parts is a challenging

PT

process. A multiphase and multiscale model is strongly beneficial in achieving this goal, as such a model can play an important role in connecting the material response at the

CE

macroscopic scale with the microstructural properties. In the present study, a computationally efficient rate independent crystal plasticity finite element (CPFE) model was used to simulate

AC

the bulging test of a three-phase (ferrite, martensite, and retained austenite) quenched and partitioned Q&P980 advanced high strength steel (AHSS) sheet. The CPFE model was developed to capture the mechanical properties of the steel phases 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 fraction in the steel sheet. In addition to random texture distribution assumption for each grain, two methods were considered in this study for the initial volume fraction distribution of the phases in steel sheet to introduce the effect of inhomogeneity 2

ACCEPTED MANUSCRIPT

based on phase distributions. The objective of developing the proposed model is to investigate the material behavior of the AHSS from the uniaxial microscale tension tests to large deformation applications through validation with the macroscale deformation based on the bulge test. The comparison between the multiphase CPFE model, experiments, and two phenomenological models (isotropic von Mises and anisotropic Hill’48 yield criteria) was conducted for validation of the proposed model. Also, the numerical results were further

CR IP T

evaluated by comparing the predicted against experimental forming limit curve (FLC). The results showed that the trend of the forming limit curve predicted by the multiphase CPFE model is in good agreement with experimental results.

AN US

Keywords: Crystal plasticity, finite element method, multiphase advanced high strength steel, Nakajima bulging test, forming limit curve

1. Introduction

The increased requirement by the federal government in the past few years on light weighting

M

of vehicle structures and improving the crashworthiness of automobiles have prompted the application of advanced high strength steels (AHSS). The development of third generation

ED

advanced high strength steels (3GAHSS) was proposed to achieve the combination of ultrahigh-strength and good ductility properties required for the manufacturing of lightweight

PT

cars with improved fuel economy and passenger safety (Matlock and Speer, 2009).

CE

The proposed quenching and partitioning (Q&P) heat treatment process could result in steel microstructures with retained austenite and controlled amounts of martensite. In general, the

AC

Q&P steels exhibit an ultrahigh strength in the range of 1000–1400 MPa and a uniform elongation in the range of 10–20%. The Q&P steel gains its enhanced strength from martensite laths, and its improved ductility from transformation induced plasticity (TRIP). That transformation of the austenite to martensite for Q&P steel contributes to the overall increase in the plastic strain, and entails a significant strain hardening due to the increase in the volume fraction of martensite phase (Wang and Speer, 2013).

The formability of the cold rolled Q&P980 sheet, which has three phases (ferrite, martensite 3

ACCEPTED MANUSCRIPT

and retained austenite) and a tensile strength of 980 MPa, is superior to most other highstrength steels with the same tensile strength; hence, they are particularly well suited for automotive structural and safety parts such as cross members, longitudinal beams, B-pillar reinforcements, sills, and bumper reinforcements (Wang and Speer, 2013). Understanding such sheet metal responses undergoing many deformation modes is important for industrial metal forming processes, particularly for further improvement of product design as well as

CR IP T

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 AHSS steels. These simple empirical constitutive models do not properly account for the complex microstructure of AHSS sheets, and the evolution of their

AN US

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,

M

applying crystal plasticity finite element method (CPFEM) is highly recommended for the reason that it bridges the gap between the polycrystalline texture and macroscopic mechanical

ED

properties, and it permits for a more profound consideration of metal anisotropy in the stamping process simulation (Kraska et al., 2009). Large scale use of such promising

PT

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

CE

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

AC

methods to coarse-grain (or homogenize) the plastic response from the microstructural scale to a composite response valid at a material point of the continuum are indispensable (Tjahjanto et al., 2007).

In general, crystal plasticity finite element models can be classified as large and small scale depending upon the scale of volume that they represent at a material point. 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 4

ACCEPTED MANUSCRIPT

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).

Large scale 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

CR IP T

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

AN US

model body centered cubic crystals (BCC) such as single phase steels (Mapar, et al., 2016, Kim et al., 2012; Tikhovskiy et al., 2008; Raabe et al., 2005). In addition, based on the fact that localized plastic deformation is strongly influenced by the microstructure (Asaro and Needleman, 1985) the CPFEM has been widely used to predict the localized necking and forming limit curves (FLC) for sheet metals based on a single phase assumption (Lévesque et

M

al., 2015; Mohammadi et al., 2014; Chiba et al., 2013; Serenelli et al., 2011; Wang et al.,

ED

2011; Inal et al., 2005).

In the past few years, industrial demand for applying AHSS to manufacture lightweight parts

PT

has placed new demands on the development of more sophisticated CPFE models for metals with more than one phase. Many of the new developed models are based on small scale

CE

CPFEM of dual phase. In small scale CPFEM simulations, the size of the finite element mesh is equal to or smaller than the grain size. With such sub-grain resolutions, heterogeneity of

AC

deformation within a single crystal is taken into account (Chen et al., 2014). The failure of dual phase steels based on small scale was the dominant approach for most of researchers (Berisha et al., 2015; Tasan et al., 2014; Choi et al., 2013; Kadkhodapour et al., 2011). Using limited number of grains in the representative volume element (RVE), and the high computational cost, makes small scale CPFEM unsuitable for simulating the large deformation of engineering components.

In addition to above studies, there are a few recent publications in which the small scale 5

ACCEPTED MANUSCRIPT

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 and 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

CR IP T

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.

AN US

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., 2016). 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

M

grains in a polycrystalline material under tension test.

Large scale application of CPFEM to sheet metal forming processes with more than one

ED

phase are scarce in the literature. In a recent study by Tjahjanto et al. (2015), the authors investigated the analysis of deep drawing of dual-phase steel (DP800) with only BCC

PT

structure (ferrite and martensite) using previously developed homogenization scheme called

CE

relaxed grain cluster (RGC) and material simulation platform DAMASK (Roters et al., 2012).

The main objective of this study is to introduce a multiscale and multiphase CPFE model to

AC

study the mechanisms of deformation in the AHSS steel. The model idealizes the constituent phases of the steel based on crystal plasticity constitutive relations, and it is based on crystallographic slip which not only accounts for the initial elastic–plastic anisotropy, but also for its gradual change during forming due to the texture evolution. Moreover, the model was implemented into the commercial finite element code (ABAQUS/Explicit) as a user material subroutine, VUMAT, which can distinguish the different phases of the steel based on their slip systems, hardening parameters, phase distribution, and volume fraction for each phase. The capability of the model was investigated based on experimental results, and compared 6

ACCEPTED MANUSCRIPT

with two phenomenological models which in general do not update the texture during forming, i.e. gradual evolution of the anisotropy during the forming process is not accounted for.

The layout of the manuscript is as follows. Section 2 briefly describes the experimental procedures; Section 3 describes the microstructure-based CPFE model; Section 4 explains the

CR IP T

phenomenological models; Section 5 describes the simulation procedure and the assumptions used. Section 6 compares the various models’ predictions with experiments in detail, and discusses their implications. The conclusions of this work are summarized in Section 7.

2. Experiments

AN US

To verify the accuracy of the proposed multiphase CPFE model, numerically generated results were validated against various experiments conducted with Q&P980 sheets. The anisotropic and strain hardening characteristics of the steel were determined using uniaxial tension tests performed at different angles measured with respect to the rolling direction. Experimental results from a hemispherical Nakajima punch test were also compared with

M

CPFEM results to verify the accuracy of the model implemented as an Abaqus/Explicit

2.1 Steel microstructure

ED

VUMAT.

A commercial Q&P980 steel sheet from Baosteel was investigated in the as-received

PT

condition. The chemical compositions of the steel sheet are detailed in Table 1. The crystallographic texture of the steel was analyzed based on electron backscatter diffraction

CE

(EBSD) measurements. Figure 1 (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

AC

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 2 (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 7

ACCEPTED MANUSCRIPT

the sheet and parallel to the rolling plane Figure 2(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.

2.2 Uniaxial tensile test

CR IP T

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. A 3D system digital image correlation (DIC) technique ARAMIS 5M was used to monitor the

AN US

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

M

of 10 mm/min, yielding an average strain rate of ~0.002 s-1, and the DIC cameras were set to

ED

capture images at a fixed rate of 2 fps.

The mechanical properties of the material 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

PT

of five times. The tensile test was conducted at a constant speed for both the yield strength and r-values (width to thickness plastic strain increment ratio in simple tension) determination. The uniaxial tensile

CE

stress–strain curves for the Q&P980 at room temperature with different alignment to the rolling direction are shown in

Figure 3. This figure also shows the Q&P980 sheet’s inherent microstructure inhomogeneity as bands

AC

appearing on its surface detected by DIC along diagonal (DD) and transverse (TD) directions at 4% strain. The measured mechanical properties of Q&P980 by the uniaxial tensile tests are summarized in Table 2. As can be seen in

Figure 3, 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.

8

ACCEPTED MANUSCRIPT

2.3 Hemispherical punch test - Nakajima Test The out-of-plane hemispherical dome stretching test (Nakajima test) is widely used to obtain forming limit curves (FLC) for sheet metals. 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 tests provide important

CR IP T

information about the material behavior that is relevant to the stamping process.

The Nakajima bulging test was conducted using Interlaken Servo-Press 150 formability test system. As shown in Figure 4(a), different blanks with varying width ranging from 20 mm up to 180 mm, and a 1mm gauge thickness were stretched by the hydraulic displacement of the

AN US

Nakajima punch at a constant tool velocity of 1 mm/s. A constant blankholder 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 4 (a) and (b). A double layer of Teflon sheet (0.1 mm sheet thickness) were placed between the blank and the hemispherical punch to prevent premature fracture of the specimen and to

M

limit the effects of friction, particularly at the pole (International Standard ISO 12004-2,

ED

2008).

A 3D DIC system similar to the one used in the uniaxial tensile test was employed to measure

PT

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

CE

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

AC

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 depicts dimensions of the Nakajima test specimens, and the deformed samples after failure.

2.4 High Energy X-ray diffraction (HEXRD) test 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 to 9

ACCEPTED MANUSCRIPT

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). This paper can be referred to for more details on the calibration and determination

CR IP T

of phase properties of the multiphase Q&P steels with in-situ HEXRD.

3. Crystal Plasticity Model 3.1 Crystal plasticity constitutive laws

To model the mechanical behavior of the multiphase Q&P980 steel sheet under complex

AN US

loading, the computationally efficient crystal plasticity model developed by Zamiri and Pourboghrat (2010) were 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)

M

slip systems for ferrite, martensite and the new martensite which is transformed from preexisting austenite were added. The hardening evolution of each phase, which represents the

ED

slip resistance of slip systems were calculated from measured lattice strain distributions obtained from high energy X-ray diffraction (HEXRD). These curves were subsequently

PT

fitted to the extended Voce hardening equation with four parameters in Eq. (11). Furthermore, the three-dimensional CPFE model was used to predict the localized deformation after

CE

necking which results in out of plane deformation in the sheet.

AC

In the proposed model, for computational efficiency and due to the low volume fraction of the newly formed martensite, 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 Xray diffraction (HEXRD).

The formulation of the crystal plasticity model is briefly described below to show the 10

ACCEPTED MANUSCRIPT

fundamental relationship between slip systems, shear rates, and hardening parameters as a function of the rate of plastic deformation.

A velocity gradient in the plastic deformation with respect to the material coordinate system

LP  DP  WP

CR IP T

can be decomposed into a rate of deformation D P and a spin tensor W P as:

(1)

An elasto-plastic 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.

AN US

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):

(2)

ED

M

trail 1 trail 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)  

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,

PT

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

AC

CE

equation for the plastic rate of deformation:

DP  

f (σ, q) σ

(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: 11

ACCEPTED MANUSCRIPT

f (σ, q) 

| σ : P |

 y

1

(4)

The constraints of problem (2) can be combined and replaced by an equivalent single constraint defined as:

where  is the slip system,

N

 N    | σ : P |    ln  exp   1         1 m  y       1

CR IP T

f (σ, q) 

(5)

is the number of slip systems,  y is the critical resolved

ρ is a parameter to determine the closeness to the condition,

AN US

shear stress (CRSS),

| σ : 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:

(6)

M

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

ED

where n is a unit vector representing the normal to the slip plane, and m is a unit vector

N

D P     P

(7)

 1

CE

PT

denoting the slip direction. The plastic deformation gradient can be expressed as:

AC

 where  are the accumulated slip rates. The spin tensor, which represents the material

rotation due to slip, can be expressed as:

N

Ω     ω P

 1

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

(8)

ACCEPTED MANUSCRIPT

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

(9)

Using Eq. (3), (5), and (7), it can be shown that during the plastic deformation of a single

CR IP T

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 )

(10)

AN US

  

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

M

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

ED

be equal to one.

PT

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

AC

CE

(Tome et al., 1984) was adopted



 0       1 

  ()   0   1  1   1  exp   

(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 13

ACCEPTED MANUSCRIPT

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

CR IP T

hardening rate for the linear hardening part of extended Voce-law. The accumulated total shear strain on all slip systems  , can be defined as:

  dt 

(12)

AN US

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

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

M

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

ED

undergoes the same deformation and deformation rate. Therefore, the deformation gradient in Eq. (7) in all grains are the same and equal to the overall deformation gradient.

PT

By assuming that grains have the same volume V c for simplicity, the homogenized stress,

AC

CE

σ , can be calculated from the crystal stress σ c :

σ 

Vc V

σ

c



c

1 M

σ

c

(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. (13) can be rewritten in the form:

 1 σ   fi  i 1  Mi 14

σ c

i c

  

(14)

ACCEPTED MANUSCRIPT

where

fi

is the volume fraction of each phase. Eq. (14) is then used in the explicit

algorithm to calculate the macroscopic stress tensor from the crystal stress tensors for all phases.

For the material texture representation, the three Euler angles that are defined based on

CR IP T

Bunge convention to define the orientation of any crystal with respect to the material coordinate system is used (Bunge, 1982). Based on the Bunge system, the orientation of a crystal can be defined by the three Euler angles 1 ,  , 2 as,

sin 1 sin  2  cos 1 sin  2 cos 

sin  2 sin  

 sin 1 sin  2  cos 1 cos  2 cos 

cos  2 sin  

 cos 1 sin 

cos 

AN US

 cos 1 cos 2  sin 1 sin 2 cos  Q(1 ,  ,  2 )    cos 1 cos  2  sin 1 sin  2 cos    sin 1 sin 

 

3.2 Explicit numerical formulation

M

(15)

For the numerical formulation, the predictor–corrector scheme (Simo and Hughes 1998; Park

ED

and Chung, 2012) was employed for the derivation of the equations and numerical

PT

implementation of the rate-independent crystal plasticity model.

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

CE

ε 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

AC

kn10  n ),

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

(16)

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: 15

ACCEPTED MANUSCRIPT

N    σ kn 10 : P     ln  exp    k 0  1    0   1  m   y ( n )       1

(17)

Otherwise, the (n+1)-th step is elasto-plastic. Then, the accumulated slip rate,  n1 , and the is satisfied. nk1  Tol



N    σ kn 1 : P      ln  exp    k  1       1  m   y ( n )     1

(18)

(19)

M

k n 1

AN US

where,

CR IP T

Cauchy stress, σ n1 , should be iteratively updated until the consistency condition in Eq. (18)

the solution accuracy.

ED

Here, Tol is the tolerance, which is numerically an infinitesimal value chosen to guarantee

PT

The Lagrangian multiplier, nk1 , the accumulated slip rate,  kn1 , and Cauchy stress, σ kn1 ,

CE

are iteratively updated until the consistency condition in Eq. (18) is satisfied. The increment of the Lagrangian multiplier is calculated by

AC

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

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

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

16

(20)

ACCEPTED MANUSCRIPT

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

k

k

k

k

(21)

Here,

k

       n 1 

CR IP T

 N sign(σ : P )    σ : P       P exp  1  k k   y   m     1 y      f (σ )             σ:P N   σ n 1  σ  n 1     m exp  1     y  m    1     

AN US

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

k

  n 1

M

where

k

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

(22)

(23)

k

ED

 σ : P    σ : P        exp  1  k k  2    m    y     f (σ )    y                     σ : P N     y n 1   y  n 1   m exp  m     1      1 y      n 1

CE

PT

(24)

AC

Also, from Eq. (11),

  y   

k

k

                       1 1  exp   0     0  1  1    exp   0     1  1   1  n 1 n1  

The following relationship is derived under the associated flow rule:

17

(25)

ACCEPTED MANUSCRIPT

d   d 

f (σ)  y

(26)

Therefore,

k

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

(27)

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

(28)

and, k

AN US

k

CR IP T

k

which is derived from the following linear elastic constitutive law:

(29)

ED

M

f (σ)   dσ  Ce  dε e  Ce  (dε  dε p )  Ce   dε  d   σ  

PT

Then, the accumulated slip rate is updated by

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

CE

(30)

AC

and the plastic slip rate is updated by

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

where

18

f (σ)  

(31)

ACCEPTED MANUSCRIPT

k

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

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

σ

k n 1

 σ

k n 1

σ

k n 1

(32)

(33)

AN US

σ

k 1 n 1

CR IP T

 sign(σ : P )    σ : P        exp  1    y    m   y f (σ )             σ:P N    m exp    1       m   y  1  n 1   

4. Phenomenological Models

For comparison, phenomenological models based on von Mises and Hill’48 yield functions and isotropic hardening rule were additionally considered in this work. The mechanical

M

properties for von Mises and Hill’48 for the Q&P980 steel have been determined from uniaxial tensile data along 0, 45◦ and 90◦ degrees with respect to the rolling direction. The

ED

uniaxial tension data along the rolling direction was considered as the reference state for

  K ( 0   )n

(34)

CE

PT

characterization of the Swift type isotropic hardening law (Swift, 1952):

where  and  are the effective stress and equivalent plastic strain, respectively, while K,

AC

 0 and n 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 rolling direction. The following material parameters have been obtained in this manner: K = 1413.64 MPa,  0 = 2.1192×10-3

and

n = 0.099.

The quadratic orthogonal yield stress function proposed by Hill (1948) is defined by:

19

ACCEPTED MANUSCRIPT

F ( yy   zz )2  G( zz   xx )2  H ( xx   yy )2  2L yz2  2M  zx2  2 N xy2   2

(35)

where x, y and z are the rolling, transverse and thickness directions of the sheet. Here, F, G, H, L, M and N are anisotropic coefficients, which were obtained using three r-values (r0, r45, and r90):

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

CR IP T

F

(36)

The other coefficients, L and M, were assumed as 1.5 for simplicity. The Hill’48 material constants for Q&P980 are listed in Error! Reference source not found.. The three

AN US

dimensional von Mises and Hill’48 continuum based yield criteria provided by Abaqus/Explicit code were used for the finite element simulations.

5. Numerical Simulation and Procedure Assumptions for Implementation of Multiphase Crystal Plasticity Model

M

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

ED

calibrated based on the weight fraction of each phase in the 1 mm thick Q&P980 sheets used

PT

in the Nakajima test with different widths.

The experimental microstructural observations of the Q&P980 steel sheet under consideration

CE

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

AC

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, since the newly transformed martensite will have low volume fraction compared to many 20

ACCEPTED MANUSCRIPT

TRIP steels. This minor effect was also reported by other researchers (Srivastava et al., 2015).

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), while the following 12 slip systems ({111} <110> slip systems) were considered for the FCC crystal (retained austenite). In general, the martensite phase has a body-centered tetragonal (BCT)

CR IP T

lattice structure, however, an approximation by BCC is possible and has been commonly used in the modeling of plastic deformation of martensite (Tjahjanto et al., 2015).

For the complex microstructure of the Q&P980 steel sheet and based on the EBSD data in Figure 1, a random crystallographic orientation was assumed for grains in an aggregate. The

AN US

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 crystal plasticity finite element simulations is often neither practical nor scientifically rewarding, since it provides more texture information than is usually necessary for predicting plastic anisotropy

M

of formed parts and the evolution of texture (Raabe and Roters, 2004).

ED

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

PT

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.

CE

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 FEA integration points have fixed

AC

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 21

ACCEPTED MANUSCRIPT

all integrations points have equal volume fractions of phases. Figure 6 shows examples of phase distributions in 180 mm and 20 mm Nakajima test blank sizes with 140 grains per integration point. In Figure 6(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

CR IP T

phase. Also, 695,170 grain orientations were assumed for austenite and the new martensite representing a volume fraction of 5% for each of them. As the distribution curve in Figure 6(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

AN US

elements, and by as many as 88 grain orientations in 2 other elements. In contrast, the other distribution curve in Figure 6(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

M

and inhomogeneity in terms of different phase distributions. A similar scheme was adopted for all the Nakajima test blank widths in this study. Figure 6(b) shows an example of different

ED

phase distributions assumed for the 20 mm wide blank.

PT

6. Model Verification and Discussions

Uniaxial and bulging test data for Q&P980 steel were used to calibrate and validate the CPFE

CE

model. The uniaxial tension test data obtained based on in-situ High Energy X-ray Diffraction (HEXRD) was used for the calibration of the hardening parameters for different phases. For

AC

validation purposes, these calibrated parameters were then adopted for the simulation of Nakajima bulging tests with different sheet widths.

6.1 Uniaxial tension The experimental data from the in-situ high energy X-ray diffraction were used for the calibration of the 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 22

ACCEPTED MANUSCRIPT

with 50 mm gauge length. 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 grains (orientations) assumed per integration point for each element.

In order to accurately predict individual phase parameters for the steel, the fit involved using experimental stress-strain data for each phase, as well as the global stress strain curve for the

CR IP T

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 code 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

AN US

parameters for each phase was carried out until the simulation and experimental curves obtained from the in-situ high energy X-ray diffraction matched with each other.

After roughly fitting the hardening parameters for each phase to its corresponding experimental stress-strain curve, another scheme was applied using LS-OPT software to

M

further optimize those parameters. In this approach, we tried to find the most accurate parameters for the hardening equation by only fitting two phases and finding the parameters

ED

for the remaining phases such that it provides the best fit to the macroscopic stress stain curve. In this case, the fitted parameters for austenite and the new martensite were kept similar to

PT

the fitted values in the first attempt since together they constituted only 10% of the steel.

CE

The final hardening parameters which were fitted to each phase using the above calibration method are listed in Table 3 . Error! Reference source not found. shows the comparison

AC

between the experimental stress-strain curve for individual phases and the best fitted curves with 60 grain orientations assumed per individual phase. Error! Reference source not found. 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 Error! Reference source not found.. 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 23

ACCEPTED MANUSCRIPT

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).

6.2. Numerical simulations of the Nakajima hemispherical dome stretching tests Numerical simulations of the Nakajima type hemispherical dome stretching tests were

CR IP T

performed with the calibrated CPFE model for the 1mm thickness Q&P980 steel to evaluate the effect of strain conditions on the forming processes. Figure shows the finite element model of the die, punch, blank and blank holder used in the Nakajima test. In contrast with many CPFEM 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

AN US

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.

M

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

ED

tooling and the deformable blanks. Since 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

PT

frictional parts in the simulation (Oberg et al., 2004). Furthermore, since the strain-rate independent material model was applied in finite element simulations, and in order to

CE

enhance the calculation speed, mass-scaling technique was used by setting the minimum step increment time to 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

AC

User’s Manual, 2000).

6.2.1 Punch force-displacement analysis Using the calibrated material parameters, finite element simulations were performed until the strain localization (necking) occurred. In order to optimize the element size in the Nakajima hemispherical dome stretching test, the effect of the element size on the punch forcedisplacement was investigated. As can be seen from Figure 10, the simulation results did not 24

ACCEPTED MANUSCRIPT

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 (a) -11 (d) show the comparison between punch force-displacement curves obtained from multiphase CPFE model, Hill’48 and von-Mises models with experimental curves for

CR IP T

four different specimen widths. Using four 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.

AN US

In the narrow blank with a 20 mm width in Figure (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 correctly predicting the localized failure point at the punch depth of 25 mm. In the case of 20 grains per gauss point, however, the punch force-displacement curve drops a little earlier compared with 60 and 140

M

grains, which perfectly match each other. The von-Mises and Hill’48 yield criteria show a little higher force than the CPFE model, and an earlier localization of around ~24 mm in the

ED

case of Hill’48, and at ~22 mm in the case of von Mises. The CPFE simulation results are in better agreement with the experimental data partly because the deformation mode is very

PT

close to the uniaxial tension, which was used along the rolling direction as the reference state to calibrate the crystal plasticity material parameters.

CE

In the case of intermediate blank width of 60 mm, Figure (b), the deformation mode is in between the uniaxial tension and plane strain tension. As can be seen, the simulated punch

AC

force-displacement 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 Hill’48 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 (c), is almost similar to the plane strain case. The CPFE results in this case again show very good match with the experimental 25

ACCEPTED MANUSCRIPT

force-displacement curve up to a displacement of about 24 mm, while the curves for the von Mises and Hill’48 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 (d), where the deformation mode is close to the

CR IP T

balanced biaxial tension, the von Mises and Hill’48 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

AN US

about 26 mm.

As mentioned in Section 6.1 and Section 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

M

model agreed very well with experiments in the Nakajima bulging test as shown in Figures 12-17. Also, in all cases, the simulated punch force based on the CPFE model matched well

ED

with experiments, while the phenomenological models over predicted the punch force in the 60 mm and 100 mm wide blanks as shown in Figures 11 (b) and (c), respectively. In general,

PT

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

CE

and Hill’48 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

AC

shown in Figure 20, 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 compared to the non-quadratic yield functions in general.

It can be noted that strain localizations observed from DIC measurements (Figures 12-17) occurred at earlier punch strokes compared with the onset of necking from the punch force versus punch stroke curves in Figure 11. As shown in Figures 12-17, the onset of strain 26

ACCEPTED MANUSCRIPT

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

CR IP T

hardening rule shown in Eq. (11) converges to the linear hardening rule as the exponential 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.

AN US

6.2.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.

M

The red dots shown in DIC images in Figures 12-17 are the fracture location.

Figures 12-14 show the comparison of predicted and measured strain distributions and strain

ED

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

PT

CPFE model with different grain orientations and phase distributions can correctly predict the evolution of strains and their distributions, with increasing applied punch displacement

CE

(Figure ). In contrast, the von-Mises and Hill’48 results show a different strain distribution.

AC

In Figure 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 , at a punch displacement of ~23.5 mm (1 mm away from fracture).

27

ACCEPTED MANUSCRIPT

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 15 ). 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 prediction in terms of strain values and distributions for the 100 mm blank width at a punch displacement of ~23 mm (Figure ). Also, the CPFE simulation based on the

observed in DIC, Figure

CR IP T

random phase distribution successfully captures the significant strain localization which was (e).

It is interesting to note that although the CPFE model predicts a higher force-displacement at the necking point, Figure (d), compared with experimental and phenomenological models, it predicts a more realistic strain distribution compared with Hill’48 model which had a closer

AN US

fit in terms of punch force displacement at the necking point Figure 17. 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 (d) shows that the experimental force-displacement curve for the 180 mm blank is lower than the CPFE model’s prediction.

agreement

with

the

M

Overall, the strain analysis based on the multiphase CPFE model shows a very good experimental

results

(obtained

from

DIC)

compared

with

ED

phenomenological models. The CPFE model was able to accurately predict strain values and

PT

localized necking locations and directions compared with experimental DIC results.

CE

6.2.3 Identification of necking strains for FLC The Nakajima test provides the data necessary to construct a forming limit curve (FLC) for a

AC

sheet metal. As explained in Section 2.3., 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 Figures 1217 show this area.

The FLC for Q&P980 multiphase steel was constructed using DIC, as well as with CPFE 28

ACCEPTED MANUSCRIPT

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 ), the time dependent linear best fitting (LBF) method was employed to determine the major and minor strains at the onset of instability.

CR IP T

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

AN US

a straight fitting line through the area of stable deformation (see Stable Line Fit, Figure ), and the second is a gliding straight line (see Unstable Line Fit, Figure ) 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

forming limit curve (FLC).

M

necking starts and the corresponding major and minor strains are taken to construct the

ED

Figures 11 (a)-(d) show that compared with experimental curves, the necking strains predicted by von Mises and Hill’48 models for various modes of deformation are inconsistent.

PT

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. Figures 12-15 also

CE

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

AC

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 ). 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= 29

ACCEPTED MANUSCRIPT

(𝜀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:

(37)

CR IP T

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

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

AN US

In Figure , 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.

M

The experimental and numerical FLC calculated using the LBF method are plotted in Figure . Forming limit strains numerically obtained by random and fixed phases distribution

ED

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

(a) and (b). However, the forming limit

PT

strains obtained numerically for other specimen widths were higher compared with experimental FLC. Compared with the random phase distribution, the forming limit strains

CE

predicted with the fixed phase distribution were slightly higher.

AC

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 6.2.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. 30

ACCEPTED MANUSCRIPT

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 distributions is very close and in good agreement with the

CR IP T

experimental strain path (Figure ). 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 that the strain path from the simulation with the fixed phase distribution is

AN US

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

(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

M

phases, compared with the stronger martensite.

ED

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

PT

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

CE

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

AC

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. 7. Conclusions A multiscale and multiphase CPFE model which accounts for the grain rotation and the evolution of anisotropy has been presented and applied to large deformation in Nakajima bulging test. A maximum of 140 grains with random orientation distribution had been 31

ACCEPTED MANUSCRIPT

specified at each integration point to represent the anisotropic behavior and texture evolution of a steel sheet during the plastic deformation. Two different methods were considered to represent the inhomogeneous distribution of phases in steel, and it was found that the random phase distribution shows more accurate results than the fixed phase distribution. The CPFE simulation results for various deformation modes, including uniaxial tension, plane strain tension and balanced biaxial tension showed good fit with the experimental data in terms of

CR IP T

punch force-displacement. In addition, the accurate prediction of strains, strain distribution, and strains at the onset of necking, offers a promising method for developing forming limit curve (FLC) for the multiphase steel, based on their constituent phase properties, and their inhomogeneous distribution. Such microstructure-sensitive models will be highly desirable

AN US

for the stamping simulation of the third-generation advance high strength steels (3GAHSS).

Acknowledgements

This material is based upon work supported by the Department of Energy under Cooperative

Partnership LLC (USAMP).

M

Agreement Number DOE DE-EE000597, with United States Automotive Materials

This report was prepared as an account of work sponsored by an agency of the United States

ED

Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or

PT

responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned

CE

rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its

AC

endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. References

ABAQUS User’s Manual, 2000. version 6.14, Hibbit, Karlsson & Sorenson, Providence, RI.

Asaro, R.J., Needleman, A., 1985. Overview no. 42 Texture development and strain hardening in rate dependent polycrystals. Acta Metall. 33, 923–953. doi:10.1016/00016160(85)90188-9 32

ACCEPTED MANUSCRIPT

ASTM Int., 2009. Standard Test Methods for Tension Testing of Metallic Materials 1. Astm 1–27. doi:10.1520/E0008 Berisha, B., Raemy, C., Becker, C., Gorji, M., Hora, P., 2015. Multiscale modeling of failure initiation in a ferritic–pearlitic steel. Acta Mater. 100, 191–201. doi:10.1016/j.actamat.2015.08.035 Bunge, H.-J., 1982. Texture analysis in materials science. Butterworth.

CR IP T

Chen, P., Ghassemi-Armaki, H., Kumar, S., Bower, A., Bhat, S., Sadagopan, S., 2014. Microscale-calibrated modeling of the deformation response of dual-phase steels. Acta Mater. 65, 133–149. doi:10.1016/j.actamat.2013.11.036

AN US

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. doi:10.1016/j.commatsci.2013.04.025 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. doi:10.1016/j.ijplas.2012.11.013

M

Grujicic, M., Batchu, S., 2002. Crystal plasticity analysis of earing in deep-drawn OFHC copper cups. J. Mater. Sci. 37, 753–764. doi:10.1023/A:1013839914584

ED

Guan, Y., Pourboghrat, F., Barlat, F., 2006. Finite element modeling of tube hydroforming of polycrystalline aluminum alloy extrusions. Int. J. Plast. 22, 2366–2393. doi:10.1016/j.ijplas.2006.04.003

PT

Hill, R., 1948. A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proc. R. Soc. Lon. Ser.-A 193, 281–297.

CE

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. doi:10.4028/www.scientific.net/KEM.549.397

AC

Hu, X., Choi, K.S., Sun, X., Ren, Y., Wang, Y., 2016. 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. doi:10.1007/s11661-016-3373-2 Inal, K., Neale, K.W., Aboutajeddine, a., 2005. Forming limit comparisons for FCC and BCC sheets. Int. J. Plast. 21, 1255–1266. doi:10.1016/j.ijplas.2004.08.001 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. 33

ACCEPTED MANUSCRIPT

Kadkhodapour, J., Butz, A., Ziaei-Rad, S., Schmauder, S., 2011. A micro mechanical study on failure initiation of dual phase steels under tension using single crystal plasticity model. Int. J. Plast. 27, 1103–1125. doi:10.1016/j.ijplas.2010.12.001 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. doi:10.1016/j.ijsolstr.2012.02.027

CR IP T

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. doi:10.1016/j.commatsci.2009.03.025 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. doi:10.1016/j.ijplas.2015.10.012

AN US

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. doi:10.1016/j.ijmecsci.2008.11.004

M

Mapar, A., Ghassemi-Armaki, H., Pourboghrat, F., Kumar, K.S., 2017. A differential exponential hardening law for non-Schmid crystal plasticity finite element modeling of ferrite single crystals. Int. J. Plasticity 91 (2017) 268-299.

ED

Matlock, D.K., Speer, J.G., 2009. Third generation of AHSS: microstructure design concepts, Microstructure and Texture in Steels. doi:10.1007/978-1-84882-454-6_11

PT

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. doi:10.1016/j.commatsci.2014.01.017

CE

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. doi:10.1017/CBO9781107415324.004

AC

Park, T., Chung, K., 2012. Non-associated flow rule with symmetric stiffness modulus for isotropic-kinematic hardening and its application for earing in circular cup drawing. Int. J. Solids Struct. 49, 3582–3593. doi:10.1016/j.ijsolstr.2012.02.015 Raabe, D., Roters, F., 2004. Using texture components in crystal plasticity finite element simulations. Int. J. Plast. 20, 339–361. doi:10.1016/S0749-6419(03)00092-5 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. doi:10.1016/j.commatsci.2004.12.072 Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D.D., Bieler, T.R., Raabe, D., 2010. 34

ACCEPTED MANUSCRIPT

Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications. Acta Mater. 58, 1152–1211. doi:10.1016/j.actamat.2009.10.058 Roters, F., Eisenlohr, P., Kords, C., Tjahjanto, D.D., Diehl, M., Raabe, D., 2012. DAMASK: The dusseldorf advanced material simulation kit for studying crystal plasticity using an fe based or a spectral numerical solver. Procedia IUTAM 3, 3–10. doi:10.1016/j.piutam.2012.03.001

CR IP T

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. doi:10.1016/j.ijsolstr.2010.12.013 Simo, J.C., Hughes, T.J.R., 1998. Computational inelasticity. In: Interdisciplinary Applied Mathematics. Springer-Verlag, Berlin.

AN US

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. doi:10.1016/j.jmps.2015.01.014 Stander, N., Roux, W., Goel, T., Eggleston, T., Craig, K., 2010. LS-OPT ® User’s Manual - A Design Optimization and Probabilistic Analysis Tool.

M

Swift, H.W., 1952. Plastic instability under plane stress. J. Mech. Phys. Solids 1, 1–18. doi:10.1016/0022-5096(52)90002-1

PT

ED

Tasan, C.C., Hoefnagels, J.P.M., Diehl, M., Yan, D., Roters, F., Raabe, D., 2014. 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. doi:10.1016/j.ijplas.2014.06.004 Taylor, G. I. 1938 Plastic strain in metals, Journal of the Institute of Metals 62, 307–325.

AC

CE

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. doi:10.1016/j.msea.2007.11.063 Tjahjanto, D., Roters, F., Eisenlohr, P., 2007. Iso-work-rate weighted-taylor homogenization scheme for multiphase steels assisted by transformation-induced plasticity effect. Steel Research International. 78(10-11), 777-783. Tjahjanto, D.D., Eisenlohr, P., Roters, F., 2015. Multiscale deep drawing analysis of dualphase steels using grain cluster-based RGC scheme. Model. Simul. Mater. Sci. Eng. 23, 45005. doi:10.1088/0965-0393/23/4/045005

Tome, C., Canova, G.R., Kocks, U.F., Christodoulou, N., Jonas, J.J., 1984. The relation 35

ACCEPTED MANUSCRIPT

between macroscopic and microscopic strain hardening in F.C.C. polycrystals. Acta Metall. 32, 1637–1653. doi:10.1016/0001-6160(84)90222-0 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. doi:10.1007/s12289-010-1012-9

CR IP T

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. doi:10.1016/j.ijsolstr.2010.12.004 Wang, L., Speer, J.G., 2013. Quenching and Partitioning Steel Heat Treatment. Metallogr. Microstruct. Anal. 2, 268–281. doi:10.1007/s13632-013-0082-8

AN US

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. doi:10.1016/j.ijplas.2006.10.012 Zamiri, A.R., Pourboghrat, F., 2010. A novel yield function for single crystals based on combined constraints optimization. Int. J. Plast. 26, 731–746. doi:10.1016/j.ijplas.2009.10.004

AC

CE

PT

ED

M

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. doi:10.1016/j.actamat.2003.03.001

Table 1 Chemical composition of Q&P980 sheet used in this study 36

ACCEPTED MANUSCRIPT

Alloying Element

Mn

P

Al

Si

C

% Weight

1.8

0.02

0.05

1.5

0.2

CR IP T

Table 2 Mechanical properties of Q&P980 Material Elasticity Constants C11 (GPa)

C12 (GPa)

(GPa) 269.23

115.38

YS (MPa)

76.92

UTS (MPa)

AN US

Direction

C44

RD (0o) DD (45o)

814.5

1019.3

0.8721

778.4

1001.6

0.9465

772.3

1005.8

0.9357

M

TD (90o)

r-value

Table 3 Anisotropic coefficients of the Hill (1948) function Hill’48 Coefficients

G

H

L

M

N

0.53415

0.46584

1.5

1.5

1.4928

ED

F 0.49785

PT

Values

Table 3 Material parameters for the CPFE model composed of four phases in Q&P980 steel

CE

(unit: MPa)

Phase

348

130

1100

117

Martensite

487

55

1994

185

Austenite

600

2

1600

550

New Martensite

496

39

2200

1788

AC

Ferrite

37

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

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

38

ACCEPTED MANUSCRIPT

(b)

AN US

(c)

CR IP T

(a)

Figure 2: 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

AC

CE

PT

ED

M

piece showing the size of the martensite band (dark regions) in the microstructure.

Figure 3: 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).

39

(a)

CR IP T

ACCEPTED MANUSCRIPT

(b)

Figure 4: Schematic view of the hemispherical dome stretching Nakajima test (a) Dimensional characteristics of the specimens used in the Nakajima tests, (b) Principles and

ED

M

AN US

dimensions of the Nakajima test with a hemispherical punch.

AC

CE

PT

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

40

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

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

41

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 7: Comparison of the experimental strain–stress curve with the CPFE model with 60

AC

CE

PT

ED

M

grain orientations assumed for each constituent phase and the Q&P980.

Figure 8: Comparison of the experimental strain–stress curve with the CPFE model in terms of the number of grain orientations assumed per integration point. 42

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

Figure 9: Schematic diagram for Finite element model of a Nakajima test.

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

43

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

(b)

44

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(c)

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

(b)

(c)

(d)

(f)

(g)

ED

(e)

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

Figure 12: Comparison between the numerically simulated and experimental results for rolling direction strain (  y ) of a 20 mm wide sheet at a punch depth of 20 mm. (a) CPFE with

PT

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

AC

CE

distribution (e) von Mises (f) Hill’48 (g) DIC.

46

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

(c)

(d)

AC

CE

PT

ED

M

(a)

(e) (f) (g) Figure 13: Comparison between the numerically simulated and experimental results for the rolling (  y ) and transvers (  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) Hill’48 (g) DIC. 47

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

(c)

(d)

AC

CE

PT

ED

M

(a)

(e) (f) (g) Figure 14: Comparison between the numerically simulated and experimental results for the rolling (  y ) and transvers (  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) Hill’48 (g) DIC. 48

(b)

M

(d)

(c)

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(e)

Figure 15: Comparison between the simulated and experimental results with blank width of

ED

60mm for rolling direction strain (  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

AC

CE

PT

random distribution (c) von Mises (d) Hill’48 (e) DIC.

49

AN US

CR IP T

ACCEPTED MANUSCRIPT

(b)

(c)

AC

CE

PT

ED

M

(a)

(d) (e) Figure 16: Comparison between the simulated and experimental results of the rolling (  y ) and transvers (  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) Hill’48 (e) DIC. 50

(b)

(c)

AC

CE

PT

ED

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(d) (e) Figure 17: Comparison between the simulated and experimental results for the rolling (  y ) and transvers (  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) Hill’48 (e) DIC. 51

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 18: Determination of instable necking strains using the linear best fitting (LBF)

AC

CE

PT

ED

M

method.

(a)

(b)

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

52

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a)

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

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 21: Comparison between experimental and simulated strain paths in balanced biaxial

AC

CE

PT

ED

M

case.

54