3:7:
.. 2..
.g-

 

 

 

 

"i. .I.|.:...}.o.
.3, . L A 4:55.
. K i f .12.... .1. . 55-.

 

 

 

fiwfifiw

. .
%
4‘

mum.
”amt...

:9. 1...

baaaé
Katyammwhsuul...

.! 5‘)-

I‘tlt. 3 3| I
’71.
6!!

. . ..

1.3!. )l.
.f! 3.1..
9;. 29.15.
K :65 3?}

«55.21....

#3135.

$545

3.
2.. x r! a}!
:15; .5: . .. {II-ii!
3.3...“fixwtil3. (.95
... .. . . 3.
1L. .§.d§...hn.m
I. t
2
a
:3:
firix. ...
if...

 

:1:
.{4 ,r...‘ f :‘\..s.
!I. .13?-

3\ 5&-
(I

 

 

 

1:35
7’3 7 LIBRARY
, on Michigan State
_ University

 

 

This is to certify that the
dissertation entitled

COMPUTATIONALLY EFFICIENT CRYSTAL PLASTICITY
MODELS FOR POLYCRYSTALLINE MATERIALS

presented by

Amir Reza Zamiri

has been accepted towards fulfillment
of the requirements for the

PhD degree in Mechanical Erflneering

 

 

(“I 'aWJM

 

Major Professor’s Signature

1 1/05/2008

 

Date

MSU is an Affirmative Action/Equal Opportunity Employer

-.—.-.--a-.-._-n---n-o—u-I--o-¢-c-n-o-n-I-I-l--l-0-1-.---o-I-o-n-o-l-l-0-1-0-0-.---o-I-a-l-o-a-I-I-n-u—¢-I-

i
i PLACE IN RETURN BOX to remove this checkout from your record.
‘i TO AVOID FINES return on or before date due.

MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 KzlProleccsPresICIRClDateDue.indd

COMPUTATIONALLY EFFICIENT CRYSTAL PLASTICITY
MODELS FOR POLYCRYSTALLINE MATERIALS
By

Amir Reza Zamiri

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY

Mechanical Engineering

2008

ABSTRACT

COMPUTATIONALLY EFFICIENT CRYSTAL PLASTICITY
MODELS FOR POLYCRYSTALLINE MATERIALS

By
Amir Reza Zamiri
Crystal plasticity models have been successfully used to take into account the effect
of microstructure in modeling and designing crystalline materials. Up to now, several
computational methods for crystal plasticity models have been proposed. The main
objectives of these computational methods have been to overcome the problem
arising from the non-uniqueness of active slip systems during the plastic deformation
of a crystal and to increase the efficiency and computational speed. The problem with
most current models is that either they are not efficient for all strain paths or are very
expensive. In this thesis, three new computational methods are developed for crystal
plasticity which are believed to be the most efficient and the fastest crystal plasticity
models currently available. In the first model, the current single crystal yield function
is modified to decrease the degree of nonlinearity of the model and increase its
computational speed. The second model, which is a temperature and microstructure
sensitive, is a bridge between computational mechanics and dislocation theory. The
results of this model have been very impressive. The third model, the so called
combined constraints model, is based on developing a new yield function for a single
crystal to overcome the problem arising from the current proposed highly nonlinear

yield functions for a single crystal. These models were implemented into

ABAQUS/Explicit finite element code and two cases of uniaxial tension and tube
hydroforming were simulated using these models. Finally, the combined constraints
model was used to investigate the plasticity induced surface roughness in
polycrystalline high purity niobium. The results of the newly proposed single crystal
model showed a good match with experimental data. Also, compared with other
crystal plasticity models such as the singular value decomposition multi-yield surface
crystal plasticity model and Gambin model, the proposed crystal plasticity models

were more efficient and computationally faster.

ACKNOWLEDGMENTS

I would like to express my deepest appreciation to my advisor, Dr. Farhang
Pourboghrat, for all his supports, guidance, and encouragement during my PhD study.
I would also like to thank the members of my PhD committee; Dr. Ronald Averill,
Dr. Martin Crimp, and Dr. Thomas Pence for their constructive criticism and
suggestions during the course of this work.

I would also like to acknowledge Dr. Frederic Barlat, Dr. Thomas Bieler, and Dr.
Alejandro Diaz for their valuable assistance and support during this work.

All my colleagues; Dr. Nader Abedrabbo, Dr. Yabo Guan, Dr. Michael Zampaloni,
Dr. Hairong Jiang, Dr. Telang, and Dr. Jitendra Prasad are greatly appreciated for
their help and suggestions on this research.

Most importantly, I would like to express my greatest appreciation to my parents for

their encouragement through the years.

iv

TABLE OF CONTENTS

LIST OF TABLES ..... vii

CHAPTER 1

INTRODUCTION ..... ............... . ....... ............ . ....... 1

CHAPTER 2

FRAMWORK FOR CRYSTAL PLASTICITY ...... . .............. 6
2.1. Basic Kinematics and Constitutive Equations .................................................... 6
2.2. Hardening ............................................................................................................. 9
2.3. Crystal Spin and Rotation ................................................................................. 10
2.4. Homogenization ............................................................................................... 104

CHAPTER 3 .............................................................................................. 16

MODIFIED SINGLE CRYSTAL POWER-LAW TYPE YIELD

FUNCTION (MODIFIED GAMBIN’S MODEL) . ....... . ............. . ............... 16
3.1 Basics of Single Crystal Power-Law Type Yield Function ................................ 16
3.2 Modified Single Crystal Power-Law Yield Surface ........................................... 18

4.1. Nonsmooth Multi Yield Surface Theory for Crystal Plasticity ........................ 20
4.2. Dual Mixed Model and its Applications to Micromechanics ............................ 27
4.3. Crystal Plasticity Based on Dual Mixed Model ................................................ 28
4.3.1. Power Type Dual Mixed Crystal Plasticity Model ........................................ 33
4.3.2. Exponential Type Dual Mixed Crystal Plasticity Model ............................... 38
4.3.3. Dual Mixed Crystal Plasticity Model Based on Gilman’s Dislocation
Velocity Model ...................................................................................................... 42
4.3.3.1. Dual mixed crystal plasticity model based on equation (4.56) ............. 44
4.3.3.2. Dual Mixed Crystal Plasticity Model Based on Equation (4.57) .......... 46
CHAPTER 5 .......................................................................... 50
YIELD SURFACE FOR SINGLE CRYSTALS ...................... 50
5.1. Combined Constraints Method in Constrained Optimization ......................... 50
5.2. A Single Crystal Yield Surface Based on Combined Constraints Method ...... 51
5.2.1. A New Model for Dislocation Velocity ........................................................ 53
5.2.2 Evaluation of The Effects of The Closeness Parameter p on The Pr0posed
Single Crystal Yield Function ................................................................................ 54

CHAPTER 6
SOULOTION ALGORITHM FOR PROPSED CRYSTL PLASTICITY

 

 

 

MODELS .................................................................................................... 57
CHAPTER 7 .............................................................................................. 64
EVALUATION OF THE PROPOSED MODELS BY MODELING OF
METAL FORMING PROCESSES ............................................................. 64
7.1. Modeling of The Uniaxial Tension Test 66
7.1.1. Modified Power-Law Yield Surface ............................................................. 68
7.1.2. Dual mixed crystal plasticity models ............................................................ 68
7.1.2.1. Power Type Dual Mixed Crystal Plasticity Model ............................... 69
7.1.2.2. Exponential Type Dual Mixed Crystal Plasticity Model ...................... 71
7.1.2.3. Dual Mixed Crystal Plasticity Models Based on Gilman’s Models ...... 72
7.1.2.4. Combined Constraints Crystal Plasticity Model .................................. 73
7.2. Modeling the Hydroforming Process of Aluminum Tubes .............................. 74
7.2.1. Modeling of The Tube Bulging Into The Round dies .................................... 74
7.2.2. Modeling The Hydroforming of An Aluminum Tube Into A Square Die ...... 86
CHAPTER 8 .............................................................................................. 94
3D MODELING OF THE SURFACE PROPERTIES OF
CRYSTALLIN E MATERIALS .................................................................. 94
8.1. Experimental And Modeling Procedure 96
8.2. Results And Discussion 97
CHAPTER 9
CONCLUSION AND FUTURE WORKS ................................................ 110
9.1. Concluding Remarks ....................................................................................... 110
9.2. Recommendation For Future Works .............................................................. 112
APPENDIX A ....................................................................................................... 115
REFERENCE .......................................................................................... 116

vi

LIST OF TABLES

Table 7.1

Components of m“ and s“ in crystal coordinate for fee. materials ..................... 67
Table 7.2

Dimensions of the Al6061-T4 tube ........................................................... 74
Table 7.3

CPU times for the simulation of tube bulging process using different crystal plasticity
models ............................................................................................ 85
Table 7.4

CPU times for simulation of hydroforming of tube into square die using different
crystal plasticity models ....................................................................... 93

vii

LIST OF FIGURES

Images in this thesis/dissertation are presented in color

Figure 4.1

Elastic domain in the multi surface plasticity model ....................................... 21
Figure 4.2

The schematic of the competition between elastic stored energy and the plastic work
in a plasticity problem ......................................................................... 26
Figure 4.3

Schematic of dual mixed model for a combined constrained problem ................. 28
Figure 4.4

Dislocation velocity versus resolved shear stress for LiF single crystal [Johnston and
Gilman (1959)] ................................................................................. 30
Figure 4.5

Dislocation velocity versus resolved shear stress for different single crystals
[Mechanical Metallurgy, Meyers (1962)] ................................................... 31
Figure 5.1

Effect of p on the shape of the combined constraints single crystal yield surface
(crystal orientation:{ 100}<001>) ............................................................. 56
Figure 7.1

a) Pole figure (111) with 1005 grain orientations representing the undeformed
textruded tube, and b) CODF for the undeformed Al6061-T4 aluminum extruded

tube ............................................................................................... 65
Figure 7.2

A uniaxial tensile model containing 469 84R shell elements used in these

simulations ....................................................................................... 66
Figure 7.3

The stress — strain curves obtained by modified single crystal power yield surface...68
Figure 7.4

Comparison of the stress-Strain curves obtained by simulation using power type dual
mixed crystal plasticity model and experiment ........................................... 70

viii

Figure 7.5
Mises stress distribution in a uniaxial tension model after 15% strain obtained by

power type dual mixed crystal plasticity model ............................................. 70
Figure 7.6

Stress-Strain curves obtained by power type dual mixed crystal plasticity model using
different number of grains ..................................................................... 71

Figure 7.7

Geometry of the die used in modeling of the tube bulging ................................ 76
Figure 7.8

Mesh and the geometry of the one-quarter model of the tube used in modeling of the
tube

bulglng76

Figure 7.9
Bulged aluminum tube .......................................................................... 77

Figure 7.10
Results of the tube bulging simulation using Anand’s singular value decomposition
multi yield surface crystal plasticity model ................................................. 77

Figure 7.11
Results of the tube bulging simulation using modified Gambin crystal plasticity
model ............................................................................................. 78

Figure 7.12
Results of the tube bulging simulation using power type dual mixed crystal plasticity
model .............................................................................................. 78

Figure 7.13
Results of the tube bulging simulation using exponential type dual mixed crystal
plasticity model .................................................................................. 79

Figure 7.14
Results of the tube bulging simulation using dual mixed crystal plasticity based on
Gilman model (equation 100) ................................................................. 79

Figure 7.15
Results of the tube bulging simulation using dual mixed crystal plasticity based on
Gilman model (equation 112) ................................................................. 80

Figure 7.16

Results of the tube bulging simulation using combined constraints crystal plasticity
model ............................................................................................ 80

ix

Figure 7.17
Comparisonof the hoop strain distribution along the length of the tube obtained by

proposed crystal plasticity models ............................................................ 81
Figure 7.18

Mesh and the geometry of the one-quarter of the tube with fine mesh (2970 shell
elements) used in modeling of the tube bulging ............................................ 82
Figure 7.19

Results of the simulation of the bulging of an aluminum tube using the modified
single crystal yield surface crystal plasticity model with a fine mesh (2970 shell

elements). . ................................................................................................................... 83
Figure 7.20

Comparison of the experiment with crystal plasticity simulation of tube bulging
process of Al606l-T4 tube at different forming stages .................................... 84
Figure 7.21

The geometry of one—half of the square die used in the hydroforming simulation of
aluminum tubes .................................................................................. 87
Figure 7.22

The finite element mesh used in the hydroforming simulation of one-quarter of the
aluminum tube ................................................................................... 87
Figure 7.23

The aluminum tube hydroformed into a square die ......................................... 88
Figure 7.24

Failure location in the hydroformed aluminum tube ....................................... 88
Figure 7.25

Hoop strain distribution in the hydroformed aluminum tube predicted by modified
Gambin crystal plasticity model .............................................................. 89
Figure 7.26

Hoop strain distribution predicted by SVD multi surface (Anand’s) crystal plasticity
model .............................................................................................. 89
Figure 7.27

Hoop strain distribution in the hydroformed aluminum tube predicted by combined
constraints crystal plasticity model ........................................................... 90
Figure 7.28

Hoop strain distribution in the hydroformed aluminum tube predicted by power type
dual mixed crystal plasticity model ........................................................... 90

Figure 7.29
Hoop strain distribution in the hydroformed aluminum tube predicted by exponential

type dual mixed crystal plasticity model .................................................... 91
Figure 7.30

Hoop strain distribution in hydroformed aluminum tube predicted by Gilman based
dual mixed crystal plasticity model ........................................................... 91
Figure 7.31

Experimental and numerically predicted hoop strain distribution along the length of
the hydroformed tube ........................................................................... 92
Figure 8.1

The computational model that was used in these investigations; a) the whole model,
b) a portion of the model showing the surface and subsurface grains ................... 98
Figure 8.2

ODFs of a) surface, and b) subsurface textures of high purity niobium ................. 99
Figure 8.3

A portion of the computational model (figurela) showing the orientation of the 4
surface and corresponding subsurface grains. All other parts of the computational
model, except for these 8 grains, are considered to be isotopic .......................... 100

Figure 8.4

a) Surface roughness developed after plastic deformation under uniaxial tension up to
40% strain, b) the transverse plastic deformation in different grains; (111) grain
(black) shows the least amount of transverse deformation while (001) grains (dark
gray in the middle) show the highest ...................................................... 102

Figure 8.5

The distribution of the misorientation in different surface grains, a) the contour plot
(blue: 0° rad and red =0.25° rad ), b) the magnitude of the misorientation across the
grains ............................................................................................. 103

Figure 8.6
A portion of the computational model (figure 1a) showing the surface grains
orientations and subsurface grains orientations. All other parts of the model, except

for these 4 grains, are considered to be isotropic .......................................... 104
Figure 8.7

The distribution of the through-thickness strain across the grains from one side to the
other side ....................................................................................... 105

xi

Figure 8.8

a) The contour plot of the plastic work inside the grains (blue=31 MPa and red=107
MPa), b) The distribution of the plastic work across grains with different orientations
from one side to the other side of the grains ............................................... 107

xii

CHAPTER 1
INTRODUCTION

Crystal plasticity has recently attracted many attentions due to its ability to relate the
plastic behavior of the crystalline materials to their rnicrostructures. Crystal plasticity
can be successfully used in modeling such phenomena as; microcrack initiation, crack
propagation, fatigue, creep in small scale plasticity, texture design, calculation of
damage parameters, evolutionary coefficients of yield functions, fracture criteria and
FLD diagram, etc in crystalline materials. The fundamental importance of single
crystal modeling is reflected by the abundance of theories on monocrystalline
plasticity in the literature. Such theories began to appear in the literature in the early
20th century and are presented by many investigators, including: Taylor, Schmid,
Hill, Rice, Hutchinson, Asaro, Havner, Bassani, Mandel and Kocks. Taylor first
showed that five independent slip systems are needed to accommodate an arbitrary
strain increment imposed on a crystal (Taylor, 1938). The Schmid law states that in
response to applied loads on a single crystal, slip occurs when the resolved shear
stress on a crystallographic slip system (uniquely defined by a slip direction and a slip
plane normal) exceeds a critical value that is a material property (Schmid and Boas,
1935). The continuum mechanics framework provided the necessary foundations to
develop a systematic computational methodology for performing the first detailed
numerical simulations of single crystals including their localization behavior (Peirce
et al., 1982). Early simulations were 2D plane strain based on the Asaro’s (1979)

planar double slip model and provided insights into nonuniform and localized

deformation in single crystals. They were further extended to include modeling of
polycrystals and other problems of interest. Due to an increase in computational
resources that became available, full-scale three dimensional simulations of
crystalline solids undergoing multislip (fcc, bcc) became increasingly feasible (e.g.,
Cuitino and Ortiz, 1992; Kalidindi et al., 1992). However, the complexity of the
structure of governing equations and large number of internal variables in crystal
plasticity models generally makes the 3D numerical simulations intensive and time
consuming.

In order to describe the plastic deformation by the crystallographic glide, three
successive steps must be accomplished; determination of the active slip systems,
determination of the increments of shear on the active slip systems, and overcoming
the problem arising from non-uniqueness of active slip systems due to
interdependency of the slip systems in an arbitrary strain path. Several algorithms
have been proposed to overcome the numerical problems due to interdependency of
slip systems in plastic deformation of crystalline materials. One approach has been
rate-dependent formulations without an elastic domain based on power-type creep
laws without differentiation of slip systems into active and inactive sets, see Pierce et
a1. (1982, 1983), Asaro and Needleman (1985), Mathur and Dawson (1989) and
Huang (1991). Some investigations show that algorithms based on multisurface yield
functions can be the most effective techniques in crystal plasticity, however, in the
rate—independent model under a particular mode of hardening, the constraints of the
multisurface framework are redundant. Cuitino and Ortiz (1992) used a multisurface

viscoplastic model with elastic domain to overcome this problem. Anand and Kothari

(1996) overcame this redundancy by using the singular value decomposition
technique in the inversion of the Jacobian of the active slip systems. Miehe and
Schroder (1999, 2001) suggested two methods of alternative general inverse where a
reduced space is obtained by dropping columns of the local Jacobian associated with
zero diagonal elements within a standard factorization procedure and diagonal shift of
the Jacobian of the active yield criterion functions to overcome this problem.
However, these procedures are very expensive and deficient at times. Raphanel et a1.
(2004) proposed a crystal plasticity model based on Runge-Kutta algorithms.
McGinty and McDowell (2006) introduced a semi-implicit integration scheme for
rate independent crystal plasticity. Nemat-Nasser et al. (1996), Knoekaert et al.
(2000), and McGinty and McDowell (2006) used a formulation based on the
consistency condition during plastic flow to improve the speed and stability relative
to the multiplicative decomposition formulation. The behavior, accuracy, and the
efficiency of some of the above crystal plasticity models have been examined by
several investigators; Busso and Cailletaud (2005), for example, evaluated the set of
active slip systems predicted by several models. Ling et aL (2005) compared the
stress-strain response and the speed of several models using several explicit and
implicit formulations. They also examined the overall running times for single-crystal
updates. Rousselier and Leclercq (2006) compare the overall running times of several
polycrystal models. Kuchnicki et a1. (2006) investigated the efficiency of the Cuitino
and Ortiz (1993) model both in implicit and explicit update procedures.

Montheillet et al. (1985), Van Houtte (1987), Lequeu (1987), Toth et al. (1991),

Arminjon (1991), Darrieulat and Piot (1996), and Gambin (1991, 1992, 1997)

proposed a rate-independent approach based on the concept of crystals with smooth
yield surfaces. In this formulation only one yield function is considered to calculate
the crystal spin and incremental shear strains on active slip systems. Therefore, the
slip system ambiguity problem in crystal plasticity is avoided. Gambin correlated the
degree of the non-linearity of the yield function (i.e., its exponential power) to the
stacking fault energy (SFE) of the material. In this model, the shape of a single
crystal’s yield surface, i.e. with sharp or rounded corners, can be correlated with high
and low SFE, respectively. Based on this formulation, for a material with high SFE,
such as aluminum, a highly nonlinear yield function should be used in order to
reasonably predict the texture of the material after deformation. This, however, makes
the computation to be very expensive. Using a lower degree of non-linearity for the
yield function (i.e., lower exponential power) produces incompatibility in different
grain spins in polycrystalline materials, which leads to an oscillation in the calculated
stresses.

The objective of this work is to find easier and more efficient approaches to solve a
crystal plasticity problem. Three distinct approaches have been followed in order to
develop effective numerical methods to address crystal plasticity problems. In the
first method, the current yield function for a single crystal is modified such that the
experimental data is captured accurately while increasing the computational speed. In
the second method, physical models expressing the relationship between the resolved
shear stress and dislocation velocity, the so called dislocation dynamics models, are
used to solve the crystal plasticity problem. Finally, in the third method a new yield

function is introduced in order to model the crystal deformation.

The fundamentals of crystal plasticity and governing equations are briefly introduced
in chapter 2. The mathematical details of the three methods are discussed in chapters
3, 4, and 5. In chapter 6, a solution algorithm for the proposed models will be
presented. The evaluation of the proposed crystal plasticity models will be discussed
in chapter 7, in which some preliminary examples of simulating the plastic
deformation of a fcc polycrystal subjected to uniaxial tension and tube hydroforming
will be presented. In chapter 8, the combined constraints crystal plasticity model are
used to model the plasticity induced surface roughness in high purity superconducting

niobium as an example of application of the models for bcc materials.

CHAPTER 2

FRAMWORK FOR CRYSTAL PLASTICITY

2.1. Basic Kinematics and Constitutive Equations

Two different coordinate systems will be considered; the material co-rotational (i.e.,
attached to a finite element) and the crystal coordinate system (i.e., defining the
crystal axes of a single crystal). A general formulation can be developed using the
material co-rotational system. Let m be a unit vector normal to the slip plane and s a

unit vector in the slip direction in the crystal coordinate system. Then, any slip system

at the beginning of the strain increment can be defined by the matrix lg as:
a _ a a
10 — m (’9 S (2.1)

The deformation in a single crystal can be decomposed into two parts; the elastic
deformation, which is responsible for the lattice rotation, and the plastic deformation,
which is caused solely by crystalline slip. Therefore, the total velocity gradient in

current state can be expressed as:

L=L*+L”=D+Q (2.2)

where the symmetric rate of deformation D and the antisymmetric spin tensor Q may

be decomposed into elastic and plastic part as follows:

D = D* + DP (2.3)
Q = 9* + 9” (2.4)

*
The lattice spin rate Q is defined by:
=I= ° * *T
Q = R R (2.5)

a:
where R is the rigid body rotation of the lattice.

Now with respect to the material co-rotational coordinate system:

I“ = R*IgR*T (2.6)

a
The symmetric and antisymmetric parts of the slip system are then expressed by P

a
and W , where:

a 1 a a
P =-2-(I +1 T) (2.6)

a 1 a 0!
w :5” ‘1 T) (2.7)

The plastic rate of deformation DP and spin rate up of a single crystal are

(Huang Y. (1991)):

a=1

(2.8)

szfif-w“

a=1

(2.9)

where N is the number of slip systems, and 7a is shear slip rate. D p corresponds to

the rate of deformation of a single crystal and up to the spin of a single crystal
defined in the co-rotational system.

Equation (2.4) shows two different rotations during a deformation increment; one is

the rotation of the material co-rotational frame, which is due to macroscopic spin 9 ,

and the other one is the rotation of the crystal coordinate system, which is due to the

*
microsc0pic spin, 9 . Therefore, the macroscopic and microscopic co-rotational

V v *
stress rate, 0' , and 6 respectively are defined as:

V
0' = ('5' — $26 + 69 (2.10)
V* * *

0' =O-flo+ofl , (2.11)

where G is the stress rate in the co-rotational coordinate system The co-rotational
stress rate on the material axes can be related to the co-rotational axes on the crystal

coordinate by:

§=§*_(o_o*)o+o(o_o*) (,1,

For metallic materials the elastic deformation is usually very small. Therefore, the

relation between the symmetric rate of stretching of the lattice and the Jaumann rate

3|:

V
of the Cauchy stress, 6 , is given by:

V*
o = Ce :D (2.13)

e
where C is the fourth order tensor of anisotropic elastic moduli.

2.2 Hardening

For a crystal plasticity problem one needs to know the evolution of the 7:, (critical
shear stress of slip systems) with plastic shear. Let us consider 3 a = f as the
current slip resistance on slip systema, which evolves from the critical resolved

or
shear stress 70 with plastic slip on the active slip systems. The evolution of g“ is

expressed by the following hardening equation [Asaro and Needleman (1985)]:
N l3
o a _ a o
g " flZh I7fl| (2.14)
=1

where 7'6 is the plastic slip rate on the active slip system ,6 and hat? are denoted

era
as the components of the hardening matrix. The h are known as the self-

hardening moduli while km? for (a it )6 ) are known as the latent-hardening moduli.

The hardening matrix plays an important role in the hardening model. Several
equations have been proposed for the hardening matrix. Hutchinson (1976), Chang

and Asaro (1981), Bassani and Wu (1991) and Peirce et al. (1982) proposed the

following simple law:

hafl = hfl[q + (1" q)5“fl] (no summation onfl) (2-15)

Here q is the so-called latent-hardening ratio which is the ratio of the latent-hardening
rate to the self-hardening rate of a slip system with the values in the range of 1< q

<1.4. It can be considered 1 for coplanar slip systems and 1.4 for non—Coplanar slip

fl. . . . . .
systems. h rs an evolutionary function denoting the self-hardening rate, Wthh can

be expressed as a function of either shear slips or resolved shear stress on slip

systems. Based on Kothari and Anand (1996), [1'6 was considered to evolve as:

.6 a [3
[1'6 :h 1_.£_ 'S 1_£g_
o 1'. gn 7. (2.16)

Where ho , a and 73 are slip system hardening parameters, which are considered to

be identical for all slip systems. ho denotes the initial hardening rate, 73 the

saturation value of the slip resistance, and a the exponent describing the shape of the

function (equation 2.16).

2.3 Crystal Spin and Rotation

The initial orientation of a grain is given by the knowledge of the mapping from the
laboratory frame to the crystallographic coordinate of the crystal. This mapping is a
rotation, which defines the initial orientation matrix Q of the crystal. The initial
orientation matrix may be constructed using Euler angles. In this case, the columns
and rows of the matrix may be interpreted as the crystallographic axes in the
symmetric and Bunge conventions for the Euler angles, respectively. One of the most

important parts of a crystal plasticity algorithm is to correctly calculate the grain

10

rotation due to an imposed plastic deformation. In the material co-rotational

coordinate system, the relative spin tensor is defined as:
' T
9 = W — RR (2.17)
*
where W is the spin tensor, and Q can be decomposed into the elastic part 9

and the plastic part up as:

Q=Q*+Q” (2.18)

And in the material coordinate system:

N
12:51 +QP=RRT+le“-7a (2.19)

a
It must be mentioned that W is in the material coordinate system.

*
Then the lattice spin 9 can be obtained from equation (2.19) as:
N
3|: _ _ a . 0
Q ‘ Q 2 w 7a (2.20)
a=1

The rigid body rotation matrix in an incremental form can be expressed as:
Rt0+At = R10 + Rt0+AtAt

_ ' T

_ [I + RWNR,0 At]R,0

= [I + O,+A,Ar]R,O

(2.21)

Based on the above equation, one can easily write a similar equation for updating the

grain orientation:

11

Qt0+At = Qto + AQ

= [I + Q’AflQto
(2.22)

= 1+[Q—fiwa-7")At Q,0

where Q is an orientation matrix defined in Appendix A.
The crystal orientation matrix updated by equation (2.22) may not have the property
of the rotation tensor. Sirno and Vu-Quoc (1986) suggested the following equation to

update the rotation matrix:

N
Qmw = Exp [9 " 2w... ' 7a)“ Q20 (2.23)

a=1

The exponential part of the above equation is an orthogonal tensor that can be

determined by the Rodrigues formula as:

Exp(Q*At) = I +

 

 

sin(QeAt) Q“ + 1— cos(S2"’At) (9.)2
96 (9e )2

met = J(O* :Q*)/2.

(2.24)

 

However, equation (2.23) seems to introduce some stability problems in the
integration algorithm especially in larger increment times.

Here a grain orientation updating procedure is used based on Pourboghrat and Barlat

*
(2002). Consider equation (2.18), then matrix Q can be considered as a matrix

made of three column vectors as:

12

S2 = [1‘1 1‘2 F3] (2.25)
Then the rotation axis r can be defined as:

I"2 ><I‘3 I“3 ><I‘l I‘l ><I‘2
r:—— r:———— r:—————
|T2 xr,| ’0‘ |F3 xI‘1|’°‘ |I‘1xI‘2|

The rotation angle is:

= sgn(r.Qi xAQ,)|AQ,.|

t8 5’ i
1 _ (roQi )2 (2.26)

where Q,- ’s are the column vectors in crystal orientation matrix Q at the beginning of

the increment and AQ) = Q QiAt .

Now it can be shown that:

Q: = COS (Q, + (r.Q,)(1— cos {I )r + sin {ir x Q, =
P1101 + P1202 + Pi3Qg 227)

I
Where Q,- s are the column vectors of Q after rotation at the end of the increment and

transformation matrix P is defined as:

I r12(1—cos{l)+cos{,’l rir2(l—cos{2)+r3sin§,"2 r1r3(1—cos{3)-rzsin§’3-
P= rzri(1—cos[1)-r38in§1 52(1—00352)+COS;2 66(1—00553)+’15in;3
_r3r1(1—cos{1)+rzsin{1 r3r2(1—cos{,’2)—rlsin;2 132(1-cosf3)+cosf3

(2.29)

 

 

l3

2.4 Homogenization

In this numerical calculation, the Taylor type homogenization model (Taylor (1983))
is used. According to the Taylor model, which is an upper bound model, each grain of
the polycrystal undergoes the same deformation and deformation rate. Therefore, the
velocity gradients in all grains are the same and equal to the macroscopic velocity
gradient.

A velocity gradient in the material coordinate system can be decomposed into a rate

of deformation and spin tensor as:

L = D + W (230)
Then according to Taylor’s model (Chin and Mammel (1969)), the spin tensor W

would also be the same in all grains and equal to the macroscopicW. It must be

* P
mentioned that Q and 9 would be different in different grains.
According to the theory of the grain averages in polycrystalline materials, the volume

average of a local tensorial quantity t(x) is defined as the volume average of each

tensor component and will be denoted by (t) :

1
(t) = V! t(x)dV (2.31)

Now the average stress tensor in a polycrystalline material can be expressed as a

summation of crystals c:

1 c
(6):? l o (x)dV (2.32)

14

C
Assuming that‘, (x) ’s are constant over the grain domains and grains have the

same volume, equation (2.32) can be rewritten as:
c l
_ c = 6
<6) " V 26:“ M 26:“ (2.33)

Where V6 is the grain volume, “C the grain stress tensor, and M the number of
grains.
Equation (2.33) is used in implicit and explicit algorithms to calculate the average

stress tensor from the crystal stress tensors.

15

CHAPTER 3

MODIFIED SINGLE CRYSTAL POWER-LAW TYPE YIELD
FUNCTION (MODIFIED GAMBIN’S MODEL)

3.1 Basics of Single Crystal Power-Law Type Yield Function

The idea of using a single crystal yield surface with round corners is a powerful
approach in overcoming the slip system ambiguity problem in a rate-independent
crystal plasticity model. In this formulation only one yield function is used in order to
calculate the crystal spin and incremental shear strains on active slip systems. The
idea first was introduced by Montheillet et al. (1985) and then was incrementally
improved by Van Houtte (1987), Lequeu (1987), Toth et al. (1991), Arminjon (1991),
Darrieulat and Piot (1996), and Gambin (1991, 1992, 1997). Based on this expression
of yield function the shear strain rate on any active slip system can be expressed as:

Sgnifa Oa)|7 2n—l
ya: 2' '7: (3.1)

 

 

a
where T is the resolved shear stress on the slip system which can be related to the

stress tensor in the fixed coordinate system by:
a A , a
T = O' . P (32)

where 0 is the rotated stress tensor back to the fixed coordinate system.
xlis a positive parameter which depends on the type of dislocation barriers and

defined as:

16

_ a a
.1 - ab meo (3.3)
where a is a constant, b is the Berger’s vector, Pm is the average dislocation density

on the slip system, and mg is the frequency of the passing barriers by dislocations.

a
To is the critical resolved shear stress on the active slip system.

Assuming that parameter 2. is the same for all slip systems, then by substituting for

shear strain rate from equation (3.1) into equations (3.2) and (3.3):

 

 

 

N a aZn—l
P _ Sgn(7 )|6:P . a
D "g r: I r: I P ‘3‘”
N sgn(t"”)|o°P“2n—1
smz . ° .

 

a=1 1'o 2'o

I a I -w (3.5)
Considering equations (3.4) and (3.5) and remembering the concept of the normality

rule, one can introduce a plastic potential with respect to crystal coordinate system as:

n

a 1 N o:P"2
¢(“’Tc’¢):§; £707— _1 (3.6)

C

where (D is a vector containing the three Euler angles defining the orientation of the

crystal with respect to the material co-rotational coordinate system at any deformation

a
step and it arises from the dependency of P on Euler angles.

17

3.2 Modified Single Crystal Power-Law Yield Surface
In contrast to Gambin’s yield function, which is defined in the fixed coordinate
system, we define the above crystal yield function in the crystal coordinate system.

Therefore, the plastic potential in equation (3.6) is independent of Euler angles, or in

other words Pa remains constant. These modifications reduce the computational
time and also allow using a value for n (the degree of yield function non-linearity)
that is half of that used in the Gambin’s original model without losing the stability in
computed stresses and grain’s spin rate. Based on this formulation, 3 crystal
orientation tensor, Q, is defined, which defines the orientation of the crystal

coordinate system with respect to the material coordinate system Stress tensor in

a]:
crystal frame, 0' , can be computed using equation (2.13) and stress tensor in

material coordinate system can be obtained by the following transformation,

o = QTo*Q (3.7)

A yield function with respect to the crystal coordinate system is defined as:

 

1 N o*°Pa 2’; \
¢(6’Tf):2n* g—F—O— —;| (3.8)

where g" is suggested to be:

1 N N a 5
g": WZZIP" :POI (3.9)

a=1 ,B=1

18

The bracket in the above equation emphasizes that {I should be an integer number. k

is a constant which depends on the rrmterial. For aluminum, k can be considered as 4
or 6 depending on the alloy. Gambin suggested that in order to obtain a reasonable
prediction for texture, the following approximate relation must be considered between

the power of the yield function, n, and stacking fault energy (SFE):
I"
n E — X 113—3 3 10
G b ( - )

where F is a value of SFE, G is the shear modulus, and b is the burger’s vector.

*
Based on this formulation 71 can be found as:

_ ~ —3
‘—" _" = _"""""" X10 (311)

19

CHAPTER 4
DUAL MIXED MODELS

4.1 Nonsmooth Multi Yield Surface Theory for Crystal Plasticity

According to the nonsmooth multi yield surface theory (Naghdi(1960), Koiter(1953),
Nemaat Naser (1983), Sirno and Hughes (1998)) the yield surface is composed of
several smooth yield surfaces which intersect nonsmoothly. This theory was initially
developed for soil and rock mechanics such as; Cam-Clay model, the cap model, and
Mohr-Coulomb model. For more information about these models one should look at
Loret and Prevost (1986), DiMaggio and Sandler (1971), Sandler and Rubin (1979),
Resende and Martin (1986). In plasticity, the classical Tresca yield function is an
example of nonsmooth multi yield surface.

Before going further, let us define two very important vectors, qanda, which are

frequently used in the theory of plasticity. q is a matrix containing state variables
such as back stress and the yield stress and a is a vector containing internal variables
such as effective plastic strain. For the simple case of rate independent plasticity with
isotropic hardening at room temperature, q would only have one variable, which is

the flow stress, 6" , and a would have one variable, which is the effective plastic

strain, 8 p . Now with these definitions let us consider the problem of crystal
plasticity as a nonsmooth multi yield surface problem. In a single crystal, depending
on the crystal structure, there are several slip systems that can be activated separately.
Here we assume that the deformation occurs near room temperature, and that when

the resolved shear stress on a slip plane exceeds a critical value plastic deformation

20

happens. Based on this definition a nonsmooth multi yield surface for a single crystal

can be defined mathematically as (Figure 4.1):
8E0 ={(o,q)l fa(o,q) =0 for some ae [1,2,...,N]}

(4.1)

and the elastic domain is defined as:

Ed ={(o,q)1 fa(o,q) s 0 for all are [1, 2,...,N]}

(4.2)

where N is the number of slip systems.

 

a — Constraint

Figure 4.1 Elastic domain in the multi surface plasticity model

21

With these definitions, the rate of plastic deformation is obtained by:

N

DP = Zzaaafa(o,q) (43)

(2:1

where aafa (Ga ‘1) represents differentiation of the yield surface with respect to

stress, and la are Lagrange multipliers that are called the plastic consistency
parameters in plasticity.

Now with the above definition, let us write the model in a computational format. Any
computational problem can be expressed as an optimization problem. The most
important part of an optimization problem is to define (an) appropriate objective
function(s). Elastic strain energy has been known to be a good objective function for
mechanical analysis. For the case of plasticity, a good objective function can be
considered to be the simultaneous minimization of the elastic stored energy and the
maximization of the release of energy (dissipation) due to plastic work. In other
words, an equilibrium must be achieved between the increase in the internal energy
due to elastic stored energy of the material and the release of the energy due to the
plastic work (Figure 4.2). This is known as the principle of the maximum plastic
dissipation (see Mandel (1964) and Lubliner (1984, 1986)). The elastic stored energy

W can be expressed as:

e 1
W =-2—(e—e”):C:(e—a”) (4.4)

Based on above equation the elastic stiffness matrix, C is obtained by:

22

_ 32W (8 — a”)
— 8(8—8” )2 (4'5)

 

There is also the following relationship between the elastic stress and strain:

6=VW(8-£p)=CZ(8-£p) (4.6)

One can develop a similar function as a measure of the plastic work or dissipation,

which is the so called hardening potential, as:
1 o
Il=-2—D.(a®a) (4.7)

where D is the so called matrix of generalized plastic moduli and is defined by:

: 8211(11)

D a (12 (4.8)

State variables are related to the internal variables by:
q = -VII(a) = -Da (4.9)

As explained before, for the simple case of rate independent plasticity with isotropic

hardening at room temperature q = —0—' and a = 5p and one may define the
_. _p n )1
strain hardening as 0 = k (8 ) or q = —k(a) . Then, D is a matrix with a

—1
single variable as D = kart .
In a plasticity problem, the aim is to find the stress tensor and the state variable
vector. In other words, the objective function must be defined in such a way that the

design variables are stresses. Therefore, for any of the above functions, it is possible

23

to define a complementary function such as the complementary elastic stored energy

potential, 1(a) , and the complementary hardening potential, A(q) , associated with

W(8 " 8p) and 11(0) , respectively. For this case the following relations can be

expressed:

_ 822(0)
C I z—aai" (4.10)

3. = 31(0)
as (4.11)

where 88 = 3 “'8p and

_ 82A(q)
D 1 = ‘3‘]?— (4.12)

= _3_A_(_€ll
aq (4.13)

With the above definitions, in a plasticity problem, the objective is to minimize the
elastic stored potential and at the same time maximize the plastic dissipation. Or, to
minimize the difference between the elastic stored energy and the plastic dissipation

as:

24

(6,0) = 214... [2(21— 1101)]
,qua

1 1
M' —2‘.:C'1:Z‘.—— ~01
2,qelEza[2 2Q ]

(4.14)

l 1
=M° —2:C’1:Z+— :D‘I:
M3205 2‘] (1]

Where E0 is defined in equation (4.2). One should notice that the matrix D is not

constant and is usually a function of (I.

Now the problem of crystal plasticity can be easily defined. For a given increment
and a given stress tensor, 2‘. , at the start of the increment (the design variable), and
the internal variable vector, q, the crystal plasticity problem is the argument (solution)
of the following constrained minimization problem. In the following, the first part of
the objective function expresses the elastic stored energy and the second part shows

the dissipation of the energy due to plastic deformation.

r 1

Min[-;-(o,’fgl — 2) : C‘1 : (of: — 2) +
1

(GnH’an) = ARG< 5(qn —q):D_1:(qn Tq)] t

. (4.15)
Subject to:

 

 

fa(2, q) S O for all as J

25

Incremental internal energy

 

V

 

s Incremental plastic work

 

Figure 4.2 The schematic of the competition between elastic stored energy and the plastic work
in a plasticity problem

As explained before, the components of q vector define the back stress for kinematic

hardening or the effective stress for the hardening rule. C is the anisotropic elastic

trial
stiffness matrix and D is the matrix of generalized plastic moduli. “n+1 is the trial

stress tensor obtained by assuming that all the incremental deformation in a given

increment is initially elastic. For a crystal plasticity problem with isotropic hardening,

fa(69 (l) is defined as:

o : P“
fa(034)=—| 0 L1 (416)
T, (q) '

A comparison of equations (2.8), (4.3), and (4.16) shows that for the case of crystal

a o
plasticity, A are the same as plastic shear deformations on slip system, 7“.

26

4.2 Dual Mixed Model and its Applications to Micromechanics

Due to the interdependency of shear slips 7a, solution to the equation (4.15) for a

crystal plasticity problem is usually singular or ill-conditioned, especially for the rate
independent case. Several algorithms have been used to overcome this singularity,
such as the singular value decomposition method, methods of generalized inverse of
the local Jacobian of the active yield criterion functions in reduced space, and the
diagonal shift method of the local jacobian matrix (See Kothari and Anand (1996)
and Guitino and Ortiz (1992)).

Though the preliminary results of above methods appeared to be effective, further
evaluation of them showed that either they are very expensive or unable to overcome
the problem arising from singularity or ill-conditioning in more complex strain paths.
A new method is introduced in this work which is faster and more efficient, compared
with existing methods. This new method is termed as “Dual mixed method”. The
method is based on a novel mathematical technique that allows us to bridge between
the computational mechanics, mathematics, or plasticity and the dislocation theory.
Based on this method, in any constrained optimization problem with several
constraints, instead of solving for several Iagrange multipliers directly, one may find
(an) intermediate function(s) that can be called “duality function(s)” and solve for (a)
Lagrange multiplier(s), which can be called “primary Lagrange multiplier(s)”, using
duality function(s). Then the design variables, state variables, and all Lagrange
multipliers of the problem can be obtained using the primary Lagrange multiplier(s)
and the duality function(s), indirectly. The duality function(s) can be developed using

existing empirical, experimental, mathematical, or physical based relationship

27

between design and state variables. Although this method is proposed for a crystal
plasticity problem, it can be used for any other type of constrained optimization
problem Figure 4.3 shows the schematic of the dual mixed model for solving a

constrained optimization problem

 

Min [800]
Subject to : fi(x) S 0 for i =1,...m

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Direct [ Or using dual mixed method
Solution
‘ Develop duality
Solve for 2i function(s)
l Solve for primary
Lagrange
Compute design multiplier(s);
and state
variables 7," J =1,...k

 

 

 

 

 

 

Figure 4. 3 Schematic of dual mixed model for a constrained optimization

4.3. Crystal Plasticity Based on Dual Mixed Model

The kinetics of the dislocations or the dynamic behavior of the dislocation under
external loads has been a subject of interest for a long time. A brief and good review
on dislocation dynamics can be found in Meyers’ book (1984) on mechanical

metallurgy. The first study on this subject may be related to the Gilman and

28

Johnston’s work (1957, 1959) in measuring the dislocations velocity in LiF. They
found a relationship between the velocities of dislocations and the shear stress (Figure
4.4). The following works by other investigators showed that for other materials, the
dislocation velocity is also a function of resolved shear stress. Some of these
investigations are: Stein and Low for Fe-Si; Ney et al. for Cu and Cu-Al; Schadler for
W; Chaudhuri et al. for Ge and Si; Erickson for Fe—Si; Gutmanas et al. for NaCl;
Rohde and Pitt for Ni; Pope et al. for Zn; Greenman et al. for Copper; Suzuki and
Ishii for Copper; Blish and Vreeland for Zn; Gorrnan et al. for aluminum;
Parameswaran and Weertman for Pb; and Parameswaran et al. for aluminum. Figure

(4.5) shows some results of these studies for different materials.

29

 

e VELOCITY OF (110) [ind SHEAR
l0 " WAVES = 3.6x 105 sin/soc

—-.-—— -_d

     
     
  

5
.n.
I

"' EDGE
COMPONENTS’I

r- \i
I
4 SCREW COMPONENTS

O
N

DISLOCATION VELOCITY, cm/sec
5.
m _.

YIELD STRESS
C.R.S.S

n i 114 n n ll'udun-l—I
QI LO I0 I00

APPLIED SHEAR STRESS, Kgf/mm?

 

 

 

Figure 4.4 Dislocation velocity versus resolved shear stress for LiF single crystal [Johnston and
Gilman (1959)]

30

 

2 2
( d9.) Fe-Si
e escrow)

DISLOCATION VELOCITY, m/s

 

 

 

l I l

 

 

10" 10° IO 102
RESOLVED SHEAR STRESS, 1.196

Figure 4.5 Dislocation velocity versus resolved shear stress for different single crystals
[Mechanical Metallurgy, Meyers (1984)]

31

Based on the above discussion, the theory of dislocation dynamics has been used to
develop a novel crystal plasticity model based on dual mixed model. For this purpose,
the duality function of the dual mixed crystal plasticity problem is obtained by using
the existing relations in dislocation dynamics theory, which expresses the average

velocity of dislocations as a function of resolved shear stress in a single crystal.

Assume that dislocations travel an average distance I- on slip systema. Then, the

plastic shear strain on slip system a is expressed by:
a a a—a

7" = ¢ .0 b l (4.17)
a

where ¢ is a correction factor taking into account the effect of activation of the

a a
different independent slip systems, ,0 is the dislocation density, and b is

Burgers vector.

The rate of shear strain then is:

7” = ¢apabaVa (4.18)

a
where V is the average dislocation velocity on slip system a.
The average velocities of dislocations on a slip system may be expressed generally as

(Gilman 1962):

—a a a
a a
where V0 is the limiting velocity and fp (7 ) is the average probability for

(I
dislocations to have a velocit V0 at iven instant. There have been several
y g

32

a
investigations, as explained in the first paragraph of section 4.3, to find fp (7 ) for

different materials. Based on the overall study of all these works, one may classify

a
fp (T ) in three different groups; power type average probability functions,

exponential type average probability functions, and average probability function
based on a linear combination of exponential functions (Gilman model). In the
following sections the proposed dual mixed crystal plasticity models based on these

three different dislocation velocity functions will be explained.

4.3.1. Power Type Dual Mixed Crystal Plasticity Model

Johnston and Gilman (1957 and 1959), Hahn (1962), and Greenman (1967) found the
following relationship between the dislocation velocity and resolved shear stress at a
constant temperature:

m

o (4.19)

where 7.. and V. are material constants.

Further investigations by Greenman et al. (1967), Pope et al. (1967), Gorman et al.
(1969), and Weertman et al. (1971 and 1972) showed a linear relationship between

the dislocation velocity and resolved shear stress for most of metals (Figure 4.5):

v-v 1-
. 7. (4.20)

and for the case of varying temperature:

33

v-V f— exp(—-g—)
o ’l' KT (4.21)

0

where Q is the activation energy, K is the Boltzmann constant, and T is the material
temperature in Kelvin. It must be mentioned that the goal of this work is not to
develop a temperature sensitive crystal plasticity model Nevertheless, because its
formulation is very similar to a non-temperature sensitive dual mixed crystal
plasticity model, it will be covered here but will not be tested or validated and
remains open for further researches.

Based on the above observation, we assume that for a small strain increment the

average velocity of dislocations on active slip system a can be obtained by:

a m
T
—a _ a _ a
V " Va a sgn(r ) (4.22)
y
a
where Ty is the critical resolved shear stress.
For the temperature sensitive case:
a m a
_ 3' Q
V“ =V°a — exp —— sgn(z'“)
Combination of equations (4.22) and (4.18) gives:
m
Ta
- _ a a a a a:
lfl-¢ p b Vo F Sgn(T ) (4.24)
y

And for the temperature sensitive case:

34

- _ a a a a Ta m Qa a
7a_¢ p b V0 :0,— exp[-E]Sgn(7) (425)

y

Combination of equations (4.24) and (2.8) yields:

N a m
a a a a T a a
D" = Z (I) p b v. —T,, sgn(T )P (4.26)
a=l y

and in the same manner substituting for the shear strain rate fiom equation (4.24) into

equation (2.9) yields:
m
Ta
P _ a a a: _ a a
O -221“); bv T, sgn(7 )w (427)
0:1 y
a n]
) .. ,4 4.4.4:... i
Assuming that _ ,0 o n T“ parameter is the same for all slip
y

systems, one can rewrite equations (4.26) and (4.27) as:

N n Ta "2
D p = ’12-? '7 Sgnwalpa (4.28)
2' 7y
"2
Q”- - 12%7 Sgn(Ta)Wa (4.29)
(1:17” Ty

where m=n1+n2 and n=n2+1.
According to the normality rule in plasticity, a potential f (o,q) can be considered in

such a way that:

35

D” = 280 f (o, q) (4.30)

Comparing equation (4.30) with (4.28) one can easily find the following relation for

active slip systems:

 

N n 2'“ ”2
a , = — — a P“
ef(6 (1) 2:1 7;,[73] 8811(7) (4.31)
Or
N ozP" n
f (o,q)=Z I 3'“ I +h(q,N) forNactiveslipsystems
0:1 y

(4.32)

In a similar manner one can easily derive a function for temperature sensitive case:

Pa
440» 4T): 2 I“,

“=1 )1

I Q“
ex —— +h N, ,T
p[ KT ( q )
(4.33)
for N active slip systems

And the rate of plastic deformation and spin for the temperature sensitive case are:

”2 a
DP =xiN n a:— ex ——Q—— s 2'“ P“
n 1'“ ”2 Q“ I
p _ __ __ a a
Q 4:1... 1', exp[ KTjsgnCt )w (4341))

36

With the above procedure, two duality functions for crystal plasticity problem (4.15)
were found. The first one, equation (4.32), is for constant temperature and the second
one, equation (4.33), is temperature sensitive. 1. is the primary Lagrange multiplier.
To complete the dual mixed crystal plasticity problem, a relation between the primary

Lagrange multiplier and the Lagrange multipliers of the crystal plasticity problem

(4.15), which are 7'25, must be found. Rewriting equations (4.24) and (4.25) by

substituting for 2 in these equations, the following expressions for ya can be

found:
no
. n T“ a
7“ *4}: 74 3811(7) (4.35.)
y y

and for temperature sensitive case:

n (1 n2 Qa
7a=27 7 exp —— sgn(2'“)
2'), 23 KT

(4.35b)
One may use the concept of the plastic work to find a general relationship between

the duality function, primary Lagrange multiplier, and the Lagrange multipliers of the

problem, 7a. The rate of the plastic work is expressed by:

 

N
W: D”. 6: 23—1; m=2z af 37 =o/lZ-aa—f32'“
=1 7

“:1 OT“ 80, (4.36)

37

And also the rate of the plastic work in a single crystal can be obtained using slip

rates as:

N .
e — e a
W " Z 707 (4.37)
a=l
Comparing the above equations, the following expression for the rate of slip

shear, 7a , can be found:

7“ = (Laraf (“9(1) (4.38)

Now in crystal plasticity problem shown in (4.15), instead of using several constraints
f a(09 (I) -<- O for all slip systems one can use a single duality function

f (0, (I) S 0 shown in equation (4.32) to find the primary Lagrange multiplier and

then rate of slip shears, Vs, as discussed above. This overcomes the problem

arising from ill-conditioning and computational intensity.

4. 3.2. Exponential Type Dual Mixed Crystal Plasticity Model

Rohde and Pitt (1967) found that a good correlation for the relationship between
dislocations velocity and stress could be obtained using an exponential type equation
for the average probability function. Based on theoretical calculation they found the
following expression for the average dislocations velocity:

_ kT AH 3(70-74')
V:—Kexp(——)exp[ kT ’] (4.39)

 

38

where h is the Plank’s constant, K a numerical constant, B an activation volume, AH

the enthalpy of activation, 7a the applied resolved shear stress, and 71 the long-

range internal stress.

For any slip system the above equation can be rewritten as:

_a_kT a AH B“(T“—t',“)
V ——h—K exp ——k—]—; exp kT (4.40)

 

a
Utilizing the above equation, a plot of log V vs T at constant temperature should

if it is assumed that B is not a function of

 

yield a straight line with a slope of

the applied stress. If B is a function of the applied stress, then the plot would not be a
straight line. For most metals, as can been seen from Figure 4.5, dislocations velocity
is a linear function of the applied stress and B remains almost constant over the range
of applied stresses.

In a strain increment for a constant temperature case, equation (4.40) can also be

written as:
--a _ Ia Ia a a a
V — K exp[C (It I— ty )] sgn(t' ) (4.41)

a
where 7y is the critical shear stress on slip system 0'.

Although the above equation correlates very well with the experimental data, it is not
in a suitable format for a computational algorithm When it is used in the above
format in a computational algorithm it causes stability problems. A better way of

expressing equation (4.41) is to modify as following:

39

cl l_ :I

a

For active slip systems '7 l‘ Ty Z O or T“ " 1— > O and also— 20 ,one can
y

also fit the experimental data using the following equations:

 

 

 

a:
’l"
-a _ a a _ a
V — K1 exp C1 1'“ 1 sgn(7 ) (442)
y
Or
2 K“ exp Ca-l—zf-l sgn(t’“)
2 2 Ta (4.43)
Y
Then using equation (4.18), the rate of slip on any slip system is expressed by:
It“!
° _ a a a a a a
ya—¢ p b Kl exp Cl Ta —1 Sgn(7 ) (444)
y
and the rate of plastic deformation is obtained by:
. lt‘l
p _ a a a a a a a
D —Z¢ p b K1 exp C1 ?—1 sgn(z' )P (4.45)

“=1 y

a a a a a a
Again let us assume that ’1 = ¢ ,0 b k1 0(7 try) is the same for all slip

e a — O a e 0
systems. Knowrng T - 0 . P , then equation (4.45) can be rewritten as:

DP=2Z—C— exp C“ halal—1 sgn(r“)P“

4.46
ar-l 45;? ( )

and spin can be expressed by:

40

 

N C“ a ozP“ a a
QP‘AZ‘TTCXP C l a "1 SW7 )w (4.47)
a: y y

Assuming that the associative flow rule can be used here, by comparing equations

(2.8), (4.30), and (4.46) one can find the following relationship:

N Ca a a : Pa a (1
30f (04(1) = 2767(1) C ——| Ta I-l sgn(3' )P (4.48)
a=l y y

Then a duality function for the crystal plasticity problem (4. 15) can be obtained as:

N a lo : Pal
f(o,q)=Zexp C 7—1 +hl(N,q)
“=1 y (4.49)
or active 5 z stems
f N ' l 'p sy

The relationship between the primary Lagrange multiplier/l and the Lagrange

multipliers of the crystal plasticity problem 7" s is given by equation (4.38) as:

° :19: Ca b.1114 a
7’” a exp 0 5811(7 )
r, r, (4.50)

And by following the same procedure for equation (4.43) one can obtain the

following duality function:

N a lozPal

g(o,q)=2exp D T +h,(N,q)
“=1 y (4.51)
forN active slip systems

The rate of the plastic deformation is given by:

41

D” =xii2iexp D“ JitTPl- sgn(t'“)P“

 

 

 

4.52
“:1 7;: , ( )
And the plastic spin is obtained by:
N D“ lo : P“|
p _ a a a
Q — [IX—1'“ exp D ——a sgn(t' )w (4.53)
“=1 y y

The Lagrange multipliers of the crystal plasticity problem are obtained by:

. D“ a o : P“ a

7“ = “—7“ exp D 2'“ sgn(r ) (4.54)

y y

It will be shown in the next chapters that with proper selection of functions h 1 and h;
in equations (4.49) and (4.51), the coefficients C“ and D“ are the same.

With the above duality functions and relations showing the relationship between
primary Lagrange multiplier and the Lagrange multipliers of the crystal plasticity

problem, the exponential type dual mixed crystal plasticity problem is complete.

4. 3. 3. Dual Mixed Crystal Plasticity Model Based on Gilman ’s Dislocation Velocity

Model

For any crystal structure, the dislocation velocity cannot increase without a bound as
stresses increase because there are always loss mechanisms present in any real
crystal. The dislocation velocity must saturate below the velocity of the sound wave

under a continuous stress field.

42

a a
The average probability function fp (7: ) must have a zero value when T = O ,

a
and must approach unity for large 7: . Gilman (1969) proposed a probability

function that has these limits:

fp (T) = exp(—D / T) (4.55)

where D is the characteristic drag stress. This function gives a good representation of
the data at medium to high stress levels. This function also can represent the plastic
yielding phenomenon by sudden transition from low to high value whenr= D/ 2.
However, it cannot represent the linear velocity-stress relations that are observed in
most metals. Also, due to the complicated derivatives it shows a poor form at low-
level stresses.

A better probability function that can represent the linear behavior at small stress

levels plus saturation at high stress levels is (Gihnan 1969):

fp (7, S) = 1 " eXP(-7/ S) (4.56)
where s is the coupling stress that acts across the glide plane. For small stress levels

this function becomes 77/ S . Therefore, the velocity has a linear relationship with
stress.
According to Gilman (1969) the behavior of almost any material can be represented

by a combination of equations (4.55) and (4.56) as:

17 = v: [1— exp(—t'/ s)] + V; exp(-D / T) (4.57)

43

Here two different dual mixed crystal plasticity models are developed based on

equations (4.56) and (4.57), respectively.

4. 3. 3. 1. Dual mixed crystal plasticity model based on equation (4.56 )

Let us rewrite equation (4.56) as:

Ta
V“=V: l—exp J2“; sgn(t'“) (4.58)

)’

Then using equation (4.18) the rate of shear slip on any slip system is expressed by:

- a a a 0: Ta a
y“=¢ p [7 VO l—exp —|—T-a-| sgn(7 ) (4.59)

)’

Assuming ’11 = ¢apabang1 (Tait?) and ’12 = ¢apabavng2(Ta,T;z)

are the same for all slip systems, then equation (4.59) can be written as:

 

n-l
. _ n If“ C: a__| “I
Va- 41—2,: ‘1? -/i2;-y;exp "C3 1.: sgn(z'“) (4.60)

Substituting for 7a in equation (2.8) the rate of plastic deformation is obtained as:

 

n—l
N n '7“ Ca ITaI
D” = — — — —3exp —C“— sgn(t'“)P“
Z} “‘4: r: “r: 3.; “-6“

Assuming the associative flow rule, equation (4.61) can be rewritten as:

a r“,r“ a r“,r“
DP :2, f‘( y)+/i, f“( ’) (4.62)
do do

 

 

Equation (4.61) shows two primary Lagrange multipliers. Therefore, this problem
has two duality functions as follow:

It

 

h<rflr§1=N> 2L:

:1 + hr (N, 7:) (4.63)
And
a a N “ital a
f2(T ,T,.N)=Zlexp “C3 7,..— +h2(N’“-y) (4.64)
a: y

Although it is possible to solve the dual mixed crystal plasticity model with two

duality functions shown above. To simplify the problem one may assume 1, = xi, = A

and rewrite equation (4.61) as:

 

n-l
1' “ Ta
P=AZ£E "I: :73? 63xP C: I751 Sgnwafl’a (4.65)
“y y y

Assuming associative flow rule, the duality function is obtained as:

 

 

 

 

 

 

N ozP“ n a ozP“
f(0.q)=Z T, +exp -C. T. +h(N.q)
“=1 y y ’ (4.66)
for N active slip systems

The plastic spin is then expressed by:

n-l

MIG—7— P“

 

”AZ-£2:

0:17)) Ty y

45

And the Lagrange multipliers of the crystal plasticity problem, 70 , are obtained by:

 

 

-l
n IozP“ n “ Io:P“I
' _ , , __3__ a__ a
7“ _2’ Ta z.a —Ta ——CXp —C3 Ta Sgn(T ) (4.68)
y y y y

4. 3. 3.2. Dual Mixed Crystal Plasticity Model Based on Equation (4.5 7)

Let us express the Gilman’s equation for dislocation velocity shown in equation

(4.57) as:

17'“- “1—ex —-|:—;I +V“ex —Sa sn(t'“)
"' P d P Ta g (4.69)

As mentioned before in equation (4.57), D is a parameter representing the plastic

yielding phenomenon by sudden transition from low to high value whens: D/ 2.

a
S has been introduced based on this concept and is assumed that
D = S“ = 27;”

Slip rate on any slip system is obtained using equation (4.18) by:

r“ 6*
V=¢“p“b“vs“ l—exp —|—7| +k“exp “€71- sgn(t'“)

Ty I T I (4.70)

where ka = Vf/V,“

46

Let us assume that there are three primary Lagrange multipliers that are the same for

all slip systems and can be expressed as ’11 = ¢ap abavp G1 (Ta, 7;“) ,
,1? = ¢“p“b“Vsz(r“,r;’), and 23 = ¢“p“b“VfG3(r“,r)“), The“
equation (4.70) can be expressed as:

° _ m ITaI ""1 Cf a_I_TaI+
V— “FIT; 42—156 XP -4C 7,,

Y Y )’

a a 4.71)
risk“£5-epr-C‘ ir—ialIIsgnW" )

Therefore the rate of the plastic deformation is given as:

N m 2'“ m“ a 1'“
NE; “FIQI —2,%epr-c:IT—aII+

y y y y

(4.72)

23k“—exp[—CS“ [5%] sgn(t'“)P“

Assuming associative flow rule, there are three duality functions in such a way that:

 

 

 

3g (7‘27“) 38 (7".7“) 8g (7“.1“)
DP = 1 y 2 y + 3 y
’11 do + ’12 do 43 do (4.73)
where
81(7a7 34:99 2 '7 _aIZmI +V1(N’7§) (4.74)
and

47

N 1'“
32(Ta,z-;',N)=Zlexp __Cfl—r-“l +y/2(N,t';') (4.75)
0= y

and
83(7“,r;’,N) = T(r“.r;’,N)+v/3(N.tf) (4.76)

where

 

 

N c“ S“ 1
‘I‘(7'“,7'“,N)= k“—5 t'“exp —C“ +C“S“ln +
y ; Ta I I 5 Ta 5 Ta

 

 

 

 

y

2 3
cgsa (car) (an) M

— _— +_— + . . .

1-l!-Ir"| 221W 331W
Now, the crystal plasticity problem can be solved using the above three duality
functions and equations (4.71) and (4.73). To simplify the problem it is assumed that
it is possible to find functionsG ,Gz , and G3 in such a way

that A) = A) = 23 = ’1 .Then equation (4.72) can be rewritten as:

m—l

 

N m IO:P“I 6' Io:P“I
DP T’I'Z a a . gexp f a +
a: 7'), 7y 1'), 7y
5a a O_' (4.78)
k“—ex —C 8 s n(7'“)P“
2')“ p 5 Io:P“I g

And therefore, the duality function is expressed as:

48

 

Io:P“I
+exp —Cf 1'“ +Y(Ta,T:,N) “I’ll/(Nail)
y

9
.l_l.
Q

for N active slip systems

(4.79)
The slip rate is defined as a function of the primary Lagrange multiplier as:
l Pa m-l a P“
. m 6 I I C I0 : I
70:]. 7 __3— ——LCXP —Cf—— +
r r r“ 7:“
y y y y
_ (4.80)

a
k“ E(fl—exp —C5“ 0-“ sgn(r“)

With above equations, the dual mixed crystal plasticity model is completely defined
for Gilman’s dislocation velocity model.

Before closing this chapter, it should be mentioned here that the duality functions are
not a yield criterion for a single crystal. In other words, a duality function can not be
used as a yield function for a single crystal. They can be only used to find the
direction of the plastic flow and also to“ find the Lagrange multipliers of the problem.
The other important point about the duality functions is that they are used for active
slip systems. Therefore, to find out about the current yielding situation or the current
number of active slip systems in a crystal, one should use equation (4.16) as a

yielding criterion.

49

CHAPTER 5
YIELD SURFACE FOR SINGLE CRYSTALS

5.]. Combined Constraints Method in Constrained Optimization

Another approach to overcome the problem arising from interdependency of plastic
slip shears is to use a single yield function for a single crystal. A new yield function
for single crystals based on combined constraints method is developed here.
According to Kreisselmeier and Steinhauser (1979) for problems with large numbers
of inequality constraints, it is possible to construct an equivalent constraint to replace

them. Let us consider the following optimization problem with several inequality

constraints ( f,- (x) 2- 09 i = 1...m ):

Min [g(x)]
Subject to: fi(X) 2 O for i :1me (5-1)

The following function can then be used as an equivalent single constraint:

h(x)=—%1n “zed-426)) (5.2)

Where p is a parameter that determines the closeness of the h(x) to the smallest

inequality, minIfoH. h(x) is usually called a KS—function. For any positive
value of p , h(x) is always more negative than the most negative constraint, forming

a lower band envelop to the equalities. As the value of p is increased, h(x) conforms

50

with minimum value of the constraints more closely. The value of h(x)is always

bounded by:

fn.(X)Sh(x)Sf.r.(X>-ln:)m) (5.3)

 

In a same manner one can easily extend the application of the above method for the

case when constraints f,(x) S 0, i=1...m . For this case the following function can be

used as an equivalent single constraint:

h(X) =%1n ieXPIPfr-(XU (5.4)

5.2. A Single Crystal Yield Surface Based on Combined Constraints Method
Let us consider the optimization problem of (4.15) for crystal plasticity problem

where the constraints are:

ozP“
fa(G,Q) =I—a—J—l $0 for a=1...N
T
y

(5.5)
where N is the number of slip systems.
The constraints of (5.5) can be combined using equation (5.4) as:
a
1 N Io : P I
f(o,q)=—ln ZCXP p _T‘1 (5.6)
p a=1 Ty

Here as the value of p is increased, f (o,q) conforms with maximum value of the

constraints more closely. The value of f (o,q) is always bounded by:

51

fmax (X) S f(6’ q) '<' fmax (X) + ln(N) (5.7)
p

 

Then the equivalent optimization problem to be considered with combined constraints

method is:

“Ml-4E (an:- 2)c 1:41:11

 

l
-(q.- 4):» 4:]
6n+’qn+ :ARG< >
( l 1) Subject to:
1 N I2:P“I
f(2.q)=—1n eXPP -1 50
p Zr: If

 

 

k

(5.8)
Based on the above combined constrained optimization method a new yield function

for single crystals can be proposed as:

p=|__° P"l_

m z.a (5.9)
y

f(0 q)— — plug Zexp—

where parameter m was introduced into equation (5.9) to give more flexibility to the
shape of the yield function.

Based on the normality rule, the rate of plastic deformation can be obtained by:

D“ = Aaaf(6,q) (5.10)

01'

52

 

 

isgn(6;P“)exp g I":faI_1 Pa
Dp-Za=1 g 7 m g - 7
. .3 (5.11)
mien, p I“: I_

sgn(o:P“) fl Io:P“I_

 

 

a exp a
V—2 g m g
— ' N p Imp/BI (5.12)
mZexp —- “—7——

/3=1 m 8

and plastic spin can be defined by

 

 

isgn(6;P“)exp [1 I“ f I we
QP— a=I 8 m g ,
. .3 (5.13)
mg... ,, Is}; I-

5.2.1. A New Model for Dislocation Velocity

Based on the proposed single crystal yield surface, a new equation for the velocity of
dislocations on a slip system, which are subjected to an external force, can be
presented. As mentioned before, the relationship between slip rate and the average

dislocation velocity on a slip system is expressed by equation (4.18). Substituting for

53

7a from equation (5.12) into equation (4.18) and assuming that

a a a a a a
A. = ¢ ,0 b TyF (7 97y 9T), the average velocity of dislocations on a slip

system can be expressed by:

ozP“
F(t'“,t'“,T)exp 5 [73—1—1 sgn(o:P“)

y

p o : P” 1 (5.14)

a—

 

<l

 

 

 

N
mzexp fl
fl=1 m “y

As discussed before, there are several empirical models that have been developed for
dislocation velocity. Unlike other models, equation (5.14) shows that stress state and
strain hardening on other slip systems also affect the dislocation velocity on a
particular slip system As will be discussed later, the parameter p is a function of
stacking fault energy (SFE) of the material. Therefore, equation (5.14) also shows
that there is a relationship between the velocities of dislocations in a crystal and its
SFE. To find the function F and validate the above-proposed model, several
experiments will be necessary. In this research, this model was not further studied and

was left as an open topic for future research.

5.2.2 Evaluation of The Effects of The Closeness Parameter ,0 on The Proposed

Single Crystal Yield Function

The so called closeness parameter, p , is an important parameters in single crystal

yield function shown in equation (5.9). The proposed yield function was studied in

54

order to find out how the yield function shape changes with p. For this purpose, let

us consider the yield function (5.9) for a plastic case to be:

 

 

 

1 N p 0' I P“
f(0, q) : Zln Zexp E Ta —1 :0 (5.15)
0= y

Then by manipulating equation (5.15) one can easily obtain the following

relationship:
N ~P“
e IE;_I_ _ e
gexp m 7'“ 1 —exp[m) (5.16)

For a plane stress condition (0'3 = 0 ), and in principle stress space, equation (5.16) is:

fiexp p I01B7+02P2gl_1 =exp[£]

“:1 m T“ m

(5.17)

Simplifying the problem with a plane stress condition makes it easy to evaluate the
effects of the yield function parameters on the behavior of the single crystal yield
surface.

Figure (5.1) shows the effect of p on the shape of the yield surface. The yield

surfaces are shown for the plane stress condition and scaled by the critical shear stress

70. As can be seen, by increasing the value of p , the shape of the yield surface

shows sharp vertices. It is well known that materials with high SFE, like aluminum,
show a yield surface with sharp vertices, which is related to the high mobility of the

dislocations in these materials. Therefore, one can see that there is a relationship

55

between p and SFE of a material. For materials with higher SFE, like aluminum,

higher value of p should be used.

 

02/10
0

 

 

 

_3 1 1 1 4L 1
-3 -2 -1 0 1 2 3
01/10

 

Figure 5.1 Effect of p on the combined constraints single crystal yield surface (crystal
orientation: (100)[001])

56

CHAPTER 6

SOULOTION ALGORITHM FOR PROPSED CRYSTL
PLASTICITY MODELS

Two different types of incremental formulation can be considered. According to an

additive decompositions formulation (for hypoelastic materials) the velocity gradient

is decomposed into the rate of deformation and spin tensors as L = D + W

while a multiplicative decomposition of deformation is expressed as L = FIN-1
for hyperelastic materials. In crystal plasticity one needs to calculate the grain
rotation at any increment. Therefore, crystal plasticity models are usually formulated
based on hyperelastic theory. However, this type of formulation significantly
increases the computation time. In this work the integrational algorithms for crystal
plasticity problems are developed based on the additive decomposition formulation.
Knowing that metals and alloys usually show a small elastic deformation and their
constitutive models can be expressed in rate form, this type of formulation is
completely valid.

Now let us rewrite the crystal plasticity problem in a strain increment [n, n+1], based

on the proposed models as:

rMin [ 24241)]
(on-i-l’an) = ARGi Subject to: >
f (2.0) S 0

(6.1)

 

 

where

57

2(2,q)=I§(o::::'—2)zc-‘ 2(02’i‘i’-2)+

%(q.-q)=D“=(q.-9)I ‘6'”

and f (2, q) is either one of the single crystal yield surfaces proposed in equations
(3.8) and (5.9) or is one of the proposed duality functions. One should notice that if
the function f (2, q) is a single crystal yield surface then the parameter N in the
equation denotes the total number of slip systems in a single crystal but if it is a
duality function then N is the number of the active slip systems.

Let assume that C and D are positive definite, or in other words, the objective
function Z is strictly convex in its defined domain. Although it is not proven, but let
us assume that all the proposed duality functions and single crystal yield surfaces
discussed in previous chapters are also convex. Now the solution to the above
constrained minimization problem is to solve the following associated Lagrangian

problem:

L(Z,q,)t) = X(2,(])+x\f(2,(l) (6.3)

Now the corresponding optimality condition is:

V BL — ria

82L = 8—2 : —C 1:("it+i — 6n+1)+ )‘aof (“n+1aqn+1) = 0
4 aqL : “n+1 — an + Aaqf(6n+l’qn+1): 0

81L 2 f (on+1’qn+l) S O (6.4)
[A Z 0’ Af(6n+1,qn+1) : 0

 

58

These equations are equal to:

85+1 = 8: +“aof (“n+1’qn+l)’
“n+1 : an +llaqf (an-rl’qnfl)’

f (6n+1’qn+l ) S 0’
2 2 0, (65’

if (“n+1 ’ qn+l) = 0

Also let us add the following expressions to the above equations;

an+1 = an + V“ (Au)

_ . P (6.6)
“n+1 — C ' (£n+1 _ 8n+1)

As mentioned before, vectors a and q define the plastic dissipation in a material. The

vector of state variables q for the case of isotropic hardening at constant temperature

will contain a single value for the grain flow stress, -O-g . The vector of internal

variables a will contain a single value for the grain level effective plastic strain,

P . .
83 . Then, accordrng to the above equations

As; = 2858 f (6,33)

(6.7)

The grain flow stress is a function of the critical shear stress of all slip systems in a

single crystal. Therefore, equation (6.7) can be rewritten as:

—p _ 1 2 a
A88 — 1H(Ty,7y,...,7y) (6.8)

59

—P
The evolution of 88 with deformation is needed for the dual mixed crystal plasticity

model, which has been deve10ped based on the Gilman’s dislocation velocity model

(equation (4.78)). For all other proposed crystal plasticity models, the evolution of

P . . .
83 wrth deformation rs not necessary.

Equations (6.5) and (6.6) can be solved implicitly or explicitly. There are several
computational techniques for solving the Optimization problem (6.1) implicitly or
explicitly. To increase the computational speed and efficiency, an explicit
formulation is used in this work. To solve the above crystal plasticity problems
explicitly, an integrational algorithm based on the cutting plane algorithm is used.
According to the cutting plane algorithm procedure, at a given increment, constraints
of the optimization problem are linearized iteratively to find hyper planes that cut the
solution domain at any iteration and remove the part of the solution domain which
dose not contain the solution to the problem. Therefore, by this manner the solution
domain shrinks at every iteration until the final solution to the problem is found for
the increment. Linearization of the duality function or the single crystal yield surface

gives:

f(0',,+1,q,,+1) : f(6;+1’q:1+1)+ aof(a:t+l’q:t+l) : A“ +

N
Zlaz-ilffflijl’q:1+l)Az-yf : O (6.9)

0+]

Combining equations (2.14) and (4.38) the following relationship can be written for

At“;

60

N
0! _ afl o 0
AT), — 2'2 h az'fl f (“n+1 ’ qn+1) (610)
[3:1
Also from equations (6.5) and (6.6):

A0 : _N’C I adf(o:l+l’q:l+l) (6-11)

Combining equations (6.9), (6.10), and (6.11) one obtains the following expression
for the incremental consistency parameter:
2 = f°
a a N N 8 B (6.12)
_I. : C : __I; _ Z Z hafl f . f

do do (Fl/3:1 at“ 8rd

 

 

Where f. is the value of the function at the beginning of the iteration. Using
equations (6.5), (6.6), and (6.12) together, the crystal plasticity problem can be solved
explicitly. If the orientation of a grain with respect to the co-rotational materials
coordinate system is defined by matrix Q (defined in appendix A) then a general

integration algorithm for proposed single crystal yield function crystal plasticity

a a
model can be expressed as following (For simplicity in notation, assume Ty = g )2

A a A
1. Given the known quantities (“j , g j , A8} , Qj) for a grain (crystal) at the

beginning of the increment j in the material co-rotational coordinate system, rotate the
strain tensor increment and stress tensor to the crystal coordinate system using grain

orientation matrix:
a. = Q .6 .QT.
J J J J

61

_ A T
- QjAngJ'

°_ (0)_ e. a(O)_ a (0)_
2.1nitialize:l—O, 0j+1—6j+C 'A£,gj+l ‘gjflmd d’ijn" O

(0) . . <0
3. Calculate fj+1 , usrng equatron(16). If j(+1) — , go to step 9

4. Compute the increment of the plastic consistency parameter from equation (31):

(l)
Alf-i): fn+l

j+1_ (')f(i)lcezaf(i)1 N N afjm Bf“)

($8
60-2212 87“] Bag

a=l ,B=1

 

5. Compute stresses, slip resistances, and the yield function:

 

() ()
(' 1)___ a;_() () df(oj‘+1, 8““! )
1+ 1 l e ,

af(o(_i) (1'))

j+1’ gj+1

 

(i) _ (01
[373-11-11344370, For a=1,..N

138a (0:21101 gfifoa Sgn 1_g5“:“) [q( + (-1 —q MT] IA7}6(1)I
Agj+l +1

For a =1,...N

 

 

a (i+l) a (i)

_ a (i)
gj+1—gj+1 + Agj-I-l For a =I,...N

62

(i+l)
6. If fj+1 S T01 go to step 8, otherwise go to step 7

7. Update the consistency parameter:

’10“) = ’10) +A20‘)

j+1 j+1 j+1
Set i<——i+1. Goto3

8. Update the grain orientation

9. Rotate the stresses back to the material co-rotational coordinate system and go to

step 1 for the next grain.

10. Homogenization

63

CHAPTER 7

EVALUATION OF THE PROPOSED MODELS BY MODELING
OF METAL FORMING PROCESSES

To obtain the parameters of a crystal plasticity model, only a texture analysis and a
uniaxial tensile test in an arbitrary direction are necessary. Therefore, crystal
plasticity models are powerful tools for analyzing metal forming processes. They can
be either used directly to models a metal forming process or to find the coefficients of
phenomenological models.

The proposed models discussed in previous chapters were applied to the modeling of
Al6061T4 polycrystalline sheet, involving uniaxial tension and tube hydroforming
with different number of grains. The initial texture for the extruded Al606l-T4 tube
was obtained using OIM measurement as shown in figure (7.1). Several integrational
algorithms based on an explicit procedure for crystal plasticity model based on the
multi yield surface proposed by Anand (1996), modified power type yield surface,
dual mixed method, and the proposed combined constraint yield surface were written

and implemented as VUMATs into the ABAQUS software.

(111) pole figure
(1005 grain orientations)

     

' w >4

 

 

 

 

(b)

Figure 7.1 a) pole figure (111) with 1005 grain orientatiom representing the undeformed
extruded tube, and b) CODF for the undeformed Al6061-T4 aluminum extruded tube

65

7.1. Modeling of The Uniaxial Tension Test
Figure (7.2) shows the finite element uniaxial tensile model with 469 S4R shell
element, used in these analyses. Hardening coefficients (equation (2.16)) for Al6061-

T4 alloy was obtained using a multisurface crystal plasticity model proposed by

Anand (1996). The components of ma and 5a, which are needed for equation
(2.1), and the 12 slip systems used for this fee material are shown in Table (7.1). The
elastic constants for the Al6061-T4 extruded aluminum tube were determined to be;
C11: lO8GPa, C12=6lGPa, and C44=29GPa and the following parameters were

obtained by curve fitting to the experimental uniaxial tensile data for Al6061-T4;
h0 =1285MPafL'S =172MPa ,7: = 62MPa, and a = 2.0,

These parameters are kept fixed in all proposed crystal plasticity models, then other
parameters of the proposed crystal plasticity models were obtained as will be

explained in the following sections.

 

Figure 7.2. A uniaxial tensile model containing 469 S4R shell elements used in these simulatiom

66

Table 7.1. Components of ma and S“ in crystal coordinate for f.c.c. materials.

 

a

 

a m“ S
1 _1____1_0 iii
JEJE fifix/i
2 _10_1_ 111
7235 75737;
3 O1___1_ 111
3J5 757575
4 _1-0_1_ _LLJ.
JEJE J§J§J§
5 _1__1_0 _11__1_
3J5 75756
6 O__1_____1_ _111
JEJE 75753
7 _101 _1__1_1_
JEJE fiWI/i
8 _1__L _1__1_1_
JEJE fififi
9 110 1_11
7577 157373
10 _110 __1__LL
ff [is/5J5
1‘ _1_0; __1___1__1_
J5 2 J3J§J3
‘2 0_L_L ___1__i_l_
JEJE J3J§J§

 

67

7.1.1. Modified Power-Law Yield Surface
As discussed in equations (3.38) and (3.9), for this model there are two other

parameters in the yield surface that needed to be identified. 12’ can be obtained using

equation (3.11) and k is obtained using curve fitting to the experimental data. For this

*

material n = 8 and k = 6 were found to be the best values. Figure (7.3) shows the
uniaxial results obtained by simulation using modified power-law yield surface model
using 50 grains. Two different values of k: 4 and k=6 were used in this model. As can
be seen, the results of the model using k=6 shows a good match with the experimental
data. Though with using appropriate values for k one may obtain good results using
the single crystal power yield surface, there are some oscillations in stresses, which

may be related to the imprOper rotation of the grains in this model.

400

 

350 w

300 ~

250 e

200 ~

Stress (MPa)

150 ~

100 1 —experiment
k=6 and n=8
- - - k=4 and n=8

 

50~

 

 

 

0 0.05 0.1 0.15 0.2 0.25 0.3
True strain

Figure 7.3. The stress - strain curves obtained by modified single crystal power yield surface

68

7.1.2. Dual mixed crystal plasticity models

7.1.2.1. Power Type Dual Mixed Crystal Plasticity Model

For power type dual mixed crystal plasticity model the function h(q,N) in equation

(4.32) should be determined. We assumed that h(q,N) = —l. The power of the yield
function was found to be n=2. Then n2=n-1 =1 and if one assumes that n1=0 then

m=1. And this completely matches with the linear relationship that exists between

dislocation velocity and the applied resolved shear stress in an aluminum single
crystal. Figures (7.4) and (7.5) show the results of the power type dual mixed crystal
plasticity model for uniaxial tension of A16061-T4 alloy using 50 grains and the same
hardening parameters as mentioned above. As can be seen the results of this model
show a much better match with experimental data. Figure (7.6) shows the results of
the simulation for the uniaxial tension using different number of grains. As can be

seen using more than 10 grains yields a similar result.

69

350

 

300 -

250 -

R

200 -

150 .

True Stress

100 <

 

— Experiment
50 .

——-Crystal plasticity 50grains

 

 

 

 

o 0.05 0.1 0.15 0.2 0.25 0.3
True Strain

Figure 7.4. Comparison of the stress-Strain curves obtained by simulation using power type dual
mixed crystal plasticity model and experiment

Mises Stress

+3 . 039e+02
+2 . 839e+02
$5.6392+02

+2 . Z40e+02
+2 _ D41e+02
+1 . 841e+02
+1 . 642e+02
+1 . 442e+02

 

Figure 7.5. Mises stress distribution in a uniaxial tension model after 15 % strain obtained by
power type dual mixed crystal plasticity model

70

350

 

 

 

 

 

 

 

 

 

 

 

 

300 .
250 ~
I
m 200 *
m
0
5
a 150 d
:3
h.
100 1 Crystal plasticity 50grains
109rains
809rains
50 . — - —1009rains
0 . T . , f
0 0.05 0.1 0.15 0.2 0.25 0.3
True Strain

Figure 7.6. Stress-Strain curves obtained by power type dual mixed crystal plasticity model using
different number of grains

7.1.2.2. Exponential Type Dual Mixed Crystal Plasticity Model

For the exponential type dual mixed crystal plasticity model the values of C a and

D“ are assumed be the same for all slip systems and equal to C and D, respectively.

The results of the simulations of the uniaxial tension of Al606l-T4 alloy using 50
grains showed that if h,(q,N) =—l for the model in equation (4.49) and
h2(q,N)=—exp(D) for the model in equation (4.51) then D=C. By keeping the

hardening parameters constant and the same as those obtained by Anand’s model, C
and D were found to be 5.0 for Al6061-T4. The uniaxial tension stress-strain curve

for Al6061-T4 predicted with exponential type dual mixed crystal plasticity model

71

was exactly the same as the results obtained with the power type dual mixed crystal

plasticity model, as shown in figure (7.4).

7.1.2.3. Dual Mixed Crystal Plasticity Models Based on Gilman ’3 Models

For dual mixed crystal plasticity models based on Gilman’s dislocation velocity

models, equations (4.66) and (4.79), it was assumed that h(N,q)=V(N,q)=—l.
Also it was assumed that for all slip systems C3“ =C3, C: = C , and C50” =C5.

Using the same hardening parameters obtained by Anand’s crystal plasticity model
and curve fitting the results of the simulation of the uniaxial tension to the
experimental data, the parameters of the duality functions (4.66) and (4.79) were

obtained to be as following:

For duality function (4.66); n=1, and C3 =16 .
For duality function (4.79); m=1, C 4 =16, and C5 =10.

Further investigation of the models showed that there is a relationship between the

SFE of the material and parameters C3 and C 4. In order to get results without

oscillation in the stress-strain curve, C3 and C4 should be set to the following

equation:

F _
C3:C4=EI;X10 3 (7.1)

This equation is the same as equation (3.10) proposed by Gambin. Also for duality

function (4.79) to avoid oscillations in stress results, the best value for C5 can be

obtained from the following equation at effective plastic strain, 5P = 0.002.

72

dc‘)‘ _
C5 2 '22:??? X10 2 (Unit of stress is MPa) (7.2)

where 5" is the effective flow stress which evolves by effective plastic strain.

For evolution of 6' in Al6061-T4 alloy the Voce equation was used:

6 = A — B exp(-CE'P ) (7.3)

The parameters of Voce equation; A, B, and C, were obtained by fitting the Voce
equation to the macroscopic uniaxial tension experimental data by which the values

of A=338.559 MPa, B=206.249 MPa, and C=9.822 were obtained for Al6061-T4

afloy.

7.1.2.4. Combined Constraints Crystal Plasticity Model

In addition to the strain hardening parameters, this model has two more unknown
parameters m and p. It was assumed that m=1, then the closeness parameter p for
Al6061-T4 alloy was found by curve fitting to the uniaxial tension experimental data,
which resulted in a value of 40.75. The comparison of this model with Gambin’s

investigation (1997) showed that one might obtain p for different materials using the

following expression:

2.51“
z — x10‘3
,0 Gb (7.4)

73

7.2. Modeling the Hydroforming Process of Aluminum Tubes

Hydroforming of aluminum-extruded tubes has recently been a topic of interest in
automotive industry. The problem with modeling the hydroforming of an extruded
tube, using a phenomenological yield function, is to find the appropriate values for
the anisotropy parameters. Crystal plasticity models, which only need as input the
uniaxial tension properties in the extrusion direction and the initial texture of the tube,
can be used to find these anisotropy parameters for phenomenological yield functions.
As part of verifying the current crystal plasticity model, tube bulging and tube
hydroforming into a square die experiments were conducted and deformed shape and
hoop strain distribution at the maximum pressure were measured for comparison.
Table (7.2) shows the dimensions of the Al606-T4 extruded aluminum tube used in
the tube hydroforming experiments. In all tube hydroforming analyses 50 grains were
used in crystal plasticity models with initial orientations obtained from ODF shown in

figure (7.1).

Table 7.2. Dimensions of the Al6061-T4 tube
Length (inch) OD (inch) Thickness (inch)
8.0 2.0 0.05

 

 

 

 

 

 

 

7.2.1. Modeling of The Tube Bulging Into The Round dies

Figure (7.7) shows the geometry of the die used to model the bulging of the Al6061-
T4 tubes. Due to the symmetry and to decrease the simulation time, only one-quarter

of the tube was considered in the modeling (figure 7.8). 174 shell elements were used

74

in the model as is shown in figure (7.8). Boundary and forming conditions used in the
simulation were similar to the experiment where a maximum axial feed of 3mm was
applied to each end of the tube while the tube was expanded with a maximum
pressure of 14.18 MPa. ABAQUS symmetry boundary condition command
(YSYMM) was used to apply boundary condition on the edges of tube in the
symmetry directions. The friction coefficient between the die and the tube was
assumed to be 0.02. During the simulation, both the internal pressure and the axial
feed were linearly increased to their maximum value as a function of time. Figure
(7.9) shows a bulged aluminum tube and figure (7.10) shows the results of the tube
bulging simulation using the crystal plasticity proposed by Anand et al. (1996) at
MIT. Figures (7.11) to (7.16) show the results of the simulation of the bulging of the
aluminum tube using the proposed crystal plasticity models. As can be seen form
these figures and figure (7.17), which shows the hoop strain distribution along the
length of the bulged tube, all crystal plasticity models predict almost the same hoop

strains with a good match with experimental data.

75

g’léé
ifff‘i

 

Figure 7.7. Geometry ofthe die med inmodeling ofthe tube bulging

 

Figure 7.8. Mesh and the geometry of the one-quarter model of the tube used in modeling of the
ulging

76

 

Figure 7.9. Bulged aluminum tube

     
 

Hoop strain
+2.023e—01
+1.857e—01
+1.69le-01
+1.525e—01
+1.360e—01
+1.194e—01
+1.028e-01
+8.622e—02
+6.964e—02
+5.307e—02
+3.649e—02
+1.991e—02
+3.337e—03

figure 7.10. Results of the tube bulging simulation using Anand’s singular value decomposition
multi yield surface crystal plasticity model

77

     
 
  

Hoop strain
+2.054e-01

+2:023e-02
+3.397e—03

Figure 7.11. Results of the tube bulging simulation using modified Gambin crystal plasticity
model

   

Hoop strain

+2.035e—01
+1.8686—01
+1.702e-Ul

+3.379e—03

Figure 7.12. Results of the tube bulging simulation using power type dual mixed crystal plasticity
model

78

     
 
  

Hoop strain
+2.045e—01

+2.011e—02
+3.349e—03

Figure 7.13. Results of the tube bulging simulation using exponential type dual mixed crystal
plasticity model

   

Hoop strain

+2.058e—01
+1.889e-01
+1.720e—01
+1.552e-01
+1.383e-01
+1.214e—01
+1.045e—01
+8.768e—02
+7.081e—02
+5.394e—02
+3.708e—02
+2.021e-02
+3.339e-03

Figure 7.14. Results of the tube bulging simulation using dual mixed crystal plasticity based on
Gilman model (equation 100)

79

     
 
  

Hoopsflam
+2.047e—01
+1.879e—01
+1.711e—01
+1.543e—01
+1.376e—01
+1.208e—01
+1.040e—01
+8.722e—02
+7.045e—02
+5.367e—02
+3.689e—02
+2.011e—02
+3.329e—03

Figure 7.15. Results of the tube bulging simulation using dual mixed crystal plasticity based on
Gilman model (equation 112)

   

HoopsUam

+2.029e-Ul
+1.863e—Ul
+1.696e-Ul
+1.53De—Ul
+1.364e—Ul
+1.1988-Ul
+1.031e—Ul
+8.6SDe—02
+6.987e—02
+5.324e—02
+3.662e-02
+1.999e—02
+3.361e—U3

Figure 7.16. Results of the tube bulging simulation using combined constraints crystal plasticity
model

80

 

 

 

 

 

 

 

 

 

 

 

---0---E irnent
0.25 _ I I Am
_ i . ———E)ponenfia|dualm'xed
_ i ! Powertypedualmixed
— # 0 ",0. —— Gilman dual mixedt
0.20 r """"" 7 i """ —-—Gilamdualmixed2
[ 1 ; Contained constraints
— j d" ; —ModiiiedGambin
,5 Q15 : ....................... €15: .............. 1 ...............
s ~ E '
in r .
a r 9 :
8 ' -' 1
a: 0.10 ;- -------------- r, ---------- 7 ---------- .3; -----------------
r I -' l G ’
0.05 r ——————————————— — ——————————————————————————————————
0.00 (II-v Gr 1 1 1 1 J 4 10
0 50 100 150 200

Location from one end (mm)

Figure 7.17. Comparison of the hoop strain distribution along the length of the tube obtained by
proposed crystal plasticity models

In order to evaluate the mesh dependency of the algorithm used in the above crystal
plasticity models, the simulation of the tube bulging was repeated using a very fine
mesh of 2970 shell elements as shown in figure (7.18). As the computational
algorithms used in all the above crystal plasticity models are based on the cutting
plane algorithm, mesh dependency was tested only on one crystal plasticity model.

For this purpose, the modified power type single crystal plasticity (modified

81

Gambin’s model) was used here. Figure (7.19) shows the hoop strain distribution for
this analysis. A comparison of this figure with the results shown in figure (7.10)
proves the mesh independency of the computational algorithm used in the above
crystal plasticity models and that with a rather coarse mesh one can obtain good

results.

 

Figlu'e 7.18. Mesh and the geometry of the one-quarter of the tube with fine mesh (2970 shell
elements) used in modeling ofthe tube bulging

82

     
 

Hoopsuam

+2.035e—01
+1.867e—Dl
+1.698e-Ul
+1.53Ue—Ul
+1.3628-Ul
+1.194e—Dl
+1.026e—Ul

+1:859e—02
+1.788e-03

Figure 7.19. Results of the simulation of the bulging of an aluminum tube using the modified
single crystal yield surface crystal plasticity model with a fine mesh (2970 shell elements)

Figure (7.20) shows the results of the crystal plasticity model for tube bulging of
aluminum at different times and its comparison with experiment. As can be seen,

there is a good match between the simulation and the experiment.

83

 

 

Experiment Simulation

Figure 7.20. Comparison of the experiment with crystal plasticity simulation of tube bulging
process of Al6061-T4 tube at different forming stages
Crystal plasticity models are usually very expensive, computationally. Therefore, it is
of interest to develop novel crystal plasticity models that are more efficient or less
expensive. To evaluate the computational speed of the proposed models, the total
CPU times for simulating tube bulging and hydroforming of tubes into square dies
were recorded and compared. Table (7.3) shows a comparison between CPU times to
simulate the tube bulging process with different proposed crystal plasticity models.
As can be seen from table 3, the exponential type crystal plasticity model has the best

computational speed, while the SVD multi-yield surface crystal plasticity model is the

most expensive among the presented crystal plasticity models. It is also clear that the
crystal plasticity models developed based on single crystal yield surface (i.e. modified
Gambin and combined constraints model), are faster than Anand’s model. Dual
mixed models are faster than both single crystal yield surface models and the Anand’s

model.

Table 7.3. CPU times for the simulation of tube bulging process using different

 

 

 

 

 

 

 

 

crystal plasticity models
Crystal plasticity model CPUtime (hour/minute/second)
SDV multi-yield surface (Anand’s model) 02:49:09
Modified Gambin model 02: 19:00
Combined constraints model 01:29:01
Power type dual mixed model 00:52:43
Exponential type dual mixed model 00:38:46
Gilman based dual mixed model 00:44:37

 

 

One may draw an interesting conclusion by comparing the power type crystal
plasticity models to exponential type crystal plasticity models. For example, the
Modified Gambin model (equation 3.8) which is based on a power type yield surface
for a single crystal, is more expensive than combined constraints model (equation 5.9)
which is based on an exponential type yield surface for a single crystal. The power
type dual mixed model is also more expensive compared to exponential type dual
mixed mode]. Therefore, constitutive or mathematical models expressed in an

exponential form are computationally faster than power type ones.

85

 

7.2.2. Modeling The Hydroforming of An Aluminum Tube Into A Square Die

To evaluate the accuracy of the model in predicting the deformation of the tube under
complex strain paths, the proposed crystal plasticity models were used in the
hydroforming simulation of the extruded aluminum tube into a square die. For this
purpose, one quarter of the extruded tube (see Table 7.2) was modeled using 2970
shell elements (S4R). The boundary conditions used in the simulation were similar to
those used in the tube bulging experiment where a maximum internal pressure of
22.34 MPa was gradually applied to expand the tube into the square die, while its
ends were fixed. Figure (7.21) shows the geometry of one-half of the square die and
figure (7.22) shows the mesh and the geometry of one-quarter of the initial tube used

in these simulations.

86

 

Figlu'e 7.21. The geometry of one-half of the square die used in the hydroforming simulation of
aluminum tubes

 

Figure 7.22. The finite element mesh used in the hydroforming simulation of one-quarter of the
aluminum tube

87

Figure (7.23) shows the comer area of the actual hydroformed tube at the maximum
pressure of 22.34 MPa. Figure (7.24) shows the failure location in the hydroformed
tube once pressure was increased beyond 22.34 MPa. Figures (7.25) to (7.30) show
the contour of the hoop strain predicted by different proposed crystal plasticity
models. By comparing figures (7.24) and (7.25) to (7.30), it can be concluded that the
crystal plasticity model has correctly predicted the location of strain localization and

potential site for fracture.

 

Figun'e 7.23. The aluminum tube hydroformed into a square die

 

Figure 7.24. Failun'e location in the hydroformed aluminum tube

88

  

.ZBle—Ol
.171e-01
.0626-01
.5198-02
.421e-02
.3246-02
.226e-02
.129e-02
.0316-02
.934e-02
.836e-02
3.386e-03

   

589e—03

Figure 7.25. Hoop strain distribution in the hydroformed aluminum tube predicted by modified
Gambin crystal plasticity model

 

\A'\ nu.
. u

a.

 

 

 

Figure 7.26. Hoop strain distribution predicted by SVD multi sun-face (Anand’s) crystal plasticity

model

+1.188e—
+1.112e-
+1.036e—
+9.596e—

 

 

'7

Itssa‘a-‘A . n. Au». A'enis an H”,

 

 

Figure 7.27. Hoop strain distribution in the hydroformed aluminum tube predicted by combined

+1.193e— I
+1.1178— Z
+1.04Ue- ‘

+9.64le—
+8.878e—
+8.115e—
+7.351e—
+6.588e—
+5.825e-
+5.062e—
+4.299e—
+3.536e—
+2.773e—
+2.009e—
+1.246e-
+4.831e—

—2.800e— I

 

constraints crystal plasticity model

'

.« eAQK'IA‘A-en1A\-1 A\I\L -aa

 

 

Figure 7.28. Hoop strain distribution in the hydroformed aluminum tube predicted by power

typedualnuxedcrymalphmficfiyrnodel

9O

 

+
m
.
D
d
.—I
(D
1
wwNNNNNNNNNNNHHHH

 

 

Figure 7.29. Hoop strain distribution in the hydroformed aluminum tube predicted by
exponential type dual mixed crystal plasticity model

+
Ch
.
Ln
\D
D
‘P

 

 

 

 

Figure 7.30. Hoop strain distribution in hydroformed aluminum tube predicted by Gilman based
dual mixed crystal plasticity model

91

Experimentally measuring the hoop strain in the corner area of the tube, shown in
figure (7.23), was very difficult and as shown in figure (7.31) there were some
oscillations in the measured hoop strains. However, on the average, the simulated
results obtained from the crystal plasticity models agreed very well with the

experimental data.

 

 

 
    
 

 

 

 

 

 

 

 

0.14 0 Experiment
- ----- Anand
—ModifiedGambin
g , -------- CorrbinedConstraints
0'10 — Ipi—-—PowerDualMxed
: ' — Exponential Dual Mixed
E ; - ----- - Gilman Dual Mxed
‘5 ______________ _____: ____________ = _____________ 1 ____________
g 0'06 . i 7 f
0.02 ------- ------------- --------- . ------------ —
i 40 80 120 160 2io
-0.02 ' ' '

 

Location from ends (mm)

Figure 7.31. Experimental and numerically predicted hoop strain distribution along the length of
the hydrofoer tube

Table (7.4) shows a comparison of the CPU times for above tube hydroforming

simulations by different crystal plasticity models.

92

Table 7.4. CPU times for simulation of hydroforming of tube into square die using
different crystal plasticity models

 

 

 

 

 

 

 

 

Crystal plasticity model CPUtime (hour/minute/second)
SDV multi-yield surface (Anand’s model) 25:40:05
Modified Gambin model 23:31:03
Combined constraints model 19: 13:02
Power type dual mixed model 13: 12:32
Exponential type dual mixed model 10:03:46
Gilman based dual mixed model 11:17:41

 

 

93

 

CHAPTER 8

3D MODELING OF THE SURFACE PRPERTIES OF
CRYSTALLINE MATERIALS

Electron work function (EWF) and phonon emission are two important surface
properties in materials. EWF usually is defined as the minimum energy required to
remove an electron from the interior of a solid to a position just outside the solid
[Ashcroft and Mermin, (1976)]. EWF is one of the fundamental characteristics of a
metallic solid and can be used in the studies of tribological phenomena including;
wear, friction, oxidation, deformation, phase transformations, changes in surface
composition, and surface adhesion, etc. See DeVechio et al. (1998), Zharin et al.
(1998), Bhushan et al. (2000), Zharin (2001).

Phonon emission is another important parameter, which is usually considered as a
measure of the thermal properties of a surface. Surface condition has a significant
effect on these two parameters. Several investigations have shown that surface
roughness, imperfections, dislocation density, grain boundary, plastic deformation,
and crystal orientation have significant effects on both EWF and phonon emission of
a surface [Van Sciever et al. (1996), Li et al. (2002), Romanowski et al. (1988), Haas
et al. (1977), Stranger et al. (1973), Li et a1 (2004)].

Li et al. (2004) reported that plastic deformation always decreases the EWF of
crystalline materials irrespectively of whether it is tensile or compressive. Levitin et
al. (1994) indicated that the EWF variation is probably related to the formation of
new surfaces caused by deformation. In the following works, Levitin et al. (2001)

suggested that the change in the EWF caused by plastic deformation might be also

94

related to dislocation behavior. More recently, Li et al. (2002) proposed a simple
model to show how dislocations impact EWF. In a different work, Zharin et al.
(2001) showed that the variation of EWF depended on the roughness of deformed
surfaces. It was also pr0posed that the change in the Fermi level or the band gap
during the plastic deformation may also affect the EWF [Mamor et al. (1995) and
Rueda et al. (1999)].

Despite intense research over the past few decades, a complete understanding of
mechanisms that are responsible for the EWF changes during the plastic deformation
of a crystalline material has not yet been achieved. Fore example, there are no
significant researches about the effect of the surface microstructure on the EWF and
phonon emission of the crystalline materials. Usually experimental study of these
effects is very difficult, expensive, or sometime impossible. Therefore, a
computational investigation can be a good replacement for experimental studies.

To investigate the relationship between the surface texture and large area EWF and
phonon emission, a multi-scale computational scheme coupled with a homogenization
technique, which allows for accounting for the real surface microstructure of
materials, should be used. However, this type of mathematical modeling is very
expensive or impossible and also the more important drawback is that a mathematical
framework relating the surface properties to the EWF and phonon emission have not
yet been established. Considering above drawbacks, the current investigation is only
focused on a few grains having the main components of a surface texture on a small
region on the surface of the material. These types of simulations would be very

powerful tools for studying the effects of the adjacent grain orientations and different

95

surface texture components on the surface dislocation density and surface roughness.
In this work the combined constraints crystal plasticity model is used to investigate
the effects of the surface texture and crystal orientation on the surface roughness and
dislocation density, and therefore on EWF and phonon emission of the plastically
deformed crystalline materials. In a given strain path, single crystals with the same
orientation but different mechanical properties and crystal structures can have
different plastic behaviors and therefore, different impacts on plasticity induced EWF
and phonon emission of materials. In this investigation, however, the high purity
niobium [Zamiri et al. (2006, 2007)] as an example of crystalline material to study the
surface texture effects on EWF and phonon emission. Niobium has a bee crystal
structure and has 24 {110}<111> and {112}<111> type slip systems [Baars et al.
(2008)]. Both EWF and phonon emission are important parameters in high purity
superconducting niobium which is used in fabrication of the superconducting

accelerator cavities, see Zamiri et al. (2006).

8.1. Experimental and Modeling Procedure

To find the main initial texture components, an Orientation Imaging Microscopy
(OIM) was used to measure the initial through thickness texture of the high purity
niobium sheet samples. The material used in this study was high purity
superconducting niobium (RRR) sheets made by Tokyo Denkai Company with 2 mm
thickness. The electropolishing technique used to prepare the sample surface for OIM
study was 90% H2SO4 and 10% HF at 0°C using a tungsten electrode at a potential of

15V. This information was necessary for the modeling part of our investigation.

96

The combined constraints crystal plasticity model was used to investigate the effect of
the surface texture on the plasticity induced surface roughness and dislocation density
in high purity superconducting niobium.

A computational model, as shown in figure 8.1, was used in this work to study the
behavior of different crystal orientations on surface behavior of the niobium during
the plastic deformation. To study only the effects of grain orientations, it is assumed
that all grains have the same geometry. The computational model represents a small
region of the surface down to the mid-thickness of the niobium. The surface and the
subsurface of the niobium sheet are presented by the grains with the same geometry
and all other layers below the mid-thickness of the niobium sheet are assumed to be
isotropic. This computational model is applied to simulate the uniaxial tension along

axis 1.

8.2. Results and Discussion

Texture measurement showed different through thickness texture layers in the
investigated superconducting niobium sheet. As the surface behavior is our interest in
this study, we only analyzed the texture of the surface and a layer close to the surface
(subsurface layer) of the niobium samples. The Orientation Distribution functions
(ODFs) for the surface texture and subsurface texture of the niobium are shown in
figure 8.2. The analysis of these ODFs shows that there is a strong (111)[11-2]
component in the subsurface texture layer of niobium while the surface has a very
strong (100)[001] component, and also (1 10)[1-10] and (111)[11-2] components with

less intensity. Other OIM measurements from other regions of the niobium sheet

97

showed the same results. By using this information, in the first simulation we

assumed that a small region on the

 

Surface grains

   
  
    

Subsurface grains
Isotropic
b)

Figure 8.1. The computational model that was used in these investigations; a) the whole model, b)
aportionofthemodelshowingthesurfaeeandsubsurfaeegraim

98

 

 

 

 

 

 

 

Figure 8.2. ODFs of a) surface, and b) subsun’face textures of high purity niobium

surface of the niobium is made up of four grains with two of those having (001)[100]
orientation and the other two having (110)[1-10] and (111)[11-2] orientations, as
shown in figure 8.3. The corresponding four grains in the subsurface layer were
assumed to have (111)[11-2] orientation as was measured by OIM. This type of
simulation will help us to study the effects of specific orientations, and also the
impact that adjacent grains have on the surface behavior during the plastic

deformation.

(100)[0011

mom-10]
(111)[1-12]

 
  
   

I“ ll l'oplt‘

(111)[1421".

Figure 8.3. A portion of the computational model (figurela) showing the orientation of the 4
surface and corresponding subsurface grains. All other parts of the computational model, except
for these 8 grains, are considered to be isotopic

Figure 8.4 shows the results from the uniaxial tension simulation with this model
along the 1 axis up to 0.4 strain. As can be seen, a surface roughness develops on the
surface of the niobium after the plastic deformation. Grains with (001) orientations on
the surface show a stronger tendency for through thickness and transverse
deformation. (110)[1-10] component shows the least through thickness deformation
among these three orientations but its transverse deformation is something between
(001) and (111) components. (111)[1-12] orientation shows a more moderate
deformation in all three directions. The plastic behavior of a crystal depends on the
number of slip systems that become activated in the crystal during a particular strain
path. The number of activated slip systems depends on the orientation of the crystal.
Therefore, in a given strain path, crystals with different orientations have different
plastic behavior. This different deformation behavior in different grains on the surface
of niobium leads to a surface roughness as seen in figure 8.4. A small scale surface
roughness is also clear within the grain with (001) orientation while it is not present
in other orientations. The surface roughness within a grain can be related to the
presence of different stress states in different regions of a grain and also at higher
strains to the development of a microtexture or an intracrystalline misorientation due
to a possibly nonuniform plastic deformation within a grain. An intracrystalline
misorientation leads to different deformation behaviors in different regions, which
causes a surface roughness to develop within a grain. We evaluated the amount of

intracrystalline misorientation by calculating the amount of rotation of normal to the

101

surface of the grain, (hkl), in different regions of the surface of a grain after the

plastic deformation.

   
    

(ll()l|1.l(l|
(001)[100]

(lllllll—ZI

a)

 

b)

Figure 8.4. a) Surface roughness developed after plastic deformation under uniaxial temion up
to 40% strain, b) the transverse plastic deformation in different grains; (111) grain (black) shows
the least amoumt of transverse defomration while (001) grains (dark gray in the middle) show the

highest

Figure 8.5 shows that the variation of intracrystalline misorientation is much higher
within the grain with (001)[100] compared to the (111)[1-12] orientation. Therefore,
the grain with {001} orientation shows much more inhomogeneous plastic

deformation on the surface compared to the grain with { 111} type orientation.

102

In a different simulation analysis, two grains with (111)[1-12] and (001)[100]
orientations are considered to be on the surface and the corresponding grains on the

subsurface have (111)[1-12] orientations (figure 8.6). The grains on the surface are

   
 

lllil)

(001)

a)

0.20

 

O_
...
G)

   

0.121.

0.08 -

Misorientation angle (rad)

.0
E

 

 

 

 

 

0 40 so 120 150 200
Dlstance across the grain (um)

b)

Figure 8.5. The distribution of the misorientation in different surface grains, a) the contour plot
(blue=0o rad and red =0.25° rad ), b) the magnitude of the misorientation across the grains

103

(ll|)[ll-2]

((ltllllllltu

  
  
 

tllllill~2|

Isotropic

Figure 8.6. A portion of the computational model (figure 8.1a) showing the surface grails
orientations and subsurface grains orientations. All other parts of the model, except for
these 4 grains, are considered to be isotropic

placed a distance apart from each other in the symmetric loading positions in order to
remove the effect of the adjacent grains on the deformation of the investigating grain.
Figure 8.7 shows the through thickness variation of the strain from one side of the
grain to the other side after the uniaxial tension up to 0.4 strain in direction 1. The
through thickness strain in all regions of the grain with cube or (001) orientation is
much higher than the grain with (111) orientation. Within the grain with (111)
orientation, strain is minimum in the center of the grain and it increases as it
approaches to the grain boundary. In other words, grains with (111) orientation have
more resistance to through thickness deformation compared to the surrounding

isotropic medium. The strain distribution is completely inverse in the grain with (001)

orientation in which the through thickness strain is maximum in the center and it

decreases as it reaches to the grain

 

    

—e—(001)[1001
- —e—(111)[11-2]
-o.14 —

Through-thickness straln

 

 

 

-o.22 '

 

Distance across the grain (urn)

Figure 8.7. The distribution of the through-thickness strain across the grains from one side to
the other side

boundaries. Based on these observations, it seems that the grain with (111) orientation
tends to develop a dome shape on the surface while the grain with (001) orientation
more likely develops a shape like a cavity on the surface. The effect of the different
grain orientations on the surface roughness can be also investigated as discussed
above.

Li et al. (2006) showed that dislocation density also has a significant impact on EWF.
It is well known that the dislocation density is a function of the rate of shear strain

[Dieter (1976)] and can be represented by:

V = ¢apabal7a (8.1)

105

a

where ¢ is a correction factor taking into account the effect of the activation of the
. . . a' . . . . b“ .

different independent slip systems, ,0 rs the dislocatron densrty, lS Burgers

—a
vector, and V is the average dislocation velocity on slip system a.

The rate of plastic work per unit of volume in a single crystal is expressed as:
N
. P _ e . a
W " Z 7“ 1'y (8.2)
a=1
a
The rate of critical resolved shear stress, 7 , is also a function of the slip rate as:

N
Til = [aha/3'70 (8.3)

Comparing equations (8.1), (8.2), and (8.3), it could be concluded that any increase in
plastic work can lead to an increase in dislocation density, without considering the

dynamics effects.

Replacing for 7a in equation (8.4) from equation (5.12), one can easily see that xi

in equation (5.12) is basically the rate of the plastic work per unit of volume in a
crystal. Therefore, we use the accumulated/1 as a measure of plastic work in the
different regions of grains to explore the effect of the crystal orientation on the
dislocation density and therefore on the EWF of niobium. The same simulation
procedure as explained in figure 8.3 was used again in this investigation. Figure 8.8
shows the contour of the plastic work and therefore the dislocation density from one
side to the other side of the grains with different orientations on the surface of

niobium. As is clear, (001) orientation shows the highest amount of the plastic work

106

while it is much lower in the grains with (110) and (l l l) orientations. Also it is clear
that the plastic deformation and therefore the dislocation density are not uniform in
different regions of the grains. For example, for the grain with (100) orientation the
plastic deformation is at the maximum near the center of the grain while it is at its
minimum level near the grain boundaries. Therefore, with plastic deformation, the
EWF can vary locally from one grain to another due to the presence of different

dislocation density in different grains with different orientations.

   
 

I|l()l

1001)

 

 

 

 

 

 

 

3)
100_
: —e—100
80; +110
A ” +111
E _
E 60 ~
1! .
3 _
.2 _
5 4°-
1:
203
o - all u lllllll‘i I] III] I Ilrutu
o 5 10 15 20 25 30

Distance across the grain (um)

b)
Figure 8.8. a) The contour plot of the plastic work inside the grains (blue=31 MPa and red=107
MPa), b) The distribution of the plastic work across graim with difl’erent orientations from one
side to the other side

107

The EWF can be also different in different regions of a grain due to an
inhomogeneous distribution of the dislocation density inside of a grain.

A random texture on the surface of the niobium leads to a surface roughness as
discussed above, therefore, it decreases the electron work function of the niobium
The electron work function of the material can be a very important parameter for
some applications. For example, a low electron work function leads to the electron
emission or the so-called field emission in the inner surface of the niobium
superconducting accelerator cavities, which decreases the performance of the cavities.
According to Leblanc et al. (1974) the electron work function of pure niobium highly
depends on the crystal orientation. Based on their measurements, the highest value for
the electron work function belongs to {110} orientations while {001} orientations
have the lowest value. As shown by the above simulations, with plastic deformation,
{001} orientations develop a rough surface and also a higher dislocation density
compared to the other investigated orientations. Based on these observation one can
conclude that a strong surface texture with high ratio of {110}/{001} components
would be desirable to have for a high value of electron work function in high purity
niobium.

Phonon emission is also an important phenomenon that helps the heat transfer from
the surface of a superconducting material to the coolant. Surface roughness is an
important parameter that helps the phonon emission from the surface of
superconducting niobium. Therefore, a more random surface texture containing some
intense {001} components is desirable for phonon emission in the deformed

superconducting niobium.

108

These types of micro-scale simulations can be easily used to design surface texture
for optimal texture-dependent physical properties of the surface. There are also
several other microstructural parameters that may affect the EWF and phonon
emission which were not investigated in this work For example, nonlinear strain
path, grain boundary, grain geometry, and size effects are important parameters that

can be further investigated.

109

CHAPTER 9

CONCLUSION AND FUTURE WORKS

9.1 Concluding Remarks

The critical issues in a crystal plasticity problem, especially in the rate-independent
case, have been to determine active slip systems, the amount of shear on active slip
systems, and the non-uniqueness due to multiplicity of slip systems. Several
mathematical and computational techniques have already been developed to address
these issues. However, due to persistent problems associated with calculation
robustness in the current crystal plasticity models, development of more stable and
efficient crystal plasticity models are still of substantial technical interest.

Three different crystal plasticity models were developed in this thesis work; Modified
Gambin crystal plasticity model, dual mixed crystal plasticity models, and combined
constraints crystal plasticity model.

Modified Gambin, which is a crystal plasticity model based on a power—law yield
surface for single crystals, was developed to decrease the degree of non-linearity and
to increase the flexibility of the previous single crystal power-law yield surfaces
developed by Gambin.

A novel mathematical technique, the so called “Dual Mixed Method”, was developed
for solving constrained optimization problem with many constraints. This method can
be used to bridge computational mathematics and other fields such as solid state

physics, materials science, etc.

110

Up to now, several empirical models for dislocation dynamics have been developed.
The dual mixed method was used to bridge the dislocation theory and the
computational mechanics to develop novel temperature and microstructure sensitive
crystal plasticity models, the so called “dual mixed crystal plasticity models”. Their
evaluation with several simulations showed that they are very CPU-efficient and fast,
compared with other crystal plasticity models. Also because they have been
developed using dislocation theory, they can be used to model the effect of
microstructure on deformation at different temperatures. Furthermore, they can be
used to model impact, room temperature creep, dynamic loading effects, micro-scale
viscous deformation, deformation in semi-crystalline polymers, etc.

A crystal plasticity model is defined as a constrained optimization problem with many
constraints. A single yield surface for a single crystal was found by combining these
constraints using a mathematical technique. A novel crystal plasticity model, the so
called “combined constraints crystal plasticity model”, then was developed using this
single crystal yield function. This crystal plasticity model was able to show a
relationship between SFE and deformation of crystalline materials. Also this crystal
plasticity model was used to develop a novel model for dislocation dynamics in
crystalline materials.

All the above models were implemented into ABAQUS as user materials to simulate
the uniaxial tension, and tube hydroforming of 6061-T4 aluminum alloy. Instead of
hyperelasticity formulation, which is typically used in the formulation of crystal
plasticity problems, an explicit updating procedure based on the consistency

condition during plastic flow in the current configuration was used to improve the

111

execution speed and stability of the crystal plasticity models. The finite element
simulation results showed that the proposed crystal plasticity models were much more
efficient and faster than previous crystal plasticity models such as SVD multi-yield
surface crystal plasticity and Gambin’s model. The combined constraints crystal
plasticity model was further used to model the effect of surface texture on surface
roughness and dislocation density distribution in crystalline materials. The model was
managed to successfully prediction the microtexture evolution, inhomogeneous
plastic deformation, and surface roughness as already reported in literature in

different surface grain of niobium as an example of its application to bee materials.

9.2 Recommendation For Future Works

Although, all finite element simulations performed in this thesis were at room
temperature, the dual mixed crystal plasticity models are temperature sensitive and
can take into account the effects of temperature on plastic deformation. In contrast
with other works in which only the strain hardening is considered to be temperature
sensitive, the dual mixed model is a temperature sensitive crystal plasticity model. It
is therefore proposed that further research be conducted where the dual mixed crystal
plasticity model is applied to high temperature applications.

Further studies can be conducted where the dual mixed method can be used for the
simulation of impact problems, room temperature creep, grain boundary slip transfer,

micro-scale viscous deformation, Bauschinger effect, semi-crystalline polymers, etc.

112

The combined constraints crystal plasticity model with minor modifications will have
the ability to model the cases when there is dynamic loading such as fatigue. A good
research can be conducted on this topic.

Modeling based on dislocation dynamic simulation is an approach that is widely used
to model the dislocation creation and dislocation interaction with second phases,
inclusions, grain boundary and other dislocations. A novel dislocation dynamics
model (equation 5-14) was developed in this work based on the combined constraints
method. This model can relate the velocity of dislocations to the SFE and dislocation
density on other slip systems. There are several unknown parameters in this model.
Further researches on different crystalline materials at different temperature should be
conducted to find these unknown parameters. Then this model can be used to deve10p
a new dislocation dynamics code to model materials at smaller length scale.

The crystal plasticity models proposed above are based on Taylor type
homogenization method and local based continuum theory and therefore, when length
scale or size effects are of interest they cannot be used. Examples of such cases in
which size effect is important are modeling of nanocrystalline materials and thin films
in which grain size must be taken into account. Nanocrystalline materials are novel
materials widely used in nanoelectronic devices and biotechnology. The above crystal
plasticity models, especially the dual mixed crystal plasticity models, can be further
improved by implementing the idea of non-local continuum theory to take into
account the size effect in modeling.

Another important future research to consider would be damage in materials. Further

work can be conducted to implement damage in the proposed crystal plasticity

113

models, especially in the dual mixed crystal plasticity models. Damage can be
implemented in these crystal plasticity models by using several micro-scale damage
parameters which evolve with deformation or time. These micro-scale parameters can
be a function of temperature and microstructural parameters. These crystal plasticity

models can then be used to obtain macro-scale damage parameters for the material.

114

APPENDIX A

Three Euler angles that are defined based on three different conventions, i.e. Bunge,
kocks and Roe, can define the orientation of any crystal with respect to the material
frame, see Bunge (1982) and Kocks, Tome and Wenk (1998). Based on the Bunge

system, the orientation of a crystal can be defined by the three Euler angles ¢,,<I>,¢2

as

cosgoospz—sinqsinpzoosCD singcospz+costqsin¢200sd> singsind)

Geisha): -oosqsin¢,—sinqooso2cos<1> —sin¢isin¢,+oosqqoos¢,cos<l> oosagsincb
singsincb —oosqtsin<l> oos<I>

(A-l)

Kocks defined the crystal orientation matrix by the three Euler angles ¢,‘P, G) as

—sin‘Psin¢—cos‘l’cos ¢cos® cos‘l’sin¢-sin‘l’cos ¢cos® cos¢sin®

Q(¢,‘l’,®) = sin ‘Pcos¢—cos‘l’sin ¢cos® —cos‘Pcos¢-sin‘l’sin ¢cos® sin ¢sin®
cos‘l’sin 9 sin ‘l’sinO cosO

(A-2)

And, based on the Roe definition the crystal orientation matrix is

—sinwsin ¢+ cos wcos¢cosa coswsin¢+sin wcos¢cost9 -cos¢sint9
Q(t1/, ¢, 6) = — sin wcos ¢ — cos visin ¢cos 6 cos t/Icos ¢ - sin i/Isin ¢cos 6 sin ¢sin 6
cos resin 6 sin visin 0 cos 6
(A-3)

115

REFERENCE

ABAQUS Manual, Habbit, Karlsoon & Sorensen, INC., Providence, RI, Version 6.3

Aifantis KB, Soer W.A., DE Hosson J. Th. M., Willis J .R., (2006) Iterface within
gradient plasticity: theory and experiment, Acta Materialia, 54, 5077-5085.

Anand L., and Kothari M., (1996) A computational procedure for rate independent
Crystal plasticity. J .Mech. Phys. Solids, 44 , 525-557.

Arminjon M. (1991) A regular form of the Schmid law. Application to the ambiguity
problem, Textures and Microstructures, Vol 14-18, pp. 1 121.

Ashcroft N.W., Merrnin N. D., (1976) Solid State Physics, Holt, Rinehart, and
Winston, New York.

Asaro RJ. (1970), Geometric effects in the inhomogenous deformation of ductile
single crystals, Acta Metal]. 27, 445-453.

Asaro R]. (1983), Micromechanics of crystals and polycrystals. Adv. Appl.Mech. 23,
1-1 15.

Asaro R]. (1983), Crystal plasticity, ASME J. Appl. Mech. 50, 921-934.

Asaro R.J., Rice J .R. (1977), Strain localization in ductile single crystals, J. Mech.
Phys. Solids 25, 309-338.

Asaro R.J., Needleman A. (1985) Texture development and strain hardening in rate
dependent polycrystals. Acta Metallurgica, 33, 923-953.

Ashmawi W.A., Zikry M.A., (2000) Effects of grain boundaries and dislocation
density evolution on large deformation modes in fee crystalline materials, Journal
of Computer Aided Materials Design, 7, 55-62.

Baars D., Jianh H., Bieler T.R., Zamiri A., Pourboghrat F., C. Compton, (2008)
ICOTOM Symposium, Pittsburgh.

Bassani J .L. (1993), plastic flow of crystals. Adv. Appl. Mech. 30, 191-258.

Bassani J .L., Wu T.Y. (1991), Latent hardening in single crystals 11. Analytical
characterization and predictions, Proc. Royal Soc. London, A 435, 21-41.

Bhushan B., Goldade A.V., (2000) Appl. Surf. Sci. 157, 373.

Blish II R.C., Vreeland T., (1969) J. Appl. Phys, 40, 884.

116

Bunge H.-J. (1982) Texture analysis in materials science, Butterworths.

Busso, E., Cailletaud, G. (2005) On the selection of active slip systems in crystal
plasticity. International Journal of Plasticity 21 (12), 2212—223 1.

Chang YW, Asaro RJ. (1981) An experimental study of shear localization in
aluminium- copper single crystals. Acta Metallurgica, 29, 241.

Chin, G. Y. and Mammel (1969), Generalization and equivalence of the minimum
work (Taylor) and maximum work (BishopHill) principles for crystal plasticity.
Trans. Metall. Sot. AZME 245, 1211-1214.

Chaudhuri A.R., Pate J.R., Rubin LG, (1962), J. App. Phys., 33, 2736.

Cottrell A.H., (1958), Trans. Metall. Soc, AIME, 212, 192-203.

Cuitino AM, Ortiz M. (1992) Computational modeling of single crystals. Modeling
and Simulation in Materials Science and Engineering, 1, 225-263.

Danieulat M. and Piot D. (1996) A method of generating analytical yield surfaces of
crystalline materials, Int. J. Plasticity, 12, 575-610.

DeVechio D., Bhushan B., (1998) Rev. scient. Instrum, 69, 3618.
Dieter GE, (1976) Mechanical Metallurgy, New York (NY), McGraw-Hill

DiMaggio F.L., Sandler LS. (1971) Material Models for Granular Soils, Journal of
Engineering Mechanics 97, 935-950.

Erickson J.S., (1962) J. App. Phys., 33, 2499.
Evers L.P., Brekelrnans W.A.M., Geers M.G.D, (2004) Scale dependent crystal
plasticity frame work with dislocation density grain boundary effects, International

Journal of Solids and Structures, 41 (18-19), 5209-5230

Gambin W., Barlat F., (1997) Modeling of deformation texture development based on
rate independent crystal plasticity, Int. J. Plasticity, 13, 75-85.

Gambin W. (1992) Refined analysis of elastic-plastic crystals, International Journal
for Solids and Structures, 29, 2013.

Gambin W., (1991) Plasticity of crystals with interacting slip systems, Engineering
Transaction, 39, 3-4, 303.

117

Gilman J .J ., Johnston W.G., (1957) Dislocation and mechanical properties of crystals,
Wiley, New York, p.116

Gilman J .J , (1969) Micromechanics of flow in solids, McGraw-Hill, New York,
p. 195

Gorman J.A., Wood BS, and Vreeland T., J. Appl. Phys., 40 (1969) 833
Greenman W.F., Vreeland T., Jr., and Wood D.S., J. AppL Phys, 38(1967) 3595
Gutmanas E.Y., Nadgomyi E.M., Stepanov A.V., (1963), 5, 743.

Has G. A., Thomas RE, (1977)]. Appl. Phys. 478, 86.

Havner KS, (1992), Finite plastic deformation of crystalline solids, Cambridge
university press, Cambridge.

Hahn G.T., Acta Met. 10(1962) 727.

Hill R. (1966), Generalized constitutive relations for incremental deformation of
metal crystals, J. Mech. Phys. Solids, 14, 95-102.

Hill R., Rice J .R., (1972), Constitutive analysis of elastic-plastic crystals at arbitrary
strain, J. Mech. Phys. Solids, 20, 401-413.

Huang Y. (1991) A User-Material Subroutine Incorporating Single Crystal Plasticity
in the ABAQUS Finite Element Program, Mech Report 178, Harvard University

Hutchison J .W., (1976) Elastic-plastic behavior of polycrystalline metals and
composites, Proceeding of the Royal Scociety London, Series A, 319, 247-272

Johnston W.G. and Gilman M., J. Appl. Phys., 33 (1959) 129
Kalidindi, S. R., Bronkhorst, C. A. and Anand, L. (1992) Crystallographic texture

evolution during bulk deformation processing of fee metals. J. Mech. Phys. Solids
40, 537-569.

Kim H.-K., Oh S.-I, (2003) Finite element analysis of grain-by-grain deformation by
crystal plasticity with couple stress, Int. J. plasticity, 19, 1245-1270.

Kocks U.F., Tome C.N. and Wenk HR. (1998) Texture and Anisotropy, Cambridge
University Press.

Koiter W.T., (1953), Stress-Strain relations, uniqueness and variational theorems for

elastic-plastic materials with a singular yield surface, Quarterly of applied
mathematics 11, 350-354.

118

Knoekaert R., Chastel Y., Massoni E. (2000) Rate-independent crystalline and
polycrystalline plasticity, application to fcc materials. International Journal of
Plasticity 16, 179—198.

Kreisselmeier G., and Steinhauser R., (1979), Systematic control design by
optimizing a vector Performance index, proceedings of IFAC symposium on
computer aided design of control systems, Zurich, Switzerland, pp. 113-117.

Kuchnicki S.N., Cuitin A.M., Radovitzky RA. (2006) E.cient and robust constitutive
integrators for single-crystal plasticity modeling. International Journal of Plasticity
22, 1988—2011

Leblanc R.P., Vanbrugghe B.C., Girouard RE, (1974) Can. J. Phys. 52, 1589.

Lequeu, Ph., Gilormini, P., Montheillet, F., Bacroix, B. and Jonas, J.J. (1987) Yield
Surfaces for Textured Polycrystals --1. Crystallographic Approach. Acta Metall.,
35(5), 439.

Lequeu, Ph., Gilormini, P., Montheillet, F., Bacroix, B. and Jonas, 1]., (1987) Yield
Surfaces for Textured Polycrystals--II. Analytical Approach. Acta Metall., 35(5),
1 159.

Levitin V.V., Pravda M.I., Serpetzky B.A., (1994) Solid State Commun 92, 973.

Levitin V.V., Grin O.L., Yatsenko V.K, (2001) Vacuum 63, 367.

Li W., Li D.Y., (2002) Mater Sci Technol 18, 1057.

Li W., Li D.Y., (2004) Philos Mag A 84, 3717.

Li W., Cai M., Wang Y., Yu S., (2006) Scripta Mater. 54,921.

Ling X., Horstemeyer M., Potirniche, G. (2005) On the numerical implementation of
3D rate dependent single crystal plasticity formulations. International Journal for
Numerical Methods in Engineering 63 (4), 548—568.

Loret, B., and J .H. Prevost, (1986) Accurate Numerical solutions for Drucker-Prager
Elastic-Plastic Models, Computer Methods in Applied Mechanics and Engineering
54, 259-277.

Lublinger J. (1984) A maximum-dissipation principle in generalized plasticity, Acta
Mechanica 52, 225-237

Lublinger J. (1986) Normality rules in large-deforrnation plasticity, Mechanics of
Materials 5, 29-34

119

Ma A., Roters F., Raabe D., (2006) On the consideration of interface between
dislocations and grain boundaries in crystal plasticity finite element modeling-
theory, experiment, and simulations, Acta Materialia, 54, 2181-2194

Mamor M., Finkman E., Meyer F., (1995) Mater Res Soc Syrup 356, 149.

Mandel J ., (1964) Contribution theorique a l’Etude de l’ ecrouissage et des lois de 1

ecoulement plastique, Proceedings of the 11‘h International Congress on Applied
Mechanics, 502-509.

Mandel J. (1974) Thermodynamics and plasticity. Proc. Int. Symp. Foundations of
Continuum Thermodynamics (ed. Delgado Domingos,J. J ., Nina M. N. R. and
Whitelaw J. H.), pp.283-3 Il. McMillan, London.

Mathur KK, Dawson PR. (1989) On modeling the development of crystallographic
texture in bulk forming processes. Int. J. Plasticity, 5, 67-94.

McGinty, R., McDowell, D. (2006) A semi-implicit integration scheme for rate
independent finite crystal plasticity. International Journal of Plasticity 22 (6), 996—
1025.

Meyers M.A., (1984), Mechanical metallurgy; Principles and applications, Prentice
Hall inc., New Jersey, p.304

Miehe C., Schroder J, and Schotte J (1999) Computational homogenization analysis
in finite plasticity simulation of texture development in polycrystalline materials,
Computer Methods in Applied Mechanics and Engineering, 171, 387-418.

Miehe C., and Schroder J (2001) A comparative study of stress update algorithms for
rate-independent and rate-dependent crystal plasticity, Int. J. Numer. Meth. Engng,
50, 273-298.

Montheillet F., Gilormini P. and Jonas J .J . (1985) Relation between Axial Stresses
and Texture Development during Torsion Testing: 3 Simplified Theory. Acta
Metall., 33, 4, 705.

Naghdi, P.M.,(1960), Stress-Strain relatin in plasticity and thermoplasticity, in
proceedings of the 2th symposium on naval structural mechanics, pergamon press,

London.

Needleman A., Asaro R.J., Lemonds J ., Peirce D., (1985). Finite element analysis of
crystalline solids, Comput. Meth. Appl. Mech.Eng., 52, 689—708.

120

Nemat-Nasser S., Okinaka T., Ni L. (1996) A physically-based constitutive model for
bcc crystals with application to polycrystal tantalum. Journal of the Mechanics and
Physics of Solids 46 (6), 1009—1038.

Ney H., Labusch R., Hassen P., (1977) Acta Met., 25, 1257

Parameswaran V.R. and Weertman J ., (1971) Met. Trans, 2, 1233

Parameswaran V.R., Urabe N ., and Weertman J., (1972) J. Appl. Phys., 43 2982

Peirce D, Asaro R, Needleman A. (1982) An analysis of nonuniform and localized
deformation in ductile single crystals. Acta Metallurgica, 30, 1087-1119.

Peirce D, Asaro R, Needleman A. (1983) Material rate dependence and localized
deformation in crystalline solids. Acta Metallurgica, 31, 1951-1983.

Peirce D., Shih C.F., Needleman A., (1984), A tangent modulus method for rate
ependent solids, Comput. Struct., 18, 875—887.

Pope D.P., Vreeland T., Jr., and Wood BS, (1967) J. Appl. Phys., 38, 40111

Pourboghrat F. and Barlat F. (2002) texture evolution during hydroforming of
aluminum extruded tubes, NSF-Design conference, Puetro Rico. 1958-1973.

Raphanel J .L., Ravichandran G., Leroy Y.M. (2004) Three-dimensional rate-
dependent crystal plasticity based on Runge—Kutta algorithms for update and

consistent linearization. International Journal of Solids and Structures 41, 5995—
6021.

Resende L., Martin J .B., (1986) Formulation of Drucker-Prager cap model, Journal of
Engineering Mechanics 111(7), 855-865.

Rice J .R., (1971) Inelastic constitutive relations for solids: an internal variables
theory and its application to metal plasticity. J. Mech.Phys. Solids, 19, 433—455.

Rohde R.W., Pitt CH, (1967) J. App. Phys., 38, 876

Romanowski S., (1988) Phys. Status Solidi B 145, 467.

Rousselier, G., Leclercq, S. (2006) A simplified “polycrystalline” model for
viscoplastic and damage finite element analyses. International Journal of Plasticity

22 (4), 685—712.

Rueda H., Slinkrnan J ., Chidambarrao D., (1999) Mater Res Soc Proc 568, 245.

121

Sandler I.S., Rubin D., (1979) An algorithm and a modular subroutine for the cap
model, International Journal of Numerical and Analytical Methods in
Geomechanics 3, 173-186.

Stranger R.W., Mackie W., Swanson L. W., (1973) Surf. Sci. 34, 225.

Schadler H.W., (1964) Acta Met., 12, 861

Schmid E., Boas W., (1935), Plasticity of crystals, Chapman and Hall, London.

Sirno J .C. and Hughes T.J.R. (1998) Computational Inelasticity, Interdisciplinary
Applied Mathematics, Springer, New York.

Sirno J .C. and Vu-Quoc L., (1986) A three dimentional finite strain rod model. Part
II: Computational aspects, Comp. Methods App. Mech. Eng. 58, 79-116.

Sirno J .C., Taylor R.L., (1985), Consistent tangent operators for rate independent
elasto-plasticity. Comput. Meth. Appl. Mech. Eng., 48, 101—118.

Stein DE and Low J.R., (1960) J. App. Phys., 33,129

Stroh A.N., (1957) Adv. Phys., 6, 418

Suzuki T., Ishii T., (1968) Suppl. Trans. Jap. Inst. Met., 9, 687.
Taylor G.I., (1938) Plastic strain in metals, J. Inst. Met. 62, 307-324.

Taylor G. I. (1938b) Analysis of plastic strain in a cubic crystal. Stephen Timoshenko
60'” Anniversary Volume, pp. 218-224. McMillan Co., New York.

Taylor G. I. and Elam, C. F. (1925) The plastic extension and fracture of aluminum
single crystals. Proc. Royal Sot. London A 108,28-51.

Taylor G.I., Elam CF, (1923), The distortion of an aluminum crystal during a tensile
test, Proc. Roy. Soc. London A. 102, 643—667.

Toth L.S., Van Houtte P. and Van Bael A. (1991) Analytical Representation of
Polycrystal Yield Sur-faces. in Boehler, J .P. and Khan, A.S., (eds), Anisotropy and
Localization of Plastic Deformation (Proc. Plasticity '91: 3rd Int. Symp. on
Plasticity and its Current Applications). Elsevier Applied Science, London, pp. 183-
186.

Van Houtte P. (1987) Calculation of the Yield Locus of Textured Polycrystals using
the Taylor and the Related Taylor Theories. Textures and Microstructures, 7(1), 29.

Van Sciever S.W., (1996) Helium Cryogenics, Plenum press.

122

Wu, T.-Y., Bassani, J. L. and Laird, C. (1991) Latent hardening in single crystals I.
Theory and experiments. Proc. Royal Sot. London A 435, l-19.

Zamiri A., Pourboghrat F., Jiang H., Bieler T.R., Barlat F., Brem J., Compton C.,
Grimm T.L., (2006) Mater Sci Eng A 435—436,658.

Zamiri A., Pourboghrat F., (2007) Int. J of Solids and Structures 44, 8627.

Zener c., (1948) The rnicromechanism of fracture, Fracturing of metals, American
Society for Metals, Metals Park, Ohio

Zharin A.L, Rigney D.A., (1998) Tribol. Lett., 4, 205.

Zharin A.L., (2001) Fundamentals of Triblogy and Bridge the Gap Between the
Macro and Micro/Nanoscales (Dordrecht: Kluwer Academic), p. 451.

123