Slip Continuity in Explicit Crystal Plasticity Simulations 

Using Nonlocal Continuum and Semi-discrete Approaches

by Richard Becker

ARL-TR-6307 January 2013

Approved for public release; distribution 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

ARL-TR-6307 January 2013

Slip Continuity in Explicit Crystal Plasticity Simulations 

Using Nonlocal Continuum and Semi-discrete Approaches

Richard Becker

Weapons and Materials Research Directorate, ARL

Approved for public release; distribution unlimited. 
ii

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)

January 2013

2. REPORT TYPE 3. DATES COVERED (From - To)

October 2010 to September 2012

4. TITLE AND SUBTITLE

Slip Continuity in Explicit Crystal Plasticity Simulations Using Nonlocal 

Continuum and Semi-discrete Approaches

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Richard Becker

5d. PROJECT NUMBER

FY11-WMR-017

5e. TASK NUMBER

5f. WORK UNIT NUMBER

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

U.S. Army Research Laboratory

ATTN: RDRL-WMP-C

Aberdeen Proving Ground, MD 21005

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-TR-6307

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

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Slip continuity across element boundaries in explicit finite element simulations is enforced through nonlocal penalty 

constraints applied to continuum crystal plasticity relations and in a nonlocal, semi-discrete crystal slip model. The continuity 

constraints provide additional coupling within the plastic deformation field, and a length scale effect is introduced. The 

deformation field becomes more diffuse with reduced physical size of the model region and the strength increases. The semidiscrete method produces dislocation pile-ups and slip gradients along discrete slip planes, but the simulation approach is

ultimately unsatisfactory. Consideration of the physical spacing of dislocations and the typical size of dislocation cells in

deformed metals elicits concerns about the applicability of continuum slip models with sub-micron spatial resolution.

15. SUBJECT TERMS

Crystal plasticity, nonlocal, gradient, dislocation

16. SECURITY CLASSIFICATION OF:

17. LIMITATION

OF 

ABSTRACT

UU

18. NUMBER

OF 

PAGES

48

19a. NAME OF RESPONSIBLE PERSON

Richard Becker

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

(410) 278-7980

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
iii

Contents

List of Figures v

List of Tables vi

Acknowledgment vii

1. Introduction 1

2. Continuum Model for Slip Continuity 2

2.1 Crystal Kinematics..........................................................................................................3

2.2 Nonlocal Model...............................................................................................................3

2.3 Crystal Flow Strength......................................................................................................5

2.4 Crystal Geometry and Boundary Conditions ..................................................................6

2.5 Finite Element Implementation.......................................................................................6

3. Continuum Finite Element Results 7

3.1 Single Crystal ..................................................................................................................7

3.1.1 Nonlocal Strength................................................................................................7

3.1.2 Stress State ..........................................................................................................9

3.1.3 Slip Rate ..............................................................................................................9

3.1.4 Stress Strain Response.......................................................................................10

3.1.5 Mesh Refinement ..............................................................................................11

3.1.6 Time Step Instability .........................................................................................12

3.2 Polycrystal.....................................................................................................................13

3.2.1 Slip Rate ............................................................................................................14

3.2.2 Nonlocal Strengthening.....................................................................................14

3.2.3 Polycrystal Stress-strain Behavior.....................................................................16

4. Considerations for Discreteness of Dislocations 17

4.1 Semi-discrete Dislocation Model ..................................................................................18

4.2 Semi-discrete Simulation Results..................................................................................20

4.2.1 Low Density of Nucleation Sites.......................................................................20

4.2.2 Moderate Density of Nucleation Sites...............................................................25
iv

4.2.3 Reduced Element Size.......................................................................................26

4.2.4 Stress-strain Response from Semi-discrete Model............................................28

5. Discussion 30

6. Conclusion 32

7. References 34

Distribution List 37
v

List of Figures

Figure 1. Initial configuration for the single crystal and the crystal lattice orientation. The 

bottom is fixed and the top is moved to the right. Periodic boundary conditions are 

applied coupling the left and right hand sides. ..........................................................................7

Figure 2. Distribution of the nonlocal contribution to the strength on the slip system aligned

vertically in the crystal, for the crystal thicknesses indicated, at shear strains of 0.03 and 

0.05.............................................................................................................................................8

Figure 3. Distribution of the normalized slip rate for single crystals of the indicated thickness 

and at shear strains of 0.03 and 0.05. The slip rates are normalized by the applied shear 

rate............................................................................................................................................10

Figure 4. Shear stress-shear strain response predicted for the four crystal thicknesses. ...............11

Figure 5. Comparison of nonlocal stress and the slip rate for the 20x100 and 40x200 element 

simulations of the 5-m-thick single crystal. The results are shown at a 0.05 shear strain.....12

Figure 6. Grain structure (a) and finite element mesh and (b) for the polycrystal simulations....13

Figure 7. Normalized slip rate distribution for simple shear deformation of idealized 

polycrystals with heights of 330, 33, and 3.0 μm. The color scale is the same for all three 

plots..........................................................................................................................................14

Figure 8. Nonlocal strength contributions on the three slip systems for polycrystal model 

sizes of 330, 33, and 3.3 μm. The color scales are consistent within each row......................15

Figure 9. Shear stress strain response for three different size scales of idealized polycrystals...16

Figure 10. Change in crystal lattice orientation, in degrees, at 0.025 shear strain for the 330-

and 3.3-μm high polycrystals...................................................................................................17

Figure 11. Number of dislocations passing through each element for the discrete dislocation 

simulations of three single crystal sizes indicated. The color levels for the 5- and 10-μm

crystals are the same. The center 80 μm is omitted from the 100-μm crystal to highlight 

the gradients at the top and bottom surfaces............................................................................21

Figure 12. Number of dislocations passed through each element along the 13 active slip 

planes in figure 11 for the 10-m-high crystal plotted against the through thickness 

location. The light lines are the individual results and the heavy, dark line is the average.....22

Figure 13. Number of dislocations currently within each element for the discrete dislocation 

simulations of three single crystal sizes. The color levels for the 5- and 10-μm crystals 

are the same. The center 80 m is omitted from the 100-μm crystal to highlight the 

gradients at the top and bottom surfaces..................................................................................23

Figure 14. Plots of the in-plane stress components at a strain of 0.01 for the 10-μm-high 

single crystal deforming by discrete single slip on 13 slip planes...........................................24

Figure 15. Number of dislocations within each element for discrete dislocation simulations 

of the 5- and 10-μm single crystals with a moderate density of nucleation sites. ...................25

Figure 16. Shear stress distribution for 10-μm-thick crystal simulation with a moderate 

dislocation nucleation site density. ..........................................................................................26
vi

Figure 17. Number of dislocations currently within each element for the discrete dislocation 

simulations of three single crystal sizes. The color levels for the 1- and 2-μm crystals are 

the same. The center 16 μm is omitted from the 20-μm crystal to highlight the gradients 

at the top and bottom surfaces. ................................................................................................27

Figure 18. Shear stress distribution for the 10 μm x 2 μm crystal in a calculation with an 

element size of 0.025 μm and 172 dislocation nucleation sites. The stress gradients are 

under-resolved..........................................................................................................................28

Figure 19. Shear stress-strain response for the semi-discrete crystal simulations of section 4. 

The legends and arrows show the model width, crystal aspect ratio, element length, and 

the percent of elements containing dislocation nucleation sites along with the 

corresponding nucleation site density. Curves of the same color have the same element 

size and element site nucleation percentage but different model aspect ratios........................29

List of Tables

Table 1. Average number of dislocations per element expected in the simulations of section 

3.1 for well-annealed, lightly deformed, and heavily deformed metals. .................................18
vii

Acknowledgment

I am grateful for the support of this work through the Director’s Research Initiative at the U.S. 

Army Research Laboratory FY11-WMR-017.
viii

INTENTIONALLY LEFT BLANK.
1

1. Introduction

Traditional metal plasticity models, formulated in terms of strain rates and stresses and 

incorporated in large-scale numerical analyses, provide useful solutions for a wide range of 

problems. Details of the material microstructure interactions that govern the deformation

response are assumed to occur at length scales not resolvable by the simulations and are captured 

implicitly in the constitutive relations. For example, dependence of the yield strength on grain 

size through the Hall-Petch effect can be incorporated by including grain size in the constitutive 

model without tracking the details of dislocation interactions with grain boundaries. 

In simulations with spatial resolution at or below the micron level, as in multiscale modeling, the 

length scales dictating some hardening mechanisms are correlated with gradients in the plastic 

strain field. The torsion experiments of Fleck et al. (1994) clearly demonstrate increased strength 

with decreasing size for wires 10’s of microns in diameter. The size effect was further observed 

in bending (Stölken and Evans, 1998) and indentation (Saha, et al. 2001), and in many 

subsequent studies. The strengthening is attributed to gradients in the crystal lattice orientation. 

These gradient microstructures both store energy and provide resistance to further dislocation 

motion (Lee et al., 1989; Fleck and Hutchinson, 1997; Gao et al., 1999; Baskaran et al., 2010; 

Schouwenaars et al., 2010). Applying traditional crystal plasticity models (Asaro, 1983; Peirce et 

al., 1983) to investigate the size scale effect in a multiscale framework will not be successful 

because the models are formulated in terms of traditional continuum variables of strain rate and 

stress, and there is no underlying microstructure length scale that would produce a size effect. 

Numerous studies over more than a decade have investigated ways to incorporate a length scale 

into continuum crystal plasticity models. Most focus on microstructure gradients, as it is 

recognized that both the Hall-Petch effect and the size scale results are tied to gradients. To cite 

just a few of the many examples, models have examined: continuously distributed dislocations 

(Acharya, 2001), dislocation density gradients (Arsenlis et al., 2004), gradient related state 

variables (Gurtin et al., 2007; Gerken and Dawson, 2008), and micro-polar theories (Mayeur et 

al., 2011). Many formulations introduce additional variables into the solution, which must be 

determined in concert with the standard degrees of freedom.

When using the traditional crystal plasticity model in a finite element code, the deformation of 

neighboring elements is only connected through shared nodal displacements and force 

equilibrium at the nodes. Even though dislocations associated with slip travel through the 

material, the transmission of dislocations from one element to another is not represented in the 

model, and the impact of coordinated slip on the deformation field is not captured. The indirect

result is that all finite element boundaries are infinite sources and sinks for dislocations. 

Advanced crystal plasticity formulations, such as those based on lattice orientation gradients, 
2

often include the continuity of dislocation flux across finite element boundaries as a byproduct

because lattice orientation gradients are created by an accumulation of like-signed dislocations.

The goal of this work is to determine if a simpler approach can be used to address the length 

scale dependence and other deficiencies of the classical crystal plasticity model. The starting 

point is simply enforcing dislocation flux continuity across finite element boundaries. Different 

slip conservation approaches were examined by Hirschberger et al. (2011). A choice is made to 

implement the model in an explicit finite element code using an operator-split approach, where it 

is assumed that the time steps are sufficiently small that the coupled mechanics can be applied 

consecutively over a time step rather than concurrently. In this approach, the slip in the elements 

can be treated as a normal state variable, whereas it would be most effective to treat it as a degree 

of freedom in fully coupled, implicit solutions. The code uses existing solution variables, and the 

boundary conditions are conceptually straightforward. Dislocation flux is unconstrained at free 

surfaces, zero at rigid boundaries, and intermediate at grain boundaries that are sources and sinks 

for dislocations. The flux gradients can be used to infer the evolution of dislocation gradients that 

lead to lattice orientation gradients.

Initial work to explore the feasibility and potential impact of enforcing dislocation flux between 

elements in a crystal plasticity model used a special purpose, explicit dynamic code that 

calculated and enforced an average slip rate on element boundaries in a multi-step constitutive 

evaluation (Becker, 2011). The results showed the desired effects of smoothing the deformation 

field and increasing the hardening rate with smaller sample sizes. However, the solution 

technique permitted spatially oscillatory fields, particularly near the boundaries, so an alternative 

method for enforcing slip continuity between elements was indicated.

The report begins with an outline of the nonlocal, gradient crystal plasticity model and its 

implementation, and results are presented for single crystal and polycrystal simulations. 

Consideration of details of the microstructure at micron size scales prompted exploration and 

evaluation of the semi-discrete approach. The discussion expounds on concerns about the 

adequacy of continuum crystal plasticity models to represent microstructures at sub-micron 

resolution.

2. Continuum Model for Slip Continuity

The goal is to develop an enhancement to traditional continuum crystal plasticity models that can 

be readily incorporated into existing finite element implementations and run with minimal 

additional computational overhead. This objective is facilitated by targeting explicit dynamic 

solutions where, because of the relatively small time steps, it is often possible to impose 

additional constraints through an operator-split approach in which the physics are applied 
3

sequentially within a time step rather than concurrently. In the current context, the operator split 

results in the constraint application lagging one time step behind its evaluation.

2.1 Crystal Kinematics

The crystal plasticity model and implementation used in the explicit dynamic finite element 

simulations is described in Becker (2004). It follows from the widely used kinematic framework 

given by Asaro (1983) with modifications to incorporate an equation of state for high pressure 

applications. The deformation gradient, , is notionally decomposed into elastic and plastic 

parts, 

 

and 

 

, respectively;

 

 

 

 

 

The elastic part accounts for distortion and rotation of the crystal lattice, and the plastic part 

represents slip on predefined crystal planes and directions that moves material but does not alter 

the underlying crystal lattice. 

The velocity gradient obtained from this kinematic description comprises an additive 

decomposition of an elastic part and an inelastic part associated slip on predefined slip systems:

 

 

 

 

 

 

 

 

  

 

 

 

 

 

and are, respectively, vectors of the current slip directions and slip plane normals 

associated with each of the slip systems (superscript ), and 

 

are the corresponding slip rates. 

The slip direction and normal vectors rotate and distort with the crystal lattice as the material 

deforms.

2.2 Nonlocal Model

In conventional crystal plasticity analyses using an explicit dynamic finite element code with 

uniform strain elements, the slip rates are evaluated at each time step using the velocity gradient,

the crystal strength, and any history variables associated with the element. There is no direct 

connection among the slip rates in neighboring elements. In order to enforce slip continuity 

between elements, a penalty approach is proposed. Differences in accumulated slip between 

adjacent elements will either increase or decrease the resistance to further slip on that slip 

system. This requires information from neighboring elements, making it a nonlocal method.

The basis for the model is that dislocations move along the slip direction from one element to 

another. The flux of dislocations crossing an element boundary for each slip system is

 

 

 

  

 

where 

 

is the dislocation density, 

 

is the dislocation velocity, and is the outward normal 

to the given element face. Continuity is enforced by requiring that the dislocation flux exiting 
4

through an element face equals the flux entering the adjacent element through the common face. 

The dot product in equation 3 accounts for the orientation of the slip system with respect to the 

element face. It is zero if the slip direction is parallel to the face. This provides the opportunity 

for sharp jumps in slip rates across parallel slip planes while enforcing continuity along slip 

planes.

The dislocation density, velocity, and the Burgers vector, b, are related to the continuum slip rate 

by Orowan’s equation

 

 

 

 

 

For a shared element face, and assuming that the Burgers vector is constant and that the slip 

directions are closely aligned across the interface, the continuity error in the accumulated 

dislocation flux between elements can be approximated as

 

  

 

 

 

 

 

The subscripts e and neighbor on the accumulate slip, 

 

, denote, respectively, the element of 

interest and the neighbor sharing the face. The approach taken is to penalize the error. It is not 

necessary to assume that the Burgers vectors are equal or that the slip directions in the adjacent 

elements are aligned, but these simplifications are made for expedience and coding clarity in this 

initial implementation.

Application of equation 5 in a penalty method requires a cumulative result over all element faces. 

This is complicated by the changing algebraic sign of the dot product for different element faces 

and the 50% probability that the accumulated slips are negative. The sign of the dot product is 

irrelevant to the physical problem, so the absolute value is taken. The algebraic sign of the slip is 

also unimportant, and the magnitude is obtained by multiplying by the sign of the slip. This is 

preferred to using the absolute values of the slips in cases where the sign of the slip is different in 

the neighboring element. With these adjustments, the driving force for the penalty method, 

equation 5, is summed as

 

 

 

 

 

 

 ign 

   

 

 

 

 

 

 

 

 

 is the distance between element centroids associated with the particular element face. It is

intended to provide a larger penalty if the slip difference occurs on a smaller spatial grid. 

 

in 

equation 6 is a weighted sum of unsigned slip gradients on a slip system for a given element. If 

the value is positive, slip is deficient in the element and additional slip is promoted. If the value 

is negative, further slip is impeded. It is significant that a constant gradient results in face 

contributions that sum to zero. Hence, a constant gradient is not suppressed. The excess 

dislocations associated with the gradient are assumed to be uniformly distributed.
5

A closer look at equation 6 for the special case of a slip direction aligned with a regular 

rectangular mesh reveals a serendipitous connection to published gradient theories. Assuming 

that the signs of the slips are the same and performing the summation,

 

 

 

 

 

 

 

 

Subscripts e–1 and e+1 denote elements on either side of element e. Dividing equation 7 by , 

this becomes the finite difference form for half of the second derivative of the slip along the slip 

direction. The second derivative appears in the micro force balance in several gradient 

formulations (e.g., Bittencourt et al., 2003; Gurtin et al., 2007). Hence, with a small coding 

change, the implementation can approximate a simplified version of an established gradient 

model. Because of this connection to established models, and since it gives a nonlinear length 

scale dependence, the modified relation

 

 

 

 

 

 

 

 ign 

   

 

 

 

 

 

  

is used for the nonlocal simulations.

2.3 Crystal Flow Strength

The rate dependent crystal constitutive model relates the loading on each slip system to the slip 

rate. The driving force for slip, 

 

, is the resolved shear stress on the slip system.

 

 

 

 

 

where is the Cauchy stress. A simple power law rate model is used here to determine the slip 

rate for the applied loading.

 

 

 

 

 

 

 

 

 

 

 

A low rate exponent, , is chosen to provide nearly rate-independent behavior while 

still smoothing the transition from elastic to plastic response at the slip system flow strength, 

 

. 

The flow strength is adjusted, based on the nonlocal contribution, to increase or reduce the 

resistance to continued deformation

 

 

 

 

 

 

 

Here, and are dimensionless parameters, , is some average shear modulus for the crystal, 

and is the Burgers vector. The base crystal strength, 

 

, could be a function of slip to capture 

strain hardening, but it is assumed constant in these analyses to simplify interpretation of the 

results. The hyperbolic tangent function is used to cap the influence of the nonlocal term. The 

contribution is approximately linear, 

 

, while the argument is small, and it is capped at 
6

a constant value of . The penalty parameters are chosen as and 

 . In all of the nonlocal continuum simulations: the base flow strength is 

 MPa;

the density is 2.7 g/cm3; the Burgers vector is nm; and shear modulus is 

 MPa.

2.4 Crystal Geometry and Boundary Conditions

An idealized two-dimensional crystal geometry is used for these analyses. The crystal consists of 

three slip systems set in an equilateral triangle configuration. This allows a multitude of slip 

modes, but redundant slip is unlikely. In the rate independent limit, it is not possible to have slip 

on all three systems simultaneously, and the rate sensitivity is too low to enable redundant slip 

activation in these simulations.

Constant strain quadrilateral elements with hourglass control (Flanagan and Belytschko, 1981) 

are used for all of the simulations. Elements in which all four faces are adjacent to an element of 

the same initial orientation are interior elements, and gradients are computed directly as indicated 

in equation 8. Elements with fewer than four faces contacting regions of the same orientation are 

either on grain boundaries or model boundaries. The face is flagged for these elements, and a 

parameter is checked to see whether it is treated as a zero flux boundary or a free boundary with 

no slip impedance. For the non-interior elements with zero-flux boundaries, ghost elements with 

opposite slip are assumed across the flagged faces. Grain boundaries and surfaces with applied 

boundary conditions are treated in this manner. For free surfaces and periodic boundaries, the 

ghost element across the flagged face is set with the same slip so that these interfaces do not 

contribute to the gradient. 

2.5 Finite Element Implementation

The model was implemented in the large-scale parallel, explicit finite element code ALE3D 

(2012). The crystal plasticity constitutive model existed previously (Becker, 2004) and the 

strengthening terms in equation 11 due to the nonlocal effects were a straightforward addition. 

The nonlocal computations occur outside of the material model when all history variables are at 

a consistent state. Since it is an Arbitrary Lagrange-Eulerian code that moves material through 

the mesh, many features and functions are already in place for the nonlocal calculations. Lists of 

neighboring elements and shared faces exist, as do functions accessing history variables in the 

adjacent elements. For parallel computations, the problem is subdivided into domains that reside 

on separate processors. Elements that are at the boundaries of these domains have neighboring 

elements that reside on other processors. The information from the neighboring elements 

residing on other processors is carried locally in ghost elements, so all of the data needed to 

perform the nonlocal calculations of equation 8 are available locally on each processor. A 

communication call prior to the nonlocal calculations assures that all of the information is 

current. A significant amount of bookkeeping and data juggling is required, but the computations 

in equation 8 are straightforward.
7

3. Continuum Finite Element Results

The effect of the slip continuity constraint (slip gradient) is evaluated on two configurations, 

each at multiple length scales. All simulations are two-dimensional. The first configuration is 

simple shear of a single crystal with one of the slip planes initially aligned orthogonal to the 

shear direction. This creates single slip conditions for a straightforward evaluation of the model. 

The second configuration is a polycrystal constructed from regular hexagons. The orientation of 

the crystal lattice for each grain is random. 

3.1 Single Crystal

Single crystal calculations were run at four size scales using a 20 x 100 mesh of square elements

(figure 1). Velocity boundary conditions are applied to the upper and lower surfaces to shear the 

top of the crystal to the right. Initial velocities of all interior nodes are prescribed consistent with 

simple shear to eliminate ringing as the explicit dynamic calculation starts. Periodic boundary 

conditions are applied to the lateral surfaces to mimic an infinitely wide crystal. The heights of 

the single crystals simulated were 50, 5, 1, and 0.5 μm, and the width of the simulation box was 

20% of the height in each case. Although the width is irrelevant with the periodic boundary 

conditions, multiple elements are used across the width to demonstrate that the boundary 

conditions are applied properly. Slip transmission is restrained on the upper and lower 

boundaries, and slip transmission is unimpeded on the lateral, periodic boundaries.

Figure 1. Initial configuration for the single crystal and the crystal lattice orientation. 

The bottom is fixed and the top is moved to the right. Periodic boundary 

conditions are applied coupling the left and right hand sides.

3.1.1 Nonlocal Strength 

The nonlocal strength contribution to equation 11 is shown in figure 2 for the four crystal sizes 

and at two shear strains. The largest crystal is on the left and the smallest is on the right. The top 

row shows the distribution at a shear strain of 0.03 and bottom shows the results at an average 

Lattice 

orientation

Periodic boundary 

conditions

Velocity boundary 

condition
8

shear strain of 0.05. The color scales are different for each crystal size, but the scales at each 

crystal size are the same at 0.03 and 0.05 shear strains to highlight the evolution.

Figure 2. Distribution of the nonlocal contribution to the strength on the slip system aligned vertically in the 

crystal, for the crystal thicknesses indicated, at shear strains of 0.03 and 0.05.

The magnitude of the nonlocal strength contribution and the relative depth that the distribution 

penetrates from the surfaces increases as the crystal thickness is decreased. The sharpest 

gradients are expected at the crystal surfaces where the slip transmission is impeded. For the 

50-μm crystal the strength is increased only within a few elements of the surface, corresponding 

to a few microns. The central portion of the crystal sees no gradient or strengthening effect, even 

as the strain increases from 0.03 to 0.05. The boundary layer also appears to penetrate a few 

microns in the 5-μm-thick crystal simulation. However, a smaller portion of the crystal is nearly 

free of gradient effects for this smaller crystal. At yet smaller crystal thicknesses, the gradient 

effect penetrates the full crystal thickness, and the evolution with increasing deformation is 

0.03 shear 

strain

0.05 shear 

strain

50 m thick 5 m thick 1 m thick 0.5 mm thick
9

evident. The strength at the center is elevated significantly by a strain of 0.03, and it continues to 

increase with further deformation. 

The minimum and maximum values are indicated in each of the plots. Except for the smallest 

crystal, the difference between the maximum and minimum increases as the crystal size

decreases and the difference also increases with increasing strain. The increase with strain 

indicates that the gradient is still evolving at a shear strain of 0.05. The trends are different for 

the 0.5-μm crystal because the gradient strengthening is beginning to run up against the cap set 

by the hyperbolic tangent function in equation 11. With continued deformation, the gradient 

strengthening is becoming more uniform, albeit at a higher level. 

3.1.2 Stress State

The momentum balance in the vertical direction (y-direction) dictates that the y-gradient in the ydirection stress component is balanced by the horizontal (x-direction) gradient in the shear stress. 

Since the periodic boundary conditions require that all horizontal gradients are zero, the stress in 

the y-direction should be constant. The magnitude of the y-direction stress is not determined by 

the momentum equation, just that it is constant. The calculations show a constant y-direction 

stress to five significant digits. 

More important for current purposes is the momentum balance in the horizontal direction. The 

lack of stress gradients in the horizontal direction requires the shear stress to be constant through 

the thickness. The simulations show that the shear stress is constant to five significant digits. 

While the x-direction stress must be constant in the x-direction, the symmetry conditions and

momentum equations provide no further constraints restricting its gradient in the y-direction.

3.1.3 Slip Rate

The normalized slip rates corresponding to the configurations in figure 2 are shown in figure 3. 

Again, the largest crystal is on the left and the smallest is on the right. The top row contains 

results at a shear strain of 0.03, and bottom row shows the normalized slip rate at a shear strain 

of 0.05. The plots are normalized by the applied shear strain rate so that a value of 1.0 would 

indicate a uniform shear. Since the shear stress and the reference strength in equation 11 are both

constant, the slip rate is approximately related to the gradient term through equations 10 and 11. 

Second-order factors influencing the slip rate include the change in lattice orientation and nonzero components of the x-direction and y-direction stresses that modify the resolved shear stress. 
10

Figure 3. Distribution of the normalized slip rate for single crystals of the indicated thickness and at shear strains 

of 0.03 and 0.05. The slip rates are normalized by the applied shear rate.

The slip rate follows the same trends as the nonlocal hardening contribution. The variation in slip 

rate is greater for the smaller crystals; and, except for the smallest crystal, the differences are 

greater with increased deformation. As a result of the boundary conditions, slip rates are low at 

top and bottom boundaries compared to the center regions. This results in the sigmoidal 

deformed shapes. More severe differences in slip rate result in greater deviation from a linear 

shear deformation profile. The slip rates for the 0.5-μm crystal become more uniform at the 

higher deformation because the gradient term is capped by the hyperbolic tangent function. The 

strength is more uniform, which results in a more uniform slip rate. 

3.1.4 Stress Strain Response

The shear stress-shear strain responses for the various crystal thicknesses are plotted in figure 4. 

The curves are identical through the linear elastic regime, and all yield at approximately the same 

0.03 shear 

strain

0.05 shear 

strain

50 m thick 5 m thick 1 m thick 0.5 mm thick
11

stress, approximately 34 MPa. The applied shear rate is 50 times the reference shear rate; and, 

accounting for the strain rate sensitivity, the apparent yield strength is calculated from equation

10 to be 33.99 MPa rather than the reference shear strength of 33.33 MPa. 

Figure 4. Shear stress-shear strain response predicted for the four crystal thicknesses.

The stress shows the expected trend of increasing strength at the smaller crystal sizes. This is 

directly related to the nonlocal strength in figure 2. The nonlinear dependence on specimen 

dimensions is due to dividing by the element size squared in equation 8. The strain hardening 

rate is fairly consistent with increasing strain for the three larger crystal sizes but not for the 

smallest. Stress in the 0.5-μm crystal peaks as the hyperbolic tangent function places a cap on the 

nonlocal strength contribution. This is also consistent with the results presented in figures 2 and 

3. A final observation from figure 4 is the kink that is most evident in the larger two specimens 

near a strain of 0.02. This marks the transition from single slip at lower strains to slip on two slip 

systems at larger strains. As the crystal lattice rotates and stresses build in the x and y directions, 

the crystals are able to accommodate the deformation more easily with multiple active systems. 

Due to the angle of the slip plane, slip constraints at the boundaries are not as severe for the 

second slip system, so the strain hardening rate is reduced. 

3.1.5 Mesh Refinement

The effect of mesh resolution on the solution is investigated by rerunning the 5-μm-thick 

simulation using twice as many elements in each direction. The results from the 40x200 mesh are 

shown along side of the 20x100 mesh results in figure 5. Other than the expected differences in 

smoothness of the fields, the nonlocal stress and slip rate distributions do not appear to be 

influenced significantly by halving the mesh size. 
12

Figure 5. Comparison of nonlocal stress and the slip rate for the 20x100 and 40x200 element simulations of 

the 5-m-thick single crystal. The results are shown at a 0.05 shear strain.

For a more quantitative assessment, the applied shear stress computed at the crystal boundaries is 

2.3% lower for the finer mesh. Part of the difference may be due to the lower energy solution 

expected with an increased number of degrees of freedom, and the remainder can be attributed to 

discretization error associated with the nonlocal computations and the operator split algorithm.

The time step was also a factor of two lower in the fine mesh calculation due to the dependence 

of the Courant stable time step on the mesh size. The smaller time step improves the accuracy of 

the operator split integration.

3.1.6 Time Step Instability

The simulation of the 0.5-μm-thick crystal experienced numerical instabilities when the strain 

increment per step was too great. This is thought to be associated with the operator split where 

the strength increase from the slip gradients creates a driving force that is too large and overcorrects the slip rate. Specifically, with a strain increment of 1.54x10–7

per time step, the slip rate

for the 0.5-μm crystal was erratic and non-zero only in scattered, isolated elements. These 

isolated regions of slip occurred briefly and died out quickly as deformation proceeds, and 

eventually strain was incremented in the entire domain, albeit unevenly. When the time step was 

reduced by a factor of two, such that the strain increment per step was 7.43x10–8, the calculation 

was well behaved. The results in figures 2 through 4 were run with a strain increment of 

3.853x10–8to be certain that the time step was small enough to suppress the instability. The 

calculations for the larger crystal sizes experienced less gradient hardening, and they were run at 

the Courant stable time step without any additional time step controls.

Nonlocal Stress Normalized Slip Rate

20x100 40x200 20x100 40x200
13

3.2 Polycrystal

An idealized polycrystal was created by filling a rectangular region with regular hexagonal

grains (figure 6a). Simple shear boundary conditions were applied by prescribing velocities to 

the nodes on the upper and lower surfaces. Periodic boundary conditions were applied on the left 

and right surfaces. The orientation of the triangular crystal lattice in each of the grains was 

random, and the rotation angle from the reference orientation is indicated on the plot. The half 

grains at the same height on the left and right of the model region were given the same 

orientation to facilitate application of periodic boundary conditions. As with the single crystal 

simulations, the initial velocity of all interior nodes was prescribed to eliminate ringing from 

abrupt imposition of boundary conditions. 

The default inter-element slip rate condition for all elements is that any element face touching 

another grain will have restricted slip. This is imposed by assuming that a ghost element across 

the interface has equal and opposite slip in equation 8. The restricted slip condition is enforced 

on the upper and lower surfaces and on grain boundaries, including those grain boundaries on the 

periodic surfaces. The half crystals on the periodic surfaces are treated differently; the element 

across the interfaces is assumed to have the same slip. This is not a truly periodic condition, but a 

data structure identifying periodic neighboring elements is not yet available. 

Three model sizes are investigated: , , and

 . All use the same mesh configuration, scaled to give the appropriate 

dimensions. Each of the 88 hexagons was discretized by an identical mesh of 2112 quadrilateral 

elements (figure 6b). Nodes are shared along the grain boundaries, so the deformation is 

continuous throughout. 

Figure 6. Grain structure (a) and finite element mesh and (b) for the polycrystal simulations.

a) b)
14

3.2.1 Slip Rate

Slip rates normalized by the applied shear rate are shown in figure 7 for the three polycrystal 

sizes. At the larger crystal size, the strain rate localizes into well-defined bands. The majority of 

the deformation is carried by two horizontal bands with some scattered activity in the central 

region of the model. Blocks of grains appear to remain elastic while localized shear along the 

horizontal and vertical bands accommodates the deformation between neighboring blocks. The 

bands do not follow the grain boundaries, but many are associated with grain boundary triple 

points.

At the intermediate polycrystal size, the nonlocal slip constraint diffuses the deformation. The 

slip bands are still fairly well defined, but the peak strain rates are not as high and regions of 

nearly elastic behavior are smaller and less well defined. The slip rates in the 3.3-μm-thick

polycrystal are considerably more diffuse and the material near the highly constrained top and 

bottom boundaries has the lowest strain rates. Grain outlines are evident as the slip rate tends to 

be high or low at the grain boundaries, and the color contrast across the boundaries accentuates 

them.

Figure 7. Normalized slip rate distribution for simple shear deformation of idealized polycrystals with heights of 

330, 33, and 3.0 μm. The color scale is the same for all three plots.

3.2.2 Nonlocal Strengthening

The nonlocal strengthening associated with the gradients is shown in figure 8 for all three slip 

systems and the three crystal sizes. The color scale in each row is the same so that the magnitude 

of the effect of the slip systems can be compared. The scales are different for each crystal size as 

the strengthening is much greater in the smaller model region. The color scale on for the 330-mm 

polycrystal is set to a relatively low value of 0.5 MPa, and even then, the gradient contribution is 

only evident at the grain boundaries or near the most highly shear regions shown in figure 7. The 

impact on the solution is minor.

330 m high 33 m high 3.3 m high
15

Figure 8. Nonlocal strength contributions on the three slip systems for polycrystal model sizes of 330, 33, and 

3.3 μm. The color scales are consistent within each row.

The nonlocal strength distribution in the 33-μm polycrystal appears in the grain interiors as well 

as at grain boundaries. The strong interior features are associated with stress concentrations at 

grain boundary triple points, and most correspond to elevated slip activity in figure 7. Many of 

the grain boundaries show strengthening on one side and softening on the other. These 

Slip System 1 Slip System 2 Slip System 3

330 m 

crystal

33 m 

crystal

3.3 m 

crystal
16

correspond to increasing slip when approaching grain boundaries for strengthening and 

decreasing slip when approaching the boundaries for softening.

The nonlocal strengthening is quite prominent in the 3.3-μm polycrystal. As with the 

intermediate size polycrystal, the sign of the gradient effect is often flipped across the grain 

boundaries. For most grains, the gradient is strongest at the grain boundaries and decays toward 

the grain center. However, there are a few notable grains where the peak values are in the 

interiors. These correspond to locations of intersecting slip activity in figure 7. The hyperbolic 

tangent function causes the gradient effect to saturate at a level between 12 and 13 MPa. It is also 

notable that the gradient strengthening occurs on only two of the three slip systems. This reflects 

the lack of redundant slip for the idealized crystal. Only two slip systems are active at any time.

3.2.3 Polycrystal Stress-strain Behavior

The shear stress-shear strain response for the three polycrystal sizes is presented in figure 9. As 

with the single crystals, the stress is higher for the smaller polycrystals. The nonlinearity with 

length scale is also clearly evident. There are, however, two important distinctions from the 

single crystal results. The first is that the initial yield point varies with crystal size, whereas it did 

not for the single crystal simulations. This is thought to be related to the single crystals yielding

throughout simultaneously, while the polycrystal yields gradually and may build up local slip 

gradients before the macroscopic yield is evident.

Figure 9. Shear stress strain response for three different size scales 

of idealized polycrystals.

The other notable difference is that the curves are not smooth. This could result from a 

combination of the evolution of the crystal lattice orientation, evolution of the slip gradients, and 

wave propagation in the explicit dynamic calculation. The change in lattice orientation is shown 

in figure 10 for the largest and smallest size scales. In the 330-μm polycrystal, where the strain 

localization is more pronounced, lattice reorientation is also localized. The local geometric 
17

softening facilitates the shear. In contrast, due to the slip continuity and gradient constraints, the 

lattice reorientation in the 3.3-μm polycrystal is smoothed over a larger region relative to the 

grain size, and the lattice within the grains rotates nearly uniformly. A larger portion of the 

polycrystal had to realign, which takes longer and results in a greater load excursion before it 

settles into a nearly steady shear mode.

Figure 10. Change in crystal lattice orientation, in degrees, at 0.025 shear strain for the 330- and 3.3-μm high 

polycrystals.

4. Considerations for Discreteness of Dislocations

The motivation for imposing dislocation flux constraints is that dislocations are discrete entities 

that propagate from one element to the next as part of the slip process. Another aspect of the 

dislocation discreteness is their spacing, which is typically quantified in terms of the dislocation 

density. For well-annealed metals, a typical dislocation density is 106– 107cm

–2

; at a few 

percent deformation, this increases to 108– 109cm–2; and for a very heavily deformed 

polycrystal, the dislocation density is in the neighborhood of 1011cm–2(Hull and Bacon 1984). 

Table 1 lists the element areas for each of the four simulations in section 3.1 and the higher 

values of dislocation density for well-annealed, lightly deformed, and heavily deformed 

polycrystalline metals. From these values, an average number of dislocations enclosed by an

element is calculated. 

330 μm 3.3 μm
18

Table 1. Average number of dislocations per element expected in the simulations of section 3.1 for well-annealed, 

lightly deformed, and heavily deformed metals.

50 μm Crystal

Elength = 0.5 μm

Earea = 0.25 μm2

5 μm Crystal

Elength = 0.05 μm

Earea = 0.0025 μm2

1 μm Crystal

Elength = 0.01 μm

Earea = 0.0001 μm2

0.5 μm Crystal

Elength = 0.005 μm

Earea = 0.000025 μm2

107cm–2(10–1μm–2) 1/40 1/4000 10–52.5x10–6

109cm–2(101μm–2) 2.5 1/40 1/1000 1/4000

1011 cm–2(103μm–2) 250 2.5 1/10 1/40

The entries in table 1 that are less than 1.0 indicate that not every element will contain a 

dislocation. For example, 1/4000 means that one out of every 4000 elements can be expected to 

contain a dislocation. An implicit assumption in continuum crystal plasticity models is that the

dislocation content in the elements is sufficient for slip to be smooth and continuous. It is clear 

that these conditions are not met for a well-annealed metal using any of the meshes in section 

3.1, since not every element would contain even one dislocation. Using the discretization 

provided by the 50-μm-thick crystal simulation, a sufficient number of dislocations would be 

represented within each element when the crystal is heavily deformed, but not in the deformation 

leading up to that state. Considering that dislocations are usually not uniformly distributed, 

element sizes of a few microns may be necessary to assure a sufficient number of dislocations 

per element for a proper continuum crystal model representation. 

Several aspects of crystal deformation cannot be represented accurately by a continuum model if 

elements are sparsely populated by dislocations. One difficulty is that only elements containing 

dislocations can slip, and other elements must either deform elastically or nucleate additional 

dislocations. A second feature not well represented is the discreteness of slip. A dislocation that 

passes through an element creates a slip increment . For a Burgers vector length 

of b = 0.3x10–3μm and an element size of Elength = 0.5 μm, the slip must occur in increments of

0.6 x10–3rather than in the arbitrarily small increments permitted by the continuum 

representation. The discrete strain increment as a dislocation moves from one element to the next 

causes a commensurate increment in stress. This stress jump can be comparable to the yield 

strength. A final aspect of crystal plasticity, which cannot be represented properly if the element 

size is less than the dislocation spacing, is dislocation interactions leading to strengthening. 

Hardening and gradient effects are a result of interactions among dislocations. If the dislocations 

are sparse in the mesh, these interactions must be accounted for explicitly rather than implicitly 

in a hardening function.

4.1 Semi-discrete Dislocation Model

In an attempt to push continuum finite element simulations to smaller length scales where 

dislocations are sparse within the elements, a semi-discrete model was developed. It is run within 

a standard explicit-dynamic finite element framework that is described in Becker (2011). The 
19

single slip constitutive model follows the traditional formulation with three major modifications: 

(1) only elements that contain dislocations or dislocation sources can slip; (2) the slip increment 

is quantized in terms of the Burgers vector and element size; and (3) elements designated to 

contain dislocation nucleation sources have a reduced flow strength. In addition to these slip 

model modifications, the code tracks dislocations moving from one element to another, and it 

also tracks the total number of dislocations that have traversed each element.

The quantization of the slip increment is created by starting with the standard continuum slip 

increment, , and using it to create an integer representing the number of Burgers 

vector lengths the dislocation traversed during increment. 

 

 

 

 

 

 

The square of the dimensionless size scale loosely accounts for the step-by-step dislocation 

progression across the element. The slip increment used in the stress update is then calculated by

 

 

 

 

 

 

Similar integer arithmetic is used to determine when a dislocation has completely traversed an 

element and will be transferred to a neighbor. 

The crystal is assumed to generate dislocations dipoles at fixed nucleation sites. These 

dislocation nucleation sites are assigned randomly at the start of the simulation to a small 

fraction of the elements using a random number generator. The yield strength in these elements 

is also set to vary randomly between 25% and 75% of the yield stress. If the resolved shear stress 

in these elements exceeds the local yield strength, the crystal slips by moving dislocations from a 

dipole in opposite directions along the slip plane. Eventually these dislocations move to the 

neighboring elements. At this time the dislocations are transferred to the neighboring elements,

and a new dislocation dipole initiates at the nucleation site. Elements with dislocations, but not 

nucleation sites, slip when the resolved shear stress exceeds the full 50-MPa yield strength. 

Elements without dislocations deform elastically. 

The semi-discrete dislocation model is assessed through simulations using square elements with 

sizes 1.25, 0.125, and 0.025 μm. The element size is important as it is anticipated that the 

solutions will be mesh-size dependent. Calculations are run for different crystal sizes and several 

nucleation densities. In all instances, the simulations are simple shear with the single slip system 

orthogonal to the direction of shear. Periodic boundary conditions are applied as they were in

section 3.1, and dislocation transmission is prohibited on the upper and lower surfaces where the 

shear is applied. The yield strength is 50 MPa, the shear modulus 30 GPa, the bulk modulus 

60 GPa, and the Burgers vector 0.25x10–3μm. 
20

4.2 Semi-discrete Simulation Results

4.2.1 Low Density of Nucleation Sites

Simple shear simulations of well-annealed crystals with a low density of dislocation nuclei were 

run for model regions 50 μm wide and three model heights. The mesh was 400 elements wide. 

The region heights were 5, 10, and 100 μm, and the numbers of elements in the thickness 

direction were, respectively, 40, 80, and 800. This gives an element size of 0.125 μm is both x

and y directions for all three models. The number of dislocation nucleation sites corresponding to 

these meshes is 7, 14, and 158, respectively. This provides roughly the same number of 

dislocations nucleation sites per unit area with a site density of approximately 3.2x106cm–2. This

is in the range of the dislocation density for a well-annealed metal.

The number of dislocations that passed through the elements at a shear strain of 0.01 is presented 

in figure 11. The complete simulation regions are shown for the 5- and 10-μm-high models. Only 

the top and bottom 10 μm of the 100-μm model are shown, since most of the center section 

appears as lines connecting the upper and lower portions. The most obvious feature is the 

discrete deformation. Slip occurs only along slip systems containing the nucleation sites. The slip 

traverses the crystal vertically along lines of elements that contains the slip planes. The 

remaining elements are elastic. Each of the 7 nucleation sites produced slip in the 5-μm crystal,

and 13 of the 14 nucleation sites in the 10-μm model are evident. One of the nucleation sites in 

this crystal is near the bottom boundary and also near another nucleation site, and it does not slip.

The 100-μm crystal contains more nucleation sites and many active slip planes are evident.
21

Figure 11. Number of dislocations passing through each element for the discrete dislocation simulations of three 

single crystal sizes indicated. The color levels for the 5- and 10-μm crystals are the same. The center 

80 μm is omitted from the 100-μm crystal to highlight the gradients at the top and bottom surfaces.

An important feature of figure 11 is the slip gradient. Since dislocations cannot pass through the 

upper and lower boundaries, the slip (number of dislocations passed) at these surfaces is zero. 

The greatest number of dislocations has passed near the center of the crystal. The slip 

distributions along the 13 active slip planes of the 10-μm crystal are plotted in figure 12. The 

average of the 13 curves is plotted as a heavy black line, and the individual slip results are shown 

as the thin gray lines. All of the curves have a parabolic appearance near the boundaries. Several 

have a flat profile across the center. These are associated with slip planes that have dislocation 

nucleation sites near the surface, and they are also near more highly slipped systems. One of the 

dislocations emitted from the dipoles travels only a short distance to the boundary, and the back 

stress due to the pile up reduces the driving force for further nucleation. The other dislocation 

travels much of the way across the crystal before it meets elevated shear resistance from the pileup at the other side of the crystal. The long travel distance creates the flat center region of the 

curve. Also, the nucleation sites near the boundaries produce fewer dislocations.

100 m

10 m

5 m
22

Figure 12. Number of dislocations passed through each element along the 13 

active slip planes in figure 11 for the 10-m-high crystal plotted against 

the through thickness location. The light lines are the individual results 

and the heavy, dark line is the average.

All of the dislocations must lie between the nucleation sites where they originate and the crystal 

boundaries. The slip distribution is directly related to the current positions of the dislocations. 

The dislocation positions are shown in figure 13. As with figure 11, the entire simulation regions 

are presented for the 5- and 10-μm-thick crystals, and only the upper and lower 10 μm are shown 

for the 100-μm thick crystal. Dislocations of opposite sign originating at the dipoles are shown 

by the red and blue. Dislocations of one sign move to the top and those of the other sign move to 

the bottom. The dislocation density is greatest at the boundaries and tapers off toward the center 

of the crystal. This is the classic picture of an edge dislocation pile-up. The dark red and blue 

elements on the crystal interior indicate the location of the dipole nucleation sites where 

dislocations can accumulate before gliding toward the boundaries.
23

Figure 13. Number of dislocations currently within each element for the discrete dislocation simulations of three 

single crystal sizes. The color levels for the 5- and 10-μm crystals are the same. The center 80 m is 

omitted from the 100-μm crystal to highlight the gradients at the top and bottom surfaces.

Each element can contain multiple dislocations, and the elements adjacent to the boundaries of 

the smaller two crystal sizes contain an average of 16. That average drops to seven in the next 

element away from the boundary, demonstrating a large gradient in dislocation density at the 

crystal boundaries. The first six elements along the slip planes and near the surface all contain 

dislocations in these two simulations. Most slip planes only have a few elements without 

dislocations along their entire lengths. Hence, the dislocation density is relatively dense for the 

5- and 10-μm crystal heights. 

The situation is different for the 100-μm-high crystal. With many more nucleation sites and 

dislocations to carry the deformation, the driving force and pile-ups are much less severe for the 

100-μm crystal, and the maximum number of dislocations per element is four. Here the gradient 

is evident in the sparseness of elements containing dislocations. Near the surface most elements 

along the slip planes contain dislocations, but as the distance from the boundary increases, the 

elements containing dislocations become increasingly sparse. 

The discrete slip creates a highly nonuniform stress field, which is illustrated by the in-plane 

stress components of the 10-μm-high crystal presented in figure 14. The shear stress, Sig-XY, is 

5 m

10 m

100 m
24

low along the slip planes due to sliding, and it builds up between the slip planes. The slip planes 

effectively partition the domain into finite size elastic blocks, each undergoing simple shear with 

the shear stress on the lateral boundaries set by the slip resistance on the slip planes. The shear 

stress rises toward the interior of the blocks as the distance from the shear stress lateral boundary 

condition increases. Displacement is prescribed on the upper and lower boundaries. This imposes 

an additional constraint that keeps the boundary nodes uniformly spaced in the horizontal 

direction, and the shear stress is more or less uniform along the upper and lower surfaces.

Figure 14. Plots of the in-plane stress components at a strain of 0.01 for the 10-μm-high 

single crystal deforming by discrete single slip on 13 slip planes.

As would be expected from an elastically sheared block, the corners experience significant 

normal stress, Sig-YY, as the blocks try to shear and rotate between the rigid platens. At every 

material point, the momentum equation in the vertical direction relates the gradient of Sig-YY in 

the vertical direction to the gradient of Sig-XY in the horizontal direction. The vertical Sig-YY 

gradient adjacent to the slip planes is small near the boundary and becomes approximately 

constant (linear stress profile) in the mid-thickness of the crystal. This coincides with the 
25

horizontal shear stress gradient adjacent to the slip planes being small at the boundary and nearly 

constant in the mid-thickness region. Likewise, the horizontal gradient in the Sig-XX stress is 

related to the vertical gradient in Sig-XY, and similar gradient correlations can be made.

4.2.2 Moderate Density of Nucleation Sites

The simulations of section 4.2.1 were repeated with a tenfold increase in the number of 

dislocation nucleation sites. These were 92, 172, and 1600 for the 5-, 10-, and 100-μm-high 

single crystals. This gives a nucleation site density of approximately 3.2x107cm–2.

Plots of the dislocation locations are shown in figure 15 for the 5- and 10-μm-thick single 

crystals deformed in simple shear to a strain of 0.01. As in the simulation with a low nucleation 

site density, the dislocation dipoles split with dislocations of one sign moving toward the upper 

surface and dislocations of the opposite sign move toward the bottom surface. However, since 

there are significantly more dislocations to accommodate the strain, they do not travel as far as in 

the prior analysis. The opposite signed dislocations are not as segregated and the total number of 

dislocations in the elements at the boundaries is significantly less. The results for the 100-μm 

simulation are not shown. This calculation had even more nucleation sites, less slip, and few 

dislocations accumulated at the boundaries. Consequently, segregation of the opposite-signed 

dislocations to the top and bottom surfaces was not apparent at a strain of 0.01.

Figure 15. Number of dislocations within each element for discrete dislocation simulations of the 5- and 10-μm 

single crystals with a moderate density of nucleation sites.

The dislocations in these moderate source density calculations are much closer together than in 

the low nucleation source runs of section 4.2.1. As a result, the dislocations are too close for the 

finite element discretization to resolve the stress field satisfactorily (figure 16). So while the 

motion of the dislocations in figure 15 appears to be rational, the lack of sufficient spatial 

resolution to resolve the stress gradients renders the solution questionable. A larger element size 

would result in more elements containing dislocations, and the stress resolution would degrade 

further. Use of a smaller element size would allow more elements between dislocations to better 

resolve the field gradients.

5 m

10 m
26

Figure 16. Shear stress distribution for 10-μm-thick crystal simulation with a moderate 

dislocation nucleation site density.

4.2.3 Reduced Element Size

Another set of simulations was run to explore the effect of element size. The mesh dimensions 

used in sections 4.2.1 and 4.2.2 were reduced by a factor of 5 in each direction such that the 

crystal width is 10 μm and the heights were 1, 2, and 20 μm. The element dimension for these 

simulations was 0.025 μm. The dislocation nucleation site seeding used in the mesh of section 

4.2.2 was also used here, so the number of dislocation nucleation sites is the same as in section 

4.2.2. However, the resulting nucleation site density is 25 times greater than in section 4.2.2

(approximately 8x108cm

–2

for each of the three crystal sizes) because the mesh is a factor of 5 

smaller in each direction. An alternative configuration would have been to keep the nucleation 

site density constant, but with the higher spatial resolution, that would have resulted in fewer slip 

planes than are evident in figure 11.

The dislocation locations from these three simple shear runs are shown in figure 17 at a strain of 

0.01. As with the previous simulations, the full crystal heights are shown for the smaller two 

crystal sizes and only the top and bottom 10% of the tallest crystal is shown. Also similar to the 

previous simulations, the dipole nucleation sites create oppositely signed pairs of dislocations 

that glide up and down toward the model boundaries. Pile-ups are evident, and it appears that the 

dislocations are more distinct and separated relative to the crystal height than those in figures 13 

and 15. This is expected since the dislocation spacing is a physical distance and the model has 

effectively zoomed in by a factor of 5. 
27

Figure 17. Number of dislocations currently within each element for the discrete dislocation simulations of three 

single crystal sizes. The color levels for the 1- and 2-μm crystals are the same. The center 16 μm is 

omitted from the 20-μm crystal to highlight the gradients at the top and bottom surfaces.

The slip planes in the 1- and 2-μm-high crystals appear distinct, and the dislocation pile-ups lie 

cleanly along the planes. In contrast, several of the nucleation sites in the 20-μm-high crystal 

have horizontal locations similar to other nucleation sites, resulting in closely spaced active slip 

planes. For these closely spaced slip planes, the locations of the dislocations near the surface 

appear to be coordinated. Such coordination could be related to a lower energy location in the 

interacting stress fields. 

The shear stress field for this higher resolution simulation is shown in figure 18. As with the 

shear stress shown in figure 16, the spatial discretization is too coarse compared to the 

dislocation spacing to resolve the gradients. Many of the lowest stresses, indicated by blue in 

figure 18, are isolated. Consequently, the quality of the solutions is questionable.

1 m

2 m

20 m
28

Figure 18. Shear stress distribution for the 10 μm x 2 μm crystal in a calculation with an 

element size of 0.025 μm and 172 dislocation nucleation sites. The stress gradients 

are under-resolved.

An important piece of information on the shear stress plots is the minimum and maximum stress. 

In figure 16, the maximum is 120 MPa and the minimum is –18 Mpa, whereas these are, 

respectively, 604 MPa and ‒307 MPa in figure 18. The reason can be traced to the quantized 

shear computed through equations 12 and 13. As a dislocation is about to leave an element, it 

exerts an elastic shear strain/stress on the neighboring element. As the dislocation crosses the 

element boundary and traverses the element, the shear stress through goes to zero and reverses 

sign. The stress due to the discrete shear peaks, and it has the opposite sign as the dislocation is 

leaving the element. Analytically, the shear stress is infinite at the dislocation core, and the 

spatial discretization effectively averages the stress over the element size. As the element size 

gets smaller, the peak stress represented in the mesh will increase, but the gradient is less well 

resolved. In figure 18, the peak shear stresses are an order of magnitude higher than the stress 

required for slip. Such high stresses dominate the solution and the error in the gradient 

overwhelms the slip resistance of the crystal. This is another reason that the quality of the 

solution is questionable.

4.2.4 Stress-strain Response from Semi-discrete Model

The shear stress-strain curves from the semi-discrete dislocation simulations are shown in 

figure 19. The element size, nucleation site density, crystal width, and crystal aspect ratio are 

indicated in the legends. The results from all of the simulations in section 4.2 are presented in 

addition to a set of three runs similar to those of section 4.2.2 but with the domains and elements 

10 times larger. As discussed in earlier, the dislocations in this calculation are significantly 

under-resolved.
29

Figure 19. Shear stress-strain response for the semi-discrete crystal simulations of section 4. The legends and arrows 

show the model width, crystal aspect ratio, element length, and the percent of elements containing 

dislocation nucleation sites along with the corresponding nucleation site density. Curves of the same 

color have the same element size and element site nucleation percentage but different model aspect 

ratios.

The color groupings contain sets of three simulations where only the crystal height was varied. 

Within these sets of three curves of the same color, the results from the short crystals are 

represented by dotted lines, the intermediate height crystals are indicated by the dashed lines, and 

the tall crystals by the solid lines. In all cases, the stress increases as the crystal height decreases. 

The green curves show the results from section 4.2.1 where the dislocation nucleation sites were 

sparse. The short crystal contained only seven active slip planes, and relatively large elastic 

regions are providing the shear stress resistance. The few slip planes do little to relieve the 

stress, and the result is nearly elastic. The intermediate height crystal was analyzed in some 

detail above. The additional slip planes created smaller elastic domains, and the stress relief 

provided by the slip planes is sufficient to have an appreciable impact on the overall stress. With 

the additional nucleation sites in the tall crystal, the elastic domains are small and the shear stress 

peaks and then saturates after a small drop. 
30

The nucleation site density was increased by an order of magnitude over the green results to 

produce the red results, section 4.2.2. The mesh and model size are identical. Here there are 

sufficient dislocations to allow pervasive slip even for the shortest crystal. The stress saturates 

and is nearly the same for all three crystal heights. Using the same number of dislocation sites, 

the mesh size is then increased by a factor of 10 in both directions over the red curves to produce 

the black curves. Here too, the stress saturates and simulations for all three crystal heights give 

approximately the same shear stress. The shear stress is lower than the red curve because the 

singular shear stress near the dislocations is averaged over a large volume, and the resulting 

lower element stress near the dislocation lowers the overall shear stress. This simulation was not 

described in earlier sections. 

The calculations producing the blue curves are described in section 4.2.3. This simulation had a 

mesh a factor of 5 smaller in each direction than the red curve, and the number of dislocation 

nucleation sites was the same. Even though there appears to be a sufficient number of 

dislocations in figure 17 to produce a smooth plastic response, the shear stress-strain curves are 

very noisy. The reason is the very high stress near the dislocations that results from the finer 

spatial discretization. Even though stress jumps occur over small volumes as the dislocations 

move from one element to another, the magnitude of the jump is sufficient to affect the average 

stress.

5. Discussion

The results from the continuum slip gradient model presented in section 3 demonstrate that an 

operator split approach can be used with an explicit dynamic time integration scheme to include 

gradient effects in crystal plasticity simulations. The nonlocal gradient formulation was an 

outgrowth of a penalty approach to enforce slip continuity between neighboring elements. The 

model was implemented into a large-scale parallel code with minimal disruption to the flow of 

the calculations. There appear to be additional time step restrictions for smaller element sizes, 

but the nature of the restrictions was not explored in the current investigation.

The single crystal and polycrystal results show the anticipated diffusion of sharp deformation 

fields and the expected stress-strain trends with size scale. While having the appearance of a 

successful modeling effort, a critical comparison of the model assumptions to the physical 

configuration, which it is supposed to represent, reveals a significant disconnect. This is

particularly evident in the micron size range as the spatial resolution of the finite element grid is 

refined. An inherent assumption in the finite element crystal plasticity model is that 

microstructure features are either explicitly resolved by the grid or the microstructure is at a scale 

sufficiently smaller than the element size such that it can be approximated as a smoothly varying 

field. Similarly, a nonlocal gradient extension to the crystal model should resolve the gradient 

over several elements. As indicated above in table 1, the element size must be larger than 1 μm 
31

to have an average of one dislocation per element for a well-annealed material. At lower

resolutions, the field within the element is not smooth and the discreteness of the microstructure 

is not resolved. A continuum crystal model at this or finer resolution does not provide a valid 

representation of the physical microstructure for an annealed crystal. Consequently, simulations 

cannot be expected to capture the microstructure mechanisms at this scale accurately. 

Looking strictly at the dislocation density for heavily deformed metals, from table 1, it appears 

that the crystal model with submicron spatial resolution could represent a smeared dislocation

microstructure within each element. However, dislocations are rarely uniformly distributed in 

heavily deformed metals. They typically organize into walls, which create a cell structure. The 

cell walls have a very high dislocation density, and the cell interiors have a low dislocation 

density. Extensive analysis of heavily deformed nickel by Hughes and Hansen (2000) shows that

these cell sizes are greater than 0.1 μm after cold rolling to a 98% reduction, or a logarithmic 

strain in excess of 3.9. In order to have a smeared representation of such a cell microstructure 

within each element, the element size would have to be on the order of 1 μm or larger. 

Traditional crystal plasticity simulations using smaller element sizes will disregard a 

microstructure that could be resolved at that spatial resolution, and the computed deformation 

mechanisms may not be accurate.

As in the modeling of composites, refinement of the spatial discretization for the standard 

continuum crystal plasticity models has a limit when the element size begins to approach the 

length scale of important microstructure features. In order to have predictive microstructurebased simulations, the modeling strategy has to change with mesh refinement to match the 

physical microstructure at that refinement level. When the element size during mesh refinement 

of a composite approaches the size of the reinforcement, a homogenized continuum model is no 

longer appropriate, and the reinforcement and matrix phases should be modeled explicitly. There 

is nothing to prevent analysis using a homogenized model with elements smaller than the 

reinforcement, but the real mechanisms operating at that scale would be missed entirely. The 

calculations would not be predictive of mechanisms occurring at that resolution.

For pure metals, the limiting microstructure features for a crystal plasticity model are the 

dislocation spacing or the dislocation cell size, whichever is larger. These are typically on the 

order of 0.1 to 1.0 μm. Further mesh refinement near this range of element sizes requires some 

explicit accounting for the dislocation microstructure and a change in the underlying constitutive 

relations. The continuum crystal model, including gradients, will run with smaller elements, but 

the solution will miss the effects of the discrete dislocation interactions and any mechanisms 

specific to dense dislocation walls and cell structures. It would give a false sense of a high 

resolution simulation since the microstructure governing the behavior is not captured spatially. 

The importance of accounting for the discreteness of dislocations in a pile up is highlighted by 

Roy et al. (2008) and Baskaran, et al. (2010). In general, however, mechanisms at this length 

scale have received little attention. Modeling at the size scale of precipitate or phase 

microstructure in alloys would also require further length scale considerations.
32

Nonlocal gradient calculations implicitly rely on the validity of the representation of the 

microstructure within the elements. If the element size is of the order of important microstructure 

features such that deformation mechanisms are not represented appropriately, the gradient 

calculations would not be based on reliable element data. Hence, gradient calculations are 

restricted to the same lower resolution limits.

A semi-discrete model based on the standard crystal plasticity model was investigated as an 

expedient method of crudely accounting for the discrete nature of dislocations at finer spatial 

resolutions. Only elements containing dislocations could slip plastically, and dislocations 

produce a definite level of slip after passing through an element. While the results of this simple 

model showed dislocation pile-ups and other expected features, the specific approach was 

unsatisfactory. The overriding issue is the singular stress field from the dislocations dominating 

the solution as the mesh was refined to resolve the stress gradients. The stress field is driven by 

the quantization of slip, which provides a higher stress magnitude related to the better resolution 

of the singularity at finer spatial resolutions. There does not appear to be a range of element sizes 

that will both resolve the stress field and not suffer from the effects of the singularity. Perhaps 

other semi-discrete approaches could be successful; this one was not.

Discrete dislocation dynamics simulations (e.g., Kubin et al., 1992 and Arsenlis et al., 2007) 

explicitly account for dislocations and their interactions and provide one means for incorporating 

dislocation microstructure at finer spatial resolutions. Finite element methods have been coupled 

with the discrete dislocation simulations though several approaches (Van der Giessen and 

Needleman, 1995; Fivel et al., 1998; Yasin et al., 2001). These types of formulations may be 

employed at the finer spatial resolutions in multiscale modeling schemes. There could be a 

transition from traditional continuum crystal plasticity to such a representation when the spatial 

resolution is fine enough, as in Wallin et al. (2008).

An area that has seen little model development activity is deformation associated with

dislocation cell walls. In addition to being sources and sinks of dislocations for the cells, 

dislocations may run within some cell walls, causing slip. These walls could be treated as 

entities in resolved calculations of dislocations cells where they would interact with models of 

the cell interior which would transmit dislocations across the cells. This approach may be more 

efficient than modeling every dislocation within the cell walls. This is an area for future research.

6. Conclusion

A nonlocal crystal model based on the second gradient of crystal slip was implemented in a 

large-scale parallel finite element code, and the results show the expected trends of decreasing 

the severity of gradients and increasing strength with decreasing physical size. Although the 

model appears successful on the surface, there are concerns over the adequacy of the traditional 
33

crystal model to represent the microstructure at smaller size scales. Elements smaller than the 

dislocation spacing or the dislocation cell dimensions cannot capture features of the deformation 

at that size scale accurately. This prompted exploration of a semi-discrete crystal plasticity 

model. The semi-discrete model also showed realistic results in terms of dislocation pile-ups, but 

the stress solution at the dislocation singularity dominated the solution at smaller element sizes, 

rendering the solution noisy and very mesh dependent. 

The size scale of the microstructure of a typical pure metal, in terms of dislocation spacing and 

size of dislocation cells, is on the order of 0.1 to 1.0 μm. If the element size from a simulation is 

small enough to be within this range, the microstructure features and mechanism should be 

modeled explicitly in order to capture the lower length scale behavior. Models particular to this 

size scale are needed.
34

7. References

Acharya, A. A Model of Crystal Plasticity Based on the Theory of Continuously Distributed 

Dislocations. J. Mech. Phys. Solids 2001, 49, 761–784.

Arsenlis, A.; Cai, W.; Tang. M.; Rhee, M.; Oppelstrup, T.; Hommes, G.; Pierce, T. G.; Bulatov, 

V. V. Enabling Strain Hardening Simulations with Dislocation Dynamics. Modelling Simul. 

Mater. Sci. Eng. 2007, 15, 553–595.

Arsenlis, A.; Parks, D. M.; Becker, R.; Bulatov, V. On the Evolution of Crystallographic 

Dislocation Density in Homogeneously Deforming Crystals. J. Mech. Phys. Solids 2004, 52, 

1213–1246.

Asaro, R. J. Micromechanics of Crystals and Polycrystals. Advances in Applied Mechanics

1983, 21, 1–115.

Baskaran, R.; Akarapu, S.; Mesarovic, S. D.; Zbib, H. M. Energies and Distributions of 

Dislocations in Stacked Pipe-ups. Int. J. Solids, Struct. 2010, 47, 1144–1153.

Becker, R. Effects of Crystal Plasticity on Materials Loaded at High Pressures and Strain Rates, 

Int. J. Plast. 2004, 20, 1983–2006.

Becker, R. Dislocation Transport in Continuum Crystal Plasticity Simulations; ARL-MR-0799; 

U.S. Army Research Laboratory: Aberdeen Proving Ground, MD, 2011.

Bittencourt, E.; Needleman, A.; Gurtin, M. E.; Van der Giessen, E. A Comparison of Nonlocal 

Continuum and Discrete Dislocation Plasticity Predictions. J. Mech. Phys. Solids 2005, 51,

281–310.

Fivel, M. C.; Robertson, C. F.; Canova, G. R.; Boulanger, L. Three-Dimensional Modeling of 

Indent-Induced Plastic Zone at a Mesoscale. Acta Mater. 1998, 46, 6183–6194.

Flanagan, D. P.; Belytschko, T. A Uniform Strain Hexahedron and Quadrilateral with 

Orthogonal Hourglass Control. Int. J. Num. Meth. Engr. 1981, 17, 679–706.

Fleck, N. A.; Hutchinson, J. W. Strain Gradient Plasticity. Adv. Appl. Mech. 1997, 33, 295–361.

Fleck, N. A.; Muller, G. M.; Ashby, M. F.; Hutchinson, J. W. Strain gradient Plasticity: Theory 

and Experiment. Acta Metall. Mater. 1994, 42, 475–487.

Gao, H.; Huang, Y.; Nix, W. D.; Hutchinson, J. W. Mechanism-Based Strain Gradient 

Plasticity-I. Theory. J. Mech. Phys. Solids 1999, 47, 1239–1263.
35

Gerken, J. M.; Dawson, P. R. A Finite Element Formulation to Solve a Non-Local Constitutive 

Model with Stresses and Strains due to Slip Gradients. Comp. Meth. Appl. Mech. Engr.

2008, 197, 1343–1361.

Gurtin, M. E.; Anand, L.; Lele, S. P. Gradient Single-Crystal Plasticity with Free Energy 

Dependent on Dislocation Densities. J. Mech. Phys, Solids 2007, 55, 1853–1878.

Hirschberger, C. B.; Peerlings, R.H.J.; Brekelmans, W.A.M.; Geers, M.G.D. On the Role of 

Dislocation Conservation in Single-Slip Crystal Plasticity. Modelling and Simul. Mater. Sci. 

Eng. 2011, 19, 085002, 1–24.

Hughes, D. A.; Hansen, N. Microstructure and Strength of Nickel at Large Strains. Acta Mater.

2000, 48, 2985–3004.

Hull, D.; Bacon, D. J. Introduction to Dislocations; 3rd Edition, Pergamon Press, Elmsford, NY.

Kubin, L. P.; Canova, G.; Condat, M.; Devincre, B.; Pontikis, V.; Brechet, Y. Dislocation 

Microstructures and Plastic Flow: A 3D Simulation. Diffusion and Defect Data-Solid State 

Data B 1992, 23–24, 455–472.

Lee, T. C.; Robertson, I. M.; Birnbaum, H. K. Prediction of Slip Transfer Mechanisms Across 

Grain Boundaries. Scripta Metall. 1989, 23, 799–803.

Mayeur, J. R.; McDowell, D. L.; Bammann, D. J. Dislocation-Based Micropolar Single Crystal 

Plasticity: Comparison of Multi- and Single Criterion Theories. J. Mech. Phys. Solids 2011, 

59, 398–422.

ALE3D: An Arbitrary Lagrange-Eulerian 2D and 3D Code System, V-4.16 (2012), A. L. 

Nichols, ed., Lawrence Livermore National Laboratory.

Peirce, D.; Asaro, R. J.; Needleman, A. Material Rate Dependence and Localized Deformation 

in Crystalline Solids. Acta Metall. 1983, 31, 1951–1976.

Roy, A.; Peerlings, R.H.J.; Geers, M.G.D.; Kasyanyuk, Y. Continuum Modeling of Dislocation 

Interactions: Why Discreteness Matters? Mater. Sci. Eng. A 2008, 486, 653–661.

Saha, R.; Xue, Z.; Huang, Y.; Nix, W. D. Indentation of a Soft Metal Film on a Hard Substrate: 

Strain Gradient Hardening Effects. J. Mech. Phys. Solids 2001, 49, 1997–2014.

Schouwenaars, R.; Seefeldt, M.; Van Houte, P. The Stress Field of an Array of Parallel 

Dislocation Pile-Ups: Implications for Grain Boundary Hardening and Excess Dislocation 

Distributions. Acta Mater. 2010, 58, 4344–4353.

Stölken, J. S.; Evans, A. G. A Microbend Test Method for Measuring the Plasticity Length 

Scale. Acta Mater. 1998, 46, 5109–5115.
36

Van der Giessen, E.; Needleman, A. Discrete Dislocation Plasticity-A Simple Planar Model. 

Modelling Simul. Mater. Sci. Eng. 1995, 3, 689–735.

Wallin, M.; Curtin, W. A.; Ristinmaa, M.; Needleman, A. Multi-Scale Plasticity Modeling: 

Coupled Discrete Dislocation and Continuum Crystal Plasticity. J. Mech. Phys. Solids 2008, 

56, 3167–3180.

Yasin, H.; Zbib, H. M.; Khaleel, M. A. Size and Boundary Effects in Discrete Dislocation 

Dynamics: Coupling with Continuum Finite Element. Mater. Sci. Eng. A 2001, 309–310,

294–299.
37

NO. OF 

COPIES ORGANIZATION 

1 ADMNSTR

ELEC DEFNS TECHL INFO CTR

ATTN DTIC OCP

8725 JOHN J KINGMAN RD STE 0944

FT BELVOIR VA 22060-6218

2 HCS US ARMY RSRCH LAB

ATTN RDRL CIH C J KNAP

ATTN RDRL CIH C P CHUNG

BLDG 310E

ABERDEEN PROVING GROUND MD 

21005

1 HC US ARMY RSRCH LAB

ATTN RDRL WML B B RICE

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMM B B LOVE

BLDG 393

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMM J BEATTY

4600 DEER CREEK LOOP

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMM R DOWDING

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMP B C HOPPEL

BLDG 390

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMP B 

S SATAPATHY

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMP C B LEAVY

BLDG 393

ABERDEEN PROVING GROUND MD 

21005

NO. OF 

COPIES ORGANIZATION 

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMP C C WILLIAMS

ABERDEEN PROVING GROUND MD 

21005

7 HCS US ARMY RSRCH LAB

ATTN RDRL WMP C J CLAYTON

ATTN RDRL WMP C R BECKER 

(4 COPIES)

ATTN RDRL WMP C S BILYK

ATTN RDRL WMP C T W BJERKE

BLDG 390

ABERDEEN PROVING GROUND MD 

21005

1 HC US ARMY RSRCH LAB

1 PDF ATTN RDRL WMP D R DONEY

(1 ELEC)

ATTN RDRL WMP E P SWOBODA 

(1 HC)

BLDG 393

ABERDEEN PROVING GROUND MD 

21005

2 PDFS US ARMY RSRCH LAB

ATTN RDRL WM P J BAKER

ATTN RDRL WMP 

S E SCHOENFELD

BLDG 309

ABERDEEN PROVING GROUND MD 

21005

1 PDF US ARMY RSRCH LAB

ATTN RDRL WM B FORCH

BLDG 4600

ABERDEEN PROVING GROUND MD 

21005-5066

1 PDF US ARMY RSRCH LAB

ATTN RDRL WM J MCCAULEY

BLDG 4600

ABERDEEN PROVING GROUND MD 

21005-5069

1 PDF US ARMY RSRCH LAB

ATTN RDRL WMP C D DANDEKAR

BLDG 390

ABERDEEN PROVING GROUND MD 

21005-5069
38

NO. OF 

COPIES ORGANIZATION 

2 HCS US ARMY RSRCH LAB

1 PDF ATTN IMAL HRA 

MAIL & RECORDS MGMT

ATTN RDRL CIO LL TECHL LIB

ATTN RDRL SER L A WICKENDEN

(1 PDF)

ADELPHI MD 20783-1197