# Modeling Directional Hardening and Intrinsic Size Effects Using a
Dislocation Density-Based Strain Gradient Plasticity Framework

Anirban Patra, Namit Pai, Parhitosh Sharma
Department of Metallurgical Engineering and Materials Science
Indian Institute of Technology Bombay
Mumbai, India - 400076
{Anirban<EMAIL_ADDRESS>

###### Abstract

This work proposes a dislocation density-based strain gradient $J_{2}$
plasticity framework that models the strength contribution due to
Geometrically Necessary Dislocations (GNDs) using a lower order, Taylor
hardening backstress model. An anisotropy factor is introduced to
phenomenologically represent the differential hardening between grains in this
$J_{2}$ plasticity framework. An implicit numerical algorithm is implemented
for the time integration of the finite deformation plasticity model. The
framework is first used to predict directional hardening due to the GND-
induced backstress during cyclic loading. Deformation contours are studied to
understand the substructure attributes contributing to directional hardening.
The framework is then used to predict the intrinsic, grain size-dependent
strengthening of polygrain ensembles. Model predictions of simulations with
different grain sizes are shown to agree with the Hall-Petch effect and also
with Ashby’s model of hardening due to GNDs in polygrain ensembles.

_K_ eywords strain gradient $\cdot$ backstress $\cdot$ $J_{2}$ plasticity
$\cdot$ directional hardening $\cdot$ Hall-Petch effect

##  1 Introduction

Plastic deformation in metallic systems is an inherently length scale-
dependent phenomenon. This is generally attributed to the formation of
Geometrically Necessary Dislocations (GNDs), which accommodate the strain
gradients arising due to discontinuities and heterogeneities in the
microstructure and the specimen during incompatible plastic deformation [1,
2]. GNDs govern the development of long range backstress and ensuing length
scale-dependent strengthening phenomena. These effects are accentuated as the
size decreases, with effects on the mechanical property most evident at micron
and sub-micron length scales. Physical manifestations of these length
scale/size effects are evident in the Hall-Petch effect [3, 4], torsion of
thin wires [5], bending of thin foils [6], micropillar compression [7],
nanoindentation of crystalline materials [8, 9] and constrained shear of films
[10, 11]. Further, it should be noted that the size effects on mechanical
properties have generally been classified into two categories: intrinsic and
extrinsic effects. While the former owe their strength to the microstructural
attributes of the material, such as grain size, the latter owe their strength
to the specimen size itself [12, 13].

Modeling frameworks have attempted to capture these length scale-dependent
effects via the development of non-local models, which account for the spatial
gradient of strain/stress and their contribution to hardening during plastic
deformation. In this regard, the Nye tensor was proposed as a function of the
gradient of plastic strain [1]. This concept has been explored and
mathematically interpreted in various forms to account for the development of
strain gradients, GNDs and associated hardening [14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Recently, grain size-dependent
strengthening has also predicted using strain gradient crystal plasticity
models [32, 33, 34]. There has been extensive research in this field over the
last three decades. Summarizing the features of these modeling studies is
perhaps beyond the scope of the present work and the reader is referred to
[35, 36] for a comprehensive review of the research done in this field.

The focus of the present work is on modeling the development of backstress due
to GNDs and their effect on the size-dependent mechanical properties. In this
regard, various backstress formulations have been proposed [17, 37, 38, 39,
40, 41, 24] that generally account for hardening contributions due to higher
order gradient terms, based on thermodynamic considerations. More recently,
Sangid and co-authors [42, 43] have proposed a lower order, GND density based
backstress model for kinematic hardening in a crystal plasticity framework,
such that the slip system-level backstress scales as the square root of the
GND density (similar to a Taylor hardening model [44]). While these recent
studies [42, 43] primarily focused on studying the GND and backstress
evolution during cyclic loading, size effects on the mechanical properties
were not studied extensively.

In the present work, we explore this concept further and propose a $J_{2}$
plasticity extension of the lower order, Taylor hardening backstress model
given by Kapoor et al. [42] in a dislocation density based framework. We also
introduce an anisotropy factor to phenomenologically represent the
differential hardening between grains with different crystallographic
orientations, in our otherwise isotropic $J_{2}$ plasticity model. Simulations
are first performed to study the directional hardening effects due to the GND-
induced backstress. Simulations are then performed to demonstrate the grain
size-dependent hardening of polygrain ensembles with varying grain sizes.
Based on the analysis of our model predictions, we establish correlations
between the GND density, backstress evolution and grain size-dependent
deformation.

##  2 Model Description

###  2.1 Finite Deformation Kinematics

This finite deformation framework is based on the multiplicative decomposition
of the deformation gradient into the elastic and plastic parts [45]:

$\mathbf{F}=\mathbf{F}^{e}\cdot\mathbf{F}^{p}$ (1)

where, $\mathbf{F}^{p}$ relates the reference configuration to a stress-free,
intermediate configuration and accounts for shear due to plastic deformation,
while $\mathbf{F}^{e}$ relates the intermediate configuration to the current,
deformed deformation and accounts for the elastic deformation.

The plastic deformation gradient is related to the velocity gradient as,
$\mathbf{\dot{F}}^{p}=\mathbf{L}^{p}\cdot\mathbf{F}^{p}$, where the velocity
gradient, $\mathbf{L}^{p}$, is given by [46]:

$\mathbf{L}^{p}=\sqrt{\frac{3}{2}}\dot{\bar{\epsilon}}^{p}\mathbf{N}^{p}$ (2)

Here, $\dot{\bar{\epsilon}}^{p}$ is the effective plastic strain rate and
$\mathbf{N}^{p}$ is the direction of the plastic flow, given by the following
relation:

$\mathbf{N}^{p}=\sqrt{\frac{3}{2}}\frac{\mathbf{S}-\mathbf{\chi}}{\bar{\sigma^{*}}}$
(3)

where, $\mathbf{S}$ is the deviatoric stress tensor, $\mathbf{\chi}$ is the
backstress tensor, and $\bar{\sigma^{*}}$ is the modified effective stress,
defined as

$\bar{\sigma^{*}}=\sqrt{\frac{3}{2}(\mathbf{S}-\mathbf{\chi}):(\mathbf{S}-\mathbf{\chi})}$
(4)

###  2.2 Kinetics and Substructure Evolution

The effective plastic strain rate is modeled using a Kocks-type, thermally
activated flow rule [47]:

$\dot{\bar{\epsilon}}^{p}=\dot{\bar{\epsilon}}^{p}_{0}\exp\left(\frac{-\Delta
F_{g}}{kT}\left(1-\left(\frac{\bar{\sigma^{*}}-S_{a}}{S_{t}}\right)^{p}\right)^{q}\right);\bar{\sigma}>S_{a}$
(5)

where, $\dot{\bar{\epsilon}}^{p}_{0}$ is the reference strain rate, $\Delta
F_{g}$ is the activation energy for dislocation glide, $S_{a}$ is the athermal
slip resistance, $S_{t}$ is the thermal slip resistance, and $p$ and $q$ are
parameters used to model the shape of the activation enthalpy curve. The
athermal slip resistance, $S_{a}$, is given by:

$S_{a}=M(\tau_{0}+k_{IH}Gb\sqrt{\rho_{SSD}})$ (6)

where, $\tau_{0}$ is threshold slip resistance, $k_{IH}$ is the Taylor
hardening coefficient [44] associated with isotropic hardening, $G$ is the
shear modulus, $b$ is the Burgers vector magnitude, and $\rho_{SSD}$ is the
Statistically Stored Dislocation (SSD) density. Here, we introduce an
anisotropy factor, $M$, to phenomenologically represent the anisotropy in the
yield stress of crystalline materials due to their crystallographic
orientations. This term may be considered to be representative of the Taylor
factor, or the inverse of the Schmid factor, in this macroplasticity framework
and is randomly assigned for each grain. As a first order approximation, $M$
is assumed to be constant and its evolution with plastic deformation is
neglected. We also note that the anisotropy factor may be expected to evolve
with deformation (cf. [48]).

The evolution of SSDs is modeled using a Kocks-Mecking-Estrin type formulation
[49, 50], albeit with the additional consideration for GNDs in the generation
of SSDs [37]:

$\dot{\rho}_{SSD}=\frac{k_{mul}}{b}\sqrt{\rho_{SSD}+\rho_{GND}}\dot{\bar{\epsilon}}^{p}-k_{rec}\rho_{SSD}\dot{\bar{\epsilon}}^{p}$
(7)

where, $k_{mul}$ is the dislocation multiplication rate constant, and
$k_{rec}$ is the recovery constant. The first term in Equation (7) represents
the multiplication of SSDs at existing dislocations, while the second term
represents the rate of annihilation of SSDs due to recovery processes.

###  2.3 Geometrically Necessary Dislocation Density and Backstress Tensor

In this finite deformation framework, we adopt the definition of the Nye
tensor proposed by Dai [51] as

$\mathbf{\Lambda}=-(\mathbf{\nabla}\times\mathbf{F}^{pT})^{T}$ (8)

Further, the Geometrically Necessary Dislocation (GND) density is defined as

$\rho_{GND}=\frac{1}{b}||\mathbf{\Lambda}||$ (9)

where, $||||$ denotes the L2 norm of the respective tensor quantity.

In the present work, we implement Equation (8) in its rate form, such that

$\mathbf{\dot{\Lambda}}=-(\mathbf{\nabla}\times\mathbf{\dot{F}}^{pT})^{T};\dot{\rho}_{GND}=\frac{1}{b}||\mathbf{\dot{\Lambda}}||$
(10)

Further, we propose the backstress tensor and its rate form as

$\mathbf{\chi}=k_{KH}Gb\sqrt{\rho_{GND}}\mathbf{N}^{p}\Leftrightarrow\mathbf{\dot{\chi}}=k_{KH}Gb\frac{\dot{\rho}_{GND}}{2\sqrt{\rho_{GND}}}\mathbf{N}^{p}$
(11)

where, $k_{KH}$ is the Taylor hardening coefficient associated with kinetic
hardening due to GNDs.

This form of the backstress tensor is a $J_{2}$ plasticity extension of the
Taylor hardening contribution of the GND density to the slip system-level
directional hardening proposed for a crystal plasticity framework [42].
Physically, this model represents the development of backstress due to GNDs
along the direction of plastic deformation, $\mathbf{N}^{p}$. The rate (or the
incremental) form of this backstress model ensures that the backstress
contribution due to the pile up of GND density along a certain loading
direction is remnant in the history-dependent backstress tensor, even after
the direction of loading is changed, for example, during cyclic loading. In
the context of prior macroplasticity frameworks, our formulation is a
deviation from the mechanism-based strain gradient plasticity models that
proposed the hardening due to the total density of SSDs and GNDs (cf. [52, 15]
and its variants). Essentially, our model delineates the contribution of the
GNDs towards the development of backstress in the direction of plastic
deformation, while the SSDs are assumed to contribute primarily to isotropic
hardening. Note, however, that higher order strain gradient frameworks do
separate the strengthening contribution of SSDs and GNDs into isotropic
hardening and kinematic hardening, respectively [37, 53, 39], although their
physical implications are different.

###  2.4 Implicit Time Integration and Numerical Implementation

We propose an implicit algorithm for the time integration of the constitutive
model. This algorithm is inspired by the implicit slip rate integration
algorithms [54, 55, 56, 57], which were given for crystal plasticity
formulations. Here we adapt the same for a $J_{2}$ plasticity formulation.

Given the deformation gradient, $\mathbf{F}$, at any time step, this algorithm
numerically estimates the effective plastic strain rate,
$\dot{\bar{\epsilon}}^{p}$, via a Newton-Raphson algorithm. At any given time
step, a function may be formulated such that the function,
$f(\dot{\bar{\epsilon}}^{p})$, given by:

$f(\dot{\bar{\epsilon}}^{p})_{i+1}=f(\dot{\bar{\epsilon}}^{p})_{i}+\frac{\partial{f(\dot{\bar{\epsilon}}^{p})}}{\partial{\dot{\bar{\epsilon}}^{p}}}\Delta\dot{\bar{\epsilon}}^{p}$
(12)

has to be minimized iteratively, i.e, $f\rightarrow 0$. Here, the subscript,
$i$, denotes the number of iterations at any given time step. In order to
satisfy this equality,

$f(\dot{\bar{\epsilon}}^{p})_{i}=-\frac{\partial{f(\dot{\bar{\epsilon}}^{p})}}{\partial{\dot{\bar{\epsilon}}^{p}}}\Delta\dot{\bar{\epsilon}}^{p}\Leftrightarrow\Delta\dot{\bar{\epsilon}}^{p}=-\frac{f(\dot{\bar{\epsilon}}^{p})_{i}}{\frac{\partial{f(\dot{\bar{\epsilon}}^{p})}}{\partial{\dot{\bar{\epsilon}}^{p}}}}$
(13)

Thus,
$\dot{\bar{\epsilon}}^{p}_{i+1}=\dot{\bar{\epsilon}}^{p}_{i}+\Delta\dot{\bar{\epsilon}}^{p}$
provides the converged value of the effective plastic strain rate at any given
time step.

This method relies on accurate calculation of the partial derivatives
associated with the
$\frac{\partial{f(\dot{\bar{\epsilon}}^{p})}}{\partial{\dot{\bar{\epsilon}}^{p}}}$
term. Taking derivatives and expanding the individual terms gives an
expression of the form:

$\begin{split}\frac{\partial{f(\dot{\bar{\epsilon}}^{p})}}{\partial{\dot{\bar{\epsilon}}^{p}}}=1-\dot{\bar{\epsilon}}^{p}\left(q\frac{\Delta{F_{g}}}{kT}\left(1-\left(\frac{\bar{\sigma^{*}}-S_{a}}{S_{t}}\right)^{p}\right)^{q-1}\right)\left(p\left(\frac{\bar{\sigma^{*}}-S_{a}}{S_{t}}\right)^{p-1}\right)\\\
\left(\frac{1}{S_{t}}\left(\frac{\partial{\bar{\sigma^{*}}}}{\partial{\dot{\bar{\epsilon}}^{p}}}-\frac{\partial{S_{a}}}{\partial{\dot{\bar{\epsilon}}^{p}}}\right)-\left(\frac{\bar{\sigma^{*}}-S_{a}}{S_{t}^{2}}\right)\frac{\partial{S_{t}}}{\partial{\dot{\bar{\epsilon}}^{p}}}\right)\end{split}$
(14)

The above expression contains derivatives of various terms with respect to
$\dot{\bar{\epsilon}}^{p}$, which are described in the following. We begin
with $\frac{\partial{\bar{\sigma^{*}}}}{\partial{\dot{\bar{\epsilon}}^{p}}}$,
which may be written as [55]

$\frac{\partial{\bar{\sigma^{*}}}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{\partial{\bar{\sigma^{*}}}}{\partial{\mathbf{\sigma}}}:\frac{\partial{\mathbf{\sigma}}}{\partial{\mathbf{\sigma}^{pk2}}}:\frac{\partial{\mathbf{\sigma}^{pk2}}}{\partial{\mathbf{E}^{e}}}:\frac{\partial{\mathbf{E}^{e}}}{\partial{\mathbf{F}^{p}}}:\frac{\partial{\mathbf{F}^{p}}}{\partial{\mathbf{L}^{p}}}:\frac{\partial{\mathbf{L}^{p}}}{\partial{\dot{\bar{\epsilon}}^{p}}}$
(15)

Here, $\mathbf{\sigma}$ is the Cauchy stress, $\mathbf{\sigma}^{pk2}$ is the
second Piola-Kirchoff stress, and $\mathbf{E}^{e}$ is the elastic Green
strain. Following McGinty [55], this expression may be simplified as

$\frac{\partial{\bar{\sigma^{*}}}}{\partial{\dot{\bar{\epsilon}}^{p}}}\approx-\sqrt{\frac{3}{2}}\mathbf{N}^{p}:\mathbf{C}:\sqrt{\frac{3}{2}}\mathbf{N}^{p}$
(16)

where, $\mathbf{C}$ is the fourth rank elasticity tensor. Further, the
derivative of the athermal stress with respect to the effective plastic strain
rate is given as

$\frac{\partial{S_{a}}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{k_{IH}Gb}{2\sqrt{\rho_{SSD}}}\frac{\partial{\rho_{SSD}}}{\partial{\dot{\bar{\epsilon}}^{p}}}$
(17)

where,

$\frac{\partial{\rho_{SSD}}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{\left(\frac{k_{mul}}{b}\sqrt{\rho_{SSD}+\rho_{GND}}-k_{dyn}\rho_{SSD}+\frac{k_{mul}}{2b\sqrt{\rho_{SSD}+\rho_{GND}}}\frac{\partial\rho_{GND}}{\partial{\dot{\bar{\epsilon}}^{p}}}\right)\Delta
t}{1-\frac{k_{mul}}{2b\sqrt{\rho_{SSD}+\rho_{GND}}}\dot{\bar{\epsilon}}^{p}\Delta
t+k_{dyn}\dot{\bar{\epsilon}}^{p}\Delta t}$ (18)

and,

$\frac{\partial{\rho_{GND}}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{1}{b}||\frac{\partial\mathbf{\dot{\Lambda}}}{\partial{\dot{\bar{\epsilon}}^{p}}}||{\Delta
t}$ (19)

Here,

$\frac{\partial\mathbf{\dot{\Lambda}}}{\partial{\dot{\bar{\epsilon}}^{p}}}=-\left(\mathbf{\nabla}\times\frac{\partial\mathbf{\dot{F}}^{p}}{\partial{\dot{\bar{\epsilon}}^{p}}}^{T}\right)^{T}$
(20)

Calculating the derivative of $\mathbf{\dot{F}^{p}}$ with respect to
$\dot{\bar{\epsilon}}^{p}$ is not straightforward. Following McGinty [55], we
write $\mathbf{F^{p}}$ at time step $t+\Delta t$ as:

$\mathbf{F}^{p}_{t+\Delta t}=exp(\mathbf{L}^{p}_{0}\Delta
t)\cdot\mathbf{F}^{p}_{t}$ (21)

where, $\mathbf{L^{p}_{0}}$ is the plastic velocity gradient in the
intermediate configuration. This expression may be expanded using the Cayley
Hamilton theorem [58, 55], i.e,

$exp(\mathbf{L}^{p}_{0}\Delta
t)=\mathbf{I}+\frac{sin\phi}{\phi}\mathbf{L}^{p}_{0}\Delta
t+\frac{1-cos\phi}{\phi^{2}}(\mathbf{L}^{p}_{0}\cdot\mathbf{L}^{p}_{0}){\Delta
t}^{2}$ (22)

where,

$\phi=\left(\sqrt{\frac{1}{2}\mathbf{L}^{p}_{0}:\mathbf{L}^{p}_{0}}\right)\Delta
t$ (23)

In order to calculate $\mathbf{\dot{F}}^{p}$, we take the finite difference of
Equation (21) with respect to $\Delta t$ such that

$\mathbf{\dot{F}}^{p}=\frac{\mathbf{F}^{p}_{t+\Delta
t}-\mathbf{F}^{p}_{t}}{\Delta t}=\frac{1}{\Delta
t}\left(exp(\mathbf{L}^{p}_{0}\Delta
t)-\mathbf{I}\right)\cdot\mathbf{F}^{p}_{t}$ (24)

Using the expression for $exp(\mathbf{L}^{p}_{0}\Delta t)$ given above and
after some mathematical manipulation, we arrive at

$\begin{split}\frac{\partial\mathbf{\dot{F}}^{p}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{1}{\Delta
t}\left(\frac{\partial(exp(\mathbf{L}^{p}_{0}\Delta
t))}{\partial{\dot{\bar{\epsilon}}^{p}}}\cdot\mathbf{L}^{p-1}_{0}\cdot\mathbf{\dot{F}}^{p}+(exp(\mathbf{L}^{p}_{0}\Delta
t)-\mathbf{I})\cdot\frac{\partial\mathbf{L}^{p-1}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}\cdot\mathbf{\dot{F}}^{p}\right)\cdot\\\
\left(\mathbf{I}-\frac{1}{\Delta t}(exp(\mathbf{L}^{p}_{0}\Delta
t)-\mathbf{I})\cdot\mathbf{L}^{p-1}_{0}\right)^{-1}\end{split}$ (25)

Further, the individual derivatives in the above expression are given as:

$\begin{split}\frac{\partial(exp(\mathbf{L}^{p}_{0}\Delta
t))}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{sin\phi}{\phi}\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}\Delta
t+\frac{1-cos\phi}{\phi^{2}}\left(\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}\cdot\mathbf{L}^{p}_{0}+\mathbf{L}^{p}_{0}\cdot\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}\right)\Delta
t^{2}\\\
+\left(\frac{cos\phi}{\phi}-\frac{sin\phi}{\phi^{2}}\right)\frac{\partial\phi}{\partial{\dot{\bar{\epsilon}}^{p}}}(\mathbf{L}^{p}_{0}\Delta
t)+\left(\frac{\sin\phi}{\phi^{2}}-\frac{1-cos\phi}{2\phi^{3}}\right)\frac{\partial\phi}{\partial{\dot{\bar{\epsilon}}^{p}}}\mathbf{L}^{p}_{0}\cdot\mathbf{L}^{p}_{0}\Delta
t^{2}\end{split}$ (26)
$\frac{\partial\phi}{\partial{\dot{\bar{\epsilon}}^{p}}}=\frac{\Delta
t}{2\sqrt{\frac{1}{2}\mathbf{L}^{p}_{0}:\mathbf{L}^{p}_{0}}}\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}:\mathbf{L}^{p}_{0}$
(27)
$\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}=\sqrt{\frac{3}{2}}\mathbf{N}^{p}$
(28)
$\frac{\partial\mathbf{L}^{p-1}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}=-\mathbf{L}^{p-1}_{0}\cdot\frac{\partial\mathbf{L}^{p}_{0}}{\partial{\dot{\bar{\epsilon}}^{p}}}\cdot\mathbf{L}^{p-1}_{0}$
(29)

Here, Equation (29) is derived using linear algebra concepts [59]. The above
set of equations provide all the terms required for the implicit Newton-
Raphson algorithm. We solve Equation (13) iteratively until the change in the
effective plastic strain rate, $\Delta\dot{\bar{\epsilon}}^{p}$, is below a
prescribed tolerance.

The constitutive model has been implemented as a user-defined material model
and interfaced with the open source finite element (FE) code, Multiphysics
Object-Oriented Simulation Environment (MOOSE) [60]. MOOSE provides the
material model with an increment of the deformation gradient at an integration
point, based on the applied global boundary conditions. The material model
performs an implicit update of the stress, internal state variables and the
consistent tangent stiffness tensor at the Gauss point, which is then passed
back to MOOSE to check for global convergence. While the FE framework
considers the displacement variables as degrees of freedom, the spatial
gradient terms are implemented using the auxiliary variable interface in
MOOSE. Note that these auxiliary variables, computed at the Gauss points, are
not used in the computation of the global Jacobian, hence reducing computation
costs significantly. In our implementation, $\mathbf{\dot{F}}^{p}$ is stored
as an auxiliary variable and its spatial gradient calculated using standard
finite element shape functions of desired order for determining the rate of
Nye tensor (cf. Equation (10)).

##  3 Model Parameters and Simulation Details

The model parameters chosen in this work are given in Table 1. In this example
application, we have not calibrated our model to the mechanical response of
any specific material and have used parameters that provide a representative
stress-strain response and flow strength expected of metallic systems.

The elastic constants and Burgers vector magnitude are chosen to be
representative of a fcc metal. The Kocks-type thermally activated flow rule
has been extensively used in crystal plasticity frameworks to model the
temperature- and strain rate-dependent response of various metals and alloys
[61, 62, 63] and our choice of flow rule parameters is inspired from these
earlier studies. The Taylor hardening coefficient associated with isotropic
hardening due to SSDs is chosen to be $0.2$ [44]. Since the focus of this work
is to highlight the contribution of GNDs to the strengthening, we have
intentionally chosen a low value of the initial SSD density as $1.0\times
10^{5}mm^{-2}$. Further, the dislocation multiplication constant for SSDs is
also chosen to have a low value of $0.01$, such that the SSD density does not
increase significantly during deformation.

As a first order approximation, the initial GND density has been assumed to be
zero. This may considered as representative of an annealed material, which has
no residual deformation. In the present work, $k_{KH}$ is the most important
parameter governing the development of backstress and its contribution to
size-dependent strengthening. In the following sections, we have demonstrated
model predictions with two different values for this parameter: $1.2$ and
$0.8$. In Section 4.2.3, we have also provided a rationale for choosing
appropriate values for this parameter.

Table 1: Constitutive model parameters. Parameter | Value | Meaning
---|---|---
$E$ | 128 GPa | Young’s modulus
$\nu$ | $0.34$ | Poisson’s ratio
$G$ | 47.76 GPa | Shear modulus
$b$ | $2.56\times 10^{-10}m$ | Burgers vector magnitude
$\dot{\bar{\epsilon}}^{p}_{0}$ | $1.0\times 10^{-2}s^{-1}$ | Reference strain rate
$\Delta F_{g}$ | $4G{b}^{3}$ | Activation energy for dislocation glide
$S_{t}$ | 400 MPa | Thermal slip resistance
$\tau_{0}$ | 50 MPa | Threshold slip resistance
$k_{IH}$ | 0.2 | Isotropic hardening coefficient due to SSDs
$k_{KH}$ | 1.2 and 0.8 | Kinematic hardening coefficient due to GNDs
$k_{mul}$ | 0.01 | Dislocation multiplication rate constant
$k_{rec}$ | 500 | Dislocation recovery constant
$\rho_{SSD}^{0}$ | $1.0\times 10^{5}mm^{-2}$ | Initial SSD density
$\rho_{GND}^{0}$ | $0$ | Initial GND density

We have performed simulations of polygrain ensembles with varying grain sizes.
The initial microstructures were instantiated using the open source Voronoi
tessellation software, Neper [64], for a 2D simulation domain of $100\times
100\mu m$. The microstructures were meshed using 2D six node, triangular
finite elements with quadratic interpolation, and a mesh size of $\approx
1.5\mu m$. There were approximately $10,000-20,000$ finite elements in each of
these microstructures. The microstructures generated by Neper were converted
to the exodus file format (required for MOOSE) using the commercial meshing
software, Trelis [65]. Further, the anisotropy factor, $M$, was randomly
assigned for each grain, such that it ranges from $1/0.27$ to $1/0.49$. Here,
$0.27$ and $0.49$ are representative of the lower and upper bounds of the
Schmid factor generally observed in cubic crystals deforming by octahedral
slip. As mentioned earlier, $M$ is not allowed to evolve during the
simulation. The primary purpose of this parameter is to introduce
intergranular heterogeneity in the $J_{2}$ plasticity simulations.

Axi-symmetric boundary conditions were used such that the left edge of the
simulation domain was constrained in the x-direction, while the bottom edge
was constrained in the y-direction. The corner node common to both these edges
was constrained in all degrees of freedom to prevent rigid body motion.
Displacement-controlled uniaxial loading was applied on the top edge at a
nominal strain rate of $1.0\times 10^{-4}s^{-1}$.

A schematic of the loading and boundary conditions along with a representative
microstructure is shown in Figure 1.

Figure 1: Representative microstructure with an average grain size of $8\mu m$
in the simulation domain of $100\times 100\mu m^{2}$, along with a schematic
of the loading and boundary conditions.

##  4 Results and Discussion

###  4.1 Directional Hardening Effects

We first present model predictions of directional hardening. We have used the
microstructure with $8\mu m$ average grain size for these simulations. The
microstructure was cyclically loaded in tension and compression up to $0.01$
strain amplitude for $500s$ at a nominal strain rate of $1.0\times
10^{-4}s^{-1}$. In order to demonstrate the effect of $k_{KH}$ on directional
hardening, we performed three simulations for the same initial microstructure.
While the first two simulations were performed with $k_{KH}=1.2$ and
$k_{KH}=0.8$, the third simulation was performed with $k_{KH}=0$ and using a
local version of the $J_{2}$ plasticity model. Essentially, all the strain
gradient-dependent terms, $\rho_{GND}$ and $\mathbf{\chi}$ were assumed to be
zero in this third simulation.

####  4.1.1 Aggregate Properties

Figure 2 (a) and (b) show the aggregate response, plotted in terms of
$\sigma_{22}$ versus $\epsilon_{22}$, and the effective stress,
$\bar{\sigma}$, versus effective cyclic strain, $\bar{\epsilon}_{cyclic}$, for
simulations with different values of $k_{KH}$. It can be seen that the flow
stress increases with increasing $k_{KH}$ in the 1st cycle of loading. For
example, $\bar{\sigma}$ is $494MPa$ for $k_{KH}=1.2$, $487MPa$ for
$k_{KH}=0.8$, and $475MPa$ for $k_{KH}=0$ after
$\bar{\epsilon}_{cyclic}=0.01$. Subsequent to tension, when the
microstructures were loaded in compression, the flow stress for the simulation
with $k_{KH}=0$ is almost the same as that in tension, indicating only
isotropic hardening. For $k_{KH}=0.8$ and $1.2$, the flow stresses in
compression are initially lower, followed by subsequent hardening. As
discussed earlier, model parameters related to isotropic hardening due to SSDs
were intentionally kept low in our simulations. Hence, the hardening observed
in our simulations is primarily due to the GND-dominated kinematic hardening.
The net hardening due to GND density may appear to be relatively small. For
example, the hardening after $\bar{\epsilon}_{cyclic}=0.01$ is $\approx 19MPa$
using $k_{KH}=1.2$ and $\approx 12MPa$ using $k_{KH}=0.8$, as compared to the
simulations with $k_{KH}=0$. However, it will be shown in the later sections
that this size-dependent strength contribution is consistent with the values
expected of fcc metals. Further, the value of the parameter, $k_{KH}$, may
simply be increased, if a higher directional hardening is desired.

Figure 2: Plot of (a) $\sigma_{22}$ versus $\epsilon_{22}$, (b) $\bar{\sigma}$
versus $\bar{\epsilon}_{cyclic}$, (c) average $\rho_{GND}$ and $\rho_{SSD}$
versus $\bar{\epsilon}_{cyclic}$, and (d) $\bar{\chi}$ versus
$\bar{\epsilon}_{cyclic}$ for simulations with different values of $k_{KH}$.
In (c), solid lines are used to denote $\rho_{GND}$, while dotted lines are
used to denote $\rho_{SSD}$. Note that the same color scheme has been used to
represent the different values of $k_{KH}$ in all the plots.

Figure 2 (c) shows the evolution of the average $\rho_{GND}$ and $\rho_{SSD}$
as a function of $\bar{\epsilon}_{cyclic}$ for the three cases shown in Figure
2 (a). The GND density is shown using solid lines, while the SSD density is
shown using dotted lines. It can be seen that the GND density is slightly
higher for $k_{KH}=1.2$, as compared to that for $k_{KH}=0.8$. While the
strain gradient itself is expected to be similar in both cases, some
additional GNDs may have developed to accommodate the stress concentrations
due to higher backstress. This will be evident in the deformation contours
discussed in the next section. Further, the SSD density is similar for both
cases, while it is 5-7 times higher than the case with $k_{KH}=0$ (no strain
gradient). Since the model parameters for SSD evolution have been
intentionally kept low, we do not see significant evolution in the SSD density
for the case with no strain gradient. For $k_{KH}=1.2$ and $0.8$, there is
some evolution of SSDs, which may be primarily attributed to the
multiplication of SSDs at GND segments, modeled using the 1st term in Equation
(7).

Before moving forward, we would like to point out that the sign of the GND
density is not considered in our $J_{2}$ plasticity framework. For example,
GNDs with opposite sign may be expected to develop during reverse loading and
hence contribute to a reduction in the net GND density. However, the GND
density is computed from the norm of the Nye tensor in our model (cf. Equation
(10)) and hence there is no way to account for the sign of the GNDs in our
framework. Nonetheless, the backstress tensor accounts for the directionality
of the hardening caused due the GND density. This is evident from the
evolution of the effective backstress, $\bar{\chi}$, with
$\bar{\epsilon}_{cyclic}$, in Figure 2(d), which shows that the effective
backstress decreases once plastic deformation commences during compression
loading, and then increases again during tension loading. This is physically
representative of the hardening contribution due to backstress changing signs,
once the loading direction is reversed. It should also be noted that the
backstress scales almost linearly with the parameter $k_{KH}$. For example,
the average value of $\bar{\chi}$ is $18.68MPa$ after
$\bar{\epsilon}_{cyclic}=0.01$ using $k_{KH}=1.2$, while it is $11.67MPa$
using $k_{KH}=0.8$, which gives a scaling of $\approx 2/3$. The GND density
and backstress are of course zero using $k_{KH}=0$.

In summary, the results presented in this section show the capability of the
model to predict directional hardening and backstress evolution under cyclic
loading conditions. While the GND density is not signed in this $J_{2}$
plasticity framework, the backstress tensor appropriately accounts for the
directional hardening effects due to the GND density.

####  4.1.2 Deformation and Substructure Contours

In order to highlight the ability of the model to predict heterogeneous
deformation, we present deformation contours from the cyclic loading
simulations in this section. For this, we choose a small set of grains from
the initial microstructure shown in Figure 1. These grains are highlighted on
the left side of Figure 3. The contours of effective strain, $\bar{\epsilon}$,
at different stages of cyclic loading are shown for the three cases in Figure
3. Essentially, we have plotted the strain contours at the end of tension
loading during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.01$), the end of
compression loading during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.03$), and
at the end of tension loading during the 2nd cycle
($\bar{\epsilon}_{cyclic}=0.05$). The corresponding contours of $\sigma_{22}$
are shown in Figure 4. The effective strain is expected to be similar at these
three stages of deformation.

Firstly, it is interesting to note that the model predicts heterogeneous
deformation in the different grains using this $J_{2}$ plasticity framework.
This can be clearly seen from the $\bar{\epsilon}$ and $\sigma_{22}$ contours,
where certain grains have higher strains, while other grains have higher
stresses. Generally, such effects can be captured using crystal plasticity
models. In our model, this is attributed to the introduction of the anisotropy
factor, $M$ (cf. Equation 6), which is used to phenomenologically represent
the differential hardening between grains. This also promotes strain and
stress concentrations at the interfaces between grains, for example in the
region indicated by A in Figures 3 and 4. It should also be noted that such
effects are observed in all cases, with and without strain gradients.

The effect of strain gradients on directional hardening can be clearly
observed in the strain contours for the different cases at
$\bar{\epsilon}_{cyclic}=0.05$. This is highlighted with circular markers in
Figure 3. In the absence of strain gradients and directional hardening
($k_{KH}=0$), the only heterogeneity in deformation may arise from the
differential hardening between grains, due to the anisotropy factor, $M$.
However, in the presence of strain gradients and with increasing values of
$k_{KH}$, the strain heterogeneities increase in the same regions. These
strain heterogeneities also increase with applied deformation, with early
precursors of strain localization being observed at
$\bar{\epsilon}_{cyclic}=0.05$.

Figure 3: Contours of effective strain, $\bar{\epsilon}$, at the end of
tension loading during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.01$), the end
of compression loading during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.03$),
and at the end of tension loading during the 2nd cycle
($\bar{\epsilon}_{cyclic}=0.05$) for simulations with different values of
$k_{KH}$. The set of grains chosen for these contours is shown on the left
side.

The $\sigma_{22}$ contours in Figure 4 show that with increasing $k_{KH}$, the
stresses are not fully reversed when the loading is reversed from tension to
compression and vice-versa. For example, the maximum and minimum values of
$\sigma_{22}$ are $\approx 550MPa$ and $\approx 450MPa$ at the end of the
tension loading for $k_{KH}=1.2$, while the same at the end of the compression
cycle are $\approx-440MPa$ and $\approx-530MPa$, respectively. For $k_{KH}=0$
and no strain gradients, the corresponding absolute values are identical
between tension and compression. Interestingly, the grains with higher
relative stress in tension have lower relative stress in compression, and
vice-versa.

Figure 4: Contours of $\sigma_{22}$ at the end of tension loading during the
1st cycle ($\bar{\epsilon}_{cyclic}=0.01$), the end of compression loading
during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.03$), and at the end of
tension loading during the 2nd cycle ($\bar{\epsilon}_{cyclic}=0.05$) for
simulations with different values of $k_{KH}$. The set of grains chosen for
these contours is shown on the left side of Figure 3.

Figure 5 shows contours of the GND density, $\rho_{GND}$, at different stages
of applied deformation for $k_{KH}=1.2$ and $k_{KH}=0.8$, while Figure 6 shows
contours of the effective backstress, $\bar{\chi}$, for the same. As can be
seen, the GND density mostly develops at the grain interfaces and triple
junctions, which are regions of stress concentration and incompatible
deformation between grains. Further, the GND density intensifies at these
regions with increasing applied deformation. As mentioned earlier, our $J_{2}$
model does not account for the sign of the GND density and although the GND
density may increase with applied deformation, their net effect may be to
reduce the backstress. The same can be seen from the effective backstress
contours in Figure 6, where the effective backstress generally reduces during
the compression cycle in all the grains. In fact, the maximum and minimum
values of $\bar{\chi}$ at the end of tension loading during the 2nd cycle
($\bar{\epsilon}_{cyclic}=0.05$) are lower than the corresponding values after
the end of tension loading during the 1st cycle
($\bar{\epsilon}_{cyclic}=0.01$). This is primarily due to the fact that the
overall backstress hardening during compression loading negates the same
developed during tension loading in the 1st cycle.

Figure 5: Contours of $\rho_{GND}$ at the end of tension loading during the
1st cycle ($\bar{\epsilon}_{cyclic}=0.01$), the end of compression loading
during the 1st cycle ($\bar{\epsilon}_{cyclic}=0.03$), and at the end of
tension loading during the 2nd cycle ($\bar{\epsilon}_{cyclic}=0.05$) for
simulations with different values of $k_{KH}$. The set of grains chosen for
these contours is shown on the left side of Figure 3. Figure 6: Contours of
$\bar{\chi}$ at the end of tension loading during the 1st cycle
($\bar{\epsilon}_{cyclic}=0.01$), the end of compression loading during the
1st cycle ($\bar{\epsilon}_{cyclic}=0.03$), and at the end of tension loading
during the 2nd cycle ($\bar{\epsilon}_{cyclic}=0.05$) for simulations with
different values of $k_{KH}$. The set of grains chosen for these contours is
shown on the left side of Figure 3.

The results presented in this section demonstrate the ability of the framework
to simulate directional hardening due to GNDs. In the following section, we
focus on establishing size effects using our model.

###  4.2 Grain Size Effects

In order to establish the grain size-dependent strengthening relations, we
have performed simulations with five different mean grain sizes: $4\mu m$,
$8\mu m$, $20\mu m$, $30\mu m$, and $40\mu m$. Essentially, we span across an
order of magnitude in terms of the mean grain size. Note that while this
represents the mean grain size provided to Neper [64], there is some variation
in the actual grain sizes obtained during the tessellations. Further, the
domain size was kept constant at $100\times 100\mu m$ for all simulations. All
simulations were performed till $0.04$ applied strain for the material loaded
in uniaxial tension. We first present the aggregate properties, followed by
the deformation contours and finally estimate the grain size-dependent scaling
relations.

####  4.2.1 Aggregate Properties

Figure 7(a) shows the aggregate response, plotted in terms of the effective
stress, $\bar{\sigma}$, versus effective strain, $\bar{\epsilon}$, for
simulations with different average grain sizes. In order to highlight the
grain size-dependent response, the plot of $\bar{\sigma}$ versus
$\bar{\epsilon}^{p}$ given by $\approx\bar{\epsilon}-\bar{\sigma}/E$ is shown
in Figure 7(b). The corresponding average GND and SSD densities are plotted as
a function of $\bar{\epsilon}$ in Figure 7(c), and the average effective
backstress, $\bar{\chi}$, is plotted as a function of $\bar{\epsilon}$ in
Figure 7(d).

Figure 7: Plot of (a) $\bar{\sigma}$ versus $\bar{\epsilon}$, (b)
$\bar{\sigma}$ versus
$\bar{\epsilon}^{p}\approx\bar{\epsilon}-\bar{\sigma}/E$, (c) average
$\rho_{GND}$ and $\rho_{SSD}$ versus $\bar{\epsilon}$, and (d) $\bar{\chi}$
versus $\bar{\epsilon}$ for simulations with different grain sizes. In (c),
solid lines are used to denote $\rho_{GND}$, while dotted lines are used to
denote $\rho_{SSD}$. Note that the same color scheme has been used to
represent the different grain sizes in all the plots.

It can be seen from Figure 7(a) and (b) that the flow stress of the material
increases as the grain size decreases. These results indicate that the model
is able to qualitatively capture the grain size-dependent strengthening [3,
66]. For example, the flow stress at the beginning of plastic deformation is
$\approx 480MPa$ for grain size $4\mu m$, while it is $\approx 465MPa$ for
grain size $40\mu m$. After 0.04 applied strain, the flow stress is $\approx
530MPa$ for grain size of $4\mu m$, while it is $\approx 482MPa$ for grain
size of $40\mu m$. Further, Figure 7(c) and (d) show that the GND density and
backstress initially increase at a high rate at the beginning of plastic
deformation and these values tend to saturate during the later stages of
deformation. An effective backstress of $\approx$ 50 MPa developed after 0.04
applied strain for grain size $4\mu m$, while the same for grain size $40\mu
m$ is less than 20 MPa. Also note that the rate of increase of GND density and
backstress is higher for the lower grain sizes and vice-versa. As mentioned in
Section 3, a low value was used for the dislocation multiplication constant,
$k_{mul}$, in our simulations, which resulted in a relatively lower average
SSD density as compared to the GND density. Increasing the value of this
parameter may contribute to increase in SSD density and increased strain
hardening during plastic deformation. These ensemble average quantities are
used to derive the size-dependent scaling relations in Section 4.2.3.

####  4.2.2 Deformation and Substructure Contours

Figure 8 shows the grain maps and the corresponding contours of
$\bar{\epsilon}$ and $\bar{\sigma}$ after 0.04 applied strain for the
different grain sizes. The corresponding contours of $\rho_{SSD}$,
$\rho_{GND}$ and $\bar{\chi}$ are shown in Figure 9. The material was strained
along the vertical direction in all cases. All deformation contours shown in
this section are from simulations using $k_{KH}=1.2$.

The effective strain contours show clear development of strain gradients at
the interfaces between grains. For example, see regions marked by the white
arrows in the strain and stress contour for grain size $40\mu m$ in Figure 8.
The effect of strain gradients can be observed in the development of GNDs and
backstress in the same regions highlighted in Figure 9. Due to the increase in
GND density and backstress in the vicinity of grain interfaces, higher plastic
deformation occurs and the SSD density increases in these regions as well.
Generally speaking, the effective strain is higher in grains with a lower
effective stress and vice-versa. For example, see grains marked 1 and 2 in
Figure 8. As mentioned earlier, our model does not consider crystallographic
orientation-dependent deformation and hardening. Rather, the anisotropy
factor, $M$, is used to introduce intergranular heterogeneity in our
simulations. While a crystal plasticity model may be more adept at capturing
these near boundary gradient zones [53, 67], our $J_{2}$ plasticity model is
still able to qualitatively simulate heterogeneous deformation near these
grain interfaces that may give rise to grain size-dependent hardening.

Figure 8: Grain maps and contours of effective strain, $\bar{\epsilon}$, and
effective stress, $\bar{\sigma}$, after 0.04 applied strain for simulations
with different grain sizes. Note that the scales are different for the
variables for each grain size. Figure 9: Contours of SSD density,
$\rho_{SSD}$, GND density, $\rho_{GND}$ and effective backstress,
$\bar{\chi}$, after 0.04 applied strain for simulations with different grain
sizes. Note that the scales are different for the different variables for each
grain size.

It can be clearly seen from Figures 8 and 9 that the density of strain
localizations, SSD density, GND density and backstress increase with decrease
in grain size. For example, the strain contours in Figure 8 show precursors of
deformation bands for grain sizes of $8\mu m$ and $4\mu m$ (see regions
highlighted with ellipses). At higher applied strain, these regions may show
significant strain localization. The GND density and backstress increase with
decreasing grain size, as the relative fraction of grain interfaces to grain
interiors increases. For example, it can be seen from Figure 9 that
$\bar{\chi}$ is 5-10 times higher near the grain interfaces as compared to the
grain interior for the simulation with grain size $40\mu m$. Further, it can
also be observed that the density of GNDs is higher than that of the SSDs by
up to an order of magnitude near the grain interfaces, even for the
simulations with high grain sizes.

Overall, it can be concluded from the contours presented in this section that
our model is able to qualitatively capture deformation traits due to the
intergranular heterogeneity via the introduction of the anisotropy factor,
$M$, in our $J_{2}$ plasticity simulations. Further, the observed trends in
the development of $\bar{\epsilon}$, $\bar{\sigma}$, $\rho_{SSD}$,
$\rho_{GND}$ and $\bar{\chi}$ for simulations with different grain sizes show
clear traits of grain size-dependent deformation.

####  4.2.3 Grain Size-Dependent Scaling Relations

The Hall-Petch equation generally used to represent intrinsic grain size-
dependent deformation is of the form [3]:

$\sigma=\sigma_{0}+\sigma_{size-dependent}=\sigma_{0}+\frac{k_{HP}}{\sqrt{D}}$
(30)

where, $\sigma$ is the flow stress, $\sigma_{0}$ is the strength contribution
due to grain size-independent mechanisms, $\sigma_{size-dependent}$ is the
grain size-dependent strength contribution, $k_{HP}$ is the Hall-Petch
coefficient and $D$ is the mean grain size. $k_{HP}$ is a constant that varies
for metals with different crystal structures [66]. Based on analysis of
relevant experimental data, Cordero et al. [66] have also shown that $k_{HP}$
is a function of the applied strain, shear modulus and Burgers vector
magnitude, i.e.,

$k_{HP}\propto G\sqrt{b\epsilon}$ (31)

where, $\epsilon$ is the macroscopic strain. This is based on Asbhy’s model
for hardening due to GNDs in polygrain ensembles [2], where it was proposed
that $\rho_{GND}\propto\epsilon/bD$. Note that variants of this scaling
relation has also been proposed in [68, 32, 33]. We try to establish these
scaling relations from our model predictions here.

In our model, the strength contribution due to GNDs may be attributed solely
to the backstress term. By inverting the flow rule given in Equation (5), we
arrive at the following equation for the flow stress:

$\bar{\sigma}^{*}=M(\tau_{0}+k_{IH}Gb\sqrt{\rho_{SSD}})+S_{t}\left(1-\left(\frac{kT}{\Delta
F_{g}}log\left(\frac{\dot{\bar{\epsilon}}^{p}_{0}}{\dot{\bar{\epsilon}}^{p}}\right)\right)^{1/q}\right)^{1/p}$
(32)

For the small strains and uniaxial loading conditions considered in our
simulations, we further approximate Equation (32) to a 1D form as

$\bar{\sigma}\approx\bar{\sigma}^{*}+\bar{\chi}\approx
M(\tau_{0}+k_{IH}Gb\sqrt{\rho_{SSD}})+S_{t}\left(1-\left(\frac{kT}{\Delta
F_{g}}log\left(\frac{\dot{\bar{\epsilon}}^{p}_{0}}{\dot{\bar{\epsilon}}^{p}}\right)\right)^{1/q}\right)^{1/p}+k_{KH}Gb\sqrt{\rho_{GND}}$
(33)

Here, the first two terms on the RHS represent the size-independent strength
contributions and are similar to the $\sigma_{0}$ term in Equation (30), while
the third term on the RHS is representative of the size-dependent strength
contribution due to GNDs, i.e.,

$\sigma_{size-dependent}=k_{KH}Gb\sqrt{\rho_{GND}}$ (34)

We first verify that the average backstress, $\bar{\chi}$, obtained from our
simulations is comparable to the value of $\sigma_{size-dependent}$ obtained
using the average GND density (Equation (34)). Further, we present results
using two different values of $k_{KH}$, i.e., 1.2 and 0.8. The same initial
microstructures, albeit with random instantiations of $M$, were used. These
values after 0.04 applied strain are presented in Table 2. While there are
some small differences between $\sigma_{size-dependent}$ and $\bar{\chi}$,
these values are comparable for all cases. $\bar{\chi}$ is higher than
$\sigma_{size-dependent}$ by $\approx 6-8\%$ for the lower grain sizes, where
stress concentrations may have developed due to the high values of backstress
near the triple junctions or grain interfaces (cf. Figure 9). Under such
multi-axial stresses, the 1D approximation used in Equation (33) may no longer
be appropriate. Nonetheless, this analysis provides verification that the
strength contribution due to the backstress term indeed scales as the square
root of the GND density in our model.

Table 2: Calculation of $\sigma_{size-dependent}$ using $\rho_{GND}$ and comparison with simulated values of $\bar{\chi}$ after 0.04 applied strain. $k_{KH}$ | Grain size ($\mu m$) |  $\rho_{GND}$ ($mm^{-2}$) |  $\sigma_{size-dependent}$ (MPa) |  $\bar{\chi}$ (MPa)
---|---|---|---|---
1.2 | 3.93 | $9.50\times 10^{6}$ | 45.03 | 48.81
| 7.93 | $6.31\times 10^{6}$ | 36.69 | 38.42
| 20.02 | $3.53\times 10^{6}$ | 27.46 | 27.37
| 30.21 | $2.53\times 10^{6}$ | 23.26 | 23.05
| 39.94 | $1.81\times 10^{6}$ | 19.67 | 18.65
0.8 | 3.93 | $8.68\times 10^{6}$ | 28.70 | 30.00
| 7.93 | $5.33\times 10^{6}$ | 22.50 | 22.84
| 20.02 | $2.82\times 10^{6}$ | 16.37 | 15.32
| 30.21 | $1.92\times 10^{6}$ | 13.50 | 12.85
| 39.94 | $1.38\times 10^{6}$ | 11.43 | 10.17

Finally, we estimate the grain size-dependent scaling relations obtained from
our simulations. For this analysis, we have fit the effective backstress,
$\bar{\chi}$, to the mean grain size, $D$, using the Hall-Petch equation
(Equation (30)) at different values of applied strain. The fit to the
simulated $\bar{\chi}$ versus $D$ data and the corresponding Hall-Petch
coefficients at different applied strains are shown in Figure 10(a) for
$k_{KH}=1.2$ and in Figure 10(b) for $k_{KH}=0.8$. Based on the obtained fit,
$k_{HP}$ ranges from $50MPa\sqrt{\mu m}$ to $108MPa\sqrt{\mu m}$ for the
different applied strains using $k_{KH}=1.2$. Detailed analysis of
experimental data for various fcc metals [66] has shown that the Hall-Petch
coefficient ranges between $70-110MPa\sqrt{\mu m}$, except for Ni, which has
$k_{HP}\approx 230MPa\sqrt{\mu m}$. We have chosen $k_{KH}$ in our model such
that the obtained Hall-Petch coefficient is in the range of experimentally
observed values. Further, it should be noted that this model parameter
provides almost a linear scaling of the strength contribution, manifested in
terms of the backstress. For example, by reducing $k_{KH}$ from 1.2 to 0.8,
the Hall-Petch coefficients at the corresponding strains scale by a factor of
$\approx 2/3$. This can be seen from the range of Hall-Petch coefficients that
lie between $32-66MPa\sqrt{\mu m}$, obtained using $k_{KH}=0.8$ in Figure
10(b).

Figure 10: Plot of $\bar{\chi}$ versus $D$, along with the fit to the Hall-
Petch equation at different applied strains for (a) $k_{KH}=1.2$, and (b)
$k_{KH}=0.8$. The legend indicates the applied strain as well as the best fit
Hall-Petch equation for each set of data.

It is worth noting that the value of the $k_{KH}$ chosen to represent
hardening due to GNDs is much higher than the corresponding value of $k_{IH}$
used to represent hardening due to SSDs. Note that there is an additional
anisotropy factor, $M$, present in the hardening contribution due to SSDs.
Based on the random instantiations, $M$ was generally found to be $\approx
2.5-2.6$ in our simulations. On dividing $k_{KH}$ by $M$, we get
$k_{KH}/M\approx 0.45-0.48$. This lies at the upper bound of strength
coefficient values generally used to represent Taylor hardening due to
dislocations [69]. It should also be noted that this parameter is expected to
be mesh-dependent and a smaller value of $k_{KH}$ may be needed to obtain the
same strength contribution using a finer mesh, where a higher GND density may
be expected [43].

As mentioned earlier, Cordero et al. [66] have shown that $k_{HP}\propto
G\sqrt{b\epsilon}$, by an extension of Ashby’s model for GNDs in polygrain
ensembles [2]. We have plotted the fitted values of $k_{HP}/G\sqrt{b}$
obtained from our simulations as a function of $\sqrt{\bar{\epsilon}}$ in
Figure 11 to verify this. The best fit straight line going through these
points has also been plotted for $k_{KH}=1.2$ and $k_{KH}=0.8$. We have also
plotted the available experimental data for Cu [70, 66] using green symbols
and that for Al [71, 66] using blue symbols in the same figure. It can be seen
that the data points follow a linear trend for a given value of $k_{KH}$.
Thus, our model predictions are shown to agree with the Ashby model of
strengthening due to GNDs [2, 66]. Further, while our model parameters were
not calibrated to any particular material, the values of $k_{HP}/G\sqrt{b}$,
as a function of $\sqrt{\bar{\epsilon}}$, predicted by our model lie in the
same range as the experimental data for Cu and Al.

In summary, the following grain size-dependent scaling relation proposed in
the literature has also been established from our model predictions:

$\sigma_{size-dependent}\propto G\sqrt{\frac{b\bar{\epsilon}}{D}}$ (35)

Figure 11: Plot of $k_{HP}/G\sqrt{b}$ versus $\sqrt{\bar{\epsilon}}$ as
obtained from Figure 10(a) and (b), along with the best fit straight line
passing through these points and comparison with experimental data for Cu and
Al. The experimental data points for Cu are taken from [70, 66] and the
experimental data points for Al are taken from [71, 66].

###  4.3 Discussion

Having presented the model and its predictions, we emphasize on the key
differences between our modeling framework and a conventional $J_{2}$
plasticity model here. Our framework accounts for both isotropic and kinematic
hardening, which may contribute to the expansion and shifting of the yield
surface, respectively. However, the isotropic hardening itself varies from
grain-to-grain and is a function of the anisotropy factor, $M$, thus leading
to individual and distinct yield in each of the constituent grains. Moreover,
we have introduced dislocation density-based constitutive models for modeling
the evolution of substructure during plastic deformation. Kinematic hardening,
which contributes to the shifting of the yield surface, is assumed to be
dominated by the strengthening due to GNDs. As has been shown in our results,
this is able to predict both directional hardening and also grain size-
dependent strengthening. While the former can also be predicted by a
conventional $J_{2}$ plasticity framework, the latter cannot. Moreover, the
deformation-induced microstructure evolution predicted by our model is
generally expected from crystal plasticity simulations. Further, it should be
noted that intrinsic grain size-dependent strengthening has generally been
simulated using crystal plasticity frameworks in the literature [32, 34, 33].
Our $J_{2}$ plasticity simulations are expected to be significantly cheaper,
in terms of computational costs, as compared to crystal plasticity
simulations.

Another discussion point worth mentioning here is that full field Crystal
Plasticity Finite Element (CPFE) models are generally more computationally
expensive than Fast Fourier Transform (FFT)-based spectral crystal plasticity
models [72]. This is due to the large number of degrees of freedom involved
and may limit the consideration of large domains for prediction of grain size
effects in crystal plasticity finite element simulations. In fact, some of the
above mentioned crystal plasticity studies [33, 34] are indeed FFT-based. Our
model may provide an alternative for prediction of grain size strengthening
using finite element simulations, which perhaps offer more convenience for
consideration of realistic geometries and imposing various realistic boundary
conditions (as compared to FFT simulations).

Prediction of anisotropic hardening and grain size-dependent strengthening
using this computationally efficient framework (as compared to a crystal
plasticity framework) is one of the main contributions of our work. In future
work, the anisotropy factor, $M$, may be allowed to evolve with deformation
[48] in order to phenomenologically represent grain rotations during
polycrystalline deformation more accurately. Orientation-dependent elasticity
could also be introduced to further represent anisotropic deformation in this
framework.

##  5 Conclusions

We have proposed a dislocation density-based strain gradient $J_{2}$
plasticity model for simulating the grain size-dependent hardening response of
metallic systems. In this framework, a lower order, Taylor hardening
backstress model is used to represent the hardening due to GNDs. Further, an
anisotropy factor, $M$, is proposed to phenomenologically represent the
differential hardening between grains. An implicit numerical algorithm has
been implemented for the time integration of the finite deformation plasticity
model.

  * •

Our model predicted increased directional hardening during cyclic loading due
to the GND-induced backstress, with higher values of the kinematic hardening
parameter. The contribution of the backstress to strain heterogeneities also
increased with higher values of this parameter.

  * •

Deformation contours also showed the development of strain gradients, GND
density and backstress in the vicinity of grain interfaces. These effects
intensified as the grain size was decreased in our simulations and are due to
the introduction of the anisotropy factor, $M$, in our model.

  * •

Aggregate properties predicted by our simulations were used to determine the
strength scaling relations. A Hall-Petch type, grain size-dependent strength
contribution [3] was predicted by our model.

  * •

Our model predictions also agree with Ashby’s model of strengthening due to
GNDs in polygrain ensembles [2, 66]. Model predictions of $k_{HP}/G\sqrt{b}$,
as a function of $\sqrt{\bar{\epsilon}}$, were found to lie in the same range
as that for fcc metals from available experimental data in the literature.

  * •

Based on the analysis of grain size-dependent strengthening contributions, a
rationale is also provided for choosing appropriate values of the kinematic
hardening parameter to capture the Hall-Petch effect.

## Acknowledgments

The authors gratefully acknowledge funding received from the Department of
Science and Technology (DST) - Science and Engineering Research Board (SERB)
for this research under grant number: CRG/2020/000593.

## References

  * [1] John F Nye.  Some geometrical relations in dislocated crystals.  Acta metallurgica, 1(2):153–162, 1953.
  * [2] MF Ashby.  The deformation of plastically non-homogeneous materials.  The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics, 21(170):399–424, 1970.
  * [3] EO Hall.  The deformation and ageing of mild steel: Iii discussion of results.  Proceedings of the Physical Society. Section B, 64(9):747, 1951\.
  * [4] NJ Petch.  The cleavage strength of polycrystals.  Journal of the Iron and Steel Institute, 174:25–28, 1953.
  * [5] NA Fleck, GM Muller, Mike F Ashby, and John W Hutchinson.  Strain gradient plasticity: theory and experiment.  Acta Metallurgica et materialia, 42(2):475–487, 1994.
  * [6] JS Stölken and AG Evans.  A microbend test method for measuring the plasticity length scale.  Acta Materialia, 46(14):5109–5115, 1998.
  * [7] Michael D Uchic, Dennis M Dimiduk, Jeffrey N Florando, and William D Nix.  Sample dimensions influence strength and crystal plasticity.  Science, 305(5686):986–989, 2004.
  * [8] AA Elmustafa and DS Stone.  Nanoindentation and the indentation size effect: Kinetics of deformation and strain gradient plasticity.  Journal of the Mechanics and Physics of Solids, 51(2):357–381, 2003\.
  * [9] Minhua Zhao, William S Slaughter, Ming Li, and Scott X Mao.  Material-length-scale-controlled nanoindentation size effects due to strain-gradient plasticity.  Acta Materialia, 51(15):4461–4469, 2003.
  * [10] Yang Mu, Ke Chen, and WJ Meng.  Thickness dependence of flow stress of cu thin films in confined shear plastic flow.  MRS Communications, 4(3):129–133, 2014.
  * [11] Yang Mu, Xiaoman Zhang, John W Hutchinson, and Wen Jin Meng.  Measuring critical stress for shear failure of interfacial regions in coating/interlayer/substrate systems through a micro-pillar testing protocol.  Journal of Materials Research, 32(8):1421–1431, 2017.
  * [12] Mark GD Geers, WAM Brekelmans, and PJM Janssen.  Size effects in miniaturized polycrystalline fcc samples: Strengthening versus weakening.  International Journal of Solids and Structures, 43(24):7304–7321, 2006.
  * [13] Julia R Greer and Jeff Th M De Hosson.  Plasticity in small-sized metallic systems: Intrinsic versus extrinsic size effect.  Progress in Materials Science, 56(6):654–724, 2011.
  * [14] A Arsenlis and DM Parks.  Crystallographic aspects of geometrically-necessary and statistically-stored dislocation density.  Acta materialia, 47(5):1597–1611, 1999.
  * [15] H Gao, Y Huang, WD Nix, and JW Hutchinson.  Mechanism-based strain gradient plasticity—i. theory.  Journal of the Mechanics and Physics of Solids, 47(6):1239–1263, 1999.
  * [16] A Acharya and JL17611250963 Bassani.  Lattice incompatibility and a gradient theory of crystal plasticity.  Journal of the Mechanics and Physics of Solids, 48(8):1565–1595, 2000.
  * [17] Morton E Gurtin.  A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations.  Journal of the Mechanics and Physics of Solids, 50(1):5–32, 2002\.
  * [18] Amit Acharya and Anish Roy.  Size effects and idealized dislocation microstructure at small scales: predictions of a phenomenological model of mesoscopic field dislocation mechanics: Part i.  Journal of the Mechanics and Physics of Solids, 54(8):1687–1710, 2006.
  * [19] MGD Geers, WAM Brekelmans, and CJ Bayley.  Second-order crystal plasticity: internal stress effects and cyclic loading.  Modelling and Simulation in Materials Science and Engineering, 15(1):S133, 2006.
  * [20] Chung-Souk Han, Anxin Ma, Franz Roters, and Dierk Raabe.  A finite element approach with patch projection for strain gradient plasticity formulations.  International Journal of Plasticity, 23(4):690–710, 2007.
  * [21] Jeffrey W Kysar, Yong X Gan, Timothy L Morse, Xi Chen, and Milton E Jones.  High strain gradient plasticity associated with wedge indentation into face-centered cubic single crystals: geometrically necessary dislocation densities.  Journal of the Mechanics and Physics of Solids, 55(7):1554–1573, 2007.
  * [22] PJ Guruprasad and AA Benzerga.  A phenomenological model of size-dependent hardening in crystal plasticity.  Philosophical Magazine, 88(30-32):3585–3601, 2008.
  * [23] SP Lele and L Anand.  A small-deformation strain-gradient theory for isotropic viscoplastic materials.  Philosophical Magazine, 88(30-32):3655–3689, 2008.
  * [24] Jason R Mayeur, David L McDowell, and Douglas J Bammann.  Dislocation-based micropolar single crystal plasticity: Comparison of multi-and single criterion theories.  Journal of the Mechanics and Physics of Solids, 59(2):398–422, 2011\.
  * [25] FPE Dunne, R Kiwanuka, and AJ Wilkinson.  Crystal plasticity analysis of micro-deformation, lattice rotation and geometrically necessary dislocation density.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 468(2145):2509–2531, 2012.
  * [26] Matthieu Mazière and Samuel Forest.  Strain gradient plasticity modeling and finite element simulation of lüders band formation and propagation.  Continuum Mechanics and Thermodynamics, 27(1):83–104, 2015.
  * [27] Emilio Martínez-Pañeda and Christian F Niordson.  On fracture in finite strain gradient plasticity.  International Journal of Plasticity, 80:154–167, 2016.
  * [28] Murat Demiral, Kai Nowag, Anish Roy, Rudy Ghisleni, Johann Michler, and Vadim V Silberschmidt.  Enhanced gradient crystal-plasticity study of size effects in a $\beta$-titanium alloy.  Modelling and Simulation in Materials Science and Engineering, 25(3):035013, 2017.
  * [29] Peter Gudmundson and Carl FO Dahlberg.  Isotropic strain gradient plasticity model based on self-energies of dislocations and the taylor model for plastic dissipation.  International Journal of Plasticity, 121:1–20, 2019.
  * [30] Mitsutoshi Kuroda and Alan Needleman.  A simple model for size effects in constrained shear.  Extreme Mechanics Letters, 33:100581, 2019.
  * [31] G Zhou, W Jeong, ER Homer, DT Fullwood, MG Lee, JH Kim, H Lim, H Zbib, and RH Wagoner.  A predictive strain-gradient model with no undetermined constants or length scales.  Journal of the Mechanics and Physics of Solids, 145:104178, 2020\.
  * [32] Fengwei Sun, Edward D Meade, and Noel P O’Dowd.  Strain gradient crystal plasticity modelling of size effects in a hierarchical martensitic steel using the voronoi tessellation method.  International Journal of Plasticity, 119:215–229, 2019.
  * [33] S Haouala, S Lucarini, J LLorca, and J Segurado.  Simulation of the hall-petch effect in fcc polycrystals by means of strain gradient crystal plasticity and fft homogenization.  Journal of the Mechanics and Physics of Solids, 134:103755, 2020\.
  * [34] Stéphane Berbenni, Vincent Taupin, and Ricardo A Lebensohn.  A fast fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals.  Journal of the Mechanics and Physics of Solids, 135:103808, 2020\.
  * [35] JR Mayeur and DL McDowell.  A comparison of gurtin type and micropolar theories of generalized single crystal plasticity.  International Journal of Plasticity, 57:29–51, 2014.
  * [36] George Z Voyiadjis and Yooseob Song.  Strain gradient continuum plasticity theories: theoretical, numerical and experimental investigations.  International Journal of Plasticity, 121:21–75, 2019.
  * [37] LP Evers, WAM Brekelmans, and MGD Geers.  Non-local crystal plasticity model with intrinsic ssd and gnd effects.  Journal of the Mechanics and Physics of Solids, 52(10):2379–2401, 2004.
  * [38] S Yefimov and E Van der Giessen.  Size effects in single crystal thin films: nonlocal crystal plasticity simulations.  European Journal of Mechanics-A/Solids, 24(2):183–193, 2005.
  * [39] CJ Bayley, WAM Brekelmans, and MGD Geers.  A comparison of dislocation induced back stress formulations in strain gradient crystal plasticity.  International Journal of Solids and Structures, 43(24):7268–7286, 2006.
  * [40] Mitsutoshi Kuroda and Viggo Tvergaard.  Studies of scale dependent crystal viscoplasticity models.  Journal of the Mechanics and Physics of Solids, 54(9):1789–1810, 2006.
  * [41] Samuel Forest.  Some links between cosserat, strain gradient crystal plasticity and the statistical theory of dislocations.  Philosophical Magazine, 88(30-32):3549–3563, 2008.
  * [42] Kartik Kapoor, Yung Suk Jeremy Yoo, Todd A Book, Josh P Kacher, and Michael D Sangid.  Incorporating grain-level residual stresses and validating a crystal plasticity model of a two-phase ti-6al-4 v alloy produced via additive manufacturing.  Journal of the Mechanics and Physics of Solids, 121:447–462, 2018\.
  * [43] Ritwik Bandyopadhyay, Sven E Gustafson, Kartik Kapoor, Diwakar Naragani, Darren C Pagan, and Michael D Sangid.  Comparative assessment of backstress models using high-energy x-ray diffraction microscopy experiments and crystal plasticity finite element simulations.  International Journal of Plasticity, 136:102887, 2021.
  * [44] Geoffrey Ingram Taylor.  The mechanism of plastic deformation of crystals. part i.—theoretical.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 145(855):362–387, 1934.
  * [45] Akhtar S Khan and Sujian Huang.  Continuum theory of plasticity.  John Wiley & Sons, 1995.
  * [46] Gustavo Weber and Lallit Anand.  Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids.  Computer Methods in Applied Mechanics and Engineering, 79(2):173–202, 1990.
  * [47] Ulrich Fred Kocks, A Argon, and M Ashby.  Thermodynamics and kinetics of slip.  1975\.
  * [48] Sumeet Mishra, Manasij Yadava, Kaustubh N Kulkarni, and NP Gurao.  A new phenomenological approach for modeling strain hardening behavior of face centered cubic materials.  Acta Materialia, 178:99–113, 2019.
  * [49] Yuri Estrin.  Dislocation-density-related constitutive modeling.  Unified constitutive laws of plastic deformation, 1:69–106, 1996\.
  * [50] UF Kocks and H Mecking.  Physics and phenomenology of strain hardening: the fcc case.  Progress in materials science, 48(3):171–273, 2003.
  * [51] Hong Dai.  Geometrically-necessary dislocation density in continuum plasticity theory, FEM implementation and applications.  PhD thesis, Massachusetts Institute of Technology, 1997.
  * [52] William D Nix and Huajian Gao.  Indentation size effects in crystalline materials: a law for strain gradient plasticity.  Journal of the Mechanics and Physics of Solids, 46(3):411–425, 1998\.
  * [53] LP Evers, WAM Brekelmans, and MGD Geers.  Scale dependent crystal plasticity framework with dislocation density and grain boundary effects.  International Journal of solids and structures, 41(18-19):5209–5230, 2004.
  * [54] Alberto M Cuitino and Michael Ortiz.  Computational modelling of single crystals.  Modelling and Simulation in Materials Science and Engineering, 1(3):225, 1993.
  * [55] Robert Davis McGinty.  Multiscale representation of polycrystalline inelasticity.  PhD thesis, Georgia institute of technology, 2001.
  * [56] RD McGinty and DL McDowell.  A semi-implicit integration scheme for rate independent finite crystal plasticity.  International Journal of Plasticity, 22(6):996–1025, 2006.
  * [57] Xianwu Ling, MF Horstemeyer, and GP Potirniche.  On the numerical implementation of 3d rate-dependent single crystal plasticity formulations.  International Journal for Numerical Methods in Engineering, 63(4):548–568, 2005.
  * [58] Morton L Curtis.  Algebraic preliminaries.  In Abstract Linear Algebra, pages 1–7. Springer, 1990.
  * [59] Kaare Brandt Petersen, Michael Syskind Pedersen, et al.  The matrix cookbook.  Technical University of Denmark, 7(15):510, 2008.
  * [60] Cody J Permann, Derek R Gaston, David Andrš, Robert W Carlsen, Fande Kong, Alexander D Lindsay, Jason M Miller, John W Peterson, Andrew E Slaughter, Roy H Stogner, et al.  Moose: Enabling massively parallel multiphysics simulation.  SoftwareX, 11:100430, 2020.
  * [61] M Kothari and L Anand.  Elasto-viscoplastic constitutive equations for polycrystalline metals: application to tantalum.  Journal of the Mechanics and Physics of Solids, 46(1):51–83, 1998\.
  * [62] Devraj Ranjan, Sankar Narayanan, Kai Kadau, and Anirban Patra.  Crystal plasticity modeling of non-schmid yield behavior: from ni3al single crystals to ni-based superalloys.  Modelling and Simulation in Materials Science and Engineering, 29(5):055005, 2021.
  * [63] Karo Sedighiani, Konstantina Traka, Franz Roters, Dierk Raabe, Jilt Sietsma, and Martin Diehl.  Determination and analysis of the constitutive parameters of temperature-dependent dislocation-density-based crystal plasticity models.  Mechanics of Materials, page 104117, 2021.
  * [64] Romain Quey, PR Dawson, and Fabrice Barbe.  Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing.  Computer Methods in Applied Mechanics and Engineering, 200(17-20):1729–1745, 2011.
  * [65] Csimsoft Trelis.  Advanced meshing for challenging simulations (2013-2016).
  * [66] Zachary C Cordero, Braden E Knight, and Christopher A Schuh.  Six decades of the hall–petch effect–a survey of grain-size strengthening studies on pure metals.  International Materials Reviews, 61(8):495–512, 2016.
  * [67] SK Mishra, P Pant, K Narasimhan, AD Rollett, and I Samajdar.  On the widths of orientation gradient zones adjacent to grain boundaries.  Scripta materialia, 61(3):273–276, 2009.
  * [68] Jaafar A El-Awady.  Unravelling the physics of size-dependent dislocation-mediated plasticity.  Nature communications, 6(1):1–9, 2015.
  * [69] Ali Argon.  Strengthening mechanisms in crystal plasticity, volume 4.  Oxford University Press on Demand, 2008.
  * [70] Noboru Ono and Seiichi Karashima.  Grain size dependence of flow stress in copper polycrystals.  Scripta Metallurgica, 16(4):381–384, 1982.
  * [71] Niels Hansen.  The effect of grain size and strain on the tensile flow stress of aluminium at room temperature.  Acta Metallurgica, 25(8):863–869, 1977.
  * [72] Philip Eisenlohr, Martin Diehl, Ricardo A Lebensohn, and Franz Roters.  A spectral method solution to crystal elasto-viscoplasticity at finite strains.  International Journal of Plasticity, 46:37–53, 2013.