PHYSICS-BASED CRYSTAL PLASTICITY MODELING OF SINGLE CRYSTAL NIOBIUM
By
Tias Maiti

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

ABSTRACT
PHYSICS-BASED CRYSTAL PLASTICITY MODELING OF SINGLE CRYSTAL NIOBIUM
By
Tias Maiti
Crystal plasticity models based on thermally activated dislocation kinetics has been successful
in predicting the deformation behavior of crystalline materials, particularly in face-centered cubic
(fcc) metals. In body-centered cubic (bcc) metals success has been limited owing to ill-defined slip
planes. The flow stress of a bcc metal is strongly dependent on temperature and orientation due
to the non-planar splitting of a/2⟨1 1 1⟩ screw dislocations. As a consequence of this, bcc metals
show two unique deformation characteristics: (a) thermally-activated glide of screw dislocations the motion of screw components with their non-planar core structure at the atomistic level occurs
even at low stress through the nucleation (assisted by thermal activation) and lateral propagation
of dislocation kink pairs; (b) break-down of the Schmid Law, where dislocation slip is driven only
by the resolved shear stress.
Since the split dislocation core has to constrict for a kink pair formation (and propagation),
the non-planarity of bcc screw dislocation cores entails an influence of (shear) stress components
acting on planes other than the primary glide plane on their mobility. Another consequence of the
asymmetric core splitting on the glide plane is a direction-sensitive slip resistance, which is termed
twinning/atwinning sense of shear and should be taken into account when developing constitutive
models.
Modeling thermally-activated flow including the above-mentioned non-Schmid effects in bcc
metals has been the subject of much work, starting in the 1980s and gaining increased interest in
recent times. The majority of these works focus on single crystal deformation of commonly used

metals such as Iron (Fe), Molybdenum (Mo), and Tungsten (W), while very few published studies
address deformation behavior in Niobium (Nb). Most of the work on Nb revolves around fitting
parameters of phenomenological descriptions, which do not capture adequately the macroscopic
multi-stage hardening behavior and evolution of crystallographic texture from a physical point of
view. Therefore, we aim to develop a physics-based crystal plasticity model that can capture these
effects as a function of grain orientations, microstructure parameters, and temperature.
To achieve this goal, first, a new dilatational constitutive model is developed for simulating the
deformation of non-linear geometries (foams or geometries with free surfaces) using the spectral
method. The model has been used to mimic the void-growth behavior of a biaxially loaded plate
with a circular inclusion. The results show that the proposed formulation provides a much better
description of void-like behavior compared to the pure elastic behavior of voids. Using the developed dilatational framework, periodic boundary conditions arising from the spectral solver has
been relaxed to study the tensile deformation behavior of dogbone-shaped Nb single crystals.
Second, a dislocation density-based constitutive model with storage and recovery laws derived
from Discrete Dislocation Dynamics (DDD) is implemented to model multi-stage strain hardening.
The influence of pre-deformed dislocation content, dislocation interaction strengths and mean free
path on stage II hardening is then simulated and compared with in-situ tensile experiments.

I dedicate this thesis to my parents, who have always loved me unconditionally and whose good
examples have taught me to work hard for the things that I aspire to achieve.

iv

ACKNOWLEDGMENTS

First, I would like to express my sincere gratitude to my advisor Dr. Philip Eisenlohr for the continuous support of my Ph.D. study and related research, for his patience, motivation, and immense
knowledge. His guidance helped me in all the time of research and writing of this thesis.
Besides my advisor, I would like to thank the rest of my thesis committee: Prof. Thomas
Bieler, Prof. Donald Morelli, and Prof. Roozbeh Dargazany, for their insightful comments and
encouragement, but also for the hard questions that gave me incentive to widen my research from
various perspectives.
My sincere thanks also go to Prof. Carl Boehlert, and Prof. Martin Crimp, who provided me
an opportunity to access to their laboratory and research facilities. Without the valuable support, it
would not be possible to conduct this research.
I am also very grateful to Dr. Pratheek Shanthraj, Max Planck Institut für Eisenforschung
GmbH, for his scientific advice and many insightful discussions and suggestions. He was my
primary resource for getting my numerical questions answered and was instrumental in helping
me develop constitutive models.
I thank my fellow colleagues Tridip Das, Aritra Chakraborty, Chen Zhang, Dr. Ajith Chakkedath, and Zhuowen Zhao for the stimulating discussions and all the fun we have had in the last four
years.
Last but not the least, I would like to thank my parents for the constant source of inspiration
throughout writing this thesis and my life in general.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

CHAPTER 1 INTRODUCTION AND RESEARCH BACKGROUND . . .
1.1 Superconducting Radio Frequency Cavity . . . . . . . . . . . . . .
1.2 Material for SRF cavity application . . . . . . . . . . . . . . . . . .
1.3 Challenges associated with using Nb as material for SRF application
1.4 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

1
2
4
6
8

CHAPTER 2 PLASTIC DEFORMATION IN BODY CENTERED CUBIC METALS
2.1 Slip systems in BCC crystals . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Orientation dependent strain hardening in BCC single crystals . . . . . . . .
2.3 Non-Schmid effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Kink-pair assisted dislocation motion . . . . . . . . . . . . . . . . . . . . . .
2.5 Temperature dependence of flow stress . . . . . . . . . . . . . . . . . . . . .
2.6 Overview of deformation in Nb single crystals . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

9
10
12
15
18
20
23

CHAPTER 3
3.1
3.2
3.3
3.4

MODELING NON-COMPACT GEOMETRIES USING SPECTRAL METHOD
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Constitutive model for void like regions . . . . . . . . . . . . . . . . . . . . . . . 29
3.1.1 Solution to mechanical boundary value problem . . . . . . . . . . . . . . . 31
Comparison between a dilatational and soft elastic void . . . . . . . . . . . . . . . 32
Uniaxial tension of dogbone-shaped sample of oligocrystalline aluminum . . . . . 35
Selection of parameters for void material description and its associated sensitivity on the simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

CHAPTER 4 CONSTITUTIVE MODELS . . . . . . . . . . . . . . . . . . . . .
4.1 Phenomenological Model . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Dislocation density based model . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Influence of elastic anisotropy on dislocation interaction strengths
4.2.2 Dislocation storage and recovery framework . . . . . . . . . . . .
4.2.3 Dislocation multiplication . . . . . . . . . . . . . . . . . . . . .
4.2.4 Dislocation annihilation . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

41
42
42
47
49
52
55

CHAPTER 5 RESULTS AND DISCUSSION . . . . . . . . . . . . . . . . . . .
5.1 Calibration of Crystal Plasticity constitutive model . . . . . . . . . . . .
5.2 Phenomenological constitutive model results . . . . . . . . . . . . . . . .
5.3 Dislocation density based crystal plasticity model results . . . . . . . . .
5.3.1 Slip activity maps as a function of orientation . . . . . . . . . . .
5.3.2 Influence of slip family on spatially resolved deformation behavior
5.3.3 Model sensitivity to initial dislocation density . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

57
57
60
62
63
64
66

vi

5.3.4
5.3.5

Influence of dislocation interaction coefficients on strain hardening . . . . . 71
Inverse optimization of dislocation mean free path coefficient . . . . . . . . 73

CHAPTER 6 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1 Recommendations for future work . . . . . . . . . . . . . . . . . . . . . . . . . . 78
CHAPTER 7 ADDITIONAL PROJECTS . . . . . . . . . . . . . . . . . . . . . . . .
7.1 Influence of geometric reconstruction accuracy on grain-averaged aggregate mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 Spatial reconstruction method . . . . . . . . . . . . . . . . . . . . . .
7.1.3 Volume statistics of real and artificial grain structures . . . . . . . . . .
7.1.4 Comparison of 3D reconstruction of grain aggregates in presence/absence
of grain volume information . . . . . . . . . . . . . . . . . . . . . . .
7.1.5 Grain averaged lattice strains of reconstructed structures . . . . . . . .
7.1.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Cumulative volume fraction of log-normal grain size distributions . . . . . . .
7.3 Grain size normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4 Generation of synthetic grain structures . . . . . . . . . . . . . . . . . . . . .

. . 80
.
.
.
.

.
.
.
.

80
80
82
83

.
.
.
.
.
.

.
.
.
.
.
.

86
89
94
96
96
97

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

vii

LIST OF TABLES

Table 1.1: Superconducting properties of some promising materials for SRF applications. .

5

Table 3.1: Material parameters; elastic constants Ci j , reference strainrate ε̇0 , stress exponent n, initial and saturation flowstress τ0 and g∞ , hardening parameters h0
and a, and Taylor factor M = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Table 4.1: Material variables for Peirce et al. [1] model; elastic constants Ci j , reference
ξ

shear rate γ̇0 , stress exponent n, initial and saturation slip resistance τ0 and
ξ
g∞ , hardening parameters h0 , qξ β , and a. . . . . . . . . . . . . . . . . . . . . . 43
Table 4.2: Dislocation interaction matrix in BCC lattice. . . . . . . . . . . . . . . . . . . . 46
Table 4.3: Types of possible interaction in BCC lattice [2] . . . . . . . . . . . . . . . . . . 46
Table 4.4: Modified dislocation interaction matrix in BCC lattice taking into account
asymmetric junction configurations [3]. . . . . . . . . . . . . . . . . . . . . . . 48
Table 4.5: Types of possible interaction in BCC lattice with asymmetric junction configurations [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Table 4.6: Material variables for dislocation density based model . . . . . . . . . . . . . . 56
Table 5.1: Comparison of CPFFT computation time between simplified reduced geometry of gauge section and a proper dogbone shaped geometry. Material described by crystal plasticity law is shown in blue and the dilatational buffer
layer is shown in light gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Table 5.2: Constitutive material parameters of pure Nb used for the crystal plasticity simulations of uniaxial tension. The asterisk (*) indicates fitting parameters. . . . . . 63
Table 5.3: Variability in constitutive parameter values resulting from minimizing the deviation in single crystal stress–strain response for eight distinct tensile directions. 75
Table 7.1: Material parameters; elastic constants Cab , reference shear rate γ̇0 , stress exponent n, initial and saturation slip resistance τ0 and g∞ , hardening parameters
h0 , qξ β , and a. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Table 7.2: DREAM.3D pipeline for generating synthetic grain structures. . . . . . . . . . . 97

viii

LIST OF FIGURES

Figure 1.1: A typical accelerating cavity geometry, showing particle beam and fundamental fields of the RF cavity. . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

Figure 1.2: Mechanical fabrication path of an SRF cavity. . . . . . . . . . . . . . . . . . .

7

Figure 1.3: Surface defects in a well prepared SRF cavity cell. . . . . . . . . . . . . . . . .

7

Figure 2.1: Potential slip systems in BCC lattice . . . . . . . . . . . . . . . . . . . . . . . 10
Figure 2.2: Areas of activation of the 110 and 112 planes in the standard triangle [4] . . . . 11
Figure 2.3: Shear stress – shear strain curves as a function of orientation in α − Fe (top)
and Nb (bottom), according to [5, 6]. . . . . . . . . . . . . . . . . . . . . . . . 12
Figure 2.4: Schematic illustration of Taylors angles χ and φ . The loading direction is
marked in the standard triangle by the angles χ and φ [7] . . . . . . . . . . . . 13
Figure 2.5: Influence of the orientation of the load on the shear stress-strain curves with
slip only in {1 1 0} system. (Left) effect of angle φ , (Right) effect of angle χ [7]. 15
Figure 2.6: Evolution of lattice reorientation with plastic deformation depending on the
initial orientation [8, 9, 10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Figure 2.7: Illustration showing core configuration of screw dislocations in a BCC lattice . 17
Figure 2.8: Intrinsic lattice resistance (Peierls-Nabarro Barrier) to slip in BCC materials. . . 19
Figure 2.9: Temperature dependence of flow stress, slip mechanisms (schematics) in
pure BCC metals (Mo) [11]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Figure 2.10: Surface topography of deformed single crystal Nb [12] . . . . . . . . . . . . . 24
Figure 2.11: In situ uniaxial tensile test results reported by Baars [12]. Left : Inverse pole
figure showing rotation of tensile axis as a function of orientation for as-cast
samples. Center: axial stress–strain curves for the as-cast samples with same
orientations in the left inverse pole figure map. Right : axial stress–strain
curves for the annealed samples with the same orientations. . . . . . . . . . . . 24
Figure 3.1: Multiplicative decomposition of deformation gradient [13] . . . . . . . . . . . 27
Figure 3.2: Schematic diagram showing calculation of S as a function of F. . . . . . . . . . 28
ix

Figure 3.3: Schematics of VON M ISES plasticity model in principle stress space. The
state of stress at a material point is represented by the vector s (black), which
can be additively decomposed into the pressure vector p (red) and the deviatoric stress vector s′ (green). The light green surface represents the yield
surface of the material point. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Figure 3.4: Comparison of mechanical field maps of a (fully periodic) two-dimensional plastic
plate containing a circular inclusion of area fraction 0.8 % and loaded biaxially to
an average deformation of ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.4. Void volume is described by an
isotropic elastic material model (left, “elastic”) or by an isotropic phenomenological
plasticity model with dilatational capabilities (right, “dilatational”). The stiffness of
the inclusion in the elastic case is set to be 1000× softer than the plate. . . . . . . . . 34

Figure 3.5: Influence of external boundary conditions on in-plane VON M ISES strain fields and
lattice reorientation. Simulated geometries with fully periodic gauge section of dogbone (first row), quasi-periodic, i.e. periodic along the thickness direction (second
row), and free surfaces normal to x and z (third and last row) compared to experimental results measured by DIC and EBSD (fourth row, from [14]). Gray (green)
semi-transparent volume in the simulated domain column uses the proposed dilatational (isotropically elastic) material model; regions of constant color reflect individual grains. Inverse pole figure (IPF, last column) coloring of lab direction z is
mapped on the undeformed configuration, except for the measured result (fourth
row). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Figure 3.6: Computational cost (top right) and void representation quality (bottom) of the
elastic and the dilatational material model for the case of biaxial expansion of
a hole in a plastic plate (top left) up to ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.2. Strain evolution
of median interface traction (lower left, collected from voxels immediately
inside and outside of interface, i.e. blue and red ring at top left) indicates
the void is reaching the flow stress level (≈ M g) of the surrounding plate in
case of the elastic and strongest dilatational material (plastic strength ratio
of 100 ); softer dilatational cases (darker red) saturate much earlier at flow
stress levels far below the plate flow stress (plastic strength ratio of 10−1 and
10−2 ). Median interface tractions collected at ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.2 (bottom
right) confirm marginal influence of elastic stiffness on void representation
quality for the dilatational material. . . . . . . . . . . . . . . . . . . . . . . . 40
Figure 4.1: A schematic figure showing dislocation reaction between segments on active
slip-system ξ (blue) and a forest system β (red). . . . . . . . . . . . . . . . . . 52
Figure 5.1: Left : Optimization function minimizing the area between two curves. The gray
line represents the standard curve, which the dashed black line represents the trial
curve. Right : Flow diagram of the optimization process used by Chakraborty and
Eisenlohr [15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

x

Figure 5.2: Comparison of stress components of simplified reduced geometry and gauge
length section of dogbone geometry in terms of kernel density estimation
maps. The top row shows the correlation of the normal stress components
and the bottom row shows the correlation of shear stress components. . . . . . . 60
Figure 5.3: Comparison of experimentally observed axial stress-curves with the ones
predicted using phenomenological crystal plasticity model as function of orientation. Solid lines indicate the measured stress-strain curves while the predicted ones are shown in dashed lines. . . . . . . . . . . . . . . . . . . . . . . 61
Figure 5.4: Adjustable phenomenological model parameters (initial CRSS τ0 saturation
CRSS g∞ and hardening slope h0 ) for pure Niobium obtained from stressstrain curve optimization of 8 single crystal orientations. The crosses represent a parameter for a particular orientation. The mean of the parameters for
all orientations is represented by a solid dot. . . . . . . . . . . . . . . . . . . . 62
Figure 5.5: Active slip family based on net accumulated shear Γacc as a function of orientation. The contribution of a slip family for each of the orientations is represented by a color varying from blue for pure {1 1 2} to red for pure {1 1 0}
slip. Left: active slip family maps at ε ∼ 0.002. Right: active slip family
maps at ε ∼ 0.050. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
Figure 5.6: Local ⟨0 1 0⟩ lattice reorientations maps of the deformed dogbone samples.
Rows in descending order represent the simulation results obtained assuming
{1 1 0}, {1 1 0}-{1 1 2} and {1 1 2} slip family respectively. . . . . . . . . . . . 65
Figure 5.7: Local von Mises strain εvM maps of the deformed dogbone samples. Rows in
descending order represent the simulation results obtained assuming {1 1 0},
{1 1 0}-{1 1 2} and {1 1 2} slip family respectively. . . . . . . . . . . . . . . . . 66
Figure 5.8: Local von Mises stress σvM maps of the deformed dogbone samples. Rows in
descending order represent the simulation results obtained assuming {1 1 0},
{1 1 0}-{1 1 2} and {1 1 2} slip family respectively. . . . . . . . . . . . . . . . . 67
Figure 5.9: Influence of pre-deformed dislocation content on stress-strain response (lower
left) simulated by dislocation density based model. The top row shows several instances of inhomogeneous distribution of dislocation content on individual slip systems for a chosen log-normal dislocation statistics. Each slip
system is denoted by a color key given in the lower right legend. . . . . . . . . 68

xi

Figure 5.10: Influence of pre-deformed dislocation content on lattice reorientation for
three different orientations simulated by dislocation density based model.
The left column shows the probability distribution of pre-deformed dislocation content with ρ m = 1 × 1013 m−2 and different variances. The corresponding lattice reorientation trajectories for three orientations are shown in
the right three columns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
Figure 5.11: Results predicted by the dislocation density based model with isotropic α −
Fe interaction coefficients. Left: ⟨0 1 0⟩ inverse pole figure for the 8 orientations showing the tensile axis re-orientation as a function of strain (εmax =
0.30. The initial tensile axis orientation is denoted by the blue dot. The final orientation after 30 % strain is indicated by the red dot. The black dots
in-between represents the trajectory of reorientation. Right: Corresponding
axial stress-strain (σ - ε ) predicted by the crystal plasticity model. . . . . . . . 71
Figure 5.12: Influence of elastic anisotropy on strain hardening behavior for a particular
orientation. The interaction coefficients for Fe assuming isotropic elasticity
was calculated by Queyreau et al. [2]. The dislocation interaction coefficients
assuming anisotropic elastic for Nb was taken from the works of Madec and
Kubin [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Figure 5.13: Comparison of experimentally observed (solid) and predicted (dashed) uniaxial stress–strain response for different orientations. The coefficients for
mean free path Khkl and dynamic recovery yhkl in the model by Devincre
et al. [16] were estimated using inverse optimization of the measured stress–
strain curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Figure 5.14: Variability in predicted single–crystal response under unidirectional tension
due to different distributions among the slip systems for fixed total initial
dislocation content. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
Figure 7.1: Polycrystalline aggregates (a) measured by various 3D techniques [17, 18,
19, 20] and associated grain statistics (blue curves in (b) and (c)). Sequence
of four gray curves represents log-normal grain volume distributions of different variance (Eq. (7.9)); the black curve corresponds to closest fit at variance σdistr = 1.2. Dotted line represents average grain volume. . . . . . . . . . 83
Figure 7.2: Synthetic polycrystals (a) generated by DREAM.3D [21, 22] using cube octahedra (top row) or super ellipsoids (bottom row) as grain shapes. Orange
and brown curves in (b) and (c) show the associated grain statistics for both
selections of grain shapes. Black line (optimal fit to experimental data) is
same as in Fig. 7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

xii

Figure 7.3: Ratio of largest to smallest principal moment of inertia (‘shape factor’) for
all cube octahedral (orange) and super ellipsoidal (brown) grains. Solid
lines correspond to synthetic (reference) structures shown in Fig. 7.2(a). The
fainter curves result after reconstructing the reference structures using either
grain volumes (dashed) or centroid distances (dotted) as a basis for Laguerre
tessellation weigths. In general, cube octahedral grain structures are closer
to being equiaxed than those based on super ellipsoids. . . . . . . . . . . . . . 86
Figure 7.4:

Joint probability of correct grain volume (normalized by average grain size) and
relative deviation in the reconstructed grain volume Vtess /Vref or grain shape factor κtess /κref . Grain reconstruction uses Laguerre tessellation with grain centroids
as seed positions. Seed point weights are derived from known grain volumes (upper half) or estimated from the distance to the nearest grain centroid (lower half).
Left and right panels correspond to cube hexahedral and super ellipsoidal grain
shapes, respectively, that were used in generating the underlying synthetic data sets
(Fig. 7.2(a)). Red maps correspond to direct use of grain centroids as seed positions, while green maps illustrate changes due to a seed position adjustment [23]
that considers the new centroid location of each resulting tessellated grain. The
comparisons of shape factor deviations (yellow maps) utilized only those latter reconstructed structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Figure 7.5: Deviations of normal, shear, and deviatoric normal components of grainaveraged Cauchy stress σtess in reconstructed grain structures relative to corresponding stress components σref resulting in the (cube octahedron) reference grain structures. Tessellations in upper panel employ known grain
volume information (requiv ) for Laguerre weights, while in lower panel the
nearest centroid distances (dnn ) are used instead. The seed positions were
adjusted [23] in both tessellation alternatives. . . . . . . . . . . . . . . . . . . 90
Figure 7.6: Joint probability of volume Vref of synthetic grains and the distance dnn from
each grain centroid to the nearest neighboring one. Quantities are normalized by the average grain volume ⟨Vref ⟩ and its equivalent radius requiv =

(3⟨Vref ⟩/4π )1/3 , respectively. Cumulative probability is scaled to unity within
each (vertical) bin of normalized grain volume. Centroid distance progressively deviates from expected relation (dashed line) with decreasing volume
of grain, resulting in a weak correlation between the distance to the nearest
grain centroid and the volume of a grain. . . . . . . . . . . . . . . . . . . . . . 93

Figure 7.7: Joint probability of stress magnitudes σref in reference structure and relative
deviations σtess /σref combined from both tessellated structures. Cumulative
probability is scaled to unity within each bin of stress magnitudes. Solid lines
exemplify relative deviations in stress due to absolute differences of 5 MPa
and 10 MPa. With increase in magnitude of stress, the relative deviations in
stress due to difference in grain morphologies diminishes. . . . . . . . . . . . 95

xiii

CHAPTER 1
INTRODUCTION AND RESEARCH BACKGROUND

Advancements in particle physics rely heavily on advances in particle accelerators. Originally,
particle accelerators were developed to do fundamental research in sub-atomic particle physics:
the study of the elemental forces and building blocks of matter. Today, many diverse real life
applications make use of particle accelerators, such as in life sciences, chemistry, medical diagnosis
and treatment, material characterization, and luggage interrogation for homeland security defense
applications. A few of the ongoing projects based on particle accelerator technology involves
x-ray laser oscillators, continuous wave (CW) and pulsed free electron lasers, energy recovery
linac (ERL) based light sources, electrons and ions colliders, nuclear waste transmutation, and
accelerator driven systems (ADS) for medical isotope production. Researchers are working on
developing a new type of nuclear reactors based on particle accelerators, which could solve the
nuclear waste disposal problem.
Many parameters influence the performance of the accelerator, and the physics they make accessible. A critical parameter among them is the energy of the particle beams. Particle beams of
higher energy allow for finer measurement of the fundamental properties of the particles constituting them. The technology for accelerating particles also depends on the type of particle being
accelerated. Heavy particle (proton) accelerators are mostly circular machines, with their technological challenges lying in the area of high field magnet development. The Tevatron, at Fermilab,
and the Large Hadron Collider (LHC) at Conseil Européen pour la Recherche Nucléaire (CERN)
are examples of these type of machines. Electron-positron colliders have been almost exclusively
circular machines up to the present. The synchrotron radiation which results from accelerating
1

charged particles in a circle makes circular machines impractical for further increase in energy.
For this reason, the next generation of electron-positron colliders will have to be linear. A linear
collider with the center of mass energy in the 1 TeV range would be a machine complementary to
the previously mentioned LHC. Several proposals have been put forward for the construction of
such a device, [24, 25] one of which would be a Super conducting Radio Frequency (SRF) based
machine. These SRF based machines are composed of smaller units stacked in series known as
SRF cavities or cells.

1.1 Superconducting Radio Frequency Cavity
An SRF cavity is a device in which particles are accelerated by a resonating radio wave. In
most electron-positron colliders, cavities are microwave resonators of "toroidal" shape connected
by tubes, to allow particles to accelerate through them. Fig. 1.1 shows a schematic of a typical
cylindrically symmetric cavity for lowest RF frequency mode (TM 010). The electrical and magnetic fields lines are depicted by arrows and crosses respectively. The electric field is parallel to
the beam axis and decays to zero radially upon approaching the cavity walls. Boundary conditions
demand that the electric field is normal to the surface of the cavity cell with the peak being located
near the iris or the region where the cavity is attached to the connecting tubes. The magnetic field
runs in the azimuthal direction, with the highest magnetic field located near the equator of the
cavity and diminishing towards the cavity axis.
As the particle beam traverses through the cavity, it experiences an acceleration along the axis
of the cavity due to the electric field. Since the RF fields alternate in time, the particle beam
must be in the proper phase with the fields so that the force is accelerating in nature rather than
decelerating. As the particles take a finite time to travel through the cavity, the accelerating field is

2

Cavity Equator

Magnetic Field

Particle
Beam

Cavity
Axis

Electric Field

Cavity Iris

Figure 1.1: A typical accelerating cavity geometry, showing particle beam and fundamental fields
of the RF cavity.
considered as the time average of the electric field along the particle path.
Surface currents in mega Ampere (MA) range are required to flow through the metal walls
of the cavity structure to achieve the desired accelerating field of 1 · 106 V m−1 . Such huge surface currents produce an enormous amount of heat in a normally conducting metal such as copper,
which is costly to cool. Also, the cavity structure near the particle beam has disruptive effects on it,
limiting its density and thereby restricting the accelerator’s use for modern day applications. These
detrimental effects can be overcomed by the use of superconducting cavities. When cooled below
a critical temperature Tc , many materials lose all of its Direct Current (DC) electrical resistance.
These materials are known as superconductors. Below Tc , the electrons of a conducting material
gain a small net attraction through their interaction with the surrounding lattice. As a result, the
electrons form "Cooper pairs" which move through the lattice without any resistance. The superiority of superconducting cavities over those made of normal-conducting metal is in its ability to
store a large amount of energy with much lower dissipation.

3

To determine which materials are suitable for producing high-performing SRF cavities, one
should look for certain properties. One of them being the Quality factor Q0 measuring the performance of a SRF cavity. The Quality factor is given as,

Q0 =

G
Rs

(1.1)

where G is a parameter, which depends on the cavity geometry and Rs is the average surface
resistance of the inner cavity walls (function of the accelerating gradient, Eacc ). In order to reach
the theoretical limit of current and future accelerators, increasing the quality factor and accelerating
gradient is one of the main topics of research on SRF cavities. For a large Eacc , a large superheating
field Bsh is desirable, as it is the relevant critical field for RF applications. In SRF cavities, the
objective is to exclude flux, not to pin it inside the superconductor. Pinned normal conducting
vortex cores dissipate energy strongly in RF fields.
It is energetically favorable for the vortices to be inside the superconductor above the first critical field B1c , and only an energy barrier prevents the flux penetration above this field. For an ideal
surface, the barrier disappears at Bsh , making it the ultimate limit for SRF applications. Surface
defects of size on the order of the coherence length ξ (spatial amplitude of the superconducting
wave function) can reduce the barrier to vortex penetration. Thus, the fields close to Bsh cannot
be realized without vortex entry. Therefore it is important to have a large enough ξ so that the
performance is not affected by the presence of small defects even on a well-prepared surface, such
as grain boundaries. How large of ξ is large enough, is still an open question.

1.2 Material for SRF cavity application
For an optimum cavity operation, the building material must possess the following attributes:
4

• Easy to fabricate in a way that it conforms to a complex geometry.
• High thermal conductivity or be able to be deposited on a substrate with high thermal conductivity to evade thermal breakdown
• Minimal surface defects to prevent field emissions.
• Less prone to surface contamination upon cleaning.
Table 1.1: Superconducting properties of some promising materials for SRF applications.
Material

λ (0)
[nm]

Nb3 Sn[26]
111
MgB2 [27, 28] 185
NbN[29]
375
Nb[30]
50

ξ (0)
[nm]
4.2
4.9
2.9
22

Bsh
[mT]
410
210
160
210

Tc
[K]
18
40
16
9.2

ρn
[µΩ cm]
8
0.1
144
2

The superconducting properties of four promising SRF cavity materials are listed in Table 1.1.
Nb3 Sn is a material with an enormous potential for SRF applications. It has a large Tc ∼18 K
(twice that of niobium). However, Nb3 Sn is brittle and has low thermal conductivity, which limits
it to be used in a film form. Therefore it faces challenges associated with thin films in SRF Cavity
applications. MgB2 is another alternative that was discovered recently. It has a Tc of approximately
40 K, but it is not as well developed as Nb3 Sn. The challenges here lies in obtaining a minimal
oxygen background in the coating chamber to fabricate MgB2 , as magnesium is highly reactive
with it. The δ phase of NbN also has the potential to be a SRF material, with a high Tc of 15–
17 K. NbN has a very complex phase diagram, which makes it very challenging to achieve the
superconducting delta phase during fabrication. Some success has been achieved with sputtering,

5

but with heat treatment in a nitrogen atmosphere, no reliable procedure has been developed that
could be applied to cavities.
Niobium is currently the material of choice for superconducting cavities. It has been an essential transition metal because of its uses in the manufacturing of stainless steels and superalloys for
aerospace applications [31]. The primary reason for the choice of niobium as an SRF material is
due to its highest critical transition temperature (Tc = 9.2 K) of all pure metals, and it is relatively
easy to fabricate. Also, it is readily available and has a high thermal conductivity. None of the
materials mentioned earlier can match niobium in terms of its ease of use and performance with
increasing RF fields. For a typical cavity made of Nb, the resonating frequencies fall in range
0.5–2.0 GHz operating at a temperature of ∼ 2.0 K. The quality factor for such a cavity is about
1010 –1011 range [32].

1.3 Challenges associated with using Nb as material for SRF application
High purity Niobium ingots are produced using high-vacuum multiple electron beam melting
process which results in a large-grain structure in the center and smaller columnar grains at the
periphery of the 20–30 cm diameter ingots. Sheets and tubes are produced by rolling and forging
these ingots, resulting in a polycrystalline microstructure. However, homogenizing this microstructure is difficult, since each of the few large and differently oriented grains responds in a strongly
anisotropic way to the deformation process, which causes variability in the crystallographic texture
of different batches of rolled niobium sheets [33]. Other problems associated with polycrystalline
Nb include magnetic flux pinning and phonon scattering near grain boundaries. To mitigate these
problems cavities can be made from wire cut slices of large-grain Nb ingots containing a handful
of single crystals.

6

Wire Cutting

Nb Disk

Nb Ingot

Electro Beam
Welding

Stamping or
hydroforming

Figure 1.2: Mechanical fabrication path of an SRF cavity.

Figure 1.3: Surface defects in a well prepared SRF cavity cell.
A typical work-flow of an SRF cavity fabrication process is shown in Fig. 1.2[34, 35]. The
inherent deformation anisotropy of the individual grains leads to variability in the forming process
(see Fig. 1.3[34, 36, 37]). Deformation paths involving surface working, such as spinning, introduce gradients of deformation from the surface inward, with a higher density of defects near the
surface. These defects act as hot-spots for electric field emission and subsequent thermal breakdown [38, 39]. Thus, to gain precise control of the forming process, it is essential to understand
7

the mesoscopic deformation behavior in terms of stresses necessary to activate dislocation slip on
various slip systems, and the resulting work hardening behavior from the interaction of dislocations in them. Identification of these criteria are important for predicting how crystal orientations
and flow behavior will evolve in more complex forming operations.

1.4 Objective
The primary objective of this thesis is to develop a finite strain crystal plasticity model at mesoscopic scale for predictions of strain hardening and lattice reorientation in pure Niobium over a
broad range crystal orientations. To achieve this goal, the boundary value problem associated with
mechanical testing of materials is solved computationally. A new dilatational constitutive model
is introduced into the finite strain framework (described in Chapter 3) for simulating deformation
of non-compact geometries using a F OURIER-based spectral method. The developed constitutive model is then used to mimic free-surface boundary conditions within a periodic hexahedral
cell, thus providing a way to test constitutive models for single crystal on a standard dogbone
shaped geometry. In Chapter 4, a physics based crystal plasticity model is combined with internal
state evolution laws consistent with Discrete Dislocation Dynamics (DDD). Parameters critical for
simulating multi-stage strain hardening are identified and discussed. Finally, in chapter 5, the developed finite strain framework combined with a DDD informed physics-based model makes it is
possible to successfully predict multistage hardening and crystal rotation with respect to the tensile
axis as a function of crystal orientation at finite strain levels.

8

CHAPTER 2
PLASTIC DEFORMATION IN BODY CENTERED CUBIC METALS

Deformation by slip in crystalline materials occurs by movement of dislocations in certain crystallographic planes along certain crystallographic directions. Generally, these are crystallographic
planes lattice planes with high planar density. For example, in Face Centered Cubic (FCC) crystals,
dislocations move on {1 1 1} slip planes and along ⟨1 1 0⟩ slip directions. These high density planes
together with a tightly packed direction constitute a slip system. When shear stresses in a particular
slip system exceed a threshold stress, a dislocation starts moving on those systems. This threshold
stress, in general, is constant for a particular slip system and is independent of the direction of slip
[40]. Deformation takes place on those systems which offer the least resistance to the dislocation
motion, which is observed macroscopically by the formation of slip traces (slip step) on a polished
surface after deformation. In other words, the slip system that has the maximum resolved shear
stress will be favorable for slip. This orientation dependence of resolved shear stress is known as
the Schmid law [40] and have been extensively verified for {1 1 1}⟨1 1 0⟩ systems in FCC and basal
slip systems {0 0 0 1}⟨1 1 2 0⟩ in Hexagonally Packed (HP) metals.
Niobium exists in the Body Centered Cubic (BCC) crystal structure from T = 0 K to its melting
point. Therefore, to predict the deformation behavior of pure Niobium, it is important to have
a detailed knowledge of the crystal structure of the material. The first part of this chapter will
therefore concentrate on describing the deformation behavior BCC metals in general. Then, in the
second part, the characteristics of dislocation microstructure in a BCC lattice will be detailed.

9

Figure 2.1: Potential slip systems in BCC lattice

2.1 Slip systems in BCC crystals
The question of slip systems in BCC metals has long been the subject of research since earlier
times. The direction of slip is always observed to be the ⟨1 1 1⟩ direction [41] of the BCC structure,
but no conclusive evidence of the slip plane has been reported. There are, however, three major
points of view in the literature concerning slip planes. Early studies of iron single crystals first
suggested that slip occurs on {1 1 0}, {1 1 2} as well as {1 2 3} planes, which are in the zone of
the slip direction ⟨1 1 1⟩ [42, 43, 44, 45]. Other authors then reported that the slip did not occur
on planes of low indices, but took place on planes close to the Maximum Resolved Shear Stress
Plane (MRSSP) [46]. This theory is often referred to as non-crystallographic glide. These two
ideas of slip in iron are mainly based on the observation of visible slip traces. But as shown in
Fig. 2.1 (right), the projections of the different planes available for slip along the slip direction are
quite close. The planes closest to the families {1 1 2}, and {1 2 3} are separated only by 10◦ . It is
therefore clear that the proximity and the multiplicity of the available slip planes make the interpretations of the slip traces rather delicate. Also, they are often very winding and branched. This
observation has consequently given rise to a third proposition according to which the macroscopic

10

slip traces, which at first sight coincide with a plane {1 2 3} [47], are in fact composed of short slip
traces on planes of the family {1 1 0} [48, 49, 50, 4].

P

Z
L
Y

X

P'

L'

Figure 2.2: Areas of activation of the 110 and 112 planes in the standard triangle [4]
Franciosi [4], Duesbery and Foxall [10] defined the activation zones for slip families {1 1 0} and
{1 1 2} in the standard triangle (see Fig. 2.2). For this purpose, they have compared the Schmid
factors of the two families of planes, assuming that the critical constraints on each of them are
equivalent. These authors thus show the existence of three activation domains: the Y domain,
which is the most extensive, corresponds to the activation of a {1 1 0} plane system, and the X
and Z domains coinciding with the zones of activation of a {1 1 2} system. The particular loading
directions which lie on the boundary between two different activation domains is equivalent to the
simultaneous activation of two systems, as was already the case with a single family of sliding
planes. For example, loading on LP leads to the activation of two systems of the family of planes
{1 1 0}, loading along the boundaries LL’ or PP’ induces the activation of two systems of different
natures {1 1 0} and {1 1 2}. The same authors also carried out tensile tests to verify the position
of these three domains and the mixed boundaries separating them (LL’ and PP’). Franciosi [4]

11

states that he confirmed the existence of some of these boundaries through the observation of slip
traces. The only exception was for a loading close to the PP’ border. The double mixed slip on the
expected {1 1 0} and {1 1 2} systems has not been observed in favor of a simple slip on a {1 1 0}
system.

2.2 Orientation dependent strain hardening in BCC single crystals
100

h1 1 1i

A

D

80

C

τ / MPa

E
60

E

B

F

40
B
A

20
0

C

0

0.2

0.4

0.6

0.8

1

F

D
h0 0 1i

γ

h0 1 1i

100

h1 1 1i

80

τ / MPa

8

60
1
40

3

8 10 2

3

10

20
21

0

0

0.2

0.4

0.6

γ

0.8

1

1
h0 0 1i

h0 1 1i

Figure 2.3: Shear stress – shear strain curves as a function of orientation in α − Fe (top) and Nb
(bottom), according to [5, 6].
In BCC single crystals two different types of deformation behavior are reported at strain-rate

12

independent temperatures. Fig. 2.3 shows several uniaxial tension deformation curves obtained in
pure Iron and Niobium at ambient temperature by Spitzig and Keh, and Keh and Nakada [5, 6, 10].
Most of the tensile axis within the standard stereographic triangle show well-known three-stage
stress-strain behavior (curves A, B, C, 10, and 3 in Fig. 2.3). Loading directions close to the high
symmetry orientations, such as [001], [111] and [011], show pseudo-parabolic behavior (curves D,
E, F, and 1 of Fig. 2.3). In these orientations, slip operates on several systems from the beginning
of the plastic deformation, which is responsible for a very high rate of hardening. The comparison
of the three curves D, E and F (top row Fig. 2.3) leads to a prioritization of the hardening between
the orientations of high symmetry. Thus the hardening is decreasing for the loadings in the order
{0 0 1}, {0 1 1}, and {1 1 1}. This hierarchy has also been observed by other authors [9, 8] and
depends on the number of systems involved.
[1 0 0]

[1 1 0]

[1 1 0]

MRSSP

[1 0 1]
[1 1 1]

[1 1 1]
Specimen
Axis

[0 1 1]

[0 1 0]

Specimen
Axis

[0 1 1]

[0 1 0]

[0 0 1]

[1 1 1]

[1 1 1]
[1 0 1]

[1 1 0]

[1 1 0]

[1 0 0]

Figure 2.4: Schematic illustration of Taylors angles χ and φ . The loading direction is marked in
the standard triangle by the angles χ and φ [7]
.
In addition to the type and number of slip systems, the loading direction also affects the shape

13

of the three-stage curves as shown by a number of studies on this subject. This is a substantial
effect, as the degree of hardening and the length of the stages are modified. There are few studies
in which systematic identification of the direction of loading in the standard triangle by two angles

χ and φ (see Fig. 2.4) are defined [49, 5, 7]. χ is the angle between the MRSSP and the plane
(1 0 1), while φ is the angle between the load direction and the boundary [0 0 1] - [0 1 1]. Most of
these works were done in materials with dominant {1 1 0} slip behavior. The effect of angle χ has
been somewhat discussed [49, 5, 7] in comparison with the angle φ [6, 7]. The results obtained by
Kumagai et al. [7] are given in Fig. 2.5. The length of stage I increases with the angle φ , while
the degree of hardening of stage II decreases. The rate of work-hardening is lowered from µ /600
(curve A) to µ /1000 (curve C). The results obtained by [6] are similar. The elongation of stage I
with φ is often explained by the rotation of the loading direction with the plastic deformation.
Jaoul and Gonzalez [8], Keh [9], Duesbery and Foxall [10] identified the initial and final orientations of the loading direction during their tensile tests. A synthesis of their results is proposed
in the standard triangle in Fig. 2.6. Under simple slip conditions on a plane {1 1 0}, the loading
direction turns in the slip direction, which is [1 1 1] in the example figure. Upon reaching the [0 0 1]
- [0 1 1] boundary, the activation of a second system should, in theory, turn the loading towards
the [0 1 1] axis. However, this is not always the case, since the direction of loading exceeds the
boundary (or overshoot) by sometimes 6◦ to 10◦ [8, 9, 7, 10]. The reason for these overruns is not
well understood yet. A loading along the [1 1 1] - [0 1 1] boundary, is directed towards the [0 1 1]
direction while loading on the [0 0 1] - [0 1 1] boundary, move towards the [0 1 1] axis. [8] showed
that a load near [0 0 1] would correspond to a simple slip along a plane {1 1 2}.
Kumagai et al. [7] also reported that with increase in the angle χ , the length of Stage I decreases
along with increased hardening in Stage II (see right column of Fig. 2.5 ). The strain hardening
14

80

80

60

τ / MPa

100

τ / MPa

100

B

C

A

40
20
0

60

A

B

D

E

C

40
20

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

γ

0.6

0.8

1

γ
h1 1 1i

h1 1 1i

C

E

A

B

B

C

D

A

h0 0 1i

h0 1 1i

h0 0 1i

h0 1 1i

Figure 2.5: Influence of the orientation of the load on the shear stress-strain curves with slip only
in {1 1 0} system. (Left) effect of angle φ , (Right) effect of angle χ [7].
rate increases from µ /800 for orientation A to µ /400 in curve E. Unfortunately, these results
can not be compared with any other study, because the majority of existing works vary χ and φ
simultaneously, preventing any conclusion on the role of each alone.

2.3 Non-Schmid effects
BCC metals are known for their direction dependent slip i.e. tension-compression asymmetry.
The underlying cause of this behavior can be traced back to the core configuration of screw dislocations in BCC lattice. Due to the polarization of screw dislocation core (see Fig. 2.7), BCC metals

15

Figure 2.6: Evolution of lattice reorientation with plastic deformation depending on the initial
orientation [8, 9, 10].
show a glide direction sensitive slip behavior. It means the slip resistance in one direction is not
the same in the opposite direction. This direction sensitive slip behavior is called twinning/antitwinning asymmetry in BCC metals and arises due to lattice symmetry based reasons. Neuman’s
principle [51] states that physical properties of the material have a similar kind of symmetry owned
by the crystal structure itself. There is no symmetry based reason to believe that positive and negative shear in the {1 1 1} slip direction to be equivalent to each other (except on {1 1 0} planes for
which the {1 1 0} dyad axis imposes this symmetry). Observations suggest that for most orientations where the MRSS is {2 1 1}, the CRSS for the slip in twinning direction is lower than that
of atwinning direction [52, 53, 54]. Hence, it is very likely for the BCC materials to show this
anisotropy irrespective of its electronic structure, type of bonding and in the absence of stress. In
practical terms, this twinning / anti-twining effect is observable for pure shear stresses in the direction of the burgers vector. Direct measurement of this phenomenon is difficult, as much larger
16

[1 1 2]
[0 1 1]

~b = a/2[1 1 1]

(0 1

1)

~b=6

[2 1 1]

[0 1 1]

[2 1 1]

~b=6

~b=3

(0

[1 1 0]

1

[1 1 0]

1)

(1 0 1)
(1 1

[1 2 1]

0)

[2 1 1]

[1 1 2]

[1 2 1]

~b=3

~b=3

~b=6

(1 0 1)

[1 0 1]

(1

10
)

b~=
6

[1 2 1]
~b=6

[1 1 2]

~b=6

[1 0 1]

(110)

(110)

(101)

(101)

(0 11)

(0 11)

Figure 2.7: Illustration showing core configuration of screw dislocations in a BCC lattice
non-glide stresses often masks it.
Moreover, due to a small fraction of edge components in the partials forming screw the dislocation core, the non-glide components of the applied stress affect the movement of dislocations
through the lattice. For these reasons, the Schmid Law breaks down in BCC metals. Macroscopically the violation of Schmid Law is observed in tension-compression asymmetry and movement
of screw dislocations even at low stresses. Qin and Bassani [55] suggested non-Schmid effects,
to capture the anomalies resulting from the core of screw dislocations. The influence of CRSS by
the components of applied stress other than that causing force on displacement was first observed
by Escaig [56, 57] while studying cross-slip of partials in FCC materials. The 1/6⟨2 1 1⟩ partials

17

dissociated from screw dislocations have edge components. Thus shear stresses in the direction
perpendicular to the total Burgers vector widens or narrows the splitting. This is known as Escaig
Effect in cross-slip and is not only important in FCC [58] but also in L12 intermetallic structures.
The edge components of Shockley partials of a/2⟨1 1 0⟩ dislocations in FCC are fixed, but in the
case of fractional dislocations forming the core of a/2⟨1 1 1⟩ screw dislocation in BCC, the edge
components change with applied stress. As long the vector sum of the edge components in the
polarized screw dislocation core is zero, it is possible for the core to have some edge character.
Edge components of fractional dislocations are small in the unstressed state. However, non-glide
components may induce edge components significantly which change the symmetry and width
of the dislocation core. Atomistic investigation suggests that components of applied stress acting to the burgers vector normal to the screw dislocation line heavily influence the width of the
dislocation core. The calculated value of Peierls stress decreases by a factor of four [53] if these
non-screw components are suppressed and thereby prevent relaxation of the screw dislocation core
along directions perpendicular to the dislocation line. The effect of non-glide components of the
applied stress is prominent in VIB metals (Cr, Mo, W) while in VB metals (V, Nb, Ta) it is weak
or absent. It is to note that most of these atomistic calculations were done on rigid dislocations at
zero temperature. At finite temperature, dislocations move by the nucleation and propagation of
double kinks, which is discussed in next section.

2.4 Kink-pair assisted dislocation motion
The flow stress of pure BCC metals strongly depends on temperature and strain-rate. In-situ
electron microscopy observation of dislocations in BCC metals suggests that edge dislocations
have high mobility compared to screws [59, 60, 61]. High mobility of edge dislocations enables

18

them to exit out of the crystal forming slip steps and leaving screw dislocations behind, which
in turn control the strain rate [62, 63]. The low mobility of screws explains the rapid increase
of flow stress with decreasing temperature or increasing strain rate [64, 62]. The difference in
mobility arises from the fact that for a dislocation to be mobile, it has to overcome the intrinsic
periodic lattice resistance (Peierls potential, Fig. 2.8), which is different for different dislocation
types [62, 65]. Edge dislocations have a low Peierls energy barrier compared to screws, and thus
have higher mobility [64, 62]. At temperatures higher than a critical one (knee temperature Tk ), slip
occurs at stresses lower than that required for overcoming the Peierls barrier [66]. Above Tk , Screw
dislocation of sufficient length can change plane on the simple condition that the effective stress in
the deflection plane is greater than that in the primary slip plane. Given the high multiplicity of the
possible plane for the deviation, such behavior is quite probable.
Peierls Barrier

Kink pair

Dislocation line

Figure 2.8: Intrinsic lattice resistance (Peierls-Nabarro Barrier) to slip in BCC materials.
Above the knee temperature, screw dislocations lying on low index crystallographic planes can
overcome their Peierls barrier [67] by lifting a short segment of the dislocation over the barrier
from time to time and letting the segments connecting the two adjacent valleys migrate along the
Peierls hills [68]. Shockley gave the name ’kinks’ to these segments of dislocations crossing the
Peierls barrier. The kinks experience periodic potential energy variation by its position z along
the hills. This position dependent potentials energy is referred to as the Peierls potential of the
19

second kind or kink potential. Whether the kink potential is of the same value or smaller than that
required for kink formation is determined by the bonding in crystals [62]. For the first case, the
kinks are confined to a separation of two neighboring atoms. Such confinement of kinks requires
directionality in bonding. In metals due to non-directional bonding, the kink potential is lower
than that required for kink formation. As a consequence of this, the kink width is larger compared
to interatomic distance. The activation volume of kink migrations is in the order of magnitude of
only b3 [68]. As a result, thermal activation is highly probable. Additionally, the applied resolved
shear stress can lower the barrier height, thereby assisting the generation of kink pairs even at
low applied stress and on slip systems with a low Schmid factor. Such anomaly of slip behavior
as a function of temperature was explained by Seeger and Wasserbach [11] in terms of kink-pair
configurations.

2.5 Temperature dependence of flow stress
Seeger and Wasserbach [11] showed that the flow stress dependence of temperature in BCC
metals could be classified into three temperature regimes and two stress regimes, based on the
active mechanism of slip at that particular stress and temperature (Fig. 2.9). The numerical values
in Fig. 2.9 correspond to the experiments performed by Hollang et al. [69] on ultra-pure Mo crystal.
As mentioned earlier, below the knee temperature TK the flow stress of BCC materials is controlled
by the mobility of a0 /2[1 1 1] screw dislocations given that most of the non-screw components have
been removed by pre-deformation. Experiments performed by Hollang et al. [69] on ultra pure Mo
demonstrated that the flow stress of BCC single crystals could be given as

σ = σM + σ ∗

20

(2.1)

1000

flow stress, τ [MPa]

Peierl s Barrier

800

Ki nk pai r

Annihilation

Dislocati on line

600
Junction

400
200
0

0

200

400

600

800

Temperature, T [K]
Figure 2.9: Temperature dependence of flow stress, slip mechanisms (schematics) in pure BCC
metals (Mo) [11].
where σM is the strain-rate-independent flow stress and σ ∗ is the effective stress required for
movement of a0 /2[1 1 1] screw dislocations. Above the knee temperature TK , the mobilities of
the screw and non-screw components are quite comparable, which is why the flow stress reduces
to a strain-rate-independent value σM . Below the knee temperature, screw dislocations have low
mobility due to a large Peierls barrier compared to the non-screw components. As mentioned
21

before, this barrier can be overcomed by thermally activated, and stress assisted formation of kink
pairs. At temperatures T (Ť < T < Tk ) and effective stress σ ∗ is less than a critical value σ̂ ∗ , the
flow stress is determined by the interaction of the pair of kinks of opposite sign and is known as
elastic-interaction regime. The formation enthalpy of kink pairs in this regime is given as

Hkp (σ ∗ ) = 2Hk (σ ∗ ) − 2ασ ∗1/2

(2.2)

where α 2 = a3 bγ0 /2. 2Hk denotes the formation energy of the two widely separated kinks
of opposite sign, a the kink height and γ0 the pre-factor proportional to the square of dislocation
strength. α can be determined experimentally, and from that, the kink height a can be deduced.
Seeger and Wasserbach [11] reported that for kink height a =

√

2a0 , which corresponds to the dis-

tance between Peierls valleys in {1 1 2} planes, the experimental data matches the above equation
supporting that at higher temperatures {1 1 2} slip is favored.
At effective stresses higher than a critical value (σ ∗ > σˆ∗ ), the kink - pair formation enthalpy
is determined by line tension of dislocations and is known as line tension regime. Ackermann et al.
[70] showed that approaching the critical effective stress σˆ∗ from both top and bottom, using line
tension model and elastic interaction model respectively, results in a different activation volume.
Hence, the mechanism of kink pair formation changes. This has been verified experimentally by
Hollang et al. [69] in ultra pure Mo single crystals. At lower temperatures T < Ť , extrapolation
of the above-described line tension model predicts a lower flow stress compared to what has been
observed in experiments. Seeger and Wasserbach [11] hypothesized that at T < Ť the core of
the screw dislocations assume ground state configuration i.e. the configuration with lowest line
energy. The screw dislocations in this configuration exhibit threefold symmetry of the ⟨1 1 1⟩

22

crystallographic direction and can glide on any one of the three {1 1 0} planes of that zone. On
application of finite shear stress, the symmetry of the dislocation core is broken in such a way that
the slip plane with highest Schmid factor is chosen as the slip plane for dislocation motion. This
phenomenon explains the high flow stress of BCC metals at low temperatures and observance of
{1 1 0} slip lines for orientations in the middle of standard stereographic triangle.

2.6 Overview of deformation in Nb single crystals
Some of the earlier work on the plasticity of pure Niobium (Nb) was carried out by Mitchell
et al. [71], Votava [72], Taylor and Christian [73] in the 1960s. They reported that under suitable
conditions of purity, orientation, temperature and strain rate, the stress-strain curves in shear show
a 3-stage hardening behavior similar to that observed in FCC metals and alloys. Multi-stage hardening in Nb single crystals grown by electron beam zone melting technique were also reported by
Duesbery et al. [74]. Mitchell et al. [71] reported that, the work hardening rate in stage II was equal
to about µ /600 for a strain rate of 4.4 × 10−5 s−1 . This value was later modified by Duesbery et al.
[74] to about µ /500 for the standard conditions used in their work. Annealing of these samples
caused the work hardening rate in stage II to increase from µ /700 to µ /400 (closer to that of FCC
metals µ /300) for an annealing time of more than 12 h. For an orientation lying in the center of the
standard stereographic triangle, slip on {1 1 0} planes were observed. Straight slip lines were seen
on the side faces while those on the top face were shorter and wavy. During the transition from
stage I to stage II and stage III, extremely coarse ’hills’ and ’valleys’ were formed on the top face
(See Fig. 2.10). Correlation of these ’hills’ with side faces indicated that the ’hills’ were formed in
the region of unrestricted primary slip and the ’valleys’ correspond to conjugate slip activity.
In a following article by the same group, a detailed study was done on the tensile axis reorien-

23

Figure 2.10: Surface topography of deformed single crystal Nb [12]
150

X
S QR
W

P

V

σ / MPa

T

V
P
W
R
X
Q
T
S

100

50

0

0

0.1

0.2

0.3

ε

V
P
W
R
X
Q
T
S

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

ε

Figure 2.11: In situ uniaxial tensile test results reported by Baars [12]. Left : Inverse pole figure
showing rotation of tensile axis as a function of orientation for as-cast samples. Center: axial
stress–strain curves for the as-cast samples with same orientations in the left inverse pole figure
map. Right : axial stress–strain curves for the annealed samples with the same orientations.
tation in pure Nb single crystals grown by zone melting. The extent and rate of re-orientation were
studied as a function of orientation. Most the specimens having orientations in the center right part
of the standard stereographic triangle showed trajectories that overshoot the [0 0 1]–[1 0 1] boundary. This behavior suggested that most of the orientations started with a single slip system until
the tensile axis had moved a small distance beyond the [0 0 1]–[1 0 1] boundary when conjugate slip
occurs. The maximum distance after crossing the boundary at which conjugate slip takes place
was about 6◦ for orientations in the central upper regions of the triangle. This angular distance was
less for the initial orientations belonging to lower right part of the standard stereographic triangle.
Recently, uniaxial in-situ tensile tests performed by Baars [12] on single crystal Nb showed

24

multi-stage hardening in macroscopic stress-strain curves and texture evolution similar to Duesbery and Foxall [10]. They also observed the formation of ’hills’ and ’valleys’ on the surface of
the sample along with bowed slip lines. The stress-strain behavior of eight single crystals in the ascast (left) conditions are shown in Fig. 2.11 [12], with the crystal orientation of the tensile axis is
indicated by letters in the standard stereographic triangle. The stress-strain curves suggest that different orientations deform very differently, but the trend is that orientations near the ⟨1 1 1⟩ corner
require the most stress, while those in the center of the triangle require the least stress. This general
trend persists despite an initial annealing treatment that reduced the dislocation content resulting
in 20 % lower flow stresses. The tensile axis re-orientation relative to the crystal frame suggested
that slip occurred predominantly on {1 1 2} planes for the as-cast condition, but this changed to slip
occurring mostly on {1 1 0} planes in the pre-annealed samples [75]. The strain hardening rates
increased as the tensile orientations move closer to the boundaries of the stereographic triangle,
consistent with the activation of slip on additional planes with the same ⟨1 1 1⟩ directions.

25

CHAPTER 3
MODELING NON-COMPACT GEOMETRIES USING SPECTRAL METHOD

A framework for simulating the deformation behavior of materials involves two important aspects,
thermodynamics and kinematics of the deformation process. The former describes why the deformation takes place, while the latter deals with the mechanism by which the material is deformed. In
this chapter, an emphasis is given to the kinetic aspect of deformation. As deformation comprises
of rotation and distortion, it is necessary to define two separate configurations i.e. an undeformed
configuration B0 and a deformed configuration B. Mechanical fields such deformation and stress
can be represented in either configuration. For a material that has been deformed to a finite extent, the deformation is defined by an operation, involving linear mapping of the materials points
in the undeformed configuration to the deformed one [76]. This linear mapping is known as the
deformation gradient, F. The boundary value problem solver (FEM or spectral solver using a fast
Fourier transform) determines the boundary condition in the form of an average deformation gradient F and the material point model provides the corresponding average first Piola-Kirchhoff stress
P. For a finite strain elasto-plastic deformation, the total deformation gradient F is decomposed
multiplicatively into two parts, the elastic deformation gradient Fe and the plastic deformation gradient Fp [77]. A schematic illustration of the multiplicatively decomposed deformation process is
provided in Fig. 3.1. The multiplicative decomposition introduces an intermediate configuration
so that the deformation process occurs in two steps. First, the material deformed plastically to a
stress-free state. After that, an elastic deformation is imposed to account for the final shape change
as observed in the deformed configuration.
Following this finite strain framework, Roters et al. [78] developed a material point model
26

Figure 3.1: Multiplicative decomposition of deformation gradient [13]
DAMASK [79], to capture the mechanical response of a material as a function of its microstructure.
In this material point model the elastic deformation gradient Fe is related to the second P IOLA–
!
"
K IRCHHOFF stress tensor (acting on intermediate configuration) S = C : Fe T Fe − I /2, where C

!
"
is the elastic stiffness tensor. The second P IOLA–K IRCHHOFF stress S = f F, Ḟ, η acts indirectly
as the driving force for the plastic velocity gradient Lp which is a function of microstructural
parameters η in the constitutive model. Depending on the type of material and computation cost
involved, a constitutive model can then be formulated as a function of η and S to calculate Lp .
Once Lp is calculated the plastic deformation gradient can be estimated as Fp = Lp −1 F˙p , which
is then used to calculate Fe . The interdependence of Lp and S is solved using self-consistency
(see Fig. 3.2) at each discretized point within a material static equilibrium is attained. Using the
formalism mentioned above, the partial differential equations describing the mechanical loading of
a material can be solved for static equilibrium with a boundary value solver such as Finite Element
Method (FEM). This has led to the development of a separate area of research in computational
mechanics known as Crystal Plasticity Finite Element Method (CPFEM).
27

S

Fe
Fe =
−1

F Fp

Fp

Lp
Ḟp = Lp Fp

Figure 3.2: Schematic diagram showing calculation of S as a function of F.
CPFEM has been the vehicle of choice to simulate and understand the mechanical behavior
of engineering components with structural dimensions approaching the grain size of the intrinsic
microstructure. Examples of such oligocrystalline mechanics can be found in the fields of leadfree solder joints, thin wires, cardiovascular stents, and mesostructured materials such as foams or
photonic crystals. The spatial resolutions achievable by CPFEM within reasonable computation
times are limited by the substantial numerical effort inherent in the FEM strategy to solve the
underlying partial differential equation system.
Spectral methods have emerged as an efficient [80, 81, 82] substitute for CPFEM [83, 84, 85,
86], but only recently were improved sufficiently to overcome their unfavorable convergence when
large gradients in mechanical properties are present [87, 88]. This enhanced capability allows the
mechanics of porous structures to be solved using spectral methods by replacing the open space (or
voids) with soft material [89] in the simulated periodic domain. Examples of such simulations have

28

been demonstrated in recent works [90, 91, 92, 19]. The methodology used in these works tightly
links the behavior of the spectral solution algorithm to the specifics of the constitutive behavior
of each material point, i.e., making a binary distinction between how to address void and filled
points. Furthermore, any nucleation of additional voids not yet present in the initial geometry
will require to propagate the information about such a change in local constitutive behavior to
the spectral solution algorithm. Such a rigid approach is akin to element elimination in FEM
simulations of crack propagation, where a particular constitutive response, i.e. loss of stiffness, is
treated through ad hoc modification of the geometry instead of properly at the constitutive level, as
achieved by, for instance, cohesive zone elements. Following this spirit of keeping any boundary
value problem solver as general as possible, i.e., independent of the constitutive material behavior
that is simulated, voids should only distinguish themselves through their constitutive model, which,
classically, has been vanishingly small elastic stiffness. Here, to describe void regions, a plastically
dilatational material is proposed that enables simulations of non-compact geometries such as open
or closed-cell foams, or other geometries with free surfaces, by means of established spectral
methods without any particular adaptations.
A plastic plate with a circular inclusion under biaxial extension and a dogbone-shaped tensile
sample of oligocrystalline aluminum serve to compare this dilatational material model to a lowstiffness isotropically elastic model of all void regions.

3.1 Constitutive model for void like regions
A finite strain framework incorporating two intermediate configurations is adopted, similar
to the work of Tjahjanto et al. [93] and following Shanthraj et al. [94]. The total deformation
gradient F = Fe Fi Fp at each material point is multiplicatively decomposed into isochoric, and

29

lattice preserving, plastic deformation Fp , non-isochoric, but stress-free, dilatation Fi , and elastic
deformation Fe , consecutively mapping the reference configuration into the “lattice” configuration
(Fp ), then into the “intermediate” configuration (Fi ), and finally into the deformed one (Fe ). Only
!
"
the elastic lattice distortion gives rise to stress, which takes the form Sp = C : 12 Fi T Fe T Fe − I Fi

in the lattice configuration and Si = Fi Sp Fi T / det Fi in the intermediate one. The stress Si drives
the dilatational velocity gradient Li (Si , η ) = Ḟi Fi −1 , while Sp drives the plastic velocity gradient

η of internal state variables. The total velocity gradient
Lp (Sp , η ) = Ḟp Fp −1 and the evolution η̇
follows as L = Le + Fe Li Fe −1 + Fe Fi Lp Fi −1 Fe −1 .
To capture the mechanical response of a void, an isotropic plasticity model that combines an
isochoric response due to the deviatoric stress with a dilatational response due to the hydrostatic
pressure (mean stress) is formulated. The kinetics and internal state parameterization of the model
are inspired by the phenomenological crystal plasticity model introduced by Peirce et al. [95] that
postulates a power-law relation and an internal deformation resistance, termed τ0 .
Thus, in an isotropic setting, the strain rate ε̇p connected to isochoric deformation is given as

ε̇p = ε̇0

#√

3J2
M τ0

$n

= ε̇0

%&

∥Sp′ ∥

3
2 M τ0

'n

,

(3.1)

with J2 as the second invariant of Sp ′ (deviatoric second P IOLA–K IRCHHOFF stress), stress exponent n, Taylor factor M, and ∥·∥ the F ROBENIUS norm. The associated plastic velocity gradient
Lp , acting in the lattice configuration, is then formulated as
Sp ′
Lp = ε̇p
= ε̇0
∥Sp ′ ∥

%&

3 1
2 M τ0

'n

Sp′ ∥Sp ′ ∥n−1 .

(3.2)

To mimic the dilatational response of a void region, a similar constitutive law but for the dilatational expansion rate ε̇i and the dilatational velocity gradient Li is formulated in the intermediate
30

configuration

#

$
p n
ε̇i = ε̇0
M τ0

(3.3)

#
$
I
p n
Li = ε̇i
= ε̇0
I
∥I∥
M τ0

(3.4)

where p is the hydrostatic pressure calculated from Si , and I is the 2nd order identity tensor. The
evolution of deformation resistance follows
(
(
#
$
(
τ0 (a
τ0
ġ = M ε̇p h0 ((1 − (( sign 1 −
g∞
g∞

(3.5)

where a, ε̇0 , and h0 are adjustable parameters.

The above formulated dilatational plasticity model is closely connected to the VON M ISES yield
surface shown in Fig. 3.3 in principal stress space. Each point on this cylindrical surface represents a state of stress at which the material deforms plastically. Any position vector s (σ1 , σ2 , σ3 )
to a point on the yield surface can be represented as a sum of a pressure vector and a devia!
"
toric stress vector, i.e. a component p σ1 = σ2 = σ3 = σsph along the diagonal and a component

"
!
s′ σ1 − σsph , σ2 − σsph , σ3 − σsph that is perpendicular to the pressure vector. The isochoric deformation rate Lp depends, usually, only on the deviatoric stress vector, while the non-isochoric

contribution Li is sensitive to (and follows) the pressure vector, i.e., there is a flow stress “cap”
that moves along the pressure axis with dilatation.

3.1.1

Solution to mechanical boundary value problem

To solve the mechanical boundary value problem of static equilibrium, the finite strain spectral
method outlined in [82, 88] and implemented as part of the Düsseldorf Advanced Material Simulation Kit (DAMASK) is used. The constitutive law described above was integrated into the flexible
31

pressure axis
σ 1 = σ2 = σ 3

σ3
s'
s

p

σ2

yield surface

σ1

Figure 3.3: Schematics of VON M ISES plasticity model in principle stress space. The state of stress
at a material point is represented by the vector s (black), which can be additively decomposed into
the pressure vector p (red) and the deviatoric stress vector s′ (green). The light green surface
represents the yield surface of the material point.
material point model offered by DAMASK. Since the solution fields resulting from the spectral
method are a superposition of a homogeneous and a fluctuating part, where the mean value of the
latter vanishes over the domain, any boundary conditions can only prescribe the average (homogeneous) fields.

3.2 Comparison between a dilatational and soft elastic void
Using a plastically isotropic plate with a circular void at its center as an exemplary case, the response of the dilatational material model outlined above is contrasted to an alternative description
of the void as a (relatively) soft and purely elastic inclusion. Table 7.1 lists the material parameters
used to model the plastic plate as well as the elastic and the dilatational version of the void. The
largest elastic stiffness constant C11 of the dilatational and elastic version of the void is scaled
down, respectively, one and three orders of magnitude relative to the plate surrounding it, with C12
and C44 calculated such that elastic isotropy and vanishing Poisson ratio for the void is ensured
in both versions. The flow stress of the dilatational version of the void is set to be two orders of
32

Table 3.1: Material parameters; elastic constants Ci j , reference strainrate ε̇0 , stress exponent n,
initial and saturation flowstress τ0 and g∞ , hardening parameters h0 and a, and Taylor factor M = 3.
plate

void
elastic

dilatational

C11
C12
C44

100
60
30

0.1
0
0.05

10
0
5

2C44
C11 −C12

1.5

1

1

ε̇0
n
τ0
g∞
a
h0

10−3
20
30
60
2
80

10−3
20
0.3
0.6
2
1

GPa
GPa
GPa

s−1
MPa
MPa
MPa

magnitude lower than that of the plate material. The chosen values reflect a compromise between
vanishing stress in void regions and the associated computational cost, as detailed in the supplementary material. A 512 × 512 × 1 regular grid is used to discretize the fully periodic geometry,
resulting in about 2100 grid points within the void taking up an area fraction of 0.8 %. The plate is
subjected to biaxial tensile elongation along the x and y direction (see top row of Fig. 3.4), which is
enforced by fixing eight of the components of the average deformation gradient rate and requiring
the remaining complementary first P IOLA–K IRCHHOFF stress component Pzz to vanish, i.e.
⎡
⎤
⎡
⎤
⎢1 0 0⎥
⎢
⎥
˙
⎢
⎥
⟨F⟩
⎢0 1 0⎥
=
⎢
⎥
10−3 s−1 ⎢
⎥
⎣
⎦
0 0 ∗

and

⎢∗ ∗ ∗⎥
⎢
⎥
⎢
⎥
⟨P⟩ ⎢
= ⎢∗ ∗ ∗⎥
⎥,
Pa
⎢
⎥
⎣
⎦
∗ ∗ 0

(3.6)

where ‘*’ indicates that stress (or deformation) needs to be iteratively adjusted since deformation
(or stress) is prescribed. In the simulations, the conditions of Eq. (3.6) are maintained for 400
increments of 1 s each (such that the final ⟨F⟩xx = ⟨F⟩yy = 1.4).
33

elastic
original
geometry &
applied
loading

dilatational

y
x

deformed
geometry

det(Fe )
1

1.5

det(Fi )
1

15

σhyd / 0.2 GPa
0
1

σvM / 0.2 GPa
0
1

Figure 3.4: Comparison of mechanical field maps of a (fully periodic) two-dimensional plastic plate
containing a circular inclusion of area fraction 0.8 % and loaded biaxially to an average deformation of
⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.4. Void volume is described by an isotropic elastic material model (left, “elastic”) or
by an isotropic phenomenological plasticity model with dilatational capabilities (right, “dilatational”). The
stiffness of the inclusion in the elastic case is set to be 1000× softer than the plate.

The elastic and dilatational version of the void are contrasted in the first and second row of
Fig. 3.4, showing the simulation results at the final simulation step in terms of deformed geometry,

34

determinant of deformation gradients (i.e. relative change in volume), hydrostatic stress, and VON
M ISES stress. The two cases reveal markedly different results. For the same final planar geometry
(second column), the dilatational void expands to about six times its original diameter, while the
purely elastically modeled void only about doubles its diameter. As the volume of the plastic plate
is preserved, its final thickness is about 10 % smaller when containing the elastically modeled void.
The volumetric expansion of the void is fully carried by det(Fe ) for the elastic version, whereas the
values of det(Fe ) are minute compared to det(Fi ) in the dilatational void version (third and fourth
column). Since only the hydrostatic elastic strain connected to det(Fe ), but not det(Fi ), is giving
rise to a hydrostatic stress σhyd = K εhyd , with K = 0.43 GPa, the σhyd inside the elastic void is
much larger than inside the dilatational void and even exceeds the hydrostatic stress experienced by
the surrounding plastic plate (column five). To the contrary, the dilatational void does not exhibit
any appreciable hydrostatic stress as well as deviatoric stress (column six) since it is plastically
growing in volume at a flowstress much lower than the flowstress of the plastic plate.

3.3 Uniaxial tension of dogbone-shaped sample of oligocrystalline aluminum
To demonstrate the usefulness of the proposed framework for non-linear geometries, the study
performed by Zhao et al. [14] that compares measured to simulated local deformation behavior
of a thin oligocrystalline Al dogbone tension sample is reevaluated. In that investigation, the
grain structure on the sample front and back face were established through electron backscatter
diffraction (EBSD) before and after straining to about 0.1. The local in-plane deformation field
was acquired using digital image correlation (DIC) at strain increments of 5 · 10−4 . Since most
of the about 20 grains in the thin gage section are almost columnar, the modeled geometry was
based on a through-thickness extrusion of the grain structure on the front face. Contrary to Zhao

35

et al. [14], three different boundary conditions are considered in the present study, namely “fully
periodic” (periodic copies of hexahedral dogbone gauge section), “quasi-periodic” (periodic along
thickness direction, i.e. infinite thickness), and “free surface” conditions. The geometries are
discretized by a regular grid of 273 × 112 points and either 10 (fully periodic and quasi-periodic)
or 20 points (free surface) across the thickness. A buffer layer of either the dilatational material
described above or an isotropic soft-elastic material (see Table 7.1) encases the dogbone geometry
for the quasi-periodic and free surface boundary conditions to arrive at a periodic hexahedral cell
required by the spectral solver (see left column in Fig. 3.5). The constitutive model used to describe
the behavior of the crystalline material is the same as used by Zhao et al. [14]. Uniaxial tension up
to a strain of 0.1 along the y direction, discretized into 1000 equal time increments, is enforced by
setting the average deformation gradient rate and complementary first P IOLA–K IRCHHOFF stress
to

⎡

⎤

⎢∗ 0 0⎥
⎢
⎥
˙
⎢
⎥
⟨F⟩
⎢
= ⎢0 1 0⎥
⎥
−3
−1
10 s
⎢
⎥
⎣
⎦
0 0 ∗

and

⎡

⎤

⎢0 ∗ ∗⎥
⎢
⎥
⎢
⎥
⟨P⟩ ⎢
= ⎢∗ ∗ ∗⎥
⎥.
Pa
⎢
⎥
⎣
⎦
∗ ∗ 0

(3.7)

Figure 3.5 maps the in-plane VON M ISES strain and lattice orientation (as inverse pole figure)
of the three simulated boundary conditions with dilatational voids, the experiment, and the free
surface boundary condition using soft-elastic voids (top to bottom). In the boundary condition sequence using dilatational voids (first to third row), the VON M ISES strain is most homogeneous for
the fully periodic case, which differs significantly from the measured

VON

M ISES strain. For the

quasi-periodic and free surface boundary conditions, the VON M ISES strain localizes in grains with
softer crystallographic orientations (i.e. higher Schmid factor) enabled by the non-compactness of
the simulated geometry. The locations, i.e. the (groups of) grains that deform most, are completely
36

different between the quasi-periodic and free surface case. When compared to the DIC measured
VON

M ISES strain, only the free surface boundary condition with either the dilatational or soft-

elastic voids captures the most relevant strain localization features (third and last row). Despite
this agreement, only with dilatational voids does the simulation properly capture the experimentally observed crystal lattice reorientations (last column). The extent and location of strain heterogeneity is not only influenced by microtexture as reported, for instance, by Raabe et al. [96]), but
it appears also essential to reflect the correct boundary conditions in a simulation when attempting
one-to-one correlations of experimental and simulated data.
In the present study, volumetric dilatation capability is added to an existing isotropic material
model to capture the behavior of void-like regions in F OURIER-based spectral solvers. The ability to mimic empty space is essential to address non-compact geometries within the periodicity
constraints imposed by such solvers. The model has been applied to simulate the plastic growth
of a circular void in a plate under biaxial tension and the known grain-scale deformation behavior of a tensile sample containing a limited number of grains. Both cases show that the proposed
formulation provides a more accurate description of void-like behavior compared to assuming a
soft-elastic behavior for void volumes.

3.4 Selection of parameters for void material description and its associated
sensitivity on the simulation results
To study the influence of material parameter values of the elastic and the dilatational material
model for voids, a systematic series of biaxial hole expansion simulations up to ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.2
was performed. Elastic stiffness ratios between void material and plate spanned the range 1 to
10−3 ; strength ratios between the dilatational (void) material and the plate ranged from 1 to 10−2 .

37

The computational cost was quantified by the total number of iterations required to arrive at the
solution (on a grid of 128 × 128 × 1). The median value of normal tractions across the plate–void
interface collected from all voxels in the immediate neighborhood of that interface (colored ring
in Fig. 3.6 top left) serves as quality measure of how well the chosen material description reflects
a void.
Figure 3.6 (bottom left) compares the evolution of void quality with biaxial extension among
both material descriptions for the two lowest stiffnesses, i.e., ratios of 10−2 and 10−3 . With
increasing hole expansion, the interface traction for the elastic material description asymptotically
reaches values slightly in excess of the flow stress of the surrounding plate (estimated by M g), i.e.,
the void turns from a relatively soft to an effectively hard inclusion. The level of strain at which this
transition occurs increases with reduced elastic stiffness of the void material. In contrast, interface
tractions of voids described by the dilatational material model are not (except at a plastic strength
ratio of 100 ) increasing up to the plate flow stress, but are limited by their intrinsic flow stress
value (as given by the respective plastic strength ratio). Therefore, variations in the elastic stiffness
have virtually no influence on the response of voids described by the dilatational material—except
for the numerical difficulty to integrate the solution with strain increments that exceed the elastoplastic transition strain (ratio of dilatational flow stress to stiffness) of the dilatational material (see
Fig. 3.6 bottom right). A comparison of the cost to solve the biaxial expansion using either of the
two tested void descriptions (Fig. 3.6 top right) shows that there is a marginal penalty (about a
factor of two) associated with using the dilatational material model. Nevertheless, this increased
effort is generally outweighed by the better quality, unless the void material undergoes only a small
amount of strain.

38

in-plane
von Mises
strain

simulated
domain

IPF(z)

yz
0

0.15

x

fully
periodic

quasiperiodic
(dilatational
rim)

free
surface
(dilatational
encasing)

experiment
[14]

free
surface
(elastic
encasing)
Figure 3.5: Influence of external boundary conditions on in-plane VON M ISES strain fields and lattice reorientation. Simulated geometries with fully periodic gauge section of dogbone (first row), quasi-periodic,
i.e. periodic along the thickness direction (second row), and free surfaces normal to x and z (third and last
row) compared to experimental results measured by DIC and EBSD (fourth row, from [14]). Gray (green)
semi-transparent volume in the simulated domain column uses the proposed dilatational (isotropically elastic) material model; regions of constant color reflect individual grains. Inverse pole figure (IPF, last column)
coloring of lab direction z is mapped on the undeformed configuration, except for the measured result (fourth
row).

39

103

101

Fxx = 1.2
Fyy = 1.2

elastic
0

10

dilatational

1

10-1

-1

10

plastic
strength
ratio

10-2

-2

10

10-3 -2
10

10-1

1

s
eed
exc in
ent stra
rem ion
inc ansit
tion c tr
gra asti
inte sto-pl
ela

median interface traction / M g

102

total # of iterations

104

1

10-2

10-4

stiffness contrast

(Fxx or Fyy)-1

Figure 3.6: Computational cost (top right) and void representation quality (bottom) of the elastic
and the dilatational material model for the case of biaxial expansion of a hole in a plastic plate (top
left) up to ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.2. Strain evolution of median interface traction (lower left, collected
from voxels immediately inside and outside of interface, i.e. blue and red ring at top left) indicates
the void is reaching the flow stress level (≈ M g) of the surrounding plate in case of the elastic and
strongest dilatational material (plastic strength ratio of 100 ); softer dilatational cases (darker red)
saturate much earlier at flow stress levels far below the plate flow stress (plastic strength ratio of
10−1 and 10−2 ). Median interface tractions collected at ⟨Fxx ⟩ = ⟨Fyy ⟩ = 1.2 (bottom right) confirm
marginal influence of elastic stiffness on void representation quality for the dilatational material.

40

CHAPTER 4
CONSTITUTIVE MODELS

A constitutive model relates the response of the material to an external load i.e. it relates a stress
(see stress measures [97]) to strain (see strain measures [97]). This relation is necessary to complement the static mechanical equilibrium and the strain compatibility condition. In comparison to
general constitutive models, crystal plasticity models take into account the anisotropic stress-strain
relation of the crystalline materials. The regular arrangement of atoms in crystalline materials at
solid state leads to anisotropic behavior. When stress applied to a crystalline material, the stress
is resolved on some particular slip systems ξ i.e. along a specific crystallographic direction (slip
direction) on some particular crystallographic planes (slip plane). The resolved component of the
applied stress is known as the Resolved Shear Stress (RSS) τ ξ is expressed as
ξ

ξ

τ ξ = S : m0 ⊗ n0

(4.1)

where mξ and nξ are the unit vectors of the slip directions and slip plane normal in the current
ξ

ξ

configuration, and m0 and n0 are the slip direction and slip plane normal vector in reference
configuration. Plastic deformation is then described by taking into account the total sum of shear
of all slip systems. By determining the shear rate γ̇ ξ (function of τ ξ ) with respect to the slip
ξ

ξ

system m0 ⊗ n0 , the plastic velocity gradient Lp as well evolution of plastic deformation can be
formulated.
ξ

ξ

Lp = F˙p Fp −1 = ∑ γ̇ ξ m0 ⊗ n0

(4.2)

ξ

In the next few sections of this chapter, two most commonly used ideologies to formulate shear-rate

γ̇ on these slip systems ξ are discussed.
41

4.1 Phenomenological Model
The phenomenological model is based on the assumption that when the resolved shear stress

τ exceeds a critical value τ0 (Critical Resolved Shear Stress - CRSS), plastic deformation takes
place. The resolved shear stress depends on the amount of the applied stress, the angle between
the loading direction and the plane normal, and the angle between the loading direction and the
slip direction Eq. (4.1). The exact relation is known as Schmid law [40]. In the present study,
the constitutive model used, is similar to the phenomenological description of Peirce et al. [1] for
ξ

cubic lattice symmetry. The internal state is captured through slip resistances τ0 per slip system ξ ,
which evolves with shear γ according to a relationship initially proposed by Hutchinson [98] and
later improved by Kalidindi et al. [99]

(
(
⎛
⎞
(
β (a
β ( (
(
τ (
τ
( (
ġξ = qξ β h0 ((1 − 0 (( sign ⎝1 − 0 ⎠ (γ̇ β ( .
β
β
(
g∞ (
g∞

(4.3)

The slip system interactions are controlled by the hardening matrix qξ β . The shear rate on each
system is given by

( (n
( ξ(
(τ (
ξ
γ̇ = γ̇0 (( (( sign τ ξ
( τ0ξ (

(4.4)

3
4
where γ̇0 is a reference shear rate, n the stress exponent, and τ ξ = S : mξ ⊗ nξ the resolved shear

stress. Table 7.1 summarizes the list of variables used in the model, with adjustable parameters h0 ,
ξ

ξ

a, τ0 , g∞ , and qξ β .

4.2 Dislocation density based model
Phenomenological models have been successful in predicting deformation behavior of polycrystalline materials [100, 99, 101, 102]. However, as shown in an earlier section, the phenomenological models fails to capture the intrinsic anisotropic deformation of BCC materials at the single
42

Table 4.1: Material variables for Peirce et al. [1] model; elastic constants Ci j , reference shear
ξ
ξ
rate γ̇0 , stress exponent n, initial and saturation slip resistance τ0 and g∞ , hardening parameters
h0 , qξ β , and a.
Elastic

Plastic

C11 C12 C44 γ̇0
GPa GPa GPa s−1

n

Strain hardening

ξ

ξ

τ0
g∞ a h0 qξ ξ
MPa MPa
MPa

qξ β

crystal level. Therefore, to better understand deformation in complex slip conditions, a crystal
plasticity model considering statistical properties of an ensemble of dislocations and their mutual
interactions is needed. These models are based on the idea that, for a dislocation to move it has to
overcome the energy barriers. As discussed earlier, there are two kinds of barrier: a short range one
which can be overcomed by thermal activation such as Peierls potentials, solid solution impurity
atoms, and long-range barriers like forest dislocations, grain boundaries and internal stresses of the
ξ

athermal type. The total resistance to slip is then defined as the sum of thermal (st ) and athermal
ξ

resistance (sa ) to dislocation motion.
ξ

ξ

sξ = st + sa

(4.5)

The average velocity of dislocations is given by an Arhenius type relationship,
ξ
ξ
vavg = v0 exp

5

−∆G
kB T

6

(4.6)

Where ∆G is the activation energy, v0 the reference velocity, kB the Boltzman constant and T the
temperature. This temperature dependence of the average velocity of dislocation can be substituted
ξ

back into Orowan equation (γ̇ ξ = ρ ξ bvavg ) to have a shear rate that is dependent on temperature
and applied stress.

γ̇ ξ =

⎧
⎪
⎪
⎪
⎨0

ξ

if τeff ≤ 0

⎪
;
<
⎪
⎪
⎩γ̇ ξ exp −∆G
0
kB T
43

ξ

if τeff > 0

(4.7)

( (
( (
ξ
ξ
ξ
ξ
where γ̇0 = ρ ξ bv0 is the reference strain rate and τeff = (τ ξ ( − sa the effective resolved shear
stress acting on slip system ξ . The activation energy barrier can be calculated using an approach
suggested by Kocks [103].
⎧
⎛
⎞ p ⎫q
ξ
⎨
⎬
τeff
⎝
⎠
∆G = ∆G0 1 −
ξ
⎩
⎭
s

(4.8)

t

Where ∆G0 is the activation energy for obstacle free dislocation motion. The parameters p and q
can be fitted to get an accurate description of the activation profile. Combining these ideas one can
arrive at an expression of shear rate that is temperature dependent.
⎧
⎡
⎛
⎞ p ⎫q ⎤
ξ
⎨
⎬
τ
−∆G0
ξ
⎦
γ̇ ξ = ρ ξ bv0 exp ⎣
1 − ⎝ eff ⎠
ξ
⎭
kB T ⎩
s

(4.9)

t

The primary state variables characterizing the microstructural state are the dislocation densities

ρ ξ on each slip system. Plastic flow at lower length scales, occurs as a succession of dislocation
avalanches with each avalanche triggered by the dislocation interactions [16]. When two dislocation segments moving in non-coplanar slip planes intersect each other, a junction is formed, if the
Frank criterion is satisfied. The two segments are then zipped together to form a sessile segment.
These sessile segments are destroyed when the work done by the local resolved stress τ on one
arm of the parent dislocations, of length lu , is sufficient to unzip the junction. The amount of work
done in formation of the junction per unit length is W j = τ blu , where b is the Burgers vectors of the
parent dislocations of length lu . The change in elastic energy due to the newly formed dislocation
segment can also be written as ∆E = α µ b2 , where α describes the strength of the junction formed.
Comparing work done with the change in elastic energy, α can be expressed as,

α=

τ lu
µb

44

(4.10)

Taylor [104] first proposed a relation between athermal resistance to deformation sa and the total
dislocation density ρ utilizing the dislocation interaction α .
√
sa = µ bα ρ

(4.11)

where µ is the isotropic shear modulus and b is the Burgers vector. This expression Eq. (4.11) was
latter vectorized by Franciosi et al. [105].
ξ
sa = µ b

@

∑ Aξ β ρ β

(4.12)

β

where the interaction coefficients Aξ β are dimensionless components of the so-called interaction
matrix (Table 4.2). They measure the strength of the interaction between the primary slip system
A
ξ and other slip systems β . The strengthening coefficients αξ β = Aξ β are often used to make
comparisons between experiment and theoretical estimates [106]. These interactions coefficient
are calculated using DDD simulations. Queyreau et al. [2] reported the interaction coefficients for
BCC α -Fe, and since then it has been used by several authors to model dislocation behavior in
other BCC materials [107, 108]. There are 6 distinct kinds of interaction possible for dislocations
in BCC lattice as shown with their corresponding values for α -Fe in Table 4.3. The values reported
in Table 4.3 were calculated at a reference dislocation forest density (ρ re f ) of 2 × 1013 m−2 . The
nature of these dislocation interactions are inherently elastic, as it removes a significant amount
of line energy from the crystal. This behavior has been well documented in experimental works
of Basinski and Basinski [109] in FCC metals. Basinski and Basinski [109] also showed that the
unzipping action of a junction is influenced by line tension and curvature of the dislocation line,
which in turn depends on the forest density. Hence, a logarithmic factor Eq. (4.13) needs to be

45

{1 1 0}

1
2
6
6
5
4
4
3
4
3
5
4

2
1
6
6
4
3
5
4
5
4
4
3

6
6
1
2
4
5
3
4
4
5
3
4

6
6
2
1
3
4
4
5
3
4
4
5

{1 1 2}

Table 4.2: Dislocation interaction matrix in BCC lattice.

6
6
4
3
3
4
6
6
4
3
6
6

6
6
3
4
4
3
6
6
3
4
6
6

4
3
6
6
6
6
3
4
6
6
4
3

3
4
6
6
6
6
4
3
6
6
3
4

{1 1 0}
5
4
4
3
1
2
6
6
3
4
5
4

4
3
5
4
2
1
6
6
4
5
4
3

4
5
3
4
6
6
1
2
5
4
3
4

3
4
4
5
6
6
2
1
4
3
4
5

4
5
4
3
3
4
5
4
1
2
6
6

3
4
5
4
4
5
4
3
2
1
6
6

5
4
3
4
5
4
3
4
6
6
1
2

4
3
4
5
4
3
4
5
6
6
2
1

6
6
4
3
3
4
6
6
3
4
6
6

6
6
3
4
6
6
3
4
6
6
3
4

4
3
6
6
4
3
6
6
6
6
4
3

3
4
6
6
6
6
4
3
4
3
6
6

3
6
4
6
6
4
6
3
4
6
3
6

4
6
3
6
6
3
6
4
3
6
4
6

6
3
6
4
4
6
3
6
6
4
6
3

6
4
6
3
3
6
4
6
6
3
6
4

3
6
6
4
4
6
6
3
6
4
3
6

4
6
6
3
3
6
6
4
6
3
4
6

6
3
4
6
6
4
3
6
4
6
6
3

6
4
3
6
6
3
4
6
3
6
6
4

1
5
6
6
5
6
6
3
5
6
3
6

5
1
6
6
6
5
3
6
6
5
6
3

6
6
1
5
6
3
5
6
3
6
5
6

6
6
5
1
3
6
6
5
6
3
6
5

{1 1 2}
3
4
6
6
6
6
4
3
4
3
6
6

4
3
6
6
4
3
6
6
6
6
4
3

6
6
3
4
6
6
3
4
6
6
3
4

6
6
4
3
3
4
6
6
3
4
6
6

4
3
6
6
4
3
6
6
6
6
4
3

3
4
6
6
6
6
4
3
4
3
6
6

6
6
4
3
3
4
6
6
3
4
6
6

6
6
3
4
6
6
3
4
6
6
3
4

5
6
6
3
1
6
5
6
5
3
6
6

6
5
3
6
6
1
6
5
3
5
6
6

6
3
5
6
5
6
1
6
6
6
5
3

3
6
6
5
6
5
6
1
6
6
3
5

5
6
3
6
5
3
6
6
1
6
6
5

6
5
6
3
3
5
6
6
6
1
5
6

3
6
5
6
6
6
5
3
6
5
1
6

6
3
6
5
6
6
3
5
5
6
6
1

Table 4.3: Types of possible interaction in BCC lattice [2]
ID

Interaction Type

1
2
3
4
5
6

self
coplanar
collinear
mixed asymmetric
mixed symmetric
edge junction

considered for these effects.

Interaction Coefficient
0.0068 ± 0.001
0.0068 ± 0.001
0.5300 ± 0.040
0.0390 ± 0.005
0.0630 ± 0.006
0.0360 ± 0.006

# &
$
β
1
ln b Aξ β ρf
ref
Aξ β = Aξ β # &
$
β
1
ln b Aξ β ρref
46

(4.13)

ξ

ξ

ρf is the existing dislocation density and ρref is the reference dislocation density at which the
interaction coefficients are calculated by DDD.

4.2.1

Influence of elastic anisotropy on dislocation interaction strengths

In pure BCC single crystals above the knee temperature, Tk ı.e athermal regime, the reactions between dislocation segments offers to most of the resistance to dislocation motion. Therefore, the
interaction coefficients defining the junctions strengths are critical in reproducing the macroscopic
deformation behavior. Initially these interaction coefficients were considered to be dependent only
on the crystal structure, until recently it was shown by Madec and Kubin [3] that this is not necessarily true. Two effects namely, “Poisson ratio” and “asymmetric junction configurations” influence the strength of interactions between dislocation segments. These two effects are more
prominent in BCC metals than in FCC metals. In their study, strengths of dislocation interaction
among ⟨1 1 1⟩{1 1 0} and ⟨1 1 1⟩{1 1 2} slip systems were determined using dislocation dynamics
simulations at high homologous temperatures for five BCC transition metals. They predicted that
the use of accurate material-dependent coefficients would notably improve the predictive ability of
current models for strain hardening.
Madec and Kubin [3] made use of Scattergood and Bacon’s shear moduli and Poisson’s ratios
[110] to approximate the behavior of curved dislocations in an anisotropic medium. The results
concluded that in a BCC metal (compared to FCC), the strength of junctions is mainly governed
by the orientation dependence (angle between the interacting slip planes) of the line energies. This
angle affects the magnitude of short-range interactions between dislocation segments influencing
the limit between attraction and repulsion for non-coplanar dislocation reactions. The value of
the Poisson ratio affects all the strengthening coefficients to various extents ranging from small
47

Table 4.4: Modified dislocation interaction matrix in BCC lattice taking into account asymmetric
junction configurations [3].

2
1
2
5
6
7
4
3
4
7
6
5

2
2
1
5
7
6
7
5
6
4
4
3

3
4
4
1
2
2
6
7
5
6
5
7

8
9
9
10
12
12
15
16
19
16
15
19

9
9
8
16
15
19
12
10
12
16
19
15

9
8
9
16
19
15
19
16
15
10
12
12

10
12
12
8
9
9
16
15
19
15
16
19

{1 1 2}

{1 1 0}

1
2
2
3
4
4
6
5
7
6
7
5

{1 1 0}
5
6
7
2
1
2
7
6
5
4
3
4

5
7
6
2
2
1
4
4
3
7
5
6

6
5
7
6
7
5
1
2
2
3
4
4

4
3
4
7
6
5
2
1
2
5
6
7

7
5
6
4
4
3
2
2
1
5
7
6

6
7
5
6
5
7
3
4
4
1
2
2

7
6
5
4
3
4
5
6
7
2
1
2

4
4
3
7
5
6
5
7
6
2
2
1

8
9
9
11
17
17
18
13
14
18
14
13

9
9
8
13
14
18
14
13
18
17
17
11

9
8
9
13
18
14
17
11
17
14
18
13

11
17
17
8
9
9
18
14
13
18
13
14

16
15
19
9
9
8
16
19
15
12
10
12

16
19
15
9
8
9
10
12
12
19
16
15

19
15
16
19
16
15
9
9
8
12
12
10

12
12
10
15
16
19
9
8
9
15
19
16

15
19
16
12
10
12
8
9
9
19
15
16

19
16
15
19
15
16
12
12
10
9
9
8

15
16
19
12
12
10
15
19
16
9
8
9

12
10
12
15
19
16
19
15
16
8
9
9

1
2
2
20
22
22
23
21
24
21
23
24

2
1
2
21
24
23
24
21
23
20
22
22

2
2
1
21
23
24
22
20
22
21
24
23

20
22
22
1
2
2
21
23
24
23
21
24

{1 1 2}
13
14
18
9
9
8
17
17
11
14
13
18

13
18
14
9
8
9
14
18
13
17
11
17

14
13
18
17
17
11
9
9
8
13
14
18

17
11
17
14
18
13
9
8
9
13
18
14

18
13
14
18
14
13
8
9
9
11
17
17

17
17
11
14
13
18
13
14
18
9
9
8

14
18
13
17
11
17
13
18
14
9
8
9

18
14
13
18
13
14
11
17
17
8
9
9

21
24
23
2
1
2
20
22
22
24
21
23

21
23
24
2
2
1
21
24
23
22
20
22

23
24
21
22
20
22
1
2
2
24
23
21

22
22
20
23
21
24
2
1
2
23
24
21

24
23
21
24
21
23
2
2
1
22
22
20

22
20
22
23
24
21
24
23
21
1
2
2

23
21
24
22
22
20
23
24
21
2
1
2

24
21
23
24
23
21
22
22
20
2
2
1

to substantial. Increasing Poisson’s ratio strengthens the collinear interactions, while decreasing
it, leads to less stable junctions due to higher line energies. Taking into consideration of all these
influences, [3] has proposed a modified interaction matrix (see Table 4.4). The corresponding
values of interaction strengths for Nb calculated by them with Scattergood and Bacon’s anisotropic
approximation are presented in Table 4.5.

48

Table 4.5: Types of possible interaction in BCC lattice with asymmetric junction configurations
[3].

4.2.2

1
2

self
collinear 60◦

0.20
0.82

{1 1 0}-{1 1 0}
interactions

3
4
5
6
7

coplanar
J1 0◦ (G)
J1 70.53◦ (G)
J2
J3

0.20
0.24
0.20
0.31
0.28

{1 1 0}-{1 1 2} or
{1 1 2}-{1 1 0}
interactions

Interaction Coefficient
(G - glissile junctions)

8
9
10
11
12
13
14
15
16
17
18
19

Collinear 90◦
Collinear 30◦
J4 0◦ (G)
J4 70.53◦ (G)
J5 0◦ (G)
J5 70.53◦ (G)
J6 29.50◦
J6 58.52◦
J7 58.52◦
J7 79.98◦
J8 35.26◦
J8 90◦

0.90
0.64
0.25
0.20
0.20
0.17
0.20
0.20
0.27
0.25
0.26
0.27

{1 1 2}-{1 1 2}
interactions

ID Interaction Type

20
21
22
23
24

J9
J10 28.56◦
J10 72.98◦
J11
J12

0.28
0.23
0.23
0.26
0.19

Dislocation storage and recovery framework

The dislocation storagerecovery framework is a rather simple concept which defines the evolution
of the total dislocation density in a uniform microstructure with time or strain, under the influence of mechanisms such as dislocation strengthening, storage, and dynamic recovery. Such a
formulation is particularly suitable for multiscale modeling because of its kinematic character and
its potential to reproduce stressstrain curves with all deformation stages in a continuous fashion.

49

Teodosiu et al. [111], Kocks and Mecking [112] proposed a generalized version of storage and
recovery framework that treats uniform dislocation densities on each slip system. This framework
assumes that the local microstructure mostly consists of dislocations immobilized (stored) by interaction with other dislocations. For each slip system ξ , the stored density ρ is much larger than
the instantaneous mobile density and is therefore assimilated to the total density. It is also assumed
that there are always enough mobile dislocations to sustain the imposed strain rate and compensate
for the losses by storage and annihilation. Thus, the continuous storage of mobile dislocations is
accounted for, and the mobile density does not explicitly appear in the model. The most important
part of the framework is defining the Mean Free Path (MFP) of dislocations. A segment of average
length l, moving by an incremental step δ l, in a volume V produces a shear strain δ γ = blVδ l . If it
is immobilized, the stored density is increased by δ ρ = Vl . The mean free path Λ is the average
length traveled by mobile dislocation segments during a strain increment. The storage probability
δ l/Λ

is then defined as δ l/Λ and the storage rate is δ ρ /δ γ = Λb . As the length of the dislocation segment changes when it moves, the mean free path is a virtual quantity that has to be defined for each
strain increments. Thus the net storage rate for a slip system ξ during a strain increment can be
expressed in terms of the shear strain rate γ̇ ξ by equation,
$( α (
( dγ (
sα
a
α
(
(
− yhkl ρ
( dt (
µ bKhkl
%
'
ξ
( (
1
s
( (
a
ρ̇ ξ =
− yhkl ρ ξ (γ̇ ξ (
b µ bKhkl

1
dρ α
=
dt
b

#

(4.14)

where Khkl is a mean free path coefficient and yhkl an orientation-dependent factor responsible for
dynamic recovery. The first term on the right-hand side of Eq. (4.14) is the dislocation multiplication or storage rate. The second term accounts annihilation of screw dislocations by thermally
activated cross-slip and its subsequent rearrangements. When dynamic recovery is not active i.e.,
50

in stage II, the model yields a constant strain-hardening rate θII /µ = α 2 /2KΛ. Introducing dynamic recovery yields the strain hardening rate in stage III.
The above dislocation recovery and storage framework requires parameters, which depends
on evolution of dislocation substructure with strain. For an accurate description of a dislocation
multiplication and annihilation rates, these parameters are calculated using DDD. The following
sections briefly discusses the DDD works of Madec and Kubin [113] for calculation of dislocation
mean free path. In their DDD framework, the dislocation densities stored for each slip system
were considered as the only state variable. They hypothesized that since both parent segments and
junction lines have similar interactions with incoming mobile segments, it is convenient to lump
junction densities into the densities stored in the active slip systems. Another advantage of this
convention is that, it indirectly considers higher order reactions or multi-junctions between dislocation lines [114, 113], which frequently forms from an interaction production of two dislocation
lines reacting with a third dislocation segment. This convention is applied in the following ways :
• mobile dislocations of system ξ may react with the stored forest density on system β or vice
versa (see Fig. 4.1).
• junction formation for both cases can be denoted as ξ → β or β → ξ , where the first identifier
indicates the active or gliding system.
• junction segments of type ξ → β or β → ξ are then assimilated into stored densities in slip
systems ξ or β respectively.
This convention also restricts stored densities of weakly active slip systems to a small number.
Once a junction line of type ξ → β is formed and lumped into slip system ξ , it cannot constitute a

51

forest obstacle for further incoming segments in system β , but it can react with all other forest slip
systems of ξ . Using this convection the rate of multiplication of dislocations in a slip system can
be defined.

nξ

ξ

nβ

β

Figure 4.1: A schematic figure showing dislocation reaction between segments on active slipsystem ξ (blue) and a forest system β (red).

4.2.3

Dislocation multiplication

Consider a scenario in which the dislocation lines in slip system ξ interacts with an inactive forest
dislocations ρ β of average length l¯β on slip system β . Then as per the convention defined earlier,
all the junctions formed during the glide of dislocations ξ are assigned to slip system ξ . Thus the
total amount of dislocation segments in ξ would contain two parent segments and one junction
segment (Fig. 4.1). During a shear increment d γ ξ , the area swept by the dislocation is dAξ =

52

d γ ξ V /b. Let dN ξ →β to be the number of interactions occurring in this swept area dAξ and ϕ
be the geometrical factor signifying the projected length of forest dislocations along the directions
perpendicular to the active slip plane. Then the total length of the dislocation line stored in the
volume dAξ ϕ l¯β after reaction with forrest density ρ β = dN ξ →β /ϕ dAξ is dN ξ →β l¯β . Therefore,
the number of interactions is given as,
dN ξ →β = ϕρ β dAξ =

ϕρ β d γ ξ V
b

(4.15)

All interactions cannot form stable junctions under stress. Out of the total number of interactions,
only a fraction f of them actually, forms stable junctions. The fraction f of stable junctions
A
depends on the interactions strength Aξ β between slip system ξ and β . Thus, the number of
stable junctions is then given as,

ξ →β

dNjct

=

A
ϕ Aξ β ρ β d γ ξ V

b
3A
4
V
ξ →β
dNjct = ϕ f
Aξ β ρ β d γ ξ
b

(4.16)

For each stable junction of length l¯ξ , the dislocations stored in the system is incremented by l¯ξ /V .
ξ →β

So for dNjct

interaction events, the total increase in dislocation density d ρ ξ is given as,
dρ ξ =

ξ →β
l¯ξ dNjct

V

(4.17)

Combining Eq. (4.17) with Eq. (4.16) for all slip systems, one can arrive at an expression for
dislocation multiplication rate as follows.
dρ ξ
dγ ξ

⎛

p l¯ξ
= 0 ⎝
b

Nslip

∑

β =1,β ̸=ξ

⎞
A
Aξ β ρ β ⎠

(4.18)

where Nslip is the total number of slip systems and p0 = ϕ f . l¯ξ is the average stored dislocaA
tion length which can be approximated as l¯ξ = k0 / ∑β Aξ β ρ β with k0 being a dimensionless
53

constant. Since the average length of stored dislocation segments is directly proportional to the
dislocation density, therefore the ratio should be constant [16].
ξ

ρjct
ρξ

=

ξ
l¯jct

l¯ξ

= κ0

(4.19)

This ratio is important as it will be used later to simplify the dislocation multiplication rate.
So far it was assumed that the forest dislocations segments were inactive. Relaxing this conditions results in a modification to the junction densities and the average length of stored dislocation
content. A dislocation segment moving on a forest system β interacts with active slip system ξ
and can form junctions on slip system β . Thus a fraction of stored dislocation content on ξ is
transferred to system β i.e.
ξ

β →ξ

ρ ξ = ρ0 − ρjct

(4.20)

ξ

where ρ0 is the dislocation density on the active slip system ξ for non-active forest. Similarly, the
average dislocation segment length can be modified as
⎛

ρξ

⎞

k0
⎝
⎠
3
4
β →ξ
Nslip
ξ
ξ →β
ρ + ρjct
∑β =1 Aξ β ρ β + ρjct

l¯ξ = &

(4.21)

Substituting Eq. (4.21) in Eq. (4.18), the dislocations multiplication rate in the case of active forest
systems can be derived.
dρ ξ
dγ ξ

=

⎛

ρξ

⎞⎛

1
p0 k0
⎝
⎠⎝ ∑
&
4
3
b
Nslip
ξ + ρ β →ξ
ξ
→
β
ρ
β =1
jct
∑β =1 Aξ β ρ β + ρjct

A

⎞⎛

A ξ β ρ β ⎠ ⎝1 −

β →ξ

ρjct

ρξ

⎞
⎠

(4.22)

Devincre et al. [16] proposed a simplified version of Eq. (4.22) using the parameters described
earlier.

$
#
1
dρ
p0 k0
ρforest
κ
=
=
1−
dγ
bΛ
b √ρtotal (1 + κ )3/2
n−1
54

(4.23)

where n is the number of slip systems, κ = κ0 /(1 − κ0 ) i.e. the fraction of junction density to the
fraction of non-junction density and Λ is the dislocation mean free path. Lumping all constant
coefficients into a dimensionless and orientation-dependent coefficient Khkl , the mean path Λ can
eventually be set in the simple form.
K
Λ = √ hkl
āρtotal

(4.24)

The coefficient Khkl depends on the average interaction strength and on the set of constants p0 , k0
and κ0 (or κ ).
4.2.4

Dislocation annihilation

Critical annihilation distance: Two screw dislocations of opposite sign gliding in parallel slip
planes separated by a distance smaller than a critical distance ys can mutually annihilate by crossslip. Although this quantity is widely used in many models, there is presently no accurate theoretical prediction for it, as it derives from a thermally activated core mechanism. Furthermore,
the density of screw dislocations is a small, but as yet unknown, fraction f of the total dislocation
density because screw dipoles annihilate more easily than edge dipoles. As the present model is
written in terms of total densities, Eq. (4.14) actually includes a length yhkl = f ys . The value of
yhkl is adjusted to yield the experimentally observed critical stress for the onset of the dynamic
recovery stage. An example is shown in [16] for copper at room temperature where, y{1 2 3} =
0.5 nm. From the classical work of [115] one has ys ≈ 50 nm in persistent slip bands of copper
crystals cycled at 300 K. In that case, the ratio of screw to edge or total densities is about 1/100
and the value determined for y looks reasonable.
The critical annihilation distance is also orientation-dependent, which entails an orientation
dependence of the onset of stage III in the stress-strain curves. It was shown by Kubin et al. [116]
55

that this effect is of geometrical origin and arises from the respective orientation dependencies of
the Schmid factors in the slip and cross-slip planes. The results obtained in Kubin et al. [116] are
expanded to yield a value of the quantity yhkl for any orientation [h, k, l] of the loading axis.
Table 4.6 summarizes the list of variables used in the model, with adjustable parameters Khkl ,
{1 1 0}

yhkl , ρ0

{1 1 2}

, ρ0

ξ

, and st .

Table 4.6: Material variables for dislocation density based model
Elastic
C11

C12

Plastic
C44

GPa GPa GPa

ξ

ρ0

m−2

ξ

Strain hardening

st

v0

b

∆G0

T

MPa

m sec−1

Å

J

K

56

p

q

Aξ β

Khkl

yhkl
nm

CHAPTER 5
RESULTS AND DISCUSSION

The outcome of the simulated deformation behavior using the constitutive framework discussed
in Chapters 3 and 4 is presented in this chapter. This chapter discusses details about the crystal
plasticity constitutive model calibration, sensitivity of the primary model parameters (ρ ξ , Aξ β )
and identification of relevant state parameter evolution that is crucial to capturing the physics of
the deformation process.

5.1 Calibration of Crystal Plasticity constitutive model
The mechanical response of a single crystal undergoing uniaxial tension (using “standard” dogbone samples) exhibits marked spatial gradients between the central and peripheral gauge sections.
Capturing the fidelity through full-field simulations requires a well suited constitutive description
with proper combinations of primary state variables (input parameters). These input parameters
are obtained by optimizing the simulated output with respect to a handful of experimental data.
In the present study, the crystal plasticity model is calibrated by optimizing the simulated uniaxial stress-strain curves with respect to experimental data collected by [12, 75]. The optimization
process used in this study has been discussed in elaborate details in Chakraborty and Eisenlohr
[15].
One of the problems faced during inverse optimization of constitutive models is the computational cost involved. The simulated stress-strain response is frequently calculated based on various
combinations of input parameters within a solution space until the desired accuracy is achieved.
The error in each minimization is calculated using an objective function. In this particular case, an

57

150

CPFFT
tension
simulation

reference
(experiment)

σ / MPa

100

constitutive
parameters

Εσ-ε
adjust
50
yes
0

0

0.1

0.2

0.3

0.4

obj. function
meets
tolerance

no

optimizer

0.5

ε

Figure 5.1: Left : Optimization function minimizing the area between two curves. The gray line represents
the standard curve, which the dashed black line represents the trial curve. Right : Flow diagram of the
optimization process used by Chakraborty and Eisenlohr [15].

objective function was designed to minimize the error in the area under the two stress-strain curves
(see Fig. 5.1 left). The computationally expensive step in this optimization process is the CPFFT
calculation step. Therefore it is essential to identify the boundary conditions which are most important to mimic the behavior of single crystal samples with a much smaller geometry (low mesh
resolution).
The second column of Table 5.1 shows the geometries used to simulate the deformation behavior. To reduce the CPFFT computation cost, a rectangular block of gauge section with reduced
mesh resolution is subjected to the same boundary conditions as the dogbone sample itself. A
buffer layer of soft, dilatational and low stiffness material (see Chapter 3) is added surrounding the
dogbone and the cube geometry to mimic free surface boundary conditions in the spectral solver.
The resolution of the reduced geometry is varied from 4 × 4 × 4 to 12 × 12 × 12 for different orientations. The corresponding average computation time for these geometries is reported in the last
column. The resulting stress (strain) response of each of these simplified reduced geometries is
then compared to the stress (strain) response of the gauge section of the dogbone geometry us58

Table 5.1: Comparison of CPFFT computation time between simplified reduced geometry of gauge
section and a proper dogbone shaped geometry. Material described by crystal plasticity law is
shown in blue and the dilatational buffer layer is shown in light gray.
Geometry

Resolution

Computation
time [min]

reduced

4×4×4
6×6×6
8×8×8
10 × 10 × 10
12 × 12 × 12

1.617
4.883
10.270
23.130
38.100

dogbone

35 × 160 × 26

2263.300

ing Kernel Density Estimation (KDE) plots. The KDE plots of all the components of the stress
(strain) behavior for both reduced and gauge section of dogbone geometry is shown in Fig. 5.2.
The top row shows the correlation of the normal stress components and the bottom row shows the
correlation of shear stress components. It is noted that the normal stress (strain) response of these
reduced geometries along the loading direction (top row middle column Fig. 5.2) is remarkably
similar (Pearson correlation r = 0.97) to that of the gauge section of dogbone geometry. However,
the other transverse normal components and the shear components vary from weak correlation
(Pearson correlation r = 0.5) to no-correlation (Pearson correlation r = 0.0053). Given that, the absolute magnitudes of the transverse and shear stress components are way below the flow stress of
the material, deviations in these components have small effects on the resultant slip behavior and
therefore can be ignored. Thus the reduced geometry can be used to approximate the deformation
of the gauge section of dogbone sample for faster optimization of input parameters.

59

1.5

1

0.5

0

0

0.5

1

1.5

80
60
40
20
0

2

2

σzz / MPa (dogbone geom.)

100

σyy / MPa (dogbone geom.)

σxx / MPa (dogbone geom.)

2

0

σxx / MPa (reduced geom.)

60

80

1

0.5

0.5

1

1.5

σxy / MPa (reduced geom.)

0

2

1

1.5

2

2

1.5

1

0.5

0

0.5

σzz / MPa (reduced geom.)

σyz / MPa (dogbone geom.)

1.5

0

0.5

0

100

2

σxz / MPa (dogbone geom.)

σxy / MPa (dogbone geom.)

40

1

σyy / MPa (reduced geom.)

2

0

20

1.5

0

0.5

1

1.5

2

σxz / MPa (reduced geom.)

1.5

1

0.5

0

0

0.5

1

1.5

2

σyz / MPa (reduced geom.)

Figure 5.2: Comparison of stress components of simplified reduced geometry and gauge length
section of dogbone geometry in terms of kernel density estimation maps. The top row shows the
correlation of the normal stress components and the bottom row shows the correlation of shear
stress components.

5.2 Phenomenological constitutive model results
Phenomenological models have been long used in crystal plasticity modeling to predict deformation response of materials, particularly FCC metals. As discussed in the earlier chapter, the
essence of a phenomenological model relies on the CRSS for a particular slip-system. Although
much work on single crystal plasticity of BCC metals has been devoted to α -Fe, Mo, W and Ta,
very little work has been done on Nb. Mapar et al. [117] used a rate-independent model with
exponential and dynamic hardening rules to predict strain hardening in Nb. As a first attempt, a
rate-dependent phenomenological crystal plasticity model (see Chapter 4) is used to predict the
strain hardening. The important model parameters τ0 , τsat and h0 are optimized using the frame60

work discussed in earlier section.

100

100

σ / MPa

σ / MPa

150

50

0

Q
50

0

0.1

0.2

0.3

0.4

0

0.5

150

100
R

50

0

0.1

0.2

ε

0.3

0.4

0

0.5

150

100

100

S
50

0

0.1

0.2

ε

150

σ / MPa

P

100

150

σ / MPa

150

0.3

0.4

0

0.5

0

0.1

0.2

ε

0.3

0.4

0.5

0.4

0.5

ε

150

150

50

0

0.1

0.2

0.3

ε

0.4

0.5

0

σ / MPa

100

σ / MPa

T
50

0

100

σ / MPa

σ / MPa

V
W

50

0

0.1

0.2

0.3

0.4

0.5

0

X
50

0

0.1

ε

0.2

0.3

ε

0.4

0.5

0

0

0.1

0.2

0.3

ε

Figure 5.3: Comparison of experimentally observed axial stress-curves with the ones predicted
using phenomenological crystal plasticity model as function of orientation. Solid lines indicate the
measured stress-strain curves while the predicted ones are shown in dashed lines.

Figure 5.3 shows the resulting stress-strain curves predicted by the phenomenological model
along with their experimentally measured ones. The model seems to work only for the orientations
which do not show marked changes in the hardening rates in between stage I and stage II. However,
for most of the orientations, the model fails to capture the change in hardening behavior over the
strain range. In almost all cases the rate of hardening predicted by the model matches closely
with the stage II hardening rate, but no sign of stage I deformation is observed in the predicted
behavior. The initial CRSS τ0 , saturation CRSS g∞ (for two slip families {1 1 0} and {1 1 2}) and
the hardening slope h0 obtained from the best-fitted simulations for each of these orientations are
shown in Fig. 5.4. A relatively small variation in hardening slope h0 and τ0 is observed along with
a wide variation in g∞ over the range of crystal orientation.
These observations suggests that for predicting a variable hardening rate with consistent CRSS

61

Best Fitted Values / MPa

450

300

150

0

Figure 5.4: Adjustable phenomenological model parameters (initial CRSS τ0 saturation CRSS
g∞ and hardening slope h0 ) for pure Niobium obtained from stress-strain curve optimization of
8 single crystal orientations. The crosses represent a parameter for a particular orientation. The
mean of the parameters for all orientations is represented by a solid dot.
over the orientation space one needs to treat the hardening behavior in different stages separately.
Such models with multiple-hardening laws for various stages has been investigated in the works of
[118, 119].

5.3 Dislocation density based crystal plasticity model results
The dislocation density based crystal plasticity model described in section 4.2 is used to predict deformation behavior of the dogbone geometries for the orientations reported by Baars [12].
Table 5.3 presents the underlying values of the model parameters with adjustable variables.

62

Table 5.2: Constitutive material parameters of pure Nb used for the crystal plasticity simulations
of uniaxial tension. The asterisk (*) indicates fitting parameters.
Parameter

Value

C11
C12
C44
µ

246.5 GPa
134.5 GPa
28.73 GPa
39.6 GPa

{1 1 0}

ρ0

{1 1 2}
ρ0

b
v0
∆G0
T

*
3.3 Å
−4
10 m s−1
2.72 · 10−19 J
300 K

st
p
q
Khkl
yhkl
Aξ β

8.5 MPa
0.85
1.27
*
*
Table 4.5

ξ

5.3.1

*

Slip activity maps as a function of orientation

To determine the slip-activity as a function of orientation, the implemented dislocation density
based crystal plasticity model is used assuming ⟨1 1 1⟩{1 1 0} and ⟨1 1 1⟩{1 1 2} as possible slipsystems populated with equal amounts of dislocations. Based on the net accumulated shear Γacc
at the onset of slip (ε ∼ 0.002), active slip family is mapped as a function of orientation. Fig. 5.5
shows the dominant slip-family determined a function of orientation. Similar slip behavior in
BCC metals was also reported in the works of Calnan and Clews [120]. However, it is observed
that with continued plastic deformation the active slip family changes drastically and is mostly of
⟨1 1 1⟩{1 1 2} type for majority of the orientations in the standard stereographic triangle except near
the orientations with very high Schmid factor on ⟨1 1 1⟩{1 1 0}.

63

h1 1 1i

h1 1 1i

h1 1 2i

h0 0 1i

h1 1 2i

h0 1 1i h0 0 1i

h0 1 1i

Figure 5.5: Active slip family based on net accumulated shear Γacc as a function of orientation.
The contribution of a slip family for each of the orientations is represented by a color varying from
blue for pure {1 1 2} to red for pure {1 1 0} slip. Left: active slip family maps at ε ∼ 0.002. Right:
active slip family maps at ε ∼ 0.050.
5.3.2

Influence of slip family on spatially resolved deformation behavior

The reduced geometry described earlier is computationally efficient and robust for model calibration, however, it cannot reproduce spatially varying heterogeneous deformation of the dogbone
shaped sample under tension. For this purpose, CPFFT simulations were done on a dogbone
geometry for all the eight orientations with the optimized dislocation density model’s optimized
parameter set. Three different slip conditions were considered in which the deformation on the slip
systems was restricted to (a) only {1 1 0} slip planes, (b) {1 1 0}, {1 1 2} slip planes and (c) only
{1 1 2} slip planes.
The local inverse pole figure, von Mises strain and von Mises stress field maps are shown in
Figs. 5.6 to 5.8 respectively for all the three slip conditions mentioned above. Some interesting
observations were made. In a few of the orientations (P, T ), the active slip family determined
the amount of necking in the latter stages of deformation. Moreover, marked differences were
observed in the final orientation of the necked region. The same behavior has also been detected in
the local von Mises stress fields. In the case of the other orientations the stress localization patterns
64

orientations
P

Q

R

S

T

V

W

X

{1 1 2}

{1 1 0} {1 1 2}

{1 1 0}

slip
family

Figure 5.6: Local ⟨0 1 0⟩ lattice reorientations maps of the deformed dogbone samples. Rows in
descending order represent the simulation results obtained assuming {1 1 0}, {1 1 0}-{1 1 2} and
{1 1 2} slip family respectively.
(alternating regions of high and low stress) seems to be dependent the active slip family. For some
of the orientations, stress localization is only observed only when slip is restricted to single slip
family. Overall, the final shape of the deformed geometry was found to be similar across all three
cases for most of the orientations.

65

orientations
P

Q

R

S

T

V

W

X

{1 1 2}

{1 1 0} {1 1 2}

{1 1 0}

slip
family

0

0.5

Figure 5.7: Local von Mises strain εvM maps of the deformed dogbone samples. Rows in descending order represent the simulation results obtained assuming {1 1 0}, {1 1 0}-{1 1 2} and {1 1 2}
slip family respectively.
5.3.3

Model sensitivity to initial dislocation density

The primary state variable for the dislocation density based model is the dislocation density content
on individual slip-systems. Given that it is impractical to measure dislocation content on individual
slip-systems experimentally, it becomes essential to know the bounds of the CPFFT simulations
up to which the predicted deformation response can be captured with reasonable accuracy.

66

orientations
P

Q

R

S

T

V

W

X

{1 1 2}

{1 1 0} {1 1 2}

{1 1 0}

slip
family

0

0.2 GPa

Figure 5.8: Local von Mises stress σvM maps of the deformed dogbone samples. Rows in descending order represent the simulation results obtained assuming {1 1 0}, {1 1 0}-{1 1 2} and {1 1 2} slip
family respectively.
The influence of initial conditions on deformation response was tested by distributing dislocation content inhomogeneously among individual slip-systems. The distribution of dislocation
densities followed a log-normal distribution with a mean dislocation density, ρmean and variance

σρ . The range of the distribution was varied by changing the variance. Three variations were tested
i.e σρ = 0.5, 1.0 and 2.0 around a mean dislocation density of ρmean = 2 × 1013 . For each distribu67

Dislocation Density, ρ

1014
1013
1012
1011
1010

0

3

6

9

12

# of Trial
103

[ 1-1 1] ( 0 1 1)
[-1-1 1] ( 0 1 1)
[ 1 1 1] ( 0-1 1)
[-1 1 1] ( 0-1 1)
[-1 1 1] ( 1 0 1)
[-1-1 1] ( 1 0 1)
[ 1 1 1] (-1 0 1)

σ / MPa

[ 1-1 1] (-1 0 1)
[-1 1 1] ( 1 1 0)
[-1 1-1] ( 1 1 0)

2

10

[ 1 1 1] (-1 1 0)
[ 1 1-1] (-1 1 0)
[-1 1 1] ( 2 1 1)
[ 1 1 1] (-2 1 1)
[ 1 1-1] ( 2-1 1)
[ 1-1 1] ( 2 1-1)
[ 1-1 1] ( 1 2 1)
[ 1 1-1] (-1 2 1)
[ 1 1 1] ( 1-2 1)
[-1 1 1] ( 1 2-1)

101

[ 1 1-1] ( 1 1 2)

0

0.1

0.2

0.3

0.4

0.5

[ 1-1 1] (-1 1 2)
[-1 1 1] ( 1-1 2)
[ 1 1 1] ( 1 1-2)

ε

Figure 5.9: Influence of pre-deformed dislocation content on stress-strain response (lower left)
simulated by dislocation density based model. The top row shows several instances of inhomogeneous distribution of dislocation content on individual slip systems for a chosen log-normal
dislocation statistics. Each slip system is denoted by a color key given in the lower right legend.
tion characteristic (ρ mean , σρ ), several random instances of the distributions were generated. The
resulting effect on stress-strain response and lattice-reorientation is reported.
Fig. 5.9 shows the influence of inhomogeneous dislocation content on the stress-strain response. It was observed that the flow stress of the material is sensitive to the initial dislocation
content. Fig. 5.10 shows the influence of σρ on the tensile axis re-orientation for three different
orientations. For each orientation and distribution characteristics, 10 different instances were simulated. It was observed that for a smaller inhomogeneity, the lattice re-orientation trajectories were
limited to only a few possible directions in orientations space for each of the tested orientations.
68

However, with the increase in the variance of the distribution, the tensile axis re-orientation became
more chaotic and unpredictable.
The above behavior is expected since the mobile dislocation content controls the deformation
process. For a given loading conditions, if the slip system with highest RSS is low on mobile
dislocation content, then the flow-stress has to increase to activate new sources for deformation
to continue. Due to this reason, the initial flow stress observed in Fig. 5.9 increased for certain
distribution instances where the slip system with highest Schmid factor has low mobile dislocation
content. Secondly, the distribution of initial dislocation content on non-active systems influences
the forest hardening. Thus, depending on the amount of particular forest dislocation content present
during flow, the tensile axis orientation can change due to different easy paths offered during
deformation. Hence, the less inhomogeneous a material is, the fewer number of paths available for
the lattice to reorient during deformation.

69

h1 1 1i

h1 1 1i

h1 1 1i

cumulative probability

1.00

0.75

h1 1 2i

h1 1 2i

h1 1 2i

0.50

0.25
var = 0.5
0.00 10
10

1011

1012

ρ (m-2)

1013

1014

h0 0 1i

h0 1 1i

h0 0 1i

h0 1 1i

h1 1 1i

h0 0 1i

h0 1 1i

h1 1 1i

h1 1 1i

cumulative probability

1.00

0.75

h1 1 2i

h1 1 2i

h1 1 2i

0.50

0.25

0.00 10
10

var = 1.0

1011

1012

ρ (m-2)

1013

1014

h0 0 1i

h0 1 1i

h0 0 1i

h0 1 1i

h1 1 1i

h0 0 1i

h0 1 1i

h1 1 1i

h1 1 1i

cumulative probability

1.00

0.75

h1 1 2i

h1 1 2i

h1 1 2i

0.50

var = 2.0

0.25

0.00 10
10

1011

1012

ρ (m-2)

1013

1014

h0 0 1i

h0 1 1i

h0 0 1i

h0 1 1i

h0 0 1i

h0 1 1i

Figure 5.10: Influence of pre-deformed dislocation content on lattice reorientation for three different orientations simulated by dislocation density based model. The left column shows the probability distribution of pre-deformed dislocation content with ρ m =
1 × 1013 m−2 and different variances. The corresponding lattice reorientation trajectories for three orientations are shown in the right
three columns.
70

5.3.4

Influence of dislocation interaction coefficients on strain hardening

Fig. 5.11 shows the tensile axis re-orientation and axial stress-strain response of the best-fitted
simulation results for the eight orientations simulated with isotropic Fe interaction coefficients
[2]. The tensile axis re-orientation predicted by the model is similar to that observed by Baars
[12], Jaoul and Gonzalez [8], Keh [9], Duesbery et al. [74]. The tensile axis for most of the
orientations lying in the center of standard stereographic triangle changes trajectory after crossing
the [0 0 1]–[1 0 1] line. The maximum amount of angular distance traversed by the stress-axis in the
adjacent standard triangle is determined by how close the initial orientation is to the {0 0 1}–{1 0 1}
boundary.

150

X
S QR
W

P

V

100

σ / MPa

T

V
P
W
R
X
Q
T
S

50

0

0

0.1

0.2

0.3

0.4

0.5

ε
Figure 5.11: Results predicted by the dislocation density based model with isotropic α − Fe interaction coefficients. Left: ⟨0 1 0⟩ inverse pole figure for the 8 orientations showing the tensile axis
re-orientation as a function of strain (εmax = 0.30. The initial tensile axis orientation is denoted
by the blue dot. The final orientation after 30 % strain is indicated by the red dot. The black dots
in-between represents the trajectory of reorientation. Right: Corresponding axial stress-strain (σ ε ) predicted by the crystal plasticity model.

71

The rate of hardening observed in stage II deformation of the simulated single crystal uniaxial
tension differs significantly from the experimental behavior. A negligible rate of hardening in
stage II is observed in axial stress-strain behavior for most of the orientations. The reason for
such low hardening is due to the fact that the used interaction coefficients were calculated for α Fe (BCC) assuming an isotropic elasticity. Until recently, it was believed that these interaction
coefficients are directionally symmetric and depend only on the lattice structure. However, Madec
and Kubin [3] predicted that the anisotropic elasticity has an effect on the dislocation interaction
coefficients and is, therefore, material dependent. Moreover, it was also demonstrated that some
of the dislocation junction reactions are asymmetric in nature.

150

Experiment

100
σ / MPa

Simulation
Nb interaction
strength

50

0

Simulation
Fe interaction
strength

0

0.1

0.2

0.3

0.4

0.5

ε
Figure 5.12: Influence of elastic anisotropy on strain hardening behavior for a particular orientation. The interaction coefficients for Fe assuming isotropic elasticity was calculated by Queyreau
et al. [2]. The dislocation interaction coefficients assuming anisotropic elastic for Nb was taken
from the works of Madec and Kubin [3].

Fig. 5.12 shows the influence of dislocation interaction coefficients (calculated under different
72

assumptions) on the rate of hardening in single crystal Nb. The black curve represents the experimentally observed stress-strain behavior. The red curve represents the stress-strain behavior
simulated assuming isotropically elastic α -Fe interaction coefficients, while the green curve represents the stress-strain behavior simulated using anisotropically elastic Nb interaction coefficients.
It was observed that for α -Fe interaction coefficients the rate of hardening is negligible, while for
the Nb interaction coefficients the rate of hardening increased rapidly at onset of stage II.

5.3.5

Inverse optimization of dislocation mean free path coefficient

The dislocation multiplication term in the storage and recovery framework is inversely proportional
to the MFP coefficient, which is orientation dependent. The MFP coefficient takes into account the
probability of forming stable junctions and the number of active slip systems for the given lattice
orientation. Devincre et al. [16] calculated the MFP coefficients in Al (FCC) for several symmetric
orientations. In their work, it was shown that for orientations with more active slip systems the
average MFP is shorter, while for low symmetry orientations (single glide case) the MFP is larger,
which explains a long easy glide stage in stress-strain response.
So far these coefficients have not been reported for any BCC material in literature. Hence, in
this work, the MFP coefficient is estimated using parameter optimization. It is to be noted that
the optimized set of MFP coefficients are in no way absolute values for the given material or its
orientations, but a result of parameter fitting with the primary objective being their qualitative trend
over orientation space.
Fig. 5.13 shows the uniaxial stress–strain response for eight single crystal orientations after
{1 1 0}

optimization of Khkl , yhkl , and the initial dislocation content ρ0

{1 1 2}

and ρ0

. The average

relative error in the difference between measured and simulated stress-strain curves for eight ori73

100

100

σ / MPa

σ / MPa

150

50

0

Q
50

0

0.1

0.2

0.3

0.4

0

0.5

150

100
R

50

0

0.1

0.2

ε

0.3

0.4

0

0.5

150

100

100

S
50

0

0.1

0.2

ε

150

σ / MPa

P

100

150

σ / MPa

150

0.3

0.4

0

0.5

0

0.1

0.2

ε

0.3

0.4

0.5

0.4

0.5

ε

150

150

50

0

0.1

0.2

0.3

ε

0.4

0.5

0

σ / MPa

100

σ / MPa

T
50

0

100

σ / MPa

σ / MPa

V
W

50

0

0.1

0.2

0.3

0.4

0.5

ε

0

X
50

0

0.1

0.2

0.3

ε

0.4

0.5

0

0

0.1

0.2

0.3

ε

Figure 5.13: Comparison of experimentally observed (solid) and predicted (dashed) uniaxial
stress–strain response for different orientations. The coefficients for mean free path Khkl and dynamic recovery yhkl in the model by Devincre et al. [16] were estimated using inverse optimization
of the measured stress–strain curves.
entations is less than 6 %. The experimentally observed variability associated with the different
crystallographic tension directions is reproduced. Except for a few orientations, the stage II hardening, which starts to become apparent at strain levels of around 0.2, is generally well captured. It
is observed that as the orientation moves from single glide conditions towards multiglide condition,
the coefficient decreases indicating shorter average MFP i.e. higher rate of dislocation reaction and
strain hardening.
Table 5.3 presents the underlying values of the adjustable parameters for all eight orientations.
The annihilation distance yhkl , which is governing the dynamic dislocation recovery, is not strongly
dependent on the specific orientation, but assumes values of around 5 nm each. This invariance is
understandable as the dislocation density accumulated during the straining is still noticeably below
that density for which Eq. (4.14) would reach a dynamic equilibrium, hence, would show a strong
influence of the dynamic recovery aspect.

74

Table 5.3: Variability in constitutive parameter values resulting from minimizing the deviation in
single crystal stress–strain response for eight distinct tensile directions.
Orientation

Khkl

{1 1 0}
{1 1 2}
yhkl
ρ0
ρ0
(nm) (1012 m−2 ) (1012 m−2 )

P
Q
R
S
T
V
W
X

20.00
7.25
12.18
17.45
15.52
5.56
7.16
9.58

4.29
7.25
7.76
8.34
5.52
7.16
4.56
8.58

0.62
0.12
0.33
0.09
0.83
0.08
0.12
0.04

0.12
0.05
0.11
0.02
0.15
0.07
0.07
0.11

The dislocation storage parameter Khkl exhibits a notably stronger dependence on the crystallographic tensile direction. As such, the simplified storage and recovery model employed here is
of limited use for general application in which complex (i.e. not only unidirectional) loadings can
occur. Therefore, as a future step, the modeling of dislocation storage needs to take into account
the instantaneous dislocation densities to more directly capture the kinematics and dynamics of
dislocation junction formation as was already indicated by Madec and Kubin [113].
{1 1 0}

The values of initial1 dislocation density (ρ0

{1 1 2}

and ρ0

) that are necessary to match the

yield stress level in each of the four single crystal experiments also exhibit an appreciable variability. Since it is difficult to accurately determine the true dislocation content in the material, a
numerical study was performed to gauge the influence of this uncertainty on the predicted single
crystal stress–strain response in unidirectional tension. Figure 5.14 presents for the exemplary orientation ‘V’ that even for a known (fixed) total initial dislocation content, if distributed unevenly
across the individual slip systems, the resulting strain hardening can markedly differ. In the case
shown, three of the realizations did not even exhibit the experimentally observed two-stage hard1 homogeneously distributed across

each of the twelve slip systems per slip family

75

150
Experiment
Simulation

σ / MPa

100

50

0

0

0.1

0.2

0.3

0.4

0.5

ε
Figure 5.14: Variability in predicted single–crystal response under unidirectional tension due to
different distributions among the slip systems for fixed total initial dislocation content.
ening, indicating that the precise matching of single crystal deformation is probably elusive even
with rather elaborate constitutive models.

76

CHAPTER 6
CONCLUSION

In the present study, volumetric dilatation capability is added to an existing isotropic material
model to capture the behavior of void-like regions in F OURIER-based spectral solvers. The ability to mimic empty space is essential to address non-compact geometries within the periodicity
constraints imposed by such solvers. The model has been applied to simulate the plastic growth
of a circular void in a plate under biaxial tension and the known grain-scale deformation behavior of a tensile sample containing a limited number of grains. Both cases show that the proposed
formulation provides a more accurate description of void-like behavior compared to assuming a
soft-elastic behavior for void volumes.
The finite strain dilatational framework is then used for crystal plasticity simulations of dogbone shaped Niobium single crystals under tensile loading with an emphasis on multi-stage hardening, orientation dependence, and non-Schmid behavior. A constitutive model with dislocation
storage and recovery rates based on Discrete Dislocation Dynamics is used to model strain hardening in stage II. Adjustable parameters in this model are identified based on an inverse strategy that
uses a Nelder–Mead simplex approach to minimize the deviation between measured and simulated
uniaxial single crystal tension experiments. The influence of dislocation interaction parameters and
the variability on stage II hardening associated with the initial dislocation content is numerically
studied and compared with tensile experiments.
The experimentally observed variability associated with the different crystallographic tension
directions is reproduced effectively. Except for few a orientations, the stage II hardening, which
starts to become apparent at strain levels of around 0.2, is generally well captured. The transition
77

from stage I to stage II is abrupt and non-smooth. This might be due to an absence of cross-slip in
the constitutive description.
The annihilation distance yhkl , which is governing the dynamic dislocation recovery, is not
strongly dependent on the specific orientation, but assumes values of around 5 nm. The dislocation storage parameter Khkl exhibits a notably stronger dependence on the crystallographic tensile
direction.
The values of initial dislocation density (ρ {1 1 0} and ρ {1 1 2} ) that are necessary to match the
yield stress level in the eight single crystal experiments also exhibits an appreciable variability.
Since it is difficult to accurately determine the true dislocation content in the material, a numerical
study was performed to gauge the influence of this uncertainty on the predicted single crystal
stress–strain response in unidirectional tension. In the case shown, some of the realizations did not
even exhibit the experimentally observed two-stage hardening, indicating that the precise matching
of single crystal deformation is likely elusive even with rather elaborate constitutive models.

6.1 Recommendations for future work
• The transition from stage I to stage II is abrupt for the simulated stress-strain behavior. For
a smooth transition in-between hardening stages, cross-slipping should be incorporated in
the crystal plasticity constitutive law. Cross-slipping has been reported by several authors
[121, 41, 122] to be a prominent phenomenon in BCC materials at higher temperatures.
For a dislocation to cross-slip in FCC materials, the Shockley partials have to combine to
form complete dislocation before gliding on a cross-slip plane, which is why cross-slip is
less probable in the materials with low Stacking Fault Energy (SFE). Following the same
ideology cross-slip in BCC materials should be easy since most of the BCC materials do

78

not form stable SF’s at room temperature and above. Moreover, the availability of a vast
number of possible slip planes with common ⟨1 1 1⟩ directions makes cross-slip very likely
all throughout deformation process.
• The current dislocation storage and recovery framework used to simulate the strain hardening is a simplified model and is of limited use for general application in which complex
(i.e. not only unidirectional) loadings can occur. Moreover, the orientation dependence of the
mean free path parameters necessitates calibration of individual orientations. Therefore, as a
future step, the modeling of dislocation storage needs to take into account the instantaneous
dislocation densities to more directly capture the kinematics and dynamics of dislocation
junction formation as already indicated by Madec and Kubin [113].
• To correct the relative deviation in the simulated stage II hardening observed for orientations having equal probability of activating ⟨1 1 1⟩{1 1 0} and ⟨1 1 1⟩{1 1 2} slip systems,
non-Schmid corrections suggested by Dezerald et al. [123] needs to implemented.

79

CHAPTER 7
ADDITIONAL PROJECTS

7.1 Influence of geometric reconstruction accuracy on grain-averaged aggregate mechanics
7.1.1

Introduction

The elastic and, in particular, plastic mechanical behavior of single crystals is generally anisotropic.
Therefore, the macroscopic mechanical response of polycrystalline materials depends primarily
on the crystallite orientation distribution (“texture”). When considering their microscopic behavior at the scale of individual grains, the distributions of size and shape of constituent grains become equally important since the inherent anisotropy of single crystal deformation coupled with
inter-granular interactions results in strong heterogeneity of local deformation even under uniform
macroscopic loading (see, for instance, [14, 124]). Since many critical engineering aspects are
governed not by average but by extreme behavior, for instance void/crack formation due to strain
concentration, knowledge about and predictability of deformation at sub-granular scales is needed.
Deformation models that explicitly account for crystal plasticity have matured during the last
decades and became generally quite successful in predicting macroscopic responses, i.e., at the
scale of typical engineering components [76]. Application of such computational models to predict grain-scale micromechanics, however, currently results in overall agreement but notable deviations in detail compared to measurements [92]. The verification of constitutive crystal plasticity
models requires experimental characterization of the initial and final material state, which primarily includes crystal orientation, plastic strain, lattice strain (stress), and defect structure. For the

80

simulation of grain-scale mechanics, the most critical information among these is the initial grain
structure. This structure needs to be characterized non-destructively and in three dimensions, since
the (easily observable) response on the surface depends notably on the actual subsurface grain
structure [125, 126, 127]. Therefore, experimental techniques capable of measuring the threedimensional arrangement of grains are essential.
With the development of third-generation synchrotron facilities such as the Advanced Photon Source (APS) at Argonne National Laboratory, X-rays of high energy (50 keV to 100 keV)
and sub-micron sized beams can be generated, which have penetration depths ranging from a
few millimeters to centimeters and are capable of determining properties in three dimensions, resolved at the intragranular scale. Over the years, the use of monochromatic high-energy X-rays
combined with rotation of the investigated samples has lead to the development of various bulk
characterization techniques generally known as Three-Dimensional X-ray Diffraction (3DXRD)
[128, 129, 130] or High Energy X-ray Diffraction Microscopy (HEDM) [131]. HEDM is a nondestructive way of characterizing microstructure in three dimensions [132, 133]. HEDM can be
utilized with near-field and/or far-field detector positions or a combination of both. In the nearfield configuration, high spatial resolution is achieved at the expense of angular resolution, by
employing a small sample-to-detector distance on the order of the diffracting volume (about 1 mm
to 10 mm). Orientation maps with high spatial resolution that reveal the grain morphology can
be determined by post processing of the data using computationally intensive algorithms [17]. In
the far-field configuration, the sample-to-detector distance is large (about 1 m) compared to the
diffracting volume defined by the beam dimension. Thus, the locations of diffracted intensity on
the detector are governed by (i) the incident beam energy and (ii) the lattice orientation and lattice
distortion of diffracting grains. Average values of crystal orientation, lattice strain, volume and
81

position (centroid) can be conveniently extracted for each identified grain from the acquired sequence of detector images. Given the large effort and data intensity to measure spatially resolved
grain maps using near-field HEDM compared to centroid locations using far-field HEDM, it needs
to be clarified (i) how accurately a spatial tessellation based on the latter information can approximate the true grain structure that the former more directly represents, and (ii) how sensitive the
grain-averaged mechanical behavior is to the geometric uncertainties around grain perimeters that
result from shape reconstruction based on limited information.
The present paper first compares the geometric reconstruction accuracy when using weighted
(Laguerre) Voronoi tessellation based only on grain centroids information to that resulting from
combining centroids and grain volume information. Next, the plastic response of reference grain
aggregates is simulated based on phenomenological crystal plasticity and compared to that of the
two reconstruction alternatives to elucidate the sensitivity of grain-averaged lattice strains on the
accuracy of grain morphology reconstruction.

7.1.2

Spatial reconstruction method

Standard Voronoi as well as weighted (a.k.a. Laguerre) Voronoi tessellation generates volumefilling cells based on a population of ‘seed’ points. The volume enclosed by the cell of site si =
(pi , wi ) , i = 1, . . . , Ngrains , is
;
<
!
"
Ci = x ∈ R3 | d 2 (x, pi ) − wi ≤ d 2 x, p j − w j , j ̸= i ,

(7.1)

with d (x, pi ) being the distance between position x and the (closest1 ) seed position pi . In case
of standard Voronoi tessellation, all weights wi = const = 0, while for Laguerre tessellation the
1 potentially considering

spatially periodic repetitions of all seed points

82

(a)
0.5

0.75

0.50
2.0

cumulative volume fraction

1.00

0.5

0.25

1.0

1.5

0.00
0.00

0.25

0.50

2.0

0.75

1.00

cumulative number fraction

(b)

10-1

1

101

V0 / 〈V0 〉

(c)

Figure 7.1: Polycrystalline aggregates (a) measured by various 3D techniques [17, 18, 19, 20] and
associated grain statistics (blue curves in (b) and (c)). Sequence of four gray curves represents lognormal grain volume distributions of different variance (Eq. (7.9)); the black curve corresponds to
closest fit at variance σdistr = 1.2. Dotted line represents average grain volume.
individual weights are variable. To use Laguerre tessellation as a tool to reconstruct an aggregate of grains based on limited (spatially-averaged) data, a connection between such data and the
tessellation sites (position and weight) needs to be established, as later outlined in Section 7.1.4.

7.1.3

Volume statistics of real and artificial grain structures

To compare the quality of grain structure reconstruction in presence/absence of grain volume information (to inform the Laguerre weights), synthetic grain structures that are statistically similar
to real ones are used.

83

cube
octahedra

super
ellipsoids

(a)
cumulative volume fraction

1.00

0.75

0.50

0.25

0.00
0.00

0.25

0.50

0.75

1.00

cumulative number fraction

(b)

10-1

1

101

V0 / 〈V0 〉

(c)

Figure 7.2: Synthetic polycrystals (a) generated by DREAM.3D [21, 22] using cube octahedra
(top row) or super ellipsoids (bottom row) as grain shapes. Orange and brown curves in (b) and (c)
show the associated grain statistics for both selections of grain shapes. Black line (optimal fit to
experimental data) is same as in Fig. 7.1.
Four exemplary lattice orientation datasets, which were measured by various 3D techniques
[17, 18, 19, 20], are shown in Fig. 7.1a and form the statistical basis of the present investigation.
Grain size populations were derived from each dataset by counting the number of voxels assigned
to each grain.2 Figure 7.1(b) presents the relation (blue curves) between cumulative number fraction and cumulative volume fraction for each grain size population. Represented in this normalized
2 Volumes

of surface grains resulted from combining every member in the lower half with a
random one from the upper half of their volume-sorted population.

84

way, all four grain populations are remarkably similar and can be approximated by a log-normal
distribution of grain volume (black fit line, see Section 7.2 for details). The same datasets and lognormal grain volume distributions are presented in Fig. 7.1(c) as a volume-weighted distribution
of relative grain volume (defined in 7.3), demonstrating again the close approximation of all data
by a log-normal distribution with variance σdistr = 1.2.3 Assuming that this particular distribution shape faithfully represents polycrystalline aggregates, it can serve as a reference for synthetic
datasets to identify systematic deviations relative to it.
Synthetic grain structures were generated by a series of filters using DREAM.3D [134] (see
7.4 for list of filters and [21, 22] for details of the generation process). A log-normal distribution
of grain sizes with mean µsize = 1.0 and variance σsize = 0.4 was created in the filter ‘StatsGenerator’ and matched well with experimentally observed grain size distributions (Fig. 7.1(b) and (c)).
Next, based on the above grain size distribution, two versions of four polycrystals (Fig. 7.2(a))
having grain shapes of either cube octahedra or super ellipsoids were generated with the filter ‘Establish Shape Types’ to investigate the potential influence of shape anisotropy. Invoking the same
statistical representation as in Fig. 7.1 demonstrates the close similarity of synthetic and measured
structures in terms of their grain volumes (see Fig. 7.2(b) and (c)). In terms of grain shapes, the
cube octahedral grain structures are close to equiaxed while the super ellipsoidal ones are generally
more elongated as can be inferred from the ratios of maximum to minimum principal moments of
inertia, i.e. the “eccentricity,” that are shown cumulatively in Fig. 7.3 (solid lines) from all synthetic
grains.
3 Translations

between the frames (b) and (c) in Figs. 7.1 and 7.2 make use of the fact that the
relative grain volume in (b) is simply given by the slope dFV /dFN since the average grain volume
in that frame is ⟨V ⟩ = Vtotal /Ngrains = 1/1 = 1, falling along the dotted line in (b).

85

The volume Vref and centroid4 of every constituting grain (Ngrains ≈ 500) was extracted from
each of the eight synthetic polycrystals to serve as a basis for comparing the reconstruction method
variations that are presented next.

supe

r ell

0.50

ipso
id

cube octa
hedron

0.75

0.25

0.00

equiaxed

cumulative probability

1.00

5

1

shape factor (κ)
Figure 7.3: Ratio of largest to smallest principal moment of inertia (‘shape factor’) for all cube
octahedral (orange) and super ellipsoidal (brown) grains. Solid lines correspond to synthetic (reference) structures shown in Fig. 7.2(a). The fainter curves result after reconstructing the reference
structures using either grain volumes (dashed) or centroid distances (dotted) as a basis for Laguerre
tessellation weigths. In general, cube octahedral grain structures are closer to being equiaxed than
those based on super ellipsoids.

7.1.4

Comparison of 3D reconstruction of grain aggregates in presence/absence of grain
volume information

Given that the identification of grain centroids is typically more common than determination of the
associated grain volumes in far-field HEDM, we investigate how accurately the known (synthetic)
polycrystals described in Section 7.1.3 can be reconstructed based solely on centroid information
compared to having additional explicit knowledge of grain volumes. For both cases, termed ‘absent’ and ‘present’, the centroids are a natural choice as seed points in a Laguerre tessellation.
4 spatial average of

the positions of all voxels that comprise a grain
86

The ‘present’ case has already been studied by Lyckegaard et al. [23] with the result that a fairly
accurate reconstruction can be achieved when the Laguerre weights are chosen as w = (requiv )2 =
(3Vref /4π )2/3 , i.e. as the squared equivalent sphere radius of grain volume Vref . Furthermore,
they observed a slight improvement in accuracy when correcting each seed position by the offset
between the initial seed position of each grain and the centroid resulting after a first tessellation.
In the ‘absent’ case, i.e. without explicit knowledge of the correct grain volumes, selection of
Laguerre weights is not straightforward. One possible choice, which is used in the present study, is
w = (dnn )2 with dnn being the distance between a grain centroid and its nearest neighboring one.
The reconstruction accuracy resulting from both above conditions, i.e. with and without explicit knowledge of the grain volume, is presented in Fig. 7.4 as joint probability density of correct
grain volume Vref and relative deviation Vtess /Vref (outer rows) or κtess /κref (inner rows) between
tessellated and correct volume or shape factor, respectively. In case of Laguerre weights being
derived from known grain volume information (‘present’, upper half in Fig. 7.4), most of the reconstructed grains are within about 20 % of the correct volume and show a systematic decrease of
reconstruction accuracy with decreasing grain size. The slight reduction in grain volume reconstruction accuracy between more equiaxed (cube octahedral) and less equiaxed (super ellipsoidal)
grain shape that is observed in the upper row of Fig. 7.4 is expected since the Laguerre tessellation is isotropic by design, hence, will become increasingly inaccurate with increasing structural
anisotropy. Similarly, the reconstruction accuracy in terms of grain shape factor (second row in
Fig. 7.4) is higher for the synthetic structures utilizing the more equiaxed cube octahedral grains
than for the super ellipsoidal grain shapes. The two-pass correction of seed points suggested by
Lyckegaard et al. [23] removes the systematic overprediction of tessellated volumes for the largest
grains (compare red to green probability maps). In case of Laguerre weights being derived from
87

Cube Octahedron

Laguerre
weighting

seed at
centroid

Super Ellipsoid

adjusted
seed position

seed at
centroid

adjusted
seed position

(requiv )2

Vtess / Vref

101

1

10-1

1

101

Vtess / Vref

‘absent’

κtess / κref

(dnn )2

1

κtess / κref

‘present’

1

10-1

10-1

1

101

Vref / 〈Vref 〉

10-1

1

101

Vref / 〈Vref 〉

10-1

1

101

Vref / 〈Vref 〉

10-1

1

101

Vref / 〈Vref 〉

Figure 7.4: Joint probability of correct grain volume (normalized by average grain size) and relative
deviation in the reconstructed grain volume Vtess /Vref or grain shape factor κtess /κref . Grain reconstruction
uses Laguerre tessellation with grain centroids as seed positions. Seed point weights are derived from known
grain volumes (upper half) or estimated from the distance to the nearest grain centroid (lower half). Left
and right panels correspond to cube hexahedral and super ellipsoidal grain shapes, respectively, that were
used in generating the underlying synthetic data sets (Fig. 7.2(a)). Red maps correspond to direct use of
grain centroids as seed positions, while green maps illustrate changes due to a seed position adjustment [23]
that considers the new centroid location of each resulting tessellated grain. The comparisons of shape factor
deviations (yellow maps) utilized only those latter reconstructed structures.

88

nearest neighbor distances of centroids (‘absent’, lower half in Fig. 7.4), the overall agreement
between correct and tessellated volume is substantially weaker—as will be discussed in connection with Fig. 7.6. Furthermore, the systematic deterioration of accuracy in grain volume as well
as grain shape reconstruction with decreasing grain size is much more pronounced in the ‘absent’
case compared to the ‘present’ case.
Given the above established uncertainties of grain volumes reconstructed by Laguerre tessellation with and without knowledge of correct grain volumes, it is informative to determine the
associated uncertainty of the grain-averaged lattice strains as margin of error when comparing
measurements to simulations. To answer this question, Crystal Plasticity Fast Fourier Transform
(CPFFT) simulations were performed, the results of which are reported in next section.

7.1.5

Grain averaged lattice strains of reconstructed structures

A finite strain framework is adopted in which the deformation gradient F = Fe Fp at each material point is multiplicatively decomposed into elastic Fe and plastic Fp components, thus introducing an intermediate (or ‘lattice’) configuration. The second P IOLA–K IRCHHOFF stress
!
"
!
"
S = C : Fe T Fe /2 = f Ḟ, η reflects the elastic lattice distortion (C being the fourth-order elas-

tic stiffness tensor) and drives the plastic velocity gradient Lp (S, η ) = Ḟp Fp −1 as well as the
evolution of internal state variables η (see [76] for details).
The constitutive model used in the present study is similar to the phenomenological description
of Peirce et al. [1] for cubic lattice symmetry. The internal state is captured through slip resistances

τ0 per slip system ξ , which evolve with shear γ according to a relationship initially proposed by

89

Laguerre
weighting

stress component ratios
normal

shear

deviatoric normal

(requiv

)2

σtess / σref

‘present’

σtess / σref

σtess / σref

101
1
10-1
101
1
10-1
101
1
10-1

1

101

102

103

1

σref / MPa

101

102

103

1

σref / MPa

101

102

103

σref / MPa

(dnn

)2

σtess / σref

‘absent’

σtess / σref

σtess / σref

101
1
10-1
101
1
10-1
101
1
10-1

1

101

102

103

1

σref / MPa

101

102

σref / MPa

103

1

101

102

103

σref / MPa

Figure 7.5: Deviations of normal, shear, and deviatoric normal components of grain-averaged
Cauchy stress σtess in reconstructed grain structures relative to corresponding stress components
σref resulting in the (cube octahedron) reference grain structures. Tessellations in upper panel
employ known grain volume information (requiv ) for Laguerre weights, while in lower panel the
nearest centroid distances (dnn ) are used instead. The seed positions were adjusted [23] in both
tessellation alternatives.
Hutchinson [98] and later improved by Kalidindi et al. [99]
⎛
⎞
β a( (
τ
( (
ġ˙ξ = qξ β h0 ⎝1 − 0 ⎠ (γ̇ β ( .
β
g∞

(7.2)

The slip system interactions are controlled by the hardening matrix qξ β . The shear rate on each

90

Table 7.1: Material parameters; elastic constants Cab , reference shear rate γ̇0 , stress exponent n,
initial and saturation slip resistance τ0 and g∞ , hardening parameters h0 , qξ β , and a.
Parameter

Value Unit

C11
C12
C44
γ̇0
n
τ0
g∞
a
h0
Coplanar qξ β
Non-coplanar qξ β

106.75
60.41
28.34
10−3
20
31
63
2.25
75
1
1.4

GPa
GPa
GPa
s−1
MPa
MPa
MPa

system is given by
( (n
( τξ (
( (
γ̇ ξ = γ̇0 ( ( sign τ ξ ,
( gξ (

(7.3)

∞

where γ̇0 is a reference shear rate, n the stress exponent, and τ ξ = S :
shear stress. The plastic velocity gradient
Lp = γ̇ ξ mξ ⊗ nξ

4
3
ξ
ξ
m ⊗ n the resolved

(7.4)

is additively composed from all shear rates, where mξ and nξ are unit vectors along the slip
direction and slip plane normal, respectively. Table 7.1 summarizes the adjustable parameters used
to mimic a material similar to aluminum.
For the periodic structures investigated in this study, spectral methods are an efficient alternative to the finite element method (FEM) for solving the mechanical boundary value problem of
static equilibrium [135, 136, 137, 80, 82]. Here, we employ a finite strain spectral solver developed by Eisenlohr et al. [82] and Shanthraj et al. [88] that includes a flexible material point model

91

(DAMASK) to incorporate the constitutive law described above. All (periodic) grain aggregates
are discretized on a regular 100 × 100 × 100 grid.
Due to the periodic nature of the geometry, boundary conditions for simulating the chosen
reference case of uniaxial tension along the y-direction have to be specified as volume-averaged
deformation gradient rate and complementary (first P IOLA–K IRCHHOFF) stress
⎡
⎤
⎡
⎤
⎢∗ 0 0⎥
⎢
⎥
˙
⎢
⎥
⟨F⟩
⎢
=
0 1 0⎥
⎥
10−3 s−1 ⎢
⎢
⎥
⎣
⎦
0 0 ∗

and

⎢0 ∗ ∗⎥
⎢
⎥
⎢
⎥
⟨P⟩ ⎢
= ⎢∗ ∗ ∗⎥
⎥,
Pa
⎢
⎥
⎣
⎦
∗ ∗ 0

(7.5)

with ‘*’ indicating mixed boundary conditions, where either stress (or deformation) is prescribed
so that deformation (or stress) needs to be iteratively adjusted. In the simulations, above conditions
are maintained during 100 increments of 1 s each (such that the final ⟨F⟩yy = 1.1).
To elucidate the influence of reconstruction quality on the average lattice strain predicted from
crystal plasticity simulations of grain aggregates, the four synthetic grain aggregate geometries
that are based on cube octahedral grain shapes were selected as reference. The above specified
deformation was simulated for those references as well as for the associated two (seed corrected)
reconstruction alternatives, where one tessellation utilizes grain volume information, the other
estimating it from nearest centroid distances (termed ‘present’ and ‘absent’ case).
Figure 7.5 compares for both tessellation cases (upper and lower panel) the relative deviation of
average Cauchy stress in each reconstructed grain from its reference result. The relative deviations
observed for the same stress component are very similar for both cases. For both reconstructed
structures, the (green) normal stress component along the loading direction agrees very closely
to that found in the respective reference structure. The magnitude of the other two normal stress
components as well as of the three shear stress components is in general substantially lower than
92

dnn / requiv

5

(d

nn

/r

eq
ui
v) 3

1

10-1

1

101

Vref / 〈Vref 〉
Figure 7.6: Joint probability of volume Vref of synthetic grains and the distance dnn from each grain
centroid to the nearest neighboring one. Quantities are normalized by the average grain volume
⟨Vref ⟩ and its equivalent radius requiv = (3⟨Vref ⟩/4π )1/3 , respectively. Cumulative probability is
scaled to unity within each (vertical) bin of normalized grain volume. Centroid distance progressively deviates from expected relation (dashed line) with decreasing volume of grain, resulting in
a weak correlation between the distance to the nearest grain centroid and the volume of a grain.
the normal stress observed along the loading direction. The deviations of these five components
progressively increase with decreasing (reference) stress magnitudes, i.e., they fan out toward the
left, occasionally even resulting in a reversal of sign (illustrated by open diamonds in Fig. 7.5).
The rightmost column in Fig. 7.5 shows the relative deviation of deviatoric normal components,
i.e., excluding the mean stress σsph = (σ11 + σ22 + σ33 ) /3. In contrast to the normal stress components shown in the leftmost column, the relative deviation of the deviatoric normal components
is markedly reduced and essentially independent of the magnitude of reference deviatoric stress.

93

7.1.6

Discussion

As shown in Fig. 7.4, the uncertainty in the volume of a reconstructed grain is much larger when
the tessellation is estimating the Laguerre weights from nearest centroid distances compared to
when directly using the known reference volumes. This increased uncertainty can be connected to
the weak correlation found between the volume of a grain and its distance dnn to the nearest neighboring grain centroid, which is illustrated in Fig. 7.6. For a given dnn , the potential grain volume
associated with that distance falls into a range of about a factor of three for the relatively large
grains (right end of Fig. 7.6), while for the relatively small grains (left end) this range increases to
about an order of magnitude. This increasing uncertainty is caused by the progressive deviation
with decreasing grain size of the most likely nearest neighbor distances toward values above the
natural correlation Vref 1/3 ∝ dnn (dashed line in Fig. 7.6). Since small grains are likely to be situated in a neighborhood of larger grains, the distance to the centroid of their nearest neighboring
grain is frequently unrelated to and far larger than its own equivalent radius. Such a disconnect is
less probable for a relatively large grain since, for many conceivable neighborhoods, the closest
centroid distance is primarily determined by the size of the large grain itself.
Although a sizable difference emerged in the quality of reconstruction between grain structures
tessellated with the Laguerre weights either derived from grain volumes or estimated from nearest
centroid distances (see Fig. 7.4), under the uniaxial loading applied in the present study, both tessellated grain structures show very similar deviations from the reference structure in terms of observed
lattice strains (i.e. stress). Figure 7.7 collects deviations across all six Cauchy stress components
from both tessellations in a joint probability plot. The observed increase in relative deviations (uncertainty) with decreasing (reference) stress magnitude is compatible with an absolute uncertainty

94

σtess / σref

101

1

Pa

M

10
M
Pa

5

10-1

1

101

102

103

σref / MPa
Figure 7.7: Joint probability of stress magnitudes σref in reference structure and relative deviations
σtess /σref combined from both tessellated structures. Cumulative probability is scaled to unity
within each bin of stress magnitudes. Solid lines exemplify relative deviations in stress due to
absolute differences of 5 MPa and 10 MPa. With increase in magnitude of stress, the relative
deviations in stress due to difference in grain morphologies diminishes.
in the range of 5 MPa to 10 MPa as illustrated by the solid lines in Fig. 7.7. The apparently much
closer agreement between deviatoric normal stress components (Fig. 7.5 right) compared to the
normal stress components containing hydrostatic stress (Fig. 7.5 left) can be rationalized with a
similar argument. The deviatoric stress projection results in an overall shift to absolutely larger
values for the two transverse stress components (xx and zz, red and blue in Fig. 7.5). Therefore,
the large relative deviations observed at small stress magnitudes tend to shrink when considered at
larger stress magnitudes.

95

7.2 Cumulative volume fraction of log-normal grain size distributions
Assuming individual grain volumes in a polycrystal to follow a log-normal distribution, their
cumulative number fraction FN as function of grain volume Vgrain is given by
#
B
C$
lnVgrain /µ
1
√
FN =
1 + erf
2
2σ

(7.6)

and the average grain volume is
⟨V ⟩ = µ exp

%

σ2
2

'

.

(7.7)

Therefore, the relative grain volume corresponding to a particular cumulative number fraction
follows as
%
'
2
√
Vgrain
σ
= exp
2 σ erf−1 [2FN − 1] −
.
⟨V ⟩
2

(7.8)

The relation between cumulative volume fraction FV and cumulative number fraction is established
by integration of Eq. (7.8), yielding
B
C
1
σ
−1
FV = erfc √ − erf [2FN − 1] .
2
2

(7.9)

Conveniently, Eq. (7.9) does not depend on the actual grain volume magnitudes but only on the
shape of their distribution, which is governed by the variance σdistr as illustrated by the sequence
of gray curves in Fig. 7.1.

7.3 Grain size normalization
The ‘relative grain volume’ of each grain in a population of Ngrains occupying an overall volume
Vtotal (typically counted in voxels) is defined as the grain volume, Vgrain , normalized by the average
grain volume:
Vrel =

Vgrain Vgrain
=
N
.
⟨V ⟩
Vtotal grains
96

(7.10)

The cumulative volume fraction as a function of relative grain volume can then be used to compare
different grain populations on the same scale.

7.4 Generation of synthetic grain structures
The pipeline listed in 7.2 is used to generate the synthetic structures shown in Fig. 7.2(a).
A log-normal grain size distribution with mean µsize = 1.0 and variance σsize = 0.4 created by
‘StatsGenerator’ was packed into the synthetic volume using grain shape types of either ‘cube
octahedron’ or ‘super ellipsoid’.
Table 7.2: DREAM.3D pipeline for generating synthetic grain structures.
Filter order
1
2
3
4
5

Filter name
StatsGenerator
Initialize Synthetic Volume
Establish Shape Types
Pack Primary Phases
Write Output

97

BIBLIOGRAPHY

98

BIBLIOGRAPHY

[1]

D. Peirce, R. J. Asaro, and A. Needleman. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metallurgica, 30(6):1087–1119, 1982.

[2]

S. Queyreau, G. Monnet, and B. Devincre. Slip systems interactions in alpha-iron determined by dislocation dynamics simulations. International Journal of Plasticity, 25(2):361–
377, 2009.

[3]

R. Madec and L. P. Kubin. Dislocation strengthening in FCC metals and in BCC metals at
high temperatures. Acta Materialia, 126:166–173, 2017.

[4]

P. Franciosi. Glide mechanisms in b.c.c. crystals: an investigation of the case of α -iron
through multislip and latent hardening tests. Acta Metallurgica, 31(9):1331–1342, 1983.

[5]

W. A. Spitzig and A. S. Keh. The effect of orientation and temperature on the plastic flow
properties of iron single crystals. Acta Metallurgica, 18(6):611–622, 1970.

[6]

A. S. Keh and Y. Nakada. Plasticity of iron single crystals. Canadian Journal of Physics,
45(2):1101–1120, 1967.

[7]

J. Kumagai, S. Takaki, S. Suzuki, and H. Kimura. Orientation dependence of quasi-threestage work hardening in high purity iron single crystals. Materials Science and Engineering:
A, 129(2):207–215, 1990.

[8]

P. B. Jaoul and D. Gonzalez. Deformation plastique de monocristaux de fer. Journal of the
Mechanics and Physics of Solids, 9(1):16–38, 1961.

[9]

A. S. Keh. Work hardening and deformation sub-structure in iron single crystals deformed
in tension at 298 K. Philosophical Magazine, 12(115):9–30, 1965.

[10] M. S. Duesbery and R. A. Foxall. A detailed study of the deformation of high purity niobium
single crystals. Philosophical Magazine, 20(166):719–751, 1969.
[11] A. Seeger and W. Wasserbach. Anomalous Slip – A Feature of High-Purity Body-Centred
Cubic Metals. physica status solidi (a), 189(1):27–50, 2002.
[12] D. C. Baars. Investigation of active slip systems in high purity single crystal niobium. PhD
thesis, 2013.
[13] F. Roters, P. Eisenlohr, T. R. Bieler, and D. Raabe. Crystal Plasticity Finite Element Methods

99

in Materials Science and Engineering. Wiley-VCH, Weinheim, 2010.
[14] Z. Zhao, M. Ramesh, D. Raabe, A. M. Cuitiño, and R. Radovitzky. Investigation of threedimensional aspects of grain-scale plastic surface deformation of an aluminum oligocrystal.
International Journal of Plasticity, 24:2278–2297, 2008.
[15] A. Chakraborty and P. Eisenlohr. Evaluation of an inverse methodology for estimating
constitutive parameters in face-centered cubic materials from single crystal indentations.
European Journal of Mechanics / A Solids, 66C:114–124, 2017.
[16] B. Devincre, T. Hoc, and L. Kubin. Dislocation Mean Free Paths and Strain Hardening of
Crystals. Science, 320(5884):1745–1748, 2008.
[17] W. Ludwig, P. Reischig, A. King, M. Herbig, E. M. Lauridsen, G. Johnson, T. J. Marrow, and
J. Y. Buffiere. Three-dimensional grain mapping by x-ray diffraction contrast tomography
and the use of Friedel pairs in diffraction data analysis. Review of Scientific Instruments, 80
(3):033905, 2009.
[18] W. Ludwig, A. King, P. Reischig, M. Herbig, E. M. Lauridsen, S. Schmidt, H. Proudhon,
S. Forest, P. Cloetens, S. R. Du Roscoat, J. Y. Buffière, T. J. Marrow, and H. F. Poulsen.
New opportunities for 3D materials science of polycrystalline materials at the micrometre
lengthscale by combined use of X-ray diffraction and X-ray imaging. Materials Science
and Engineering: A, 524(1-2):69–76, 2009.
[19] R. Pokharel, J. Lind, S. F. Li, P. Kenesei, R. A. Lebensohn, R. M. Suter, and A. D. Rollett.
In-situ observation of bulk 3D grain evolution during plastic deformation in polycrystalline
Cu. International Journal of Plasticity, 67:217–234, 2015.
[20] P. A. Shade, M. A. Groeber, J. C. Schuren, and M. D. Uchic. Experimental measurement of
surface strains and local lattice rotations combined with 3D microstructure reconstruction
from deformed polycrystalline ensembles at the micro-scale. Integrating Materials and
Manufacturing Innovation, 2(1):5, 2013.
[21] M. Groeber, S. Ghosh, M. D. Uchic, and D. M. Dimiduk. A framework for automated analysis and simulation of 3D polycrystalline microstructures. Part 1: Statistical characterization.
Acta Materialia, 56(6):1257–1273, 2008.
[22] M. Groeber, S. Ghosh, M. D. Uchic, and D. M. Dimiduk. A framework for automated
analysis and simulation of 3D polycrystalline microstructures. Part 2: Synthetic structure
generation. Acta Materialia, 56(6):1274–1287, 2008.
[23] A. Lyckegaard, E. M. Lauridsen, W. Ludwig, R. W. Fonda, and H. F. Poulsen. On the Use of
Laguerre Tessellations for Representations of 3D Grain Structures. Advanced Engineering

100

Materials, 13(3):165–170, 2011.
[24] H. and Padamsee, editor. Proceedings of the 1st TESLA Workshop. 1990.
[25] Proceedings of International Workshop on Next-Generation Linear Colliders. Technical
report, 1988.
[26] M. Hein. High-Temperature-Superconductor Thin Films at Microwave Frequencies.
Springer Berlin Heidelberg, Berlin, Heidelberg, 1999.
[27] Y. Wang, T. Plackowski, and A. Junod. Specific heat in the superconducting and normal state
(2-300 K, 0-16 T), and magnetic susceptibility of the 38 K superconductor MgB2: evidence
for a multicomponent gap. Physica C: Superconductivity, 355(3-4):179–193, 2001.
[28] X. X. Xi, A. V. Pogrebnyakov, S. Y. Xu, K. Chen, Y. Cui, E. C. Maertz, C. G. Zhuang, Q. Li,
D. R. Lamborn, J. M. Redwing, Z. K. Liu, A. Soukiassian, D. G. Schlom, X. J. Weng, E. C.
Dickey, Y. B. Chen, W. Tian, X. Q. Pan, S. A. Cybart, and R. C. Dynes. MgB2 thin films by
hybrid physical-chemical vapor deposition. Physica C: Superconductivity, 456(1-2):22–37,
2007.
[29] D. E. Oates, A. C. Anderson, C. C. Chin, J. S. Derov, G. Dresselhaus, and M. S. Dresselhaus.
Surface-impedance measurements of superconducting NbN films. Physical Review B, 43
(10):7655–7663, 1991.
[30] B. W. Maxfield and W. L. McLean. Superconducting Penetration Depth of Niobium. Physical Review, 139(5a):A1515–A1522, 1965.
[31] H. Stuart. Niobium: Proceedings of the International Symposium. Metallurgical Society of
AIME, Warrendale, Pa., 1984.
[32] P. Dhakal, G. Ciovati, P. Kneisel, and G. R. Myneni. Enhancement in Quality Factor of SRF
Niobium Cavities by Material Diffusion. IEEE Transactions on Applied Superconductivity,
25(3):1–4, 2015.
[33] T. R. Bieler, N. T. Wright, F. Pourboghrat, C. Compton, K. T. Hartwig, D. Baars, A. Zamiri,
S. Chandrasekaran, P. Darbandi, H. Jiang, E. Skoug, S. Balachandran, G. E. Ice, and W. Liu.
Physical and mechanical metallurgy of high purity Nb for accelerator cavities. Physical
Review Special Topics – Accelerators and Beams, 13(3):031002, 2010.
[34] FRIB SRF Cavities, 2012. 7th SRF Materials Workshop.
[35] Ingot niobium – the Choice Materials for Future SRF Cavity Applications?, 2012. 7th SRF
Materials Workshop.

101

[36] A. Mapar, D. Kang, T. R. Bieler, F. Pourboghrat, and C. C. Compton. Refining Crystal
Plasticity Finite Element Model for Deformation of Single Crystal Niobium. In 7th SRF
Materials Workshop, 2012.
[37] Overview of cavity requirements for high intensity and high energy, 2012. 7th SRF Materials
Workshop.
[38] H. Padamsee. RF Superconductivity. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim,
Germany, 2009.
[39] J. Knoblock. Field emission and thermal breakdown in superconducting niobium cavities
for accelerators. IEEE Transactions on Appiled Superconductivity, 9(2):1016–1022, 1999.
[40] E. Schmid and W. Boas. Kristallplastizität. Springer Berlin Heidelberg, Berlin, Heidelberg,
1935.
[41] G. I. Taylor and C. F. Elam. The Distortion of Iron Crystals. Proceedings of the Royal
Society A: Mathematical, Physical and Engineering Sciences, 112(761):337–361, 1926.
[42] H. Gough. The Behaviour of a Single Crystal of Formula-Iron Subjected to Alternating
Torsional Stresses. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 118(780):498–534, 1928.
[43] W. Fahrenhorst and E. Schmid. Über die plastische Dehnung von α -Eisenkristallen.
Zeitschrift für Physik, 78(5-6):383–394, 1932.
[44] J. J. Cox, G. T. Horne, and R. F. Mehl. Trans ASM, 49:118, 1956.
[45] H. Bilger. On the recovery of plastically deformed iron single crystals above room temperature. physica status solidi (b), 18(2):K135–K138, 1966.
[46] Jr F. L. Vogel and R. M. Brick. J. Metall AIME, page 700, 1953.
[47] A. Seeger. Why anomalous slip in body-centered cubic metals? Materials Science and
Engineering A, 319-321:254–260, 2001.
[48] J. R. Low Jr. and A. M. Turkalo. Slip band structure and dislocation multiplication in
silicon-iron crystals. Acta Metallurgica, 10(3):215–227, 1962.
[49] T. Taoka, S. Takeuchi, and E. Furubayashi. Slip Systems and Their Critical Shear Stress in
3% Silicon Iron. Journal of the Physical Society of Japan, 19(5):701–711, 1964.
[50] S. Takeuchi, E. Furubayashi, and T. Taoka. Orientation dependence of yield stress in 4.4%
102

silicon iron single crystals. Acta Metallurgica, 15(7):1179–1191, 1967.
[51] F. E. Neumann. Vorlesungen über die Theorie der Elastizität der festen Körper und des
Lichtäthers. B. G. Teubner-Verlag, Leipzig, 1885.
[52] J. W. Christian. Some Surprising Features of the Plastic Deformation of the Body-Centered
Cubic Metals and Alloys. Metallurgical Transactions A, 14:1237–1256, 1983.
[53] M. S. Duesbery. volume 8 of Dislocations in Solids, chapter The dislocation core and
plasticity, pages 67–173. Elsevier, Amsterdam, 1989.
[54] V. Vitek. Theory of the core structures of dislocations in bcc metals. Crystal Lattice Defects,
5:1, 1974.
[55] Q. Qin and J. L. Bassani. Non-schmid yield behavior in single crystals. Journal of the
Mechanics and Physics of Solids, 40(4):813–833, 1992.
[56] B. Escaig. Sur le glissement dévié des dislocations dans la structure cubique a faces centrées.
Journal de Physique, 29(2-3):225–239, 1968.
[57] B. Escaig. Dissociation And Mechanical Properties.Dislocation Splitting And The Plastic
Glide Process In Crystals. Le Journal de Physique Colloques, 35(C7):C7–151–C7–166,
1974.
[58] J. Bonneville and B. Escaig. Cross-slipping process and the stress-orientation dependence
in pure copper. Acta Metallurgica, 27(9):1477–1486, 1979.
[59] E-Ichi Furubayashi. Behavior of Dislocations in Fe-3% Si under Stress. Journal of the
Physical Society of Japan, 27(1):130–146, 1969.
[60] F. Louchet, L. P. Kubin, and D. Vesely. In situ deformation of b.c.c. crystals at low temperatures in a high-voltage electron microscope Dislocation mechanisms and strain-rate equation. Philosophical Magazine A, 39(4):433–454, 1979.
[61] T. Imura. In M. Meshii, editor, Proc. US-Japan Symp. Mechanical Properties of BCC Metals, pages 65–73, Warrendale, PA, 1982. TMS-AIME.
[62] A. Seeger. Progress and problems in the understanding of the dislocation relaxation processes in metals. Materials Science and Engineering: A, 370(1-2):50–66, 2004.
[63] R. A. Foxall and C. D. Statham. Dislocation arrangements in deformed single crystals
of niobium-molybdenum alloys and niobium-9 at.% rhenium. Acta Metallurgica, 18(11):
1147–1158, 1970.

103

[64] J. P. Hirth and J. Lothe. Theory of dislocations. John Wiley & Sons, New York, 1982.
[65] R. Peierls. The size of a dislocation. Proceedings of the Physical Society, 52(1):34–37,
2002.
[66] J-Y Kim, D. Jang, and J. R. Greer. Insight into the deformation behavior of niobium single
crystals under uniaxial compression and tension at the nanoscale. Scripta Materialia, 61(3):
300–303, 2009.
[67] R. Peierls. The size of a dislocation. Proceedings of the Physical Society, 52:34–37, 1940.
[68] A. Seeger. Peierls barriers, kinks, and flow stress: Recent progress. Zeitschrift für Metallkunde, 93:760–777, 2002.
[69] L. Hollang, M. Hommel, and A. Seeger. The Flow Stress of Ultra-High-Purity Molybdenum
Single Crystals. physica status solidi (a), 160(2):329–354, 1997.
[70] F. Ackermann, H. Mughrabi, and A. Seeger. Temperature- and strain-rate dependence of the
flow stress of ultrapure niobium single crystals in cyclic deformation. Acta Metallurgica,
31(9):1353–1366, 1983.
[71] T. E. Mitchell, R. A. Foxall, and P. B. Hirsch. Work-hardening in niobium single crystals.
Philosophical Magazine, 8(95):1895–1920, 1963.
[72] E. Votava. Eine neue Methode zur Herstellung verformungsfreier Einkristall-Zugproben
hochschmelzender Metalle und einige Ergebnisse über de plastische Deformation von NiobEinkristallen. physica status solidi (b), 5(2):421–434, 1964.
[73] G. Taylor and J. W. Christian. The effect of high vacuum purification on the mechanical
properties of niobium single crystals. Acta Metallurgica, 13(11):1216–1218, 1965.
[74] M. S. Duesbery, R. A. Foxall, and P. B. Hirsch. The Plasticity Of Pure Niobium Single
Crystals. Le Journal de Physique Colloques, 27(C3):C3–193–C3–204, 1966.
[75] D. Kang, D. C. Baars, T. R. Bieler, and C. Compton. Study of Slip in High Purity Single
Crystal Nb for Accelerator Cavities. In C. Antoine, S. Bousson, and G. Martinet, editors,
Proceedings of SRF 2013, pages 461–465. JACoW, 2014.
[76] F. Roters, P. Eisenlohr, L. Hantcherli, D. D. Tjahjanto, T. R. Bieler, and D. Raabe. Overview
of constitutive laws, kinematics, homogenization, and multiscale methods in crystal plasticity finite element modeling: Theory, experiments, applications. Acta Materialia, 58:1152–
1211, 2010.

104

[77] E. H. Lee. Elastic-plastic deformation at finite strains. J. Appl. Mech. ASME, 36(1):1–6,
1969.
[78] F. Roters, P. Eisenlohr, C. Kords, D. D. Tjahjanto, M. Diehl, and D. Raabe. DAMASK: The
Düsseldorf Advanced Material Simulation Kit for studying crystal plasticity using an FE
based or a spectral numerical solver. In O. Cazacu, editor, Procedia IUTAM: IUTAM Symposium on Linking Scales in Computation: From Microstructure to Macroscale Properties,
volume 3, pages 3–10, Amsterdam, 2012. Elsevier.
[79] P. Eisenlohr, F. Roters, M. Diehl, C. Kords, and P. Shanthraj. DAMASK Web page. 2014.
[80] A. Prakash and R. A. Lebensohn. Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms. Modelling and Simulation in Materials
Science and Engineering, 17:064010–64016, 2009.
[81] B. Liu, D. Raabe, F. Roters, P. Eisenlohr, and R. A. Lebensohn. Comparison of finite element
and fast Fourier transform crystal plasticity solvers for texture prediction. Modelling and
Simulation in Materials Science and Engineering, 18(8):085005, 2010.
[82] P. Eisenlohr, M. Diehl, R. A. Lebensohn, and F. Roters. A spectral method solution to crystal
elasto-viscoplasticity at finite strains. International Journal of Plasticity, 46:37–53, 2013.
[83] A. K. Kanjarla, R. A. Lebensohn, L. Balogh, and C. N. Tomé. Study of internal lattice strain
distributions in stainless steel using a full-field elasto-viscoplastic formulation based on fast
Fourier transforms. Acta Materialia, 60(6-7):3094–3106, 2012.
[84] M. Montagnat, O. Castelnau, P. D. Bons, S. H. Faria, O. Gagliardini, F. Gillet-Chaulet,
F. Grennerat, A. Griera, R. A. Lebensohn, H. Moulinec, J. Roessiger, and P. Suquet. Multiscale modeling of ice deformation behavior. Journal of Structural Geology, 2013.
[85] D. Ma, P. Eisenlohr, E. Epler, C. A. Volkert, P. Shanthraj, M. Diehl, F. Roters, and D. Raabe.
Crystal plasticity study of monocrystalline stochastic honeycombs under in-plane compression. Acta Materialia, 103:796–808, 2016.
[86] S. L. Wong, M. Madivala, U. Prahl, F. Roters, and D. Raabe. A crystal plasticity model for
twinning- and transformation-induced plasticity. Acta Materialia, 118:140–151, 2016.
[87] N. Lahellec, J. C. Michel, H. Moulinec, and P. Suquet. Analysis of inhomogeneous materials at large strains using fast Fourier transforms. In C. Miehe, editor, IUTAM Symposium
on Computational Mechanics of Solid Materials at Large Strains, volume 108 of Solid Mechanics and its Applications, pages 247–258, Dordrecht, The Netherlands, 2001. Kluwer
Academic Publishers.

105

[88] P. Shanthraj, P. Eisenlohr, M. Diehl, and F. Roters. Numerically robust spectral methods for
crystal plasticity simulations of heterogeneous materials. International Journal of Plasticity,
66:31–45, 2015.
[89] M. Schneider and M. Kabel. The Lipmann-Schwinger equation in elasticity for porous
media. In The Mechanics of Multifunctional Materials, pages 79–82, 2014.
[90] R. A. Lebensohn, J. P. Escobedo, E. K. Cerreta, D. Dennis-Koller, C. A. Bronkhorst, and
J. F. Bingert. Modeling void growth in polycrystalline materials. Acta Materialia, 61(18):
6918–6932, 2013.
[91] R. A. Lebensohn and R. Pokharel. Interpretation of Microstructural Effects on Porosity
Evolution Using a Combined Dilatational/Crystal Plasticity Computational Approach. JOM,
66(3):437–443, 2014.
[92] R. Pokharel, J. Lind, A. K. Kanjarla, R. A. Lebensohn, S. F. Li, P. Kenesei, R. M. Suter,
and A. D. Rollett. Polycrystal Plasticity: Comparison Between Grain-Scale Observations of
Deformation and Simulations. Annual Review of Condensed Matter Physics, 5(1):317–346,
2014.
[93] D. D. Tjahjanto, S. Turteltaub, and A. S. J. Suiker. Crystallographically based model for
transformation-induced plasticity in multiphase carbon steels. Continuum Mechanics and
Thermodynamics, 19(7):399–422, 2008.
[94] P. Shanthraj, B. Svendsen, L. Sharma, F. Roters, and D. Raabe. Elasto-viscoplastic phase
field modelling of anisotropic cleavage fracture. Journal of the Mechanics and Physics of
Solids, 99:19–34, 2017.
[95] D. Peirce, R. J. Asaro, and A. Needleman. Material rate dependence and localized deformation in crystalline solids. Acta Metallurgica, 31(12):1951–1976, 1983.
[96] D. Raabe, M. Sachtleber, Z. Zhao, F. Roters, and S. Zaefferer. Micromechanical and
macromechanical effects in grain scale polycrystal plasticity experimentation and simulation. Acta Materialia, 49(17):3433–3441, 2001.
[97] M. Diehl. A spectral method using fast Fourier transform to solve elastoviscoplastic mechanical boundary value problems. Master’s thesis, München, 2010.
[98] J. W. Hutchinson. Bounds and self-consistent estimates for creep of polycrystalline materials. Proceedings of the Royal Society A, 348:101–127, 1976.
[99] S. R. Kalidindi, C. A. Bronkhorst, and L. Anand. Crystallographic texture evolution in bulk
deformation processing of fcc metals. Journal of the Mechanics and Physics of Solids, 40

106

(3):537–569, 1992.
[100] S. R. Kalidindi. Modeling anisotropic strain hardening and deformation textures in low
stacking fault energy fcc metals. International Journal of Plasticity, 17(6):837–860, 2001.
[101] A. D. Rollett. Strain hardening at large strains in aluminum alloys. PhD thesis, Drexel
University, Philadelphia, 1988., 1988.
[102] A. D. Rollett, R. A. Lebensohn, M. Groeber, Y. Choi, J. Li, and G. S. Rohrer. Stress hot
spots in viscoplastic deformation of polycrystals. Modelling and Simulation in Materials
Science and Engineering, 18(7):074005, 2010.
[103] U. F. Kocks, A. S. Argon, and M. F. Ashby. Thermodynamics and Kinetics of Slip. Progress
in Materials Science, 19:1–291, 1975.
[104] G. I. Taylor. The Mechanism of Plastic Deformation of Crystals. Part I. Theoretical. Proceedings of the Royal Society A, 145(855):362–387, 1934.
[105] P. Franciosi, M. Berveiller, and A. Zaoui. Latent hardening in copper and aluminium single
crystals. Acta Metallurgica, 28(3):273–283, 1980.
[106] R. Madec, B. Devincre, and L. P. Kubin. From Dislocation Junctions to Forest Hardening.
Physical Review Letters, 89(25), 2002.
[107] D. Cereceda, A. Stukowski, M. R. Gilbert, S. Queyreau, L. Ventelon, M-C Marinica, J. M.
Perlado, and J. Marian. Assessment of interatomic potentials for atomistic analysis of static
and dynamic properties of screw dislocations in W. Journal of Physics: Condensed Matter,
25(8):085702–, 2013.
[108] D. Weygand, M. Mrovec, T. Hochrainer, and P. Gumbsch. Multiscale Simulation of Plasticity in bcc Metals. Annual Review of Materials Research, 45(1):369–390, 2015.
[109] S. J. Basinski and Z. S. Basinski. Dislocations in solids, vol. 4, chapter Plastic deformation
and work hardening, pages 261–362. North Holland, Amsterdam, 1979.
[110] S. Aubry, S. P. Fitzgerald, S. L. Dudarev, and W. Cai. Equilibrium shape of dislocation shear
loops in anisotropic ?-Fe. Modelling and Simulation in Materials Science and Engineering,
19(6):065006–, 2011.
[111] C. Teodosiu, J. L. Raphanel, and L. Tabourot. Finite element simulation of the large elastoplastic deformation of multicrystals. In A.A. Balkema, editor, Large plastic deformation :
fundamental aspects and application to metal forming, pages 153–168, Rotterdam, 1993.

107

[112] U. F. Kocks and H. Mecking. Physics and phenomenology of strain hardening: the FCC
case. Progress in Materials Science, 48(3):171–273, 2003.
[113] R. Madec and L. P. Kubin. Second-order junctions and strain hardening in bcc and fcc
crystals. Scripta Materialia, 58(9):767–770, 2008.
[114] V. V. Bulatov, L. L. Hsiung, M. Tang, A. Arsenlis, M. C. Bartelt, W. Cai, J. N. Florando,
M. Hiratani, M. Rhee, G. Hommes, T. G. Pierce, and T. D. de la Rubia. Dislocation multijunctions and strain hardening. Nature, 440(7088):1174–1178, 2006.
[115] U. Essmann and H. Mughrabi. Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities. Philosophical Magazine A, 40(6):731–756,
1979.
[116] L. P. Kubin, B. Devincre, and T. Hoc. Inhibited dynamic recovery and screw dislocation
annihilation in multiple slip of fcc single crystals. Philosophical Magazine, 86(25-26):
4023–4036, 2006.
[117] A. Mapar, D. Kang, T. R. Bieler, F. Pourboghrat, and C. C. Compton. A comparative study
of stress update algorithms for rate-independent and rate-dependent crystal plasticityNb. In
Proceedings of SRF2015, number MOPB057, 2015.
[118] A. Mapar, H. Ghassemi-Armaki, F. Pourboghrat, and K. S. Kumar. A differentialexponential hardening law for non-Schmid crystal plasticity finite element modeling of ferrite single crystals. International Journal of Plasticity, 91:268–299, 2017.
[119] J. L. Bassani and T-Y Wu. Latent Hardening in Single Crystals II. Analytical Characterization and Predictions. Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 435(1893):21–41, 1991.
[120] E. A. Calnan and C. J. B. Clews. LXV. The development of deformation textures in
metals.—Part II. Body-centred cubic metals. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42(329):616–635, 1951.
[121] T. E. Mitchell and W. A. Spitzig. Three-stage hardening in tantalum single crystals. Acta
Metallurgica, 13(11):1169–1179, 1965.
[122] D. P. Ferriss, R. M. Rose, and Wulff. J. Deformation of tantalum single crystals. Trans.
AIME Metall. Soc, 224:975–81, 1962.
[123] L. Dezerald, D. Rodney, E. Clouet, L. Ventelon, and F. Willaime. Plastic anisotropy and
dislocation trajectory in bcc metals. Nature Communications, 7:11695, 2016.

108

[124] H. Lim, J. D. Carroll, C. C. Battaile, T. E. Buchheit, B. L. Boyce, and C. R. Weinberger.
Grain-scale experimental validation of crystal plasticity finite element simulations of tantalum oligocrystals. International Journal of Plasticity, 60:1–18, 2014.
[125] A. Zeghadi, S. Forest, A-F Gourgues, and O. Bouaziz. Ensemble averaging stress–strain
fields in polycrystalline aggregates with a constrained surface microstructure – Part 2: crystal plasticity. Philosophical Magazine, 87(8-9):1425–1446, 2007.
[126] C. Zhang, H. Li, P. Eisenlohr, W. J. Liu, C. J. Boehlert, M. A. Crimp, and T. R. Bieler. Effect
of realistic 3D microstructure in crystal plasticity finite element analysis of polycrystalline
Ti-5Al-2.5Sn. International Journal of Plasticity, 69:21–35, 2015.
[127] M. Diehl, P. Shanthraj, P. Eisenlohr, and F. Roters. Neighborhood influences on stress and
strain partitioning in dual-phase microstructures. An investigation on synthetic polycrystals
with a robust spectral-based numerical method. Meccanica, 51(2):429–441, 2016.
[128] H. F. Poulsen. An introduction to three-dimensional X-ray diffraction microscopy. Journal
of Applied Crystallography, 45(6):1084–1097, 2012.
[129] J. Oddershede, S. Schmidt, H. F. Poulsen, H. O. Sørensen, J. Wright, and W. Reimers.
Determining grain resolved stresses in polycrystalline materials using three-dimensional Xray diffraction. Journal of Applied Crystallography, 43(3):539–549, 2010.
[130] J. Oddershede, B. Camin, S. Schmidt, L. P. Mikkelsen, H. O. Sørensen, U. Lienert, H. F.
Poulsen, and W. Reimers. Measuring the stress field around an evolving crack in tensile
deformed Mg AZ31 using three-dimensional X-ray diffraction. Acta Materialia, 60(8):
3570–3580, 2012.
[131] S. F. Li, J. Lind, C. M. Hefferan, R. Pokharel, U. Lienert, A. D. Rollett, and R. M. Suter.
Three-dimensional plastic response in polycrystalline copper via near-field high-energy Xray diffraction microscopy. Journal of Applied Crystallography, 45(6):1098–1108, 2012.
[132] W. Ludwig, S. Schmidt, E. M. Lauridsen, and H. F. Poulsen. X-ray diffraction contrast tomography: a novel technique for three-dimensional grain mapping of polycrystals. I. Direct
beam case. Journal of Applied Crystallography, 41(2):302–309, 2008.
[133] G. Johnson, A. King, M. G. Honnicke, J. Marrow, and W. Ludwig. X-ray diffraction contrast
tomography: a novel technique for three-dimensional grain mapping of polycrystals. II. The
combined case. Journal of Applied Crystallography, 41(2):310–318, 2008.
[134] M. A. Groeber and M. A. Jackson. DREAM.3D: A Digital Representation Environment for
the Analysis of Microstructure in 3D. Integrating Materials and Manufacturing Innovation,
3(1):5, 2014.

109

[135] H. Moulinec and P. Suquet. A fast numerical method for computing the linear and nonlinear
properties of composites. Comptes rendus de l’Académie des sciences. Série II, Mécanique,
physique, chimie, astronomie, 318:1417–1423, 1994.
[136] H. Moulinec and P. Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Computer Methods in Applied Mechanics
and Engineering, 157(1-2):69–94, 1998.
[137] R. A. Lebensohn. N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform. Acta Materialia, 49(14):2723–2737, 2001.

110