A Multiscale Gradient Theory for Single 

Crystalline Elastoviscoplasticity 

by John D. Clayton, David L. McDowell, and Douglas J. Bammann 

ARL-RP-117 February 2006 

A reprint from the International Journal of Engineering Science, 

vol. 42, pp. 427–457, 2004.

Approved for public release; distribution is unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory 

Aberdeen Proving Ground, MD 21005-5066 

ARL-RP-117 February 2006 

A Multiscale Gradient Theory for Single 

Crystalline Elastoviscoplasticity 

John D. Clayton 

Weapons and Materials Research Directorate, ARL 

David L. McDowell

Georgia Institute of Technology 

Douglas J. Bammann

Sandia National Laboratories 

A reprint from the International Journal of Engineering Science, 

vol. 42, pp. 427–457, 2004. 

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and 

maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

February 2006 

2. REPORT TYPE 

Reprint 

3. DATES COVERED (From - To) 

1 January 2003–1 January 2006 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

A Multiscale Gradient Theory for Single Crystalline Elastoviscoplasticity 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

1L1622618AH80 

5e. TASK NUMBER 

6. AUTHOR(S) 

John D. Clayton, David L. McDowell,* and Douglas J. Bammann†

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

ATTN: AMSRD-ARL-WM-TD 

Aberdeen Proving Ground, MD 21005-5066 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-RP-117 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

13. SUPPLEMENTARY NOTES 

A reprint from the International Journal of Engineering Science, vol. 42, pp. 427–457, 2004. *

George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0405. 

†

Department of Science-based Materials Modeling, Sandia National Laboratories, Livermore, CA 94551. 

14. ABSTRACT 

Explicit volume averaging procedures are used to motivate a gradient-type description of single crystalline elastoviscoplasticity. 

Upon regarding local elastic and plastic deformation gradients within the crystal as continuously differentiable fields, we arrive 

at a three-term multiplicative decomposition for the volume-averaged deformation gradient, consisting of a recoverable elastic 

term associated with the average applied stress and average lattice rotation, an inelastic term associated with the average plastic 

velocity gradient, and a (new) third term reflecting the presence of the residual microelastic deformation gradient within the 

volume and providing a representation of the kinematics of grain subdivision via formation of low-angle subgrain boundaries, 

for example. A variant of the classical Eshelby stress tensor provides the driving force for homogenized viscoplastic flow, with 

slip resistances dictated by densities of geometrically necessary and statistically stored dislocations. Distinctive features of the 

continuum model include coupling of internal elastic strain energy densities associated with residual and applied stresses, 

dependency of the single crystalline effective elastic moduli upon evolution of lattice substructure, and a characteristic length 

potentially based upon both the size of the crystal element used in volume averaging and the grain subdivision measure. 

15. SUBJECT TERMS

crystal plasticity, dislocations, disclinations, lattice defects, multiscale models 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

John D. Clayton 

a. REPORT 

UNCLASSIFIED 

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF ABSTRACT 

UL 

18. NUMBER 

OF PAGES 

38 19b. TELEPHONE NUMBER (Include area code) 

(410) 306-0975 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
A multiscale gradient theory for single crystalline

elastoviscoplasticity

J.D. Clayton a, D.L. McDowell b,*, D.J. Bammann c

a Impact Physics Branch, US Army Research Laboratory, Aberdeen Proving Ground, MD 21005, USA

b George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta,

GA 30332-0405, USA c Department of Science-based Materials Modeling, Sandia National Laboratories, Livermore, CA 94551, USA

Received 25 August 2003; accepted 25 August 2003

(Communicated by JOHN W. HUTCHINSON)

Abstract

Explicit volume averaging procedures are used to motivate a gradient-type description of single crystalline elastoviscoplasticity. Upon regarding local elastic and plastic deformation gradients within the

crystal as continuously differentiable fields, we arrive at a three-term multiplicative decomposition for the

volume-averaged deformation gradient, consisting of a recoverable elastic term associated with the average

applied stress and average lattice rotation, an inelastic term associated with the average plastic velocity

gradient, and a (new) third term reflecting the presence of the residual microelastic deformation gradient

within the volume and providing a representation of the kinematics of grain subdivision via formation of

low-angle subgrain boundaries, for example. A variant of the classical Eshelby stress tensor provides the

driving force for homogenized viscoplastic flow, with slip resistances dictated by densities of geometrically

necessary and statistically stored dislocations. Distinctive features of the continuum model include coupling

of internal elastic strain energy densities associated with residual and applied stresses, dependency of the

single crystalline effective elastic moduli upon evolution of lattice substructure, and a characteristic length

potentially based upon both the size of the crystal element used in volume averaging and the grain subdivision measure.

 2003 Elsevier Ltd. All rights reserved.

*

Corresponding author. Tel.: +1-404-8945128; fax: +1-404-8940186.

E-mail address: david.mcdowell@me.gatech.edu (D.L. McDowell).

0020-7225/$ - see front matter  2003 Elsevier Ltd. All rights reserved.

doi:10.1016/j.ijengsci.2003.08.001

International Journal of Engineering Science 42 (2004) 427–457

www.elsevier.com/locate/ijengsci
1. Introduction

A rather large number of single crystal and polycrystal plasticity models have been proposed

incorporating higher than first-order gradients (e.g., strain gradients, lattice curvature, gradientbased dislocation density measures) of deformation in the material response [1–6]. The higher

order gradients render these classes of models non-local, effectively injecting a length scale (e.g.,

normalization constant) into the formulation on dimensional grounds, the value of which is

presumably associated with a characteristic dimension of the microstructure, such as the evolving

dislocation cell size in ductile FCC crystals, for example (cf. [7]).

Higher order deformation gradients have been included within plasticity theories for a variety

of different reasons; we elaborate briefly on a few of these here. Their inclusion has permitted

resolution of numerical difficulties associated with the solution of boundary value problems of

strain softening materials in which strain localization occurs (cf. [8,9]). Gradient based approaches

have also been used to model dislocation dynamics and pattern formation [10–12] and to describe

single and periodic shear bands in single crystals (cf. [13]) and polycrystals [12]. Other recent

applications of non-local theories include characterization of stress and strain fields, without

singularities in the field variables, at dislocation cores and crack tips [14,15] and modeling the

evolution of the plastic spin tensor in macroscopic finite deformation plasticity theory [16,17].

Perhaps the most recurrently reported motivation for use of gradient theories of plasticity has

been representation of the observed trend of increasing strength with decreasing size of considered

volume or microstructural features. Often cited is the Hall–Petch relation, in which hardness

properties (i.e., yield stress and cleavage strength) increase with decreasing grain size in polycrystals (specifically an inverse square-root dependence [18,19]), a phenomenon that classical local

plasticity theory, being devoid of a material length scale, is unable to capture. Shu and Fleck [20]

and Forest et al. [21] used couple stress theories (pioneered by Cosserat and Cosserat [22,23]) to

characterize Hall–Petch behavior in bicrystals and polycrystals, respectively. Fleck et al. [24]

employed a couple stress model of strain gradient plasticity to describe an increase in flow stress

with decreasing diameter of twisted thin copper wires. Shu and Fleck [25] and Hwang et al. [26]

used variations of the same strain gradient-couple stress theory [27] to capture an observed increase in hardness with decreasing indentor size in pure metals.

Many theories have included a higher order gradient of elastic, plastic, or total deformation

in the expression for the yield stress, slip system hardening rate, flow stress, or the backstress,

with details of incorporation of gradients in the material response functions varying widely

among different models. In many cases the higher order gradients are associated, upon invocation of differential-geometric arguments, with the density of ‘‘geometrically necessary dislocations’’ (GNDs) in single crystals (cf. [28–37]) and in polycrystals (cf. [38,16,17,9,39]). These

dislocations are required to sustain the compatibility of the total deformation gradient, in contrast to the ‘‘statistically stored’’ dislocations (SSDs) [40] that accumulate under homogeneous

plastic flow.

Physical experiments have demonstrated how the processes of grain subdivision and dislocation substructure formation substantially influence slip system activity, strain hardening,

stored lattice energy, and texture evolution in single and polycrystals [41–47]. Local models

have been proposed that explicitly embed subdivision and related dislocation substructure

effects into the single crystal kinematics [47–50] and the hardening and grain interaction laws

428 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
of polycrystal plasticity theory [51,52], without explicitly considering the higher order gradients of deformation associated with geometrically necessary defects. Non-local solutions in terms

of Greens functions [14,53] have also been developed to address defect substructure within grains

and its influence on the micro-stress fields and commensurate hardening behavior. Also

noteworthy is the recent approach of Ortiz and co-workers [7,54] for modeling additional

internal degrees of freedom associated with dislocation substructure via a non-local sequential

lamination theory. This approach permits dislocation substructure development, when energetically favorable, to occur within single crystals even under uniform monotonic macroscopic

deformation, in contrast to many gradient-based approaches requiring heterogeneity of the

deformation gradient field to drive evolution of GNDs. Such lamellar dislocation structures

have been observed experimentally within ductile single crystals under monotonic loading at

large deformations [42,43,55] and are thought to influence the flow stress across a wide range

of temperatures and strain rates [44,56]. Carstensen et al. [57] viewed evolution of plastic flow

from the standpoint of a constrained energy minimization problem and remark how heterogeneous plastic flow (e.g., subdivision or localization) may result from a lack of convexity of the free

energy.

Also observed within pure ductile metals and certain alloys at large deformations and/or high

temperatures are long range internal stress fields associated with misoriented subgrain boundaries [58–60]. These internal stress fields––attributed to misorientations that develop among

neighboring subgrains––occur with a periodicity on the order of the subgrain size (cf. [43]). The

study of Gibeling and Nix [59] found that while local internal stress fields in the unloaded

configuration (i.e., residual stresses) are typically quite small compared to the average applied

stress, the applied loads can alter the arrangement of subgrain walls such that the stress fields of

adjacent dislocations do not cancel, thereby biasing the internal stress fields to a substantial

degree.

The remainder of this work is organized as follows. Kinematics and balance laws are discussed

in Sections 2 and 3, respectively, from the perspective of multiscale volume averaging. Section 4

more fully develops constitutive relations within the context of a higher order gradient, single

crystal plasticity theory. In the interest of brevity, thermal effects (i.e., temperature rates and heat

fluxes) and dynamic effects (i.e., acceleration and body forces) are often neglected.

We employ the following notation. Vector and tensor quantities are typically represented with

boldface type, while scalars and individual components of vectors and tensors are written in

italics. The index notation is often used for clarity, following the Einstein summation convention

and distinguishing between covariant (subscript) and contravariant (superscript) components.

Current configuration indices are written in lower case Latin, reference configuration indices in

upper case Latin, and intermediate configuration indices are written using Greek symbols. Juxtaposition implies summation over two repeated adjacent indices (e.g., ðABÞ

b

a 1⁄4 AacBcb). The dot

(scalar) product of vectors is represented by the symbol ‘‘Æ’’ (e.g., a  b 1⁄4 aagabbb, with gab components of the metric tensor). Angled brackets denote a dual (scalar) product (e.g., for secondrank tensors, hA; Bi 1⁄4 trðABÞ 1⁄4 AabBba, and for contra covector pairs, ha; bi 1⁄4 aaba). The colon

denotes contraction over repeated pairs of indices (e.g., A : B 1⁄4 trðATBÞ 1⁄4 AabBab and

C : E 1⁄4 CabcdEcd ). The symbol ‘‘’’ represents the tensor (outer) product (e.g., ða bÞ

ab 1⁄4 aabb).

Superposed )1, T, and ‘‘Æ’’ denote inverse, transpose, and material time derivative operations,

respectively. Additional notation is clarified as needed later in the text.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 429
2. Multiscale kinematics

A fundamental difference between the model developed in the present work and the numerous

strain gradient approaches already cited is our usage of rigorous volume averaging procedures to

characterize the kinematics of elasto(visco)plastic deformation. Continuum elements (i.e., domains for volume integration) representing a deforming crystal or region within are shown in Fig.

1. The symbols vref, ~vint, and vcur denote volume elements in reference, intermediate, and current

configurations. Dimensions of vref are assumed to adhere to

a0  ‘ref 6 LG; ð1Þ

with a0 and LG the lattice parameter and average grain diameter, respectively, and with

‘ref  ffiffiffiffiffiffi

vref p3 . We also frequently invoke the following notation for configurations at multiple length

scales: vref  bref  Bref, ~vint  ~bint  Beint, and vcur  bcur  Bcur, where global reference, intermediate, and current configurations of the entire macroscopic body are labeled Bref, Beint, and Bcur. All

configurations will be described in more detail later in the text. Local, scalar ‘‘differential’’ volume

elements within reference, intermediate, and current configurations are labeled dvref, d~vint, and

dvcur, respectively. Differential elements are required to be smaller than their associated volume

elements, but not necessarily infinitesimal in size.

The motion of points of the body from the reference to current configuration is written

xa 1⁄4 ua

ðXA; tÞ ða; A 1⁄4 1; 2; 3Þ; ð2Þ

where xa, XA, and t denote, respectively, current configuration coordinates, reference configuration coordinates, and time. The local motion ua is assumed to be continuously differentiable

within vcur. Local intermediate configuration coordinates are also assumed to be continuous within

~v. These are denoted by ~xa 1⁄4 ~xaðXA; tÞða 1⁄4 1; 2; 3Þ. Further remarks on the availability of continuous, single-valued coordinates ~xa follow later.

Basis vectors in each configuration are labeled as

Fig. 1. Configurations, differential volumes, and coordinate systems for crystal volume element: (a) reference, (b) intermediate, and (c) current.

430 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
GA 

o

oXA ; ~ga 

o

o~xa ; ga 

o

oxa : ð3Þ

Basis covectors (i.e., dual vectors) are introduced as

GA  dXA; ~ga  d~xa; ga  dxa; ð4Þ

satisfying

hGA

; GBi 1⁄4 dA

B; h~ga

; ~gbi 1⁄4 da

b; hga

; gbi 1⁄4 da

b: ð5Þ

Metric tensors are introduced on each configuration, obeying the relations

GAB 1⁄4 GA  GB; g~ab 1⁄4 ~ga  ~gb; gab 1⁄4 ga  gb: ð6Þ

The notation G  det G, g~  det ~g, and g  det g denotes determinants of metric tensors given in

(6). Since we often use volume averaging operations over the crystal element to define certain

tensorial quantities, we restrict the basis vectors (3), basis covectors (4), and metric tensors (6) to

be constant, but not necessarily Cartesian, within each element (i.e., vref, ~vint, and vcur) in each

configuration, such that covariant and partial differentiation are equivalent operations. However,

these variables are permitted to vary from volume element to volume element (if curvilinear

coordinates are useful), and also from configuration to configuration.

The configurations of the crystalline volume element shown in Fig. 1 are now defined in turn.

The reference configuration volume vref consists of the crystalline lattice as it existed prior to

application of external forces (i.e., at t 1⁄4 0), such that it is free of traction along the external

boundary sref. It may or may not contain dislocations, internal residual elastic lattice strains, or

residual plastic deformation. The current configuration volume element, vcur, is the elastoplastically-deformed crystal material, with possibly non-vanishing traction vector t per unit reference

area applied on external boundary scur.

The local deformation gradient f for material points with local coordinates XA within the

volume element is defined as the tangent mapping

f  TuX 1⁄4 oua

oXA ga GA: ð7Þ

The volume-averaged deformation gradient F for the crystal element is then defined by the motion

of its external boundary, which is equivalent to the volume-averaged local deformation gradient

upon invocation of Gausss theorem [61]:

F a

A 

1

vref Z

sref

xaðdsrefÞA 1⁄4 1

vref Z

sref

oxa

oXA dvref 1⁄4 1

vref Z

sref

f a

A dvref: ð8Þ

In Eq. (8), ðdsrefÞA  ðnrefÞAdsref is an oriented differential surface element of sref (Fig. 1), with nref

a unit covariant vector normal to sref.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 431
An average elastic stretch tensor Ve is associated with the average external stress applied to scur.

The intermediate configuration bint (Fig. 2) of the crystal element reached upon hypothetical

instantaneous elastic unloading from the current configuration via the inverse of the elastic stretch

Ve1 corresponds to null traction conditions on the external boundary of the crystal volume element vint (i.e., the traction tint 1⁄4 0 along sint), as shown in Fig. 2. The left elastic stretch tensor Ve

is determined explicitly from

Ve 1⁄4

Z

sref

x



dsref Z

sref

x



dsref1

; ð9Þ

where x are the local coordinates of the external boundary of the element corresponding to

traction-free intermediate configuration bint. Configuration bint arises from instantaneous removal

of traction along the boundary of vint, constrained in such a way that the global rotation of the

volume element, Re1

, does not occur upon stress relaxation. Upon unloading, plastic rearrangements are idealized as rate independent and inertial effects are neglected.

The total ‘‘elastic’’ rotation tensor Re is determined from the solution of the following integrodifferential equation for the elastic spin and associated initial conditions:

R_ e

Re1 1⁄4 1

vref Z

vref

r_

e

r

e1 dvref; Re

ðt 1⁄4 0Þ 1⁄4 r

e

ðt 1⁄4 0Þ 1⁄4 1; ð10Þ

where 1 is the identity map and re is the local elastic and rigid body rotation exhibited by dvref as

it is deformed to its current representation dvcur. Note that ReT 1⁄4 Re1 follows from (10) and

reT 1⁄4 re1, since volume averaging preserves the anti-symmetric property of the spin.

The intermediate configuration ~bint is defined by the net unloading procedure Fe1 1⁄4 ReTVe1

,

i.e. unloading by the external forces and subsequent rotation by the inverse (i.e., transpose) of the

average lattice rotation. After this unloading procedure, local coordinates ~xa describe the positions

of particles within the volume element ~vint. The local residual deformation gradient ~f is the tangent

mapping

Fig. 2. Configurations for crystal volume element.

432 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
~f  Tu~X 1⁄4 ou~a

oXA ~ga GA

; ð11Þ

with ~xa 1⁄4 u~aðxaðXA; tÞ; tÞ the motion for differential volume elements between bref and ~bint. The

volume-averaged residual deformation gradient Fe for the crystal element is deduced from the

intermediate configuration coordinates of its external boundary, and is equivalent to the volumeaveraged residual local deformation gradient upon invocation of Gausss theorem (compare with

(8))

F

ea

A 

1

vref Z

sref

~xa

ðdsrefÞA 1⁄4 1

vref Z

sref

o~xa

oXA dvref 1⁄4 1

vref Z

sref

~f a

A dvref: ð12Þ

We now make two additional assumptions regarding the kinematics of single crystalline elastoplasticity on a ‘‘pointwise’’ basis, for each differential volume element [62]

f 1⁄4 f

e

f

p

; ~f 1⁄4 ~f

e~fp

: ð13Þ

In Eq. (13)1, f

e 1⁄4 ve

re is the total lattice stretch ðveÞ and rotation ðreÞ for a differential volume

element deformed to the current configuration, and f

p is the remaining (plastic) deformation

gradient attributed to the history of motion of dislocations. In Eq. (13)2, ~f

e embodies the residual

elastic lattice stretch and rotation for a differential volume element contained within the externally

unloaded element ~vint, and ~f

p is the residual plastic deformation gradient arising from dislocation

motion. In order to illustrate (13), we label differential volume elements as dvref  vref, d~vint  ~vint,

and dvcur  vcur in Fig. 3 with differential elements visually enlarged relative to their surrounding

crystal volume elements for clarity. The local deformation gradients of (13) then act as mappings

between tangent spaces of local ‘‘differential configurations’’:

f : T ðdvrefÞ ! T ðdvcurÞ; f

p : T ðdvrefÞ ! T ðdvpÞ; fe : T ðdvpÞ ! T ðdvcurÞ;

~f : T ðdvrefÞ ! T ðd~vintÞ; ~f

p : T ðdvrefÞ ! T ðd~vpÞ; ~fe : T ðd~vpÞ ! T ðd~vintÞ; ð14Þ

Fig. 3. Tangent mappings and volume elements at multiple length scales.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 433
Notice that of the six tangent mappings in (14), only f and ~f are necessarily compatible deformation gradient fields, in agreement with definitions (7) and (11).

The total plastic deformation gradient F

p for the volume element emerges from the following

integro-differential equation written in terms of the plastic velocity gradient:

_

F

p

F

p1



1

vref Z

vref

_

f

p

f

p1 dvref; ð15Þ

with initial conditions F

p

ðt 1⁄4 0Þ 1⁄4 f

p

ðt 1⁄4 0Þ 1⁄4 1 if a perfect reference lattice is assumed. Notice

that if the local plastic flow is isochoric, i.e. trð_f

p

f

p1

Þ 1⁄4 0, then the volume-averaged plastic flow

of (15) is volume-preserving as well, i.e. trð_F

p

F

p1

Þ 1⁄4 0.

Next we assume a three-term multiplicative decomposition for the total deformation F:

F 1⁄4 Fe Fei

F

p

|ffl{zffl}

eF

1⁄4 Fe

Fe; ð16Þ

where the net residual deformation gradient Fe enters the decomposition as shown, in accordance

with our previous definitions of Fe and Fe. Combining Eqs. (8), (12), (13)2, and (16) leads to the

definition of Fei

:

Fei  FeF

p1 1⁄4 1

vref Z

vref

~f

e~fp

dvref F

p1

: ð17Þ

The Fei (two-point) tensor can be thought of as an indicator of residual elasticity in configuration ~bint (and corresponding residual stresses and interaction energies). The deformation tensor Fei

contains both residual elastic ð~f

e

Þ and plastic ð~f

p

; F

p

Þ contributions. It is clear from Eq. (17)2 that

Fei ! 1 as ~f

e ! 1 and ~fp ! F

p

, as would be the case for homogeneous deformation of the entire

crystalline element, such that external unloading by Ve1 relieves all local internal stresses and

plastic rearrangements do not occur upon instantaneous unloading. On the other hand, if the

elastoplastic deformation fields are heterogeneous throughout the volume element, Fei and commensurate residual stresses may not vanish. Such is the case when cellular or laminar dislocation

substructures misoriented from one other by finite lattice rotations (embodied here locally by the

orthogonal part ~re of ~f

e

) evolve and refine during large strain deformations in FCC metals at room

temperatures [43,44]. When the volume element encompasses an entire single crystal (‘ref 1⁄4 LG in

Eq. (1)), we thus conclude that Fei represents the contribution of grain subdivision processes to the

total deformation gradient F for the crystal, as originally proposed by Butler and McDowell [50].

Additionally, if heterogeneity of local deformation and stresses due to elastoplastic incompatibility (e.g., increased slip system activity and stress concentrations) are intensified in the vicinity of

misoriented grain and subgrain boundaries (cf. experimental and numerical results in [62]), then

one would expect the largest contributions to Fei to come from these regions.

The deformation tensor Fei represents in an average sense the incompatibility of the local microelastic deformation f

e1 within the volume element. If the local elastic unloading fe1 is uniform

(and hence, compatible) throughout vcur, then f

e 1⁄4 Fe

, ~f

e 1⁄4 1, and Fe 1⁄4 ~f 1⁄4 F

p 1⁄4 ~f

p 1⁄4 fp

, such that

Fei 1⁄4 1 (Fig. 3). However, Fei tells us nothing about the compatibility, or lack thereof, of the

434 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
average recoverable elastic deformation Fe1 from volume element to volume element. If we regard each vref as an entire single crystal, then the incompatibility of Fe1 measures the intergranular incompatibility between grains, while if we regard each vref as a subgrain, then the

incompatibility of Fe1 is an intragranular measure. By incompatibility of Fe1 we mean lack of

continuous, single-valued coordinates ~xa spanning the union of i local intermediate configuration

volume elements V

eint  [~v

ðiÞ

int. Thus, while the ~xa are assumed to be available and differentiable

within each local volume element ~v

ðiÞ

int, they may not be so in the global configuration Beint  V

eint. If

the ~xa are multi-valued or discontinuous in the global configuration, they are typically called

anholonomic coordinates [63].

We now appeal to the continuum theories of continuously distributed dislocations [64–66]

to characterize the incompatibility of global configuration Beint (i.e., the anholonomicity or lack

of integrability of Fe1

). The coefficients Ca

bc of the ‘‘crystal connection’’ (cf. [31,67]) acquire

the following form when referred to the current configuration:

Ca

bc 1⁄4 F eaa

oF e1a

oxb 1⁄4 F ea

a F e1ac;b 1⁄4 F e1ac F eaa;b; ð18Þ

with the subscripted comma denoting partial differentiation with respect to spatial coordinates as

indicated. In (18) we require Fe1 to be spatially differentiable to first-order; thus, for the present

discussion, local misorientations across grain and subgrain boundaries (i.e., between volume

elements vcur) are envisioned as steep gradients of lattice rotation, as opposed to actual discontinuities in the lattice arrangement. The torsion of the connection Ca

bc is written

t 1⁄4 2F ea

a F e1a1⁄2b;c ga gb gc 1⁄4 2F eaa F e1a1⁄2b j c

ga gb gc; ð19Þ

where we have replaced partial differentiation (subscripted comma) with partial covariant differentiation with respect to the Levi-Civita connection on Bcur whose Christoffel symbols stem

from the components of the metric gab of Eq. (6)3 and thus are symmetric in covariant indices. We

use the notation of vertical bars for such covariant differentiation. Additionally, the bracketed

indices are anti-symmetrized according to 2A1⁄2ab 1⁄4 Aab  Aba.

We now relate the torsion tensor t of Eq. (19) to the anholonomicity of the inverse elastic

deformation gradient field Fe1

. The incompatibility of regions within global configuration Beint is

associated with the net Burgers vector Be, defined in terms of the closure failure of the line integral

of dx~  Fe1 dx over a closed loop ~c in Beint [65,66]

Be   I

~c

dx~ 1⁄4  I

c

Fe1 dx: ð20Þ

Applying Stokes theorem (cf. [4]), we express Eq. (20) as



I

c

F e1a

a dxa 1⁄4

Z

a

F e1a

1⁄2a j b dxa ^ dxb 1⁄4

Z

a

e

abcF e1a

a;b nc da; ð21Þ

where dxa ^ dxb 1⁄4 eabcnc da 1⁄4 ffiffiffi 1 g p eabcnc da is the differential area two-form corresponding to oriented area element nda bounded by curve c (Fig. 4), which in turn is the current configuration

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 435
image of Burgers circuit ~c. The wedge product is denoted by ^, and components of the contravariant permutation tensor (i.e., alternator tensor) in the current configuration are denoted by

eabc 1⁄4 ð ffiffiffi 1 g p Þeabc, with eabc the standard permutation symbols.

The incompatibility is expressed in terms of the torsion of the crystal connection by combining

(18)–(21)1, i.e.

Bea 1⁄4  1

2

Z

a

F e1a

c t

c

ab dxa ^ dxb

: ð22Þ

By defining the operator

curlðÞ  ðoðÞ=oxÞ : e; ð23Þ

we may write for (21)2

Be 1⁄4

Z

a

curlðFe1Þ  nda 1⁄4

Z

a

Ae  nda; ð24Þ

where Ae  curlðFe1

Þ and the vector nda 1⁄4 ðna daÞga. Invoking Nansons formula ~nd~a 1⁄4

Je1n  Fe da (cf. [35]), we can transform (24) to an area integral over an intermediate configuration

region ~a

Be 1⁄4

Z

~a

~n  Je

Fe1AeT

|fflfflfflfflfflffl{zfflfflfflfflfflffl}

Aee

d~a 1⁄4

Z

~a

~n  A~e

d~a; ð25Þ

where the dot product in the integrand of (25) denotes contraction of the covariant index of ~n with

the contravariant index of Fe1

. Note that A~e only measures GNDs and does not account for

SSDs (e.g., closed dislocation loops) within vref [68,69].

Discussion regarding the length scales inherent in the line integral (20) is now in order. Two

distinct length scales are included in the integral (20) and are intimately related. One is the choice

of size of ~c: we let Dð~cÞ denote the diameter of a circle equivalent to ~c. The other is the choice of

‘ref, which defines the size of the volume element vref 1⁄4 ‘3

ref, and thus ~vint, and implicitly affects the

definition of the average elastic deformation gradient for the element, Fe

. Eq. (20) is really valid

Fig. 4. Burgers circuit in current configuration.

436 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
only for Dð~cÞ Plref, if we take ‘ref to be a characteristic diameter of the volume element ~vint

considered. Under this condition, Fe and its gradient are defined in regions within circuit ~c. If the

converse is true, such that Dð~cÞ < ‘ref, we need to consider the incompatibility of the inverse

microelastic deformation within the volume element, f

e1

, instead since Fe1 is not defined within

the volume element. Repeating the derivation leading up to (24) for the microelastic deformation

f

e1

, we arrive at an expression for a local Burgers vector ~be

~be   I

c

f

e1

dx 1⁄4

Z

a

curlðf

e1

Þ  nda: ð26Þ

Notice that ~be 1⁄4 ~bea

~ga 2 T ðdvpÞ resides in a different tangent space than the vector Be of (25). We

use barred Greek indices to denote components in this space, which arises from the local plastic

deformation f

p

, as shown in Fig. 3. Under certain conditions we can further characterize (26) in

terms of the spatial gradient of the inverse of the residual elastic deformation gradient ~f

e

. Since

both d~vint and dvcur inhabit local volume elements deformed in a compatible fashion, we can define

a local compatible tangent mapping between them: f

c1  ox~

ox : T ðdvcurÞ ! T ðd~vintÞ. If no additional

local plastic deformation occurs upon external unloading of the volume element, in the process

symbolized by Fe1

, then microscopically we have f

p 1⁄4 ~fp and dvp 1⁄4 d~vp, as shown in Fig. 5.

Furthermore, we may write

f

e 1⁄4 fc~fe

; ð27Þ

thereby decomposing the total microelastic deformation into a residual part ~f

e and a compatible

part f

c due to the applied stress [70]. Substituting (27) into (26) and invoking the compatibility

conditions f c1b

1⁄2a;b 1⁄4 0 then gives

~bea 1⁄4

Z

a

e

cabð~f e1a

b;b f c1ba þ ~f e1ab f c1b

a;b Þnc da 1⁄4

Z

a

e

cabð~f e1a

b;b f c1b:a Þnc da; ð28Þ

Fig. 5. Tangent mappings and volume elements when f

p 1⁄4 ~fp

.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 437
demonstrating that ~be of Eq. (26), representing the anholonomicity of the local elastic deformation gradient, depends only upon the spatial gradient of the (inverse of the) residual elastic lattice

deformation map ~f

e

. and not higher order gradients of the locally compatible and recoverable

elastic lattice deformation field f

c

.

3. Multiscale balance laws

Neglecting body forces, let the crystal element be subjected to arbitrary but self-equilibrating

surface traction t, measured per unit area in the reference configuration, with corresponding

outward unit normal covector nref (Fig. 2). Let sAa represent the contravariant components of the

local two-point nominal stress tensor (the transpose of the first Piola–Kirchhoff stress tensor). The

traditional (non-polar) local balances of linear and angular momentum, respectively, are written

for quasi-static conditions as

s

Aa

jA 1⁄4 0; f aASAb 1⁄4 SAaf bA; ð29Þ

with nrefASAa

1⁄4 t

a on Sref, and with the vertical bar denoting covariant differentiation with respect

to the symmetric Levi–Civita connection on bref. Notice that Eq. (29) are applicable locally, for

points (i.e., differential elements) within the crystal element. We next define the average nominal

stress tensor for the volume element as (cf. [71])

S 

1

vref Z

vref

sdvref 1⁄4 1

vref Z

sref

X tdsref; ð30Þ

where Gausss theorem has been used for Eq. (30)2, along with the assumptions of quasi-static

conditions, Eq. (29)1, and stress continuity within the volume element. The definition of the

contravariant Cauchy stress R then follows as

Rab  J 1

F a

ASAb; ð31Þ

with J  detF ffiffiffiffiffiffiffiffiffig=G p . The macroscopic rate of the deformation gradient, F_ , is given by

F_ 

1

vref Z

vref

_

f dvref 1⁄4 d

dt

1

vref Z

vref

f dvref  1⁄4 1

vref Z

sref

x_ dsref; ð32Þ

where x_ðX; tÞ 1⁄4 d

dtðxagaÞ is a compatible (material) surface velocity. Consider the following

equality [72]:

1

vref Z

vref

_

fsdvref  FS_ 1⁄4 1

vref Z

sref

ðx_  FX_ Þðnrefðs  SÞÞdsref: ð33Þ

The right-hand side of Eq. (33) is zero for linear displacement, constant traction, and certain

periodic boundary conditions specified with respect to the reference configuration. Assuming

438 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
henceforward such boundary conditions apply as representative for our crystal element, the

volume-averaged stress work rate per unit volume is found as

1

vref Z

sref

t

a

gabx_

b dsref 1⁄4 1

vref Z

vref

SAagab _

f b

A dvref 1⁄4 SAaF_ aA; ð34Þ

where we have used Gausss theorem in (34)1.

The balance of energy for a purely mechanical process, considering the localized form of (34)

in the absence of heat flux due to conduction or internal heat supply, is written as

u_ 0 1⁄4 SA

a

_

f a

A; ð35Þ

with u0 the local internal energy per unit volume in the reference configuration. The corresponding

averaged energy balance then follows from Eq. (34)2 and (35) as

U_ 0 1⁄4 SA

aF_ aA; ð36Þ

where U0  ðvrefÞ

1 R

vref u0dvref is the average internal energy per unit reference configuration

volume. Notice that Eqs. (35) and (36) are constitutive assumptions, prohibiting conjugate microstresses to second-order gradients of deformation from performing mechanical work on the

external boundaries of dvref and vref, respectively. These rather strong assumptions are in agreement with previous treatments by Teodosiu [3], Steinmann [32], Acharya and Bassani [33,73],

Menzel and Steinmann [9], Bammann [34], Bassani [74], and Regueiro et al. [39]. On the other

hand, they are contradictory to the models of Teodosiu [1], Dillon and Kratochvil [5], Naghdi

and Srinivasa [38,28], Le and Stumpf [30,31], Shizawa and Zbib [16,17], and Gurtin [13,36].

Introducing g0 as the average specific entropy per unit reference volume, the Clausius–Duhem

inequality is written

g_0 P0; ð37Þ

again in the absence of heat flux due to conduction or internal heat supply. The average Helmholtz free energy per unit volume in the reference configuration, w0, is then defined by

w0  U0  g0h; ð38Þ

with h the absolute temperature. Substituting Eqs. (36) and (38) into the Clausius–Duhem inequality (37) and assuming stationary temperature then yields a reduced form for the entropy

inequality:

SA

aF_ aA P w_

0: ð39Þ

The remaining balance laws at the level of the volume element are now presented:

q0 1⁄4 qJ; SAa

jA 1⁄4 0; F aASAb 1⁄4 SAaF bA: ð40Þ

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 439
Eq. (40)1 provides the usual relationship between current, q, and reference, q0, average mass

densities. Eq. (40)2 and (40)3 follow from volume averaging the microscopic balance laws in Eq.

(29) and by assuming that boundary conditions imposed on the crystal volume element are such

that the right-hand side of Eq. (33) is zero throughout the deformation process. The last of (40)

leads to symmetry of the Cauchy stress R, as defined in Eq. (31).

Notice that inequality (39) is written in terms of two-point tensors S and F and the free energy

per unit reference volume w0. We shall later find it useful to express (39) in terms of quantities with

all components referred to intermediate configuration ~bint (i.e., the relaxed intermediate configuration of elastoplasticity), which we view as the most convenient configuration for deducing

thermodynamic restrictions and posing constitutive assumptions (see also [16,17,39]). Expression

of (39) in the intermediate configuration is achieved here by first mapping to the current

configuration, then pulling back by the two-point elastic deformation gradient Fe

. From (31) we

have

SA

aF_ aA 1⁄4 JRab gbcF_ cAF 1a A

|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} Lba

1⁄4 JRabLba; ð41Þ

with L the spatial velocity gradient. Substituting into (39) and multiplying by the inverse of the

Jacobian invariant eJ  det Fe ffiffiffiffiffiffiffiffiffig~=G p > 0 then gives the Clausius–Duhem inequality in the current

configuration, i.e.

ðRabLba þ w0J_ 1

Þ Pw_; ð42Þ

where the free energy per unit current volume is w  J 1w0. From the symmetry of R and g1, and

using the standard identities _ J 1J 1⁄4 0 and J_ 1⁄4 JgabLba, Eq. (42) can be simplified to read

hðR  wg1

Þ;LvðgÞiP2w_; ð43Þ

with Lv the Lie derivative with respect to the spatial velocity field v 1⁄4 x_  u1, and with

hA;Bi 1⁄4 trðABÞ the dual product defined for second-rank tensors A and B. We next define the

velocity gradient Le referred to intermediate configuration ~bint as

eLa

b  F e1a

a La

bF eb

b ð44Þ

and we define the mixed-variant elastic second Piola–Kirchhoff stress Se and the Mandel stress [75] fM as

eSa

b  Je

F e1a

a RacgcbF eTbb ; Me a

b  Je

F eTa

a RacgcbF eTbb 1⁄4 Ceeav eSv

b; ð45Þ

where Ceea

v

 g~abF ea

v gabF ebb and Je  det Fe ffiffiffiffiffiffiffiffi

g=g~ p . Then from relations (39), (41) and (45)2, we can

obtain

eJ eLa

bMe Tb

a Pw_

0: ð46Þ

440 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
Multiplying (46) by eJ 1 > 0 and using the relation _ eJ 1 1⁄4 eJ 2 _

eJ 1⁄4 eJ 1 _

F

ea

AF

e 1A

a , along with

we  eJ 1

w0, then gives

eLa

bMe Tb

a  we _

F

ea

AF

e 1A

a P _

we; ð47Þ

which is the reduced entropy inequality mapped to configuration ~bint, with energetic quantities

defined on a per-unit-volume basis with respect to local element volume ~vint.

4. Crystal plasticity model

Hereafter we more fully develop a continuum crystal plasticity model, focusing on a formulation couched at a single length scale ‘ref rather than one phrased in terms of volume averages.

This is necessary for typical numerical implementations, wherein each volume element vref corresponds to the local region about an integration point in a finite element simulation, for example.

The microscopic deformation gradients f, ~f, f

e

, f

p

, ~f

e

, ~f

p

, and f

c were introduced in Section 2 for

illustrative purposes and are not calculated explicitly in this model. The same can be said for the

microscopic kinetic variables––such as t, s and u0, for example––of Section 3. Instead, we consider

here the evolution of the average deformation gradients for the crystalline volume element (e.g., F,

Fe

, Fei, and F

p

), the average stresses (e.g., S and R), and the average energies (e.g., U0 and w0). Since

the microscopic deformations and local lattice rearrangements within the volume element are no

longer explicitly tracked, we are unable to invoke Eq. (9) to determine Ve

, Eq. (10) for Re, Eq. (15)

for F

p

, Eq. (17) for Fei, or Eq. (26) to calculate the local incompatibility ~be within microscopic

subregions of the crystal element. Likewise, we are unable to calculate local residual stress fields

and elastic energies associated with ~f

e in the model forthcoming in Section 4. Instead, we must rely

upon additional constitutive assumptions and corresponding balance relations to ready our model

for implementation––these are considered in detail in what follows. Of course, the constitutive

assumptions we make hereafter are motivated by the physically-based, multiscale averaging

treatment of Sections 2 and 3.

4.1. Kinematics, balance laws, and thermodynamics: summary

The fundamental thermomechanical relations already derived in Sections 2 and 3 and applied

now to a ‘‘single-scale’’ crystal plasticity formulation are restated here for clarity and ease of

reference.

Deformation gradient:

F 1⁄4 ox

oX 1⁄4 Fe Fei

F

p

|ffl{zffl}

eF

1⁄4 VeRe

|ffl{zffl}

Fe

VeiRei

|ffl{zffl}

eFi

F

p

: ð48Þ

Dislocation density tensor (GNDs):

A~e 1⁄4 Je

Fe1ðcurlFe1Þ

T

: ð49Þ

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 441
Mechanical stresses:

fM 1⁄4 ~g1Cee

Se 1⁄4 Je

FeTRgFeT 1⁄4 eJ 1FeTFSgFeT: ð50Þ

Mass conservation:

q0 1⁄4 qJ: ð51Þ

Balance of linear momentum (quasi-static):

DivðSÞ 1⁄4 0: ð52Þ

Balance of angular momentum (quasi-static, intermediate configuration):

Se  SeT 1⁄4 0: ð53Þ

Balance of energy (mechanical case, intermediate configuration):

hfMT; Fe1

FF_ 1

Fe

|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}

eL

i 1⁄4 eJ 1

U_ 0: ð54Þ

Entropy inequality (mechanical case, intermediate configuration):

hfMT; Leihwe~1; _

FeFe1i P _we: ð55Þ

4.2. Free energy and consequences of the dissipation inequality

We make the following general assumption regarding the dependency of the Helmholtz free

energy function for the crystal volume element, referred to relaxed intermediate configuration ~bint

and neglecting temperature effects:

we 1⁄4 weðCee

; Vei; A~e; ~eSS; ~gÞ 1⁄4 weðCee

ab; V

e i

a

b ; Aeeab

; ~eSS; g~abÞ: ð56Þ

The covariant elastic strain tensor Cee

ab  F ea

a gabF ebb is included to model the change in average

energy with a change in external loads, a standard assumption in finite crystalline elastoplasticity

theories (cf. [76]). The left stretch tensor derived from Fei

, denoted by Vei, is included to reflect the

contribution to the free energy from residual microelasticity within the volume element, and is

non-negligible when the deformation within the volume element is heterogeneous (e.g., during

grain subdivision). Kratochvil [77] and Lion [78] made similar constitutive assumptions. The

elastic energy due to the average elastic curvature in the volume element attributed to GNDs is

reflected by the inclusion of A~e

, a particular assumption also suggested explicitly by Bammann

[34], Gurtin [36], and Regueiro et al. [39]. Since the GND density tensor does not include a

442 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
measure of the total length of all dislocation lines (e.g., SSDs consisting of closed dislocation

loops and dislocation lines of opposing signs), the scalar parameter ~eSS  ~b ffiffiffiffiffiffi

q~SS p is included to

model the net contribution of the elastic self-energy of the SSDs to the total free energy (cf.

[79,80]). The scalar norm of the Burgers vector (e.g., the lattice parameter) is denoted by ~b, while

the line length of SSDs per unit intermediate configuration volume is written q~SS. We regard ~eSS as

a local residual lattice strain measure due to the presence of SSDs [34].

It is emphasized that all independent variables in the free energy function are ‘‘elastic’’ variables, in contradiction to works of Naghdi and Srinivasa [38,28], Gurtin [13], and Svendsen [37],

among others, who include a dependency upon the plastic deformation gradient and/or higher

gradients of plastic deformation. We treat the plastic deformation gradient Fep as a continuum

idealization of the rigid sliding of portions of the lattice due to relative motion of elastic blocks of

material contained between slip planes. Such motions clearly do not alter the energetic properties

of material within the blocks, and therefore should not influence the Helmholtz free energy.

Reasoning behind inclusion of the metric ~g will be explained later.

Notice that each of the constitutive variables in (56) is invariant with respect to superposed

rigid body motions in the current configuration: x ! Qx þ c, with Q 2 SO3 a rigid body rotation

matrix and c a constant translation vector:

Fe ! QFe ) Cee ! Cee

; Fei ! Fei ) Vei ! Vei;

A~e 1⁄4 Je

Fe1ðcurlFe1Þ

T ! Je

ðdetQÞFe1QTðcurlðFe1QTÞÞT 1⁄4 A~e;

q~SS ! q~SS ) ~eSS ! ~eSS; ~g ! ~g:

ð57Þ

Even though Fei is unaffected by rigid transformations of the current configuration [77], we choose

not to include the rotational part of Fei

, denoted by Rei, in the free energy function because rigidbody rotations of the volume element in the unloaded configuration ~bint do not influence the

stored elastic energy.

We define the conjugate stresses to the independent state variables included in (56):

See 

owe

oCee ; Sei 

owe

oVei

; SeGN 

owe

oA~e ; eS SS 

owe

o~eSS

; Seg 

owe

o~g : ð58Þ

Notice that SeGN is a couple stress, with units of force/length, while the remaining stress measures,

See

, Sei, eS SS, and Seg possess the standard stress dimensions of force/length2. Expanding _we in (55)

using the definitions of (58) gives

eLa

bMe Tb

a  we _

F

ea

AF

e 1A

a P eSeab _

Cee

ab þ eSi

	b

a

_

V

ei

a

b þ eS GN

ab ~_

A

eab

þ eS SS~_

eSS þ eSgab ~_

gab; ð59Þ

where the dual map [81,82] is denoted by ð Þ	 and satisfies ðAa

bÞ

	 1⁄4 A	a

b . From the three-term

multiplicative decomposition of the deformation gradient, F 1⁄4 Fe

Fei

F

p

, of Eqs. (16) and (48), the

velocity gradient in the intermediate configuration, Le, can be partitioned as

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 443
Le  Fe1 FF_ 1

|ffl{zffl}

L

Fe 1⁄4 Fe1

F_ e

|fflffl{zfflffl}

eLe

þ _Fei

Fei1

|fflffl{zfflffl}

eLi

þ Fei _F

p

F

p1

Fei1

|fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}

eLp

: ð60Þ

Further algebraic manipulations give

hLe; fMTi 1⁄4 Se~g1

;

1

2

_

Cee

  þ hfMT; Lei þ Lep

i; _Fei

Fei1

|fflffl{zfflffl}

eLi

1⁄4 _Vei

Vei1 þ Vei _

Rei

ReiT

|fflffl{zfflffl}

Wei

Vei1

; ð61Þ

where the spin fWi is skew via g~abWe i

b

v 1⁄4 g~vbWe i

b

a . Using Eqs. (60) and (61) in inequality (59) then

results in

1

2

eSa

vg~1vb 

 eSeab _

Cee

ab þ ððMe Tb

v  wd e b

vÞV

ei1v

a  eSi	ba Þ _

V

ei

a

b  eSgab ~_

gab

þ ðV

ei

b

v Me Tvc V

e i1c

a ÞWe i

a

b  eS GN

ab ~_

A

eab

 eS SS~_

eSS þ ðMe Tb

a  wd e b

aÞL~pa

b P 0: ð62Þ

Assuming that (62) must hold for independent specification of the elastic strain rate and the other

rate variables [83,84] we arrive at standard relationships between the (mechanical) elastic stresses

Se, R, S, and T, and the lattice stress See conjugate to Cee

2See 1⁄4 2 owe

oCee 1⁄4 Se~g1 1⁄4 Je

Fe1RFe	 1⁄4 eJ 1FSF e e	 1⁄4 eJ 1Fe1TFe	; ð63Þ

where T  ST is the contravariant first Piola–Kirchhoff stress tensor, and the final equality follows

from (40)3. More familiar and compact constitutive equations for the Cauchy stress and first

Piola–Kirchhoff stress are then derived readily with the chain rule of differentiation

Rab 1⁄4 2Je1 owe

oCee

ab

F ea

a F ebb

|fflffl{zfflffl}

1⁄4oeCe

ab=ogab

1⁄4 2Je1 owe

ogab


 2 ow

ogab

; ð64Þ

T aA 1⁄4 2eJ F 1A

b F eba

owe

oCee

ab

F ea

b 1⁄4 eJ F

e 1A

a

ow~

oF ec

a |fflfflfflfflffl{zfflfflfflfflffl}

1⁄4oew=oF c

A

gca 1⁄4 eJ

owe

oF c

A

gca 


ow0

oF c

A

gca: ð65Þ

The final equality in (64) is rigorous only when oJe1=og 1⁄4 0 (e.g., Cartesian current coordinates),

while the final equality in (65) is ensured only when oJ

~=oF 1⁄4 0 (e.g., when the total inelastic

deformation is isochoric). We next define (cf. [85])

444 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
Pe  fMT  we~1 1⁄4 Je

Fe1 ðRg  w1Þ

|fflfflfflfflfflffl{zfflfflfflfflfflffl}

p

Fe 1⁄4 eJ 1Fe ðFTTG  w010Þ

|fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}

P

Fe1 ð66Þ

as the push-forward of the mixed-variant reference configuration Eshelby 1 [86] energy–momentum tensor P, or equivalently as the pull-back of the current configuration Eshelby stress

tensor p. In Eq. (66), 10, ~1, and 1 denote mixed-variant identity maps on configurations bref, ~bint,

and bcur, respectively. Upon substituting (63) and (66) into (62) and assuming that ~_

g 1⁄4 0, we arrive

at a reduced form of the dissipation inequality:

hPe Vei1  Sei	

; _Veiþ hVei

Pe Vei1

; fWiihSeGNT; ~_AeiheS SS; ~_

eSSiþhPe; De p þ fWp

i P0: ð67Þ

Notice that both the plastic deformation rate Depa

b  g~adð~gLep

ÞðdbÞ and the plastic spin

We pa

b  g~adð~gLep

Þ1⁄2db contribute to the dissipation, as does the spin associated with residual elasticity/subdivision fWi

. Note also the prominent role of the Eshelby stress tensor Pe as a force

conjugate in the plastic dissipation term (see also [85,87]).

4.3. Constitutive model

A more explicit form of the strictly mechanical, intermediate configuration free energy function

(56) is now proposed on physical grounds:

we 1⁄4 1

2

hEee

;ðCeeffðEeiÞÞ : Eeei þ c1lhEei; Eeei þ c2l

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

hEei; Eeii

q

þ c3lð~eSSÞ

2 þ bGNlhA~e

#; A~e#i: ð68Þ

The strains Eee and Eei follow from the decomposition of the total strain tensor Ee referred to ~bint:

Ee 

1

2 FeTðFTF  10ÞFe1 1⁄4 Eee þ Eei þ Eep

; ð69Þ

where the elastic ðEee

Þ, heterogeneity ðEeiÞ, and plastic ðEepÞ parts are given by

Eee 

1

2

ð~g1Cee  ~1Þ; Eei 

1

2

ð~1  VeiTVei1

Þ; Eep 

1

2 FeiTð~1  F

pT

F

p1

ÞFei1: ð70Þ

Notice that Eee and Eei are functions of Cee and Vei

, respectively, in agreement with the general form

for the free energy in proposition (56). The mixed-variant GND variable A~e

# is defined as

A~e

#  A~e

~g: ð71Þ

1 Eshelby [86] in fact used the negative of our P as the referential energy–momentum tensor, with the material

gradient of displacement in place of the deformation gradient. We follow here for convenience the definition and sign

convention of Le and Stumpf [87] and Le et al. [76].

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 445
Components of the fourth-rank effective elastic modulus tensor for the volume element in the

intermediate configuration are denoted by ðCeeffÞ

av

bd

, and l is a characteristic elastic constant for a

perfect reference lattice (e.g., an average reference shear modulus). Finally, c1; c2, and c3 are dimensionless scalar constants, and bGN is a (possibly) evolving scalar parameter that will be addressed more later (Eqs. (72) and (73)), with dimensions of length2. Notice the role of ~g in defining

the mixed-variant tensors Eee and A~e

#, rendering its inclusion in the general free energy function

(56) a necessity [87].

We now discuss the physical reasoning behind our choice of each term in the specific free energy

function (68). The first term, 1

2 hEee

;ðCeeffðEeiÞÞ : Eeei, is reminiscent of the quadratic form seen in

finite linear hyperelasticity, differing only in the sense that the effective elastic modulus tensor for

the volume element, referred here to the intermediate configuration ~bint, is assumed to depend

explicitly upon Eei

. This assumption is intended to reflect the influence of local microelastic rotations which arise, quite possibly, from grain subdivision and intragranular substructure development during plastic deformation (cf. Hughes et al. [43], Butler and McDowell [50]). For

example, for a single crystal consisting of subgrains exhibiting a random distribution of local

misorientations unrestricted in magnitude, the effective moduli will approach those of a random

polycrystal (e.g., elastic isotropy in many metals).

The second term, c1lhEei

; Eeei, accounts for the aforementioned amplification of internal residual microstress fields (and the corresponding internal elastic energy) at flexing subgrain

boundaries [59,60] with increases in the applied stress. The third term, c2l

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

hEei; Eeii

q

, represents a

portion of the stored energy of cold work attributed to heterogeneous elastoplasticity within the

volume element. The linear dependency of energy upon the effective ‘‘meso-incompatibility’’

strain Eei is motivated from previous solutions obtained from computational micromechanics [88]

in which low-angle and high-angle grain boundary arrangements were assigned as initial conditions to a deforming polycrystalline aggregate. Since the present model is intended for single

crystals, we presume that data for low-angle grain boundary arrangements from this previous

study are most applicable (Fig. 6(a)), and our form of the free energy term reflects an assumption

of similitude between single crystals and polycrystals with low-angle boundaries. Thus, this term

should be regarded as a rough initial approximation.

The fourth term, bGNlhA~e

#; A~e#i, represents a quadratic dependency of the free energy upon the

average density tensor of GNDs. Such an assumption is rather standard in gradient finite strain

single- and polycrystalline elastoplasticity theories in the literature (cf. [9,32,34,37,39]), as is the

linear dependence upon the total length of SSDs per unit volume, q~SS 1⁄4 ð~eSS=~bÞ

2 [34,39,37]. Additional experiments and/or numerical simulations are clearly needed to determine the parameters

c1, c2, and c3 for a particular material.

We remark that the parameter bGN provides a squared effective ‘‘radius of non-local action’’ in

our model (see also later Eqs. (90) and (93)). As bGN ! 0, the contribution of GNDs to the free

energy function becomes negligible, and the model assumes a local character. We hypothesize that

‘ref 1⁄4 ffiffiffiffiffiffi

vref p3 ––the size of the representative crystal volume element (Eq. (1))––should influence the

value of bGN in some way. In other words, the size of the local crystal volume element over which

the average dislocation density tensor A~e is calculated should have some effect upon the range of

non-local interactions. For example, if we assume that the strength of non-locality is inversely

related to scale (cf. [24]), then an obvious choice is

446 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
bGN 1⁄4 ‘

2

GN

‘GN

‘ref c4

; ð72Þ

with ‘GN a length parameter that is characteristic for a given material (cf. [39]) and with c4 another

dimensionless material constant. In other words, if we regard ‘GN and c4 as fixed parameters

(constants) for a given material, then the non-local radius will scale with the size of the volume

element according to ‘

c4=2

ref . Of course, we would now need a more extensive series of tests to

determine the two values ‘GN and c4, rather than just bGN. Interestingly, if we take ‘ref 1⁄4 LG, the

grain size, and c4 1⁄4 1, then a Hall–Petch type relation [18,19] for hardening due to GNDs is

acquired via Eq. (72).

Additionally, one may argue that since the radius of non-locality should reflect a characteristic

dimension of the microstructure, such as cellular structure size in ductile single crystals that

subdivide (cf. [7]), then the effective length ‘GN should be an evolving, rather than stationary,

parameter. In our framework, we could easily extend (72) to

Fig. 6. Residual elastic energy versus

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ð2=3ÞhEei; Eeii

q

(a) and evolution of

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

hEei; Eeii=hEee; Eeei

q

versus applied effective

strain (b).

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 447
bGN 1⁄4 ^‘

2

GN

^‘GN

‘ref !c4

; ^‘GN 1⁄4 ^‘GNðVeiÞ; ð73Þ

allowing the effective scaling factor ^‘GN to evolve with the heterogeneity parameter (i.e., the grain

subdivision metric) Vei

. If Eq. (73) is invoked, then the presence of ^‘GN in the free energy function

(68) will influence the values of thermodynamic force conjugates calculated via Eqs. (58) contributing to the dissipation inequality (67).

Consider a boundary value problem where F and its rate are prescribed incrementally as a

function of time. In order to characterize completely the state of the material from the kinematic

standpoint of Eq. (48), we need to specify the deformation gradient measures Fe

, Fei, and F

p

. In a

hyperelastic setting (cf. [89,90]), a standard approach is to formulate evolution laws for (time rates

of) F

p and Fei

, thus leaving Fe to be determined from the product FðFeiF

p

Þ

1

, assuming that the rate

of F and the driving forces for each inelastic deformation measure are known at the beginning

of each time increment. Such is our approach: the requisite evolution equations for F

p and Fei are

considered next.

The time rate of average plastic deformation gradient in the single crystal volume element is

assumed to follow from the standard kinematic framework of finite crystal plasticity [80,91–93],

i.e.,

_

F

p 1⁄4 Xn

i1⁄41

c_

i

s

i m i

 !

|fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

Lp

F

p

; ð74Þ

where c_

i

, si, and m i are the shearing rate, slip direction (a unit vector), and slip plane normal (a

unit covector) for slip system i, all defined with respect configuration bp of Fig. 2. As written in

(74), Lp  _

F

p

F

p1 is the mixed-variant average plastic velocity gradient in configuration bp. When

the single crystalline lattice is initially heterogeneous in the reference state, or when the lattice

deforms heterogeneously within the volume element, c_

i

, si, and m i are understood to be suitablydefined spatial averages of their fluctuating local counterparts. In accordance with classical crystal

plasticity theory, the representation of the slip directions and slip plane normals is assumed

identical in each of the unloaded configurations bref and bp. The slip direction unit vectors and slip

plane normal covectors are orthogonal in the unloaded configurations (i.e., hm i

;sii 1⁄4 0), and are

typically given as initial conditions in a boundary value problem. However, their representation

changes in the current configuration as a result of the elastic deformation gradient Fe (cf. [76,93])

ðs

i

Þ

a 1⁄4 F ea

a da

aðs

i

Þ

a

; ðmiÞa 1⁄4 F e	a

a da

aðmi

Þa: ð75Þ

In Eq. (75), si and mi are push-forwards of their counterparts in bp (we have used barred Greek

indices for components in bp) and are not necessarily of unit length when the elastic stretch Ve is

significant, although each pair of si and mi does remain orthogonal in configuration bcur, as is

easily verified by direct calculation. These orthogonality relations effectively prohibit dislocation

climb, a non-isochoric process. Two useful consequences arise from kinematic assumption (74)

and the slip plane normal-slip direction orthogonality relations. The first is

448 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
_

J p 1⁄4 J p

trðLpÞ 1⁄4 X

i

c_

i

hm i;s

i

i 1⁄4 0; ð76Þ

meaning that the plastic flow is isochoric. The second is a reduced form of the ‘‘purely plastic’’

dissipation in the entropy inequality (67):

hPe; De p þ fWp

i1⁄4hFei1

Pe FeiT

|fflfflfflfflfflffl{zfflfflfflfflfflffl}

P

;Lpi1⁄4hP;Lpi 1⁄4 X

i

pi

c_

i

; ð77Þ

where P is the Eshelby stress tensor mapped to ‘‘plastic’’ configuration bp and

pi  Pa



b ðs

i

Þ

b

ðmiÞa ð78Þ

is the projected Eshelby stress on system i. We could equivalently use the pull-back of the

transposed Mandel stress, Fei1fMTFeiT , in place of P in (77) and (78) since the plastic flow is

volume-preserving. By assuming a particular flow potential U for the slip rates c_

i

, we can ensure

that the dissipation in (77) is non-negative, such that the entropy inequality is satisfied automatically when the contributions of Fei

, A~e, and q~SS to the dissipation are neglected. We conveniently prescribe

U 1⁄4 Xn

i1⁄41

U^i

ðpi; Fei; SeGN; eS SSÞ 1⁄4 _k

X

i

ji

m þ 1

pi

ji





















mþ1

; ð79Þ

where each slip potential U^i is a scalar function of its arguments. In Eq. (79)2, _k > 0 (dimensions

of 1=time) and m (dimensionless) are material parameters which we assume are constant on each

slip system, while ji

ðFei; SeGN; eS SSÞ > 0 are scalar slip resistances (dimensions of force/length2) that

evolve with inelastic deformation. The slip rates are then assumed to adhere to

c_

i 1⁄4 oU

opi 1⁄4 oU^i

opi 1⁄4 _

ksgnðpiÞ pi

ji





















m

ðno sum on iÞ; ð80Þ

a power-law form reminiscent of typical flow rules proposed for classical crystal viscoplasticity

[94], except that the resolved Eshelby stress is used as a conjugate variable to the slip rates in our

theory. Substituting (80) into (77), we find that the purely plastic dissipation is unconditionally

non-negative:

X

i

pi

c_

i 1⁄4 _

k

X

i

pi

sgnðpiÞ pi

ji





















m

P 0: ð81Þ

Our thermodynamic motivation for choosing the resolved Eshelby stress as a driving force for

viscoplasticity is clear from (67) and (81). Our physical reasoning for selecting the Eshelby stress

stems its relationship to the configurational force on a dislocation [85], as opposed to the usual

resolved Cauchy stress of Schmids law [95]. We also cite the work of Le et al. [76] who used

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 449
thermodynamic and physical arguments similar to ours to motivate usage of the Eshelby stress as

a driving force for plastic deformation in a finite strain, rate-independent single crystal plasticity

model.

We propose the following dependency for the positive scalar slip resistances ji

:

ji

ðFei; SeGN; eS SSÞ 1⁄4 jiðJi1Fei	ðSeGNÞ

TFei

|fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

SGN

; ffiffiffiffiJi 

p1 eS SS

|fflfflfflffl{zfflfflfflffl}

SSS

Þ 1⁄4 jiðSGN; SSSÞ; ð82Þ

with SGN and SSS complete pull-backs of SeGN and eS SS to configuration bp, and with

Ji  det Fei ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

det ~g= det g p > 0 (g is now introduced as a local metric on the anholonomic space bp).

The thermodynamic forces SGN and SSS satisfy

ðSGNÞab 1⁄4 Ji1F

ei	a

a

owe

oA~eab F

eib

 b 1⁄4 owe

oAeab ; SSS 1⁄4 ffiffiffiffi

Ji 

p1 owe

o~eSS

1⁄4 owe

oeSS

; ð83Þ

with

Ae  Ji

Fei1

A~eTFei	; eSS  ffiffiffiffi

Ji p ~eSS ð84Þ

the densities of GNDs and SSDs, respectively, pulled back from ~bint to bp. As a more particular

constitutive assumption, we represent the hardening variable for each slip system as a linear

combination of the thermodynamic force conjugates to the densities of GNDs and SSDs partitioned to all systems, i.e.,

ji 1⁄4 ji

0 þ ffiffiffiffiffiffiffiffi

bGN



p1 hi

GNjSj

GN þ hi

SSjSj

SS; ð85Þ

where ji

0 is the initial friction or threshold stress for system i, hiGNj and hiSSj comprise n  n slip

system interaction matrices (dimensionless units) to be obtained from experimental measurements

of self- and latent hardening (cf. [96]), bGN is the length parameter for normalization,

Sj

GN  jhs

j

; SGN  m jij ð86Þ

is the scalar magnitude of the projected lattice couple stress associated with the GND tensor on

slip system j (with the dot product here denoting a contraction of indices on the cotangent space

of bp), and Sj

SS is the lattice stress associated with the SSDs on slip system j. The strain measures of

SSDs, ei

SS, are defined by an additive decomposition of the statistically stored line densities on

each slip system, qi

SS:

qSS 1⁄4 Xn

i1⁄41

qi

SS ) eSS  ~b ffiffiffiffiffiffi

qSS p 1⁄4 ~b

ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Xn

i1⁄41

qi

SS s



ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Xn

i1⁄41

ðei

SSÞ

2

s

) e

i

SS 1⁄4 ~b

ffiffiffiffiffiffi

qi

SS q

: ð87Þ

The corresponding thermodynamic force is then, from (58)4, (84), and (87),

Si

SS 

owe

oei

SS

1⁄4 owe

oeSS

oeSS

oei

SS

1⁄4 owe

o~eSS

|{z}

c3l~eSS

o~eSS

oeSS

|{z}ffiffiffi

J  i p1

ei

SS

eSS ! 1⁄4 c3l ffiffiffiffi

Ji 

p1

~e

i

SS 1⁄4 c3lJi1

e

i

SS; ð88Þ

450 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
providing a hardening contribution in (85) proportional to the square root of the dislocation line

density per unit volume, as originally suggested by Taylor [97] and since verified experimentally

numerous times (cf. [56,98]), and used in both local crystal plasticity theories (cf. [99,100]) and

gradient-based theories [34,39,37]. The GND (scalar) couple stresses Si

GN are found explicitly in

terms of the tensorial density of GNDs by inserting (58)3 into (86):

Si

GN 1⁄4 s

ia

Ji1F

ei

	a

a

owe

oA~eab

F

eib

 b g

bv

mi

v





















1⁄4 bGNlJi1 js

ia

F

ei

	a

a A~ebaF

eib

 b g

bv

mi

vj

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Ai

GN

1⁄4 bGNlJi1

Ai

GN ð89Þ

with Ai

GN a scalar measure of GNDs on slip system i. Substituting Eqs. (88) and (89) into Eq. (85)

then gives

ji 1⁄4 ji

0 þ lJi1 ffiffiffiffiffiffiffiffi

bGN p hi

GNjAj

GN 

þ c3hi

SSj

e

j

SS	

: ð90Þ

The contribution of GNDs to individual glide system resistances was previously suggested by

Acharya and Bassani [33,73], Acharya and Beaudoin [101], and Bassani [74], for rate-dependent

and/or rate-independent single crystal plasticity, similarly to the second term in our Eq. (90).

Many have used the conjugate stress to a tensorial GND measure as a contribution to a

backstress on each slip system, rather than a contribution to slip system resistance ji as in our

Eq. (90), with a variety of different ways proposed for projecting the tensorial backstress onto

individual glide systems [9,28,32,34,36,37,101]. If both a backstress and a glide stress are desirable

in order to completely characterize the hardening behavior (cf. [102]), then we can generalize the

flow rule in Eq. (80) to

c_

i 1⁄4 _

ksgnðpi  viÞ pi  vi

ji





















m

; ð91Þ

where increases in the friction stress ji > 0 now manifest only from SSDs:

ji 1⁄4 ji

0 þ ðc3lJi1

Þhi

SSj

e

j

SS; ð92Þ

and where the backstress vi

, not necessarily positive in sign, depends solely upon the density of

GNDs (cf. [34]):

vi 1⁄4 ðlJi1 ffiffiffiffiffiffiffiffi

bGN p Þhi

GNjAj

GN: ð93Þ

In Eqs. (92) and (93), hi

SSj and hiGNj are slip system interaction matrices, different in value, but not

in function, from those introduced already in (85). Different, and perhaps more inclusive, ways of

expressing the contributions of SGN to hardening on each slip system than the direct projection

method of Eq. (86) also merit further exploration (see e.g. [28,35,36]).

While the density of geometrically necessary dislocations can be calculated directly from spatial gradients of the elastic lattice deformation (Eq. (49)), the densities of statistically stored

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 451
dislocations on each slip system are modeled here as internal state variables (ISVs) equipped with

separate equations dictating their evolution. By definition, the time rate of change in dislocation

density for a given slip system, q_ i

SS, and the rate of the corresponding dimensionless lattice strain

measure e_i

SS for that slip system are related by

2e_i

SS 1⁄4 ~b

ffiffiffiffiffiffi

qi

SS

1

 q

q_ i

SS: ð94Þ

A general evolution equation for each dislocation density qi

SS is now suggested:

q_ i

SS 1⁄4 q_ iSSðfSj

SSg; fSjGNg; fpj

g; SiÞ ðj 1⁄4 1; ... ; nÞ: ð95Þ

As is clear from (95), the dislocation density rate on a particular system i generally includes

contributions from the slip rate on the j-th system, via the pi1⁄4j dependence, as well as the contributions of the slip rates of other systems, through the pi61⁄4j dependence. Likewise, the rate of qi

SS

may be influenced by the current statistically stored densities on all slip systems through the fSj

SSg

dependency [80,91]. Since interactions between SSDs and GNDs are not ruled out (e.g., trapping

or annihilation), we also allow the densities of GNDs on all slip systems to enter the evolution law

for a given system, through the fSj

GNg dependency. Finally, the dependency of q_ i

SS on

Si  Fei1

Sei

FeiT, the pull-back to bp of the force conjugate to the stretch associated with microscopic heterogeneity within the volume element, is included since the formation of intragranular

cellular structures will likely spur dislocation generation and annihilation at misoriented subgrain

boundaries [60]. Notice also that Fei will implicitly affect the density of GNDs through its presence

in the multiplicative decomposition of the total deformation gradient, Eq. (48). More specific rate

equations for the SSDs, with dislocation populations further divided into mobile and immobile

parts [4,99,100] and accounting for hardening and recovery on each system at different stages of

plastic deformation (cf. [34,39]) are also envisioned as possible implementations of (95).

With the evolution of F

p and the strain hardening prescribed by Eqs. (74)–(95), we now turn

our attention to the evolution of the Fei tensor. The following general evolution equation, referred

to configuration ~bint, is now proposed for the rate of stretching attributed to Fei

:

_

Vei

1⁄4 _VeiðPe; Sei; SeGN; eS SSÞ; ð96Þ

and a similar dependency is proposed for the rate of rotation arising from Fei

, i.e.,

fWi  _

Rei

ReiT 1⁄4 fWi

ðPe; Sei; SeGN; eS SSÞ: ð97Þ

More specific forms of (96) and (97) have not yet been developed. However, previous computational micromechanical solutions indicate that the magnitude of Vei should remain small in

comparison to the magnitudes of both the total applied strain and the volume-averaged plastic

strain in deforming ductile FCC single crystals [88]. The normalized stretch associated with Vei was

non-negligible, however, and did attain a magnitude comparable to that of the recoverable elastic

strain Cee in these calculations, as shown in Fig. 6(b). Guidance in formulating Eq. (96) can also

stem from the contribution of Vei to the stored energy of cold work, as demonstrated in the data of

452 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
Fig. 6(a) and suggested in the Helmholtz free energy function (68). We also must note that the

calculations used to generate data for Fig. 6 employed local crystal plasticity theory within

subgrains, so the influence of GNDs as defined in the present work was neglected in that study

[88].

While we currently lack any specific data for the evolution of the rotation tensor Rei

, additional

motivation for its inclusion in our crystal plasticity model is acquired by considering the decomposition of the total vorticity in the intermediate configuration ~bint:

ðLeÞskew  fW 1⁄4 ðReTR_ e

Þskew

|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}

Wee

þðReTVe1

V_ e

ReÞskew þ ð _Vei

Vei1

Þskew þ ðVeif_

Ri

ReiTVei1

Þskew

|fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

1⁄4ðeViWeieVi1Þskew

þfWp;

ð98Þ

where all skew-symmetrization is conducted with respect to the metric ~g. Assume for the moment

that plastic spin fWp

, spin due to residual stretching rate _Vei, and total spin fW are prescribed,

respectively, via the evolution equations for crystalline slip (Eqs. (74) and (80) or (91)), the

evolution equation for strain due to subdivision/heterogeneity (Eq. (96)), and the displacement

boundary conditions. Furthermore, assume that second term on the right-hand side of Eq. (98),

the elastic spin due to V_ e

, which is known through an appropriate rate form of the constitutive

relationship (63) between the applied stress and elastic strain, is negligible in comparison to the

other terms in (98) by the typical assumption of small elastic strains in engineering metals. As

mentioned previously, since Vei is expected to not exceed the same order of magnitude reached by

Ve

, neglecting the third term in (98) may also be a valid assumption. Under these assumptions, the

skew elastic spin fWe is then found as

fWe 
 fW  fWp  ðVeifWi

Vei1

Þskew: ð99Þ

Recall that Fe provides the current configuration orientations of the ‘‘average’’ slip plane normals

and slip directions in our crystal plasticity framework by virtue of Eq. (75); i.e., texture evolution

is essentially specified by the time rate of Fe (cf. [103]). With the assumption of small recoverable

elastic strains, we have Fe 
 Re

, and the time rate of elastic rotation dictates texture evolution.

Essentially, we are at liberty to completely control texture evolution through prescription of fWi

,

which in turn can be thought of as the contribution to the total spin from intragranular dislocation substructure formation (e.g., grain subdivision processes). For example, if we specify

fWi 1⁄4 ðVei1

ðfW  fWpÞVeiÞskew, then texture evolution (i.e., elastic spin) will be largely precluded

due to substructure development within the crystalline volume element. On the other hand if we

set fWi 1⁄4 0 for our evolution law (97), then crystal lattice orientations will evolve as in the classical

crystal plasticity theory. Numerical simulations [50,104] have shown that the classical theory, with

fWi 1⁄4 0, gives texture predictions that are consistently too-sharp for certain FCC polycrystals, at

least within the framework of Taylors [105] constraints. So by correlating potential evolution

Eqs. (97) with experimental texture measurements, one could conceivably formulate evolution

laws for fWi

, as speculated by Butler and McDowell [50]. Since grain subdivision and substructure

development are energetically favorable, and hence expected, for some materials even under

application of macroscopically uniform boundary conditions [7,43,98] _Vei

and/or fWi should attain

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 453
non-zero values during some point in the history of uniform deformation applied to such materials.

In closing the discussion of evolution equations, we emphasize that when completely developing Eqs. (95)–(97), which we have left quite general in this work, one should also consider the

restrictions on thermodynamic admissibility stemming from the complete dissipation inequality,

Eq. (67). While we have included terms associated with dislocation defect densities and residual

elasticity in the free energy function (68), in most metals their total cumulative contribution to the

first and second laws of thermodynamics, (54) and (55), will remain relatively small in magnitude

compared to that of the plastic dissipation due to dislocation motion (77), in agreement with

experimental measurements of the stored energy of cold working [106,107]. This observation is

fully consistent with the framework presented here (cf. Eq. (68)); the contribution of incompatibility- related tensors (i.e., Fei and defect densities) to the free energy is small when averaged over

the entire volume of the crystal(s). The relatively small magnitude of this energetic contribution is

scaled correctly by the appropriate choice of material parameters in (68). However, in our

opinion, the residual elastic energy associated with crystal defects and elastoplastic incompatibility should not be neglected, as its local density may attain much greater magnitudes in the

vicinity of local heterogeneities such as grain or phase boundaries, and its release may facilitate

void growth, fracture, shear localization, recrystallization, and/or phase transition processes in

many metals [88,108].

5. Conclusions

We have used explicit volume averaging procedures (Sections 2 and 3) to motivate a continuum

formulation of gradient crystal plasticity (Section 4). Notable features include

• A three-term multiplicative decomposition for the deformation gradient, F 1⁄4 Fe

Fei

F

p

, with Fe

associated with the average applied stress and total lattice rotation, Fei accounting for the possible development of heterogeneous dislocation substructures and residual elasticity, and Fp

accounting for the history of plastic deformation due to dislocation motion.

• Prescription of a resolved Eshelby-type stress measure as a driving force for plastic shearing

rates, in order to ensure a positive dissipation contribution from the plastic velocity gradient.

• Slip system-level strain hardening dependent upon conjugate thermodynamic forces to densities

of GNDs and SSDs.

• Allowance of Fei

, the kinematic measure of intragranular heterogeneity and possible subdivision, to influence the evolution of dislocation densities, and vice-versa.

• Allowance of a measure of residual elasticity––the stretch associated with Fei

––to account for a

fraction of the stored energy of cold work, influence the effective elastic moduli in the unloaded

configuration, and provide a residual free energy term biased by the applied stress.

Our framework offers powerful tools for addressing multiscale issues in plasticity, most notably

providing a very rigorous (i.e., mathematically precise) methodology for characterizing kinematics

and thermodynamics with changes in scale of observation. Future work should focus on determination of more precise evolution equations for Fei and the dislocation densities, determination

454 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
of the needed material parameters (e.g., hardening coefficients and length scale constants), and

exploration of boundary conditions. Micromechanical simulations of representative defect populations invoking discrete dislocation dynamics relations [109–111] are foreseen as valuable tools

for deducing evolution laws and model parameters. Extensions to other defects (e.g., disclinations,

point defects, and damage), temperature effects, and dynamic conditions are also envisioned.

Acknowledgements

J.D.C. and D.L.M. are grateful for the support of the Army Research Office (ARO Award

DAAD19-99-1-0221). D.J.B. acknowledges the support of Sandia National Laboratories. Sandia

is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for

the United States Department of Energy under Contract DE-AC04-94AL85000.

References

[1] C. Teodosiu, Rev. Roum. Sci. Tech.––Mech. Appl. 12 (1967) 961–977.

[2] C. Teodosiu, Rev. Roum. Sci. Tech.––Mech. Appl. 12 (1967) 1291–1308.

[3] C. Teodosiu, A dynamic theory of dislocations and its application to the theory of elastic–plastic continuum, in:

Proceedings of the Conference on Fundamental Aspects of Dislocation Theory, Washington, 1969, pp. 837–876.

[4] R.W. Lardner, Z. Agnew Math. Phys. 20 (1969) 514–529.

[5] O.W. Dillon, J. Kratochvil, Int. J. Solids Struct. 6 (1970) 1513–1533.

[6] O.W. Dillon, P. Perzyna, Arch. Mech. 24 (1972) 727–747.

[7] M. Ortiz, E.A. Repetto, J. Mech. Phys. Solids 47 (1999) 397–462.

[8] H.M. Zbib, E.C. Aifantis, Acta Mech. 92 (1992) 209–225.

[9] A. Menzel, P. Steinmann, J. Mech. Phys. Solids 48 (2000) 1777–1796.

[10] D.L. Holt, J. Appl. Phys. 41 (1970) 3197–3201.

[11] D.J. Bammann, E.C. Aifantis, Acta Mech. 45 (1982) 91–121.

[12] E.C. Aifantis, Int. J. Plast. 3 (1987) 211–247.

[13] M.E. Gurtin, J. Mech. Phys. Solids 48 (2000) 989–1036.

[14] A.C. Eringen, J. Appl. Phys. 56 (1984) 2675–2680.

[15] E.C. Aifantis, J. Eng. Mater. Tech. 121 (1999) 189–202.

[16] K. Shizawa, H. Zbib, Int. J. Plast. 15 (1999) 899–938.

[17] K. Shizawa, H. Zbib, J. Eng. Mater. Tech. 121 (1999) 247–253.

[18] E.O. Hall, Phys. Soc. London Proc. 64 (1951) 747–753.

[19] N.J. Petch, J. Iron, Steel Inst. 174 (1953) 25–28.

[20] J.Y. Shu, N.A. Fleck, J. Mech. Phys. Solids 47 (1999) 297–324.

[21] S. Forest, F. Barbe, G. Cailletaud, Int. J. Solids Struct. 37 (2000) 7105–7126.

[22] E. Cosserat, F. Cosserat, C.R. Acad. Sci. Paris 145 (1907) 1139–1142.

[23] E. Cosserat, F. Cosserat, Theorie des Corps Deformables, Hermann, Paris, 1909, p. vi+, 226.

[24] N.A. Fleck, G.M. Muller, M.F. Ashby, J.W. Hutchinson, Acta Metall. Mater. 42 (1994) 475–487.

[25] J.Y. Shu, N.A. Fleck, Int. J. Solids Struct. 35 (1998) 1363–1383.

[26] K.C. Hwang, H. Jiang, Y. Huang, H. Gao, N. Hu, J. Mech. Phys. Solids 50 (2002) 81–99.

[27] N.A. Fleck, J.W. Hutchinson, J. Mech. Phys. Solids 41 (1993) 1825–1857.

[28] P.M. Nagdhi, A.R. Srinivasa, Int. J. Eng. Sci. 32 (1994) 1157–1182.

[29] K.C. Le, H. Stumpf, Proc. R. Soc. London A 452 (1996) 359–371.

[30] K.C. Le, H. Stumpf, Int. J. Plast. 12 (1996) 611–627.

[31] K.C. Le, H. Stumpf, Int. J. Eng. Sci. 34 (1996) 339–358.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 455
[32] P. Steinmann, Int. J. Eng. Sci. 34 (1996) 1717–1735.

[33] A. Acharya, J.L. Bassani, J. Mech. Phys. Solids 48 (2000) 1565–1595.

[34] D.J. Bammann, Mater. Sci. Eng. A 309–310 (2001) 406–410.

[35] P. Cermelli, M.E. Gurtin, J. Mech. Phys. Solids 49 (2001) 1539–1568.

[36] M.E. Gurtin, J. Mech. Phys. Solids 50 (2002) 5–32.

[37] B. Svendsen, J. Mech. Phys. Solids 50 (2002) 1297–1329.

[38] P.M. Nagdhi, A.R. Srinivasa, Philos. Trans. R. Soc. London A 345 (1993) 425–458.

[39] R.A. Regueiro, D.J. Bammann, E.B. Marin, K. Garikipati, J. Eng. Mater. Technol. 124 (2002) 380–387.

[40] M.F. Ashby, Philos. Mag. 21 (1970) 399–424.

[41] D. Kuhlmann-Wilsdorf, N. Hansen, Scr. Met. Mater. 25 (1991) 1557–1562.

[42] B. Bay, N. Hansen, D.A. Hughes, D. Kuhlmann-Wilsdorf, Acta Metall. Mater. 40 (1992) 205–219.

[43] D.A. Hughes, Q. Liu, D.C. Chrzan, N. Hansen, Acta Mater. 45 (1997) 105–112.

[44] D. Kuhlmann-Wilsdorf, Philos. Mag. A 79 (1999) 955–1008.

[45] G.C. Butler, S.R. Stock, D.L. McDowell, Key Eng. Mater. 177–180 (2000) 165–170.

[46] N.R. Barton, P.R. Dawson, Metal. Mater. Trans. 32A (2001) 1967–1975.

[47] T. Leffers, Int. J. Plast. 17 (2001) 469–489.

[48] T. Leffers, Mater. Sci. Forum 157–162 (1994) 1815–1820.

[49] T. Leffers, Int. J. Plast. 17 (2001) 491–511.

[50] G.C. Butler, D.L. McDowell, Int. J. Plast. 14 (1998) 703–717.

[51] M.F. Horstemeyer, D.L. McDowell, Mech. Mater. 27 (1998) 145–163.

[52] M.F. Horstemeyer, D.L. McDowell, R.D. McGinty, Mod. Sim. Mater. Sci. Eng. 7 (1999) 253–273.

[53] M. Berveiller, D. Muller, J. Kratochvil, Int. J. Plast. 9 (1993) 633–652.

[54] M. Ortiz, E.A. Repetto, L. Stainier, J. Mech. Phys. Solids 48 (2000) 2077–2114.

[55] D. Hughes, N. Hansen, Mater. Sci. Technol. 7 (1991) 544–553.

[56] D. Kuhlmann-Wilsdorf, Metal. Trans. 16A (1985) 2091–2108.

[57] C. Carstensen, K. Hackl, A. Mielke, Proc. R. Soc. London A 458 (2002) 299–317.

[58] S. Takeuchi, A.S. Argon, J. Mater. Sci. 11 (1976) 1542–1566.

[59] J.C. Gibeling, W.D. Nix, Acta Metall. 28 (1980) 1743–1752.

[60] A.S. Argon, S. Takeuchi, Acta Metall. 29 (1981) 1877–1884.

[61] R. Hill, Proc. R. Soc. London A 326 (1972) 131–147.

[62] J.D. Clayton, B.M. Schroeter, D.L. McDowell, S. Graham, J. Eng. Mater. Technol. 124 (2002) 302–313.

[63] J.A. Schouten, Ricci Calculus, second ed., Springer-Verlag, Berlin, 1954.

[64] K. Kondo, On the geometrical and physical foundations of the theory of yielding, in: Proceedings of the 2nd

Japan Congress on Applied Mechanics, Japan National Committee for Theoretical and Applied Mechanics,

Science Council of Japan, Tokyo, 1953, pp. 41–47.

[65] B.A. Bilby, R. Bullough, E. Smith, Proc. R. Soc. London A 231 (1955) 263–273.

[66] E. Kroner, Arch. Rat. Mech. Anal. 4 (1960) 273–334. €

[67] G.A. Maugin, Material Inhomogeneities in Elasticity, Chapman & Hall, London, 1993.

[68] A. Arsenlis, D.M. Parks, Acta Mater. 47 (1999) 1597–1611.

[69] E. Kroner, Int. J. Solids Struct. 38 (2001) 1115–1134. €

[70] B.A. Bilby, E. Smith, Proc. R. Soc. London A 232 (1956) 481–505.

[71] S. Nemat-Nasser, Mech. Mater. 31 (1999) 493–523.

[72] R. Hill, Math. Proc. Camb. Philos. Soc. 95 (1984) 481–494.

[73] A. Acharya, J.L. Bassani, Incompatible lattice deformations and crystal plasticity, in: N. Ghoniem (Ed.), Plastic

and Fracture Instabilities in Materials, AMD 200/MD 57, ASME, NY, 1995, pp. 75–80.

[74] J.L. Bassani, J. Mech. Phys. Solids 49 (2001) 1983–1996.

[75] J. Mandel, Thermodynamics and plasticity, in: J.J. Delgado et al. (Eds.), Continuum Thermodynamics,

Macmillan, New York, 1974, pp. 283–304.

[76] K.C. Le, H. Schutte, H. Stumpf, Int. J. Plast. 14 (1998) 1109–1131.

[77] J. Kratochvil, Finite-strain theory of inelastic behavior of crystalline solids, in: A. Sawzcuk (Ed.), Foundations

of Plasticity, Warsaw, 1972, pp. 401–415.

456 J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457
[78] A. Lion, Int. J. Plast. 16 (2000) 469–494.

[79] E. Kroner, J. Math. Phys. 42 (1963) 27–37. €

[80] C. Teodosiu, F. Sidoroff, Int. J. Eng. Sci. 14 (1976) 713–723.

[81] E. Van der Giessen, F.G. Kollmann, Z. Agnew. Math. Mech. 76 (1996) 447–462.

[82] H. Stumpf, U. Hoppe, Z. Agnew Math. Mech. 77 (1997) 327–339.

[83] B.D. Coleman, W. Noll, Arch. Rat. Mech. Anal. 13 (1963) 167–178.

[84] B.D. Coleman, M.E. Gurtin, J. Chem. Phys. 17 (1967) 597–613.

[85] G.A. Maugin, Int. J. Plast. 10 (1994) 393–408.

[86] J.D. Eshelby, J. Elast. 5 (1975) 321–335.

[87] K.-Ch. Le, H. Stumpf, Acta Mech. 100 (1993) 155–170.

[88] J.D. Clayton, D.L. McDowell, Int. J. Plast. 19 (2003) 1401–1444.

[89] J.C. Simo, M. Ortiz, Comp. Meth. Appl. Mech. Eng. 49 (1985) 221–245.

[90] A. Cuitino, M. Ortiz, Eng. Computat. 9 (1992) 437–452. ~

[91] J.R. Rice, J. Mech. Phys. Solids 19 (1971) 433–455.

[92] R. Hill, K.S. Havner, J. Mech. Phys. Solids 30 (1982) 5–22.

[93] R.J. Asaro, J. Appl. Mech. 50 (1983) 921–934.

[94] J.W. Hutchinson, Proc. R. Soc. London A 348 (1976) 101–127.

[95] E. Schmid, in: Proceedings of the International Congress of Applied Mechanics, Delft, 1924, p. 342.

[96] U.F. Kocks, Metall. Trans. 1 (1970) 1121–1143.

[97] G.I. Taylor, Proc. R. Soc. London A 145 (1934) 362–415.

[98] D. Kuhlmann-Wilsdorf, Mater. Sci. Eng. A 113 (1989) 1–41.

[99] M.A. Zikry, M. Kao, J. Mech. Phys. Solids 44 (1996) 1765–1798.

[100] T. Kameda, M.A. Zikry, Scr. Mater. 38 (1998) 631–636.

[101] A. Acharya, A.J. Beaudoin, J. Mech. Phys. Solids 48 (2000) 2213–2230.

[102] D.L. McDowell, J.C. Moosbrugger, Acta Mech. 93 (1992) 73–87.

[103] R.J. Asaro, A. Needleman, Acta Metall. 33 (1985) 923–953.

[104] U.F. Kocks, M.G. Stout, A.D. Rollett, The influence of texture on strain hardening, in: Strength of Metals and

Alloys, vol. 1, 1988, pp. 25–34.

[105] G.I. Taylor, J. Inst. Metals 62 (1938) 307.

[106] G.I. Taylor, H. Quinney, Proc. R. Soc. London A 143 (1934) 307–326.

[107] P. Rosakis, A.J. Rosakis, G. Ravichandran, J. Hodowany, J. Mech. Phys. Solids 48 (2000) 581–607.

[108] S.P.A. Gill, M.G. Cornforth, A.C.F. Cocks, Int. J. Plast. 17 (2001) 669–690.

[109] H.M. Zbib, M. Rhee, J.P. Hirth, Int. J. Mech. Sci. 40 (1998) 112–127.

[110] J.L. Bassani, A. Needleman, E. Van der Giessen, Int. J. Solids Struct. 38 (2000) 833–853.

[111] A. Needleman, Acta Mater. 48 (2000) 105–124.

J.D. Clayton et al. / International Journal of Engineering Science 42 (2004) 427–457 457
NO. OF 

COPIES ORGANIZATION 

 1 DEFENSE TECHNICAL 

 (PDF INFORMATION CTR 

 ONLY) DTIC OCA 

 8725 JOHN J KINGMAN RD 

 STE 0944 

 FORT BELVOIR VA 22060-6218 

 1 US ARMY RSRCH DEV & 

 ENGRG CMD 

 SYSTEMS OF SYSTEMS 

 INTEGRATION 

 AMSRD SS T 

 6000 6TH ST STE 100 

 FORT BELVOIR VA 22060-5608 

 1 INST FOR ADVNCD TCHNLGY 

 THE UNIV OF TEXAS 

 AT AUSTIN 

 3925 W BRAKER LN 

 AUSTIN TX 78759-5316 

 1 DIRECTOR 

 US ARMY RESEARCH LAB 

 IMNE ALC IMS 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

 3 DIRECTOR 

 US ARMY RESEARCH LAB 

 AMSRD ARL CI OK TL 

 2800 POWDER MILL RD 

 ADELPHI MD 20783-1197 

ABERDEEN PROVING GROUND

 1 DIR USARL 

 AMSRD ARL CI OK TP (BLDG 4600) 
NO. OF 

COPIES ORGANIZATION 

ABERDEEN PROVING GROUND

 30 DIR USARL 

 AMSRD ARL CI HC 

 P CHUNG 

 C CORNWELL 

 AMSRD ARL WM 

 M FERMEN COKER 

 T WRIGHT 

 J MCCAULEY 

 AMSRD ARL WM MA 

 W NOTHWANG 

 AMSRD ARL WM TA 

 S SCHOENFELD 

 AMSRD ARL WM TD 

 S BILYK 

 T BJERKE 

 D CASEM 

 J CLAYTON (10 CPS) 

 D DANDEKAR 

 M GREENFIELD 

 Y HUANG 

 K IYER 

 B LOVE 

 M RAFTENBERG 

 E RAPACKI 

 M SCHEIDLER 

 S SEGLETES 

 T WEERASOORIYA 
INTENTIONALLY LEFT BLANK.