Mechanical C
Resume:
PHYSICAL REVIEW B **, 104***-****
Quantifying the early stages of plasticity through nanoscale experiments and simulations
Krystyn J. Van Vliet,1,* Ju Li,2, Ting Zhu,2 Sidney Yip,1,2 and Subra Suresh1,
1
Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
2
Department of Nuclear Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
Received 14 June 2002; published 25 March 2003
Nucleation and kinetics of defects at the atomic scale provide the most fundamental information about the
mechanical response of materials and surfaces. Recent advances in experimental and computational analyses
allow us to study this phenomenon in the context of nanoindentation and localized mechanical probing of
surfaces. Here, we present an analytical formulation of the elastic limit that predicts the location and slip
character of a homogeneously nucleated defect in crystalline metals, and extend this formulation to the atomic
scale in the form of an energy-based, local elastic stability criterion, termed the criterion. We demonstrate
that this approach can be incorporated ef ciently into computational methods such as molecular dynamics and
nite-element models. Furthermore, we validate and calibrate the criterion directly through nanoindentation
experiments and two-dimensional experimental analogs such as the bubble raft model. We outline explicitly a
compact and ef cient application of the criterion within the context of a nonlinear, interatomic potential
nite-element model IPFEM . Further, we report three-dimensional molecular dynamics simulations in several
face-centered cubic systems that elucidate the transition from the initiation to the early stages of plasticity
during nanoindentation of metals, as characterized by homogeneous and heterogeneous nucleation of up to
hundreds of dislocations. Correlation of these simulations with direct observations from nanoindentation ex-
periments provides atomistic insights into the early stages of plasticity.
DOI: 10.1103/PhysRevB.67.104105 PACS number s : 62.25. g, 02.70. c, 81.07. b, 81.40. z
I. INTRODUCTION opaque crystals, as well as the challenges associated with
A fundamental transition in the mechanical behavior of conducting suf ciently clean experiments at this length scale
materials is the shift from elastic reversible to plastic irre- in initially defect-free materials.
versible deformation. In crystalline metals at room tempera- Atomic scale simulations such as molecular dynamics
ture, this process is mediated predominantly by the nucle- MD can provide the required visualization and precise con-
trol over characterization of the evolving defect microstruc-
ation and motion of dislocations. Theoretical frameworks
ture over a very well de ned, though limited, time period.
have been established for dislocation nucleation at the
atomic scale,1 4 but currently no such framework exists for Indeed, several researchers have modeled indentation of 3D,
fcc crystals and observed dislocation activity.6 10 However,
the so-called homogeneous nucleation mechanism, where de-
the length scales in such studies are limited by computational
fects are nucleated in the crystal in the absence of other
resources, resulting in indenter tip dimensions and strain
nearby, preexisting defects. Indeed, although we have ana-
rates that are incongruent with those attainable in experi-
lytical formulations of the self-energy of a single dislocation,
ments. Even when computational analyses incorporate mul-
the conditions necessary to nucleate dislocations in the bulk
and in thin lms are unresolved.5 Currently we rely on con- tiple length scales to approach realistic dimensions, such as
in the quasicontinuum method of Tadmor et al.,11 the obser-
tinuum concepts such as critical resolved shear stress and
vations have not been suf ciently well calibrated against ex-
yield strength to quantify the onset of plasticity in metals,
periments and do not culminate in predictive capabilities.
even though this transition is initiated at the atomic scale.
Thus, it appears that although both experiments and compu-
The need for nanoscale interpretation of defect nucleation
tational modeling have the potential to explore the same phe-
is even more pressing now that a sophisticated range of ex-
nomenon, neither has been constructed in such a way that its
perimental and computational tools are available to explore
results truly address the hypotheses and quantitative capabili-
material response near the atomic level. In fact, recent devel-
ties of the other.
opments in both multiscale computational modeling and
In this work, we outline experimental, analytical, and
nanoscale experiments allow these two distinct methods of
computational analyses of nanoindentation that were de-
material analysis to explore mechanical behavior at length
signed together to elucidate the mechanisms by which local-
scales comparable to one another. Nanoscale experiments via
ized contact can induce the transition from elastic to plastic
nanoindentation, atomic force microscopy, and other tech-
deformation in an initially defect-free crystal, and the mecha-
niques can probe mechanically the elastic stability of crystals
nisms by which that plasticity proceeds at early stages, i.e.,
at subnanometer displacements and large local contact
mechanisms of incipient plasticity. This discussion builds on
strains. Results from these and other related experiments
our earlier results, reported by Li et al.12 First, we summarize
have indicated that such contact can induce elastic instabili-
experimental results of nominally sharp nanoindentation for
ties that have been deposited as evidence of homogeneous
both 3D fcc crystals and a two-dimensional 2D experimen-
dislocation nucleation. The chief obstacle to further valida-
tal analog. Next, we present an analytical description of elas-
tion of such experiments is the considerable dif culty in de-
tic stability at the atomic scale in the form of the criterion.
termining atomic scale activities in three dimensional 3D,
0163-1829/2003/67 10 /104105 15 /$20.00 67-104***-*-**** The American Physical Society
VAN VLIET, LI, ZHU, YIP, AND SURESH PHYSICAL REVIEW B 67, 104***-****
TABLE I. Experimental measurements of critical indentation load in fcc crystals.
E *a b
Element Orientation Indenter radius Pc G /2 Reference
max
nm (N GPa GPa GPa
Aluminum 110 50 30 72.6 7.2 4.36 13
111 50 8 76.0 4.7 4.6 13
133 50 15 63.6 5.7 4.46 13
NAd
113 11.3 65.8 4.0 47
111 c
Copper 50 35 123 10.7 7.6 46
Gold 110 205 7.9 117 2.5 4.6 49
111 205 31.3 81.8 3.1 6.6 49
111 70 14.5 117 6.2 6.6 8
100 205 43.8 43 2.2 2.46 49
Iridium 100 50 58.5 370 26.4 34.3 50
a
Equation 16 for E i 1100 GPa and 0.07.
i
b
Equation 15 .
c
60 80 % 111 texture.
d
Reference 48; not available.
We verify the corresponding predictions of defect nucleation rough correlation cannot be used to rationalize subsequent
and slip character with experiments and MD simulations of displacement bursts, or to determine the initial location and
2D indentation, as well as with MD and nonlinear nite el- structure of such homogeneously nucleated defects.
ement modeling FEM of cylindrical indentation in elasti- Using reliable electronic structure methods, Roundy et al.
cally anisotropic crystals. Finally, we present experimental calculated the ideal shear strengths on 111 planes of Al and
results and 3D MD simulations of nominally sharp nanoin- Cu to be 1.85 and 2.65 GPa, respectively, when the remain-
ing ve stress components are fully relaxed.14 The nominal
dentation that reveal that, subsequent to the initial elastic
instability, the development of incipient plasticity under con- difference between the above theoretical predictions and
tact loading is not due to continued, homogeneous nucleation Table I is signi cant, but direct comparison may not be
of dislocations. Rather, instabilities observed beyond the rst meaningful due to the following facts. As Roundy et al.
event are due to the formation of heterogeneous defect noted, the tangential modulus is softened signi cantly just
sources that operate at local stresses well below those re- prior to yielding of a crystal, causing the linear elasticity
quired for homogeneous nucleation. Although our experi-
mental and computational correlations focus on nanoscale
contact, the analytical and computational tools presented
herein are applicable to the general phenomenon of large
strain deformation, a process relevant to a wide range of
disciplines.
II. EXPERIMENTAL OBSERVATIONS
A. Nanoindentation
The results of load-controlled, nominally sharp nanoin-
dentation experiments on several different fcc crystals have
been reported, as summarized in Table I. Although these ex-
periments were conducted on various instruments and the
elastic moduli E of the crystals ranged from 69 to 528 GPa,
all of the indentation responses exhibited displacement ex-
cursions beyond some critical indentation load P c see Fig.
1 . In fact, as suggested by Gouldstone et al.,13 a continuum
approximation of the maximum shear stress corresponding to
P c for each crystal indicates that the rst displacement ex- FIG. 1. Schematic of indentation load-displacement response
cursion occurs when the Hertzian maximum shear stress max indicating initial defect nucleation at a critical load P c, and subse-
is on the order of the theoretical shear strength of the crystal, quent defect activity at increased loads and displacements. a Dis-
estimated as G /2 where G is the resolved shear modulus. placement control, showing sharp decreases in load dips corre-
This observation suggests that the local stresses prior to these sponding to defect activity. b Load control, showing large
discontinuities are suf cient to induce homogeneous nucle- increases in indenter displacement bursts corresponding to defect
ation in crystals of initially low defect density. However, this activity.
104105-2
QUANTIFYING THE EARLY STAGES OF PLASTICITY . . . PHYSICAL REVIEW B 67, 104***-****
to a given applied load depends directly on the dimensions of
the indented substrate. Thus, although it can be inferred from
continuum analysis that the position of nucleation corre-
sponds to that of the shear stress maximum, the applied
stress prior to nucleation cannot be measured experimentally.
For this reason, the quantitative results from nanoindentation
( P c) and bubble raft experiments nucleation position and
dislocation structure cannot be compared directly.
III. ANALYTICAL CRITERION OF ELASTIC STABILITY
Unequivocal validation of hypotheses regarding contact-
mediated instabilities requires not only precise correlation of
atomic scale motion and the far- eld load-depth response,
FIG. 2. Color online Nanoindentation normal to the terminal but also a theoretical framework that offers predictive capa-
edge of a close-packed plane of a 2D bubble raft produces subsur- bilities from which further experiments can be designed. To
face dislocation dipole nucleation. Nucleation depth z c scales with this end, we rst review the work of Hill19 and Rice20 on
indenter tip radius R, but the dislocation dipole structure and motion continuum formulations of elastic stability, and then extend
upon nucleation is invariant. The scale bar is 10 mm, or 3 nm in the
these concepts to the atomic scale in the Cauchy-Born frame-
atomistic size scale analogy of Gouldstone et al. Ref. 16 .
work for crystalline metals.
extrapolation15 that Table I uses to overestimate the local
A. Derivation of the criterion
stress at the site of homogeneous nucleation. On the other
19
Hill developed the general concept of continuum me-
hand, the stress condition at the defect nucleation site be-
dium stability in terms of acceleration waves, which in the
neath an indenter is far more complicated than the pure shear
context of this paper can be understood simply as long-
condition that Roundy et al. studied. Our MD calculations
wavelength phonons or elastic waves. He derived a tensor
show that the mean pressure at the defect nucleation site is
representing the connection between wave frequency and the
on the order of tens of GPa; therefore signi cant pressure-
corresponding wave vector k and displacement vector w: 19
induced shear hardening is possible. These two countering
factors make the comparison between the idealized elec-
k, w C i jkl w i w k k jkl, 1
tronic structure method results and the particular response jl
under indentation-induced strains a dif cult task. The reso-
where j l is the Cauchy stress and C i jkl is the isothermal
lution of this discrepancy can be mediated by a rigorous,
elastic constant, de ned as,21
localized elastic stability criterion that is applicable under
general loading conditions. 1,
FY FX X X C X
ij ij
2 i jkl ij kl
B. Bubble raft analog 2
In order to visualize readily the response of atoms to lo- where F ( X ) and F ( Y ) are the Helmholtz free energies of
calized nanoscale contact, Gouldstone et al.16 employed the homogeneous con gurations X and Y and is the Lagrang-
Bragg-Nye bubble raft model, a two-dimensional experimen- ian strain of Y with respect to X. For a given k, the value of
tal analog to close-packed crystals.17 The bubble raft pro- w that gives rise to an extremum in ( k, w) is the phonon
vides a perfectly arranged, defect-free, 2D lattice comprising eigendisplacement, and in this context is proportional to
several thousands of bubbles arranged analogously to the the vibrational frequency squared. In the long-wavelength
closest-packed plane in an fcc crystal. The surface tension of limit, k and w are immaterial, and hence we stipulate that
the solution and the bubble diameter can be designed such
that the interaction among bubbles re ects a simple descrip- k 1, w 1, 3
tion of the interatomic potential of Cu.18 By relating the
and hereafter k and w are used only to denote directions.
indentation-induced motion of millimeter-scale bubbles in a
If there exists a pair of k and w such that ( k, w) is
close-packed raft to that of nanometer-scale atoms in a close-
negative, then the elastic stability of the homogeneous me-
packed atomic plane, it was shown that homogeneous dislo-
dium is lost: negative signi es negative concavity of the
cation dipole nucleation occurs subsurface, at the location of
free energy surface, and thus the frequency of the corre-
maximum shear stress predicted by 2D continuum contact
sponding elastic wave is imaginary and its amplitude in-
mechanics. Further experiments con rmed that, although the
creases exponentially with time. This would ultimately result
depth of dislocation nucleation scales with the width of con-
in the destruction of the homogeneous medium and the cre-
tact, the dipole structure and subsequent formation of surface
ation of defect singularities within. The growth of the soft
steps are independent of length scales such as the indenter
mode wave proceeds in four stages: i linear growth of the
radius see Fig. 2 . One inherent limitation of these experi-
unstable elastic wave, well described by continuum; ii non-
ments is that the load depth response does not converge for
linear evolution at larger amplitudes, over which the wave
2D indentation. That is, the indentation depth corresponding
104105-3
VAN VLIET, LI, ZHU, YIP, AND SURESH PHYSICAL REVIEW B 67, 104***-****
pro le steepens; iii progressive steepening of the wave models such as FEM with atomistic simulations. Addition-
front until its length scale approaches atomic spacing, upon ally, as the criterion relies on the long-wavelength limit
which its description must be transferred from a continuum and thermoelastic material response, it continues to function
to an atomistic basis; iv arrest of the atomistically sharp at nite temperatures, when the phonon description loses its
wave front in a low-dimensional atomic energy landscape, at rigorous de nition. We will demonstrate that in realistic
which point a defect is nucleated. If the wave is transverse, loading con gurations such as nanoindentation, Eq. 1 cor-
the instability will likely result in the formation of a disloca- rectly predicts the critical load and the defect nucleation site.
tion or twin of slip plane normal k and Burgers vector direc- In addition, this energetic criterion reasonably predicts the
tion w; if the wave is longitudinal, a microcrack would likely character of the defect when we combine information about
result. The exact character of the resulting defect is a func- the soft mode k and w with crystallographic structure. That
tion of both the continuum-scale linear instability caused by is, we can predict whether the linear instability will lead to a
external loading and the atomic scale interactions within the
dislocation loop or a twinning embryo, specifying the slip
material. This dynamic process is studied in detail in another
plane and the Burgers vector. This allows true predictive
paper22 and is, in general, complicated. Acknowledging this
power in nite-element models, afforded by a defect nucle-
complexity, we make the following observations for most
ation criterion with no adjustable parameters for a given
monatomic crystals:
description of atomic interactions.
1 If the angle between k and w is greater than 60, i.e.,
the initial soft mode is chie y transverse, then the nal de- B. Incorporation of the criterion
fect created is likely to be a dislocation loop or a twinning into computational analyses
embryo. The structure of the nal defect should correlate
Two practical dif culties prevented the criterion, for-
strongly with whether k is the plane normal of a crystallo-
mulated at the continuum scale several decades ago, from
graphic dislocation slip system or twinning slip system.
being widely used in practice. The rst challenge is the lack
2 If the the angle between k and w is less than 30, i.e.,
of accurate constitutive models to estimate i j and C i jkl at
the initial soft mode is mostly longitudinal, then the nal
very large strain, near elastic instabilities. The parametriza-
defect created is likely to be a microcrack.
tion of the strain energy function or stress functions in six-
Rice20 derived the same stability criterion in the context
dimensional general strain space is by no means a simple
of a particular continuum phenomenon, that of shear band
task, either conceptually or numerically. Prior to the advent
formation. From this viewpoint, one considers an interface of
of modern, reliable electronic struture codes and their deriva-
normal k between a plastically deforming shear band and
tive empirical interatomic potentials, few data existed. An-
elastically deforming material, across which there is a dis-
other obstacle is the accurate calculation of min . Even
placement discontinuity w. When ( k, w) becomes nega-
tive, Rice showed that traction equilibrium across the inter- though four-dimensional minimization since k 1,
face can no longer be maintained, and the shear band thus w 1 ) is not an especially dif cult computation to perform
grows catastrophically. Although the nomenclature and ar- a single time, a reliable and ef cient algorithm is necessary
ticulation of the criterion formulated by Hill and Rice con- due to the iteration of this procedure for each material point.
trast slightly, they are entirely equivalent in content. Math- For example, in the context of MD, the calculation of min is
ematically, this criterion can be summarized as performed for each of up to 2 106 atoms at each time step.
The rst challenge can be resolved naturally in an atom-
min C i jkl w i w k k jkl, 4
min jl
istic calculation because one can conceptualize and evaluate
k 1, w 1
atomistic local stress i j and atomistic local elastic constant
and if min is negative, the speci c material point or repre- C i jkl by partitioning the virial and Ray sums,25,24 respec-
sentative volume element RVE becomes elastically un- tively, down to individual atoms. For pair and embedded-
stable.
atom type interatomic potentials, such partitioning is
A full phonon spectrum analysis of a strained crystal can
straightforward because the total virial and Ray sums are
be carried out; this has been performed for binary com-
comprised algebraically of contributions from atom pairs,
pounds such as SiC.23 Such calculations show the connection
and one can readily attribute half of every pair contribution
between Brillouin zone boundary soft phonons and the ideal
to each of the two participating atoms. More complicated
strength of these systems. Indeed, because the phonons con-
interatomic formulations and partitioning schemes26 would
stitute a complete basis set for all atomistic excitations in a
work equally well, provided that basic symmetry and sum
crystal, such analysis is necessary and suf cient to ascertain
rule requirements are satis ed. This partitioning scheme is
the stability of any crystal at T 0 . The present elastic wave
justi able insofar as the atom under consideration is a few
analysis, however, does not guarantee stability: It is a neces-
atomic spacings away from any topological defects, even if it
sary condition. Despite this lack of suf ciency, the concept
is under a large elastic strain. In other words, as long as
and implementation of the criterion are extremely power-
several nearest-neighbor levels of an atom retain their bulk
ful and cover the majority of instabilities that occur in mon-
crystalline coordination, the concept of atomistic local stress
atomic crystals. The fact that Eq. 1 is a tractable analytical
formula involving only well-known mechanical parameters and elastic constant is well posed and one can evaluate them
i j and C i jkl, which are also well studied in atomistic directly from an interatomic potential, or even from ab initio
simulations,24 makes it an ideal means to connect continuum calculations.
104105-4
QUANTIFYING THE EARLY STAGES OF PLASTICITY . . . PHYSICAL REVIEW B 67, 104***-****
general purpose nite-element package such as ABAQUS28
To overcome the second dif culty, we observe that
and ii an accurate constitutive relation be provided to the
( k, w) is a biquadratic function in w and k, and can be
nite-element code. The latter is facilitated by the model of
alternatively written as
Cauchy-Born elasticity by which the local stress state at a
material point is dictated by the movements of underlying
k, w C i jkl k j k l w i w k jlk jk l, 5
lattices.29,30 Speci cally, FEM-imposed boundary conditions
where the order of the k, w summation is swapped from Eq. cause displacement of the nodes. The stress at each node is
1 . If we assume w is a constant in Eq. 1 and minimize calculated in a subroutine comprising a small atomic lattice
against k only, then the solution for that simple quadratic for which the lattice interaction is determined by an inter-
form is straightforward: Diagonalize the 3 3 symmetric atomic potential i.e., molecular statics at T 0 and ensemble
matrix A j l C i jkl w i w k j l and choose the softest eigen- averaging at nite T ). When this lattice is strained according
vector to be k. Alternatively, if we assume that k is a con- to the FEM nodal displacement, interatomic interactions dic-
stant in Eq. 5 and minimize against w, then the solution is tate stresses on each atom according to the potential, giving a
also straightforward: Diagonalize the 3 3 symmetric ma- local nodal value of i j . Both C i jkl and i j can be calcu-
trix B i k C i jkl k j k l and choose the softest eigenvector to be lated through manipulation of Eq. 2 . The quantity ( k, w)
w. Therefore, we propose the following iterative algorithm can also be calculated and minimized for each node, accord-
to minimize ( k, w): ing to the procedure outlined above. In this manner, the de-
formation of large-scale models can be dictated by an essen-
a Choose an arbitrary w.
tially atomic scale strain-to-failure criterion; this criterion is
b Minimize Eq. 1 for a xed w, identifying an initial
parameter-free for a given interatomic potential. As the con-
value of k that minimizes ( k, w xed).
stitutive response of this FEM approach is dictated by inter-
c Minimize Eq. 5 for a xed k identi ed in the previ-
atomic potential IP, we will refer to it hereafter as IPFEM.
ous step, identifying an initial value of w that minimizes
Although IPFEM is not the rst computational approach that
( k xed, w).
allows observation of defect nucleation, this model repre-
d Repeat until mutual convergence is reached.
sents the rst framework that incorporates a predictive, local
The above uses only 3 3 symmetric matrix diagonaliza-
elastic stability criterion that is directly linked to defect
tion, which can be implemented ef ciently using the analyti-
cal, Cardano cubic equation solution.27 nucleation.
Using the approaches outlined above, a complete recipe
for criterion implementation in atomistic simulations is
IV. COMPUTATIONAL MODELING
proposed:
OF ELASTIC STABILITY
1 Select all atoms suf ciently separated from topologi-
cal defects, even if they are withstanding large elastic strain. Computational simulations such as MD are particularly
2 Evaluate the local atomistic stress and elastic constant well suited to quantify the relationship between internal
on an atom-by-atom basis by distributing the virial and Ray atomic activities and external loading sustained by a crystal:
sums.25,24 This allows ( k, w) to be properly de ned for The only assumption of material behavior in MD is that a
each atom. certain interatomic potential correctly describes the interac-
3 Minimize ( k, w) for each atom, and record min, k, tion among atoms. However, due to the computational cost
and w for each atom. min can be interpreted as the micros- associated with keeping track of large numbers of atoms, the
tiffness of the atom, and its sign can be interpreted as the physical time scales and/or length scales are often unrealis-
concavity of the free energy F. tically small. By contrast, continuum-based computational
4 Among all atoms in the con guration, identify the approaches such as FEM are, by de nition, designed to
atom with the smallest min as the softest atom for the simulate behavior over any length scale for which the mate-
particular displacement con guration. rial response can be described by continuum parameters
Despite the large number of atoms in a typical 3D MD e.g., Young s modulus E and yield strength y). The objec-
simulation of nanoindentation ( 2 106 ), the ef ciency of tive of IPFEM, then, is to incorporate the accuracy and
the above minimization procedure is suf cient to allow si- parameter-free nature of MD and related simulations into the
multaneous simulation and calculation of min . At the onset realistic length and time scales of FEM.
of homogeneous nucleation, min of the softest atom tends IPFEM is distinct from the quasicontinuum QC
method31 in that it does not transfer the basis of the material
to zero. Note that the criterion does not force w and k to
correspond to close-packed directions and planes. Rather, description from continuum-based elements to atomistic-
correspondence of w and k to 110 and 111, respectively, based MD when the elemental response approaches the elas-
in the 3D simulations discussed below is a natural outcome tic limit. Rather, IPFEM is a fully continuum nite-element
of the interatomic potential within the fcc crystal structure. model at all levels of strain, for which the nodal response of
Because instabilities predicted by the criterion are elas- crystalline RVEs is governed by an atomistic Cauchy-Born
tic in nature,21 and since all atomic information is incorpo- constitutive relation.29 In contrast, QC hybridizes the use of
rated through C i jkl and i j, this defect nucleation criterion linear and Cauchy-Born nonlinear elasticity with MD, trans-
also can be readily incorporated within nite-element analy- ferring the constitutive relations according to the level of
sis. Such implementation requires that i the large strain local strain, and thus utilizes a numerical, rather than analyti-
formalism be used e.g., the nonlinear geometry option in a cal, nucleation criterion. The elds calculable via IPFEM
104105-5
VAN VLIET, LI, ZHU, YIP, AND SURESH PHYSICAL REVIEW B 67, 104***-****
are, by design, easily integrated with the ( k, w) analysis
tool. As such, IPFEM is conceptually less general than the
quasicontinuum method and can handle only a subset of the
problems the quasicontinuum method is theoretically capable
of treating, for which the results of the two methods should
be largely the same. That is, although the quasicontinuum
method can be used to identify and validate a defect nucle-
ation criterion, to our knowledge the quasicontinuum method
has not been applied to this task. Three considerations lead
us to think that IPFEM is suf ciently important and unique
to be identi ed separately from the general quasicontinuum
method:
1 Difference in perspective: We do not regard the small
atomic lattice in question as being embedded in a larger,
inhomogeneously strained atomic con guration. Instead, we
consider the small atomic lattice to be embedded in an in -
nite, but homogeneously strained lattice, representing a ma-
terial point in the FEM calculation. This more primitive FIG. 3. Color online Generic short-ranged interbubble poten-
viewpoint enables modular design. That is, we can couple tial with equilibrium separation distance r 0 0.85. Both the poten-
FEM to either interatomic potential or to more accurate ab tial well depth and cutoff distance are set to unity and serve as
initio periodic boundary condition PBC unit-cell normalizations of energy and length, respectively. This potential
calculations.32 Although the incorporation of ab initio consti- was used in molecular dynamics simulations of the bubble raft, our
tutive descriptions is not a large advance conceptually, this model 2D system.
ability would have great practical impact with respect to re-
liability and accuracy of observations and predictions made which is plotted in Fig. 3. Here, the primary units of length
at the atomic scale. and energy are the cutoff distance r c at which the potential is
2 Ease of implementation: One may utilize any general- valued zero and potential well depth, respectively, and are
purpose FEM package such as ABAQUS to implement both scaled to 1. Consequently, a single parameter r 0, the
IPFEM, simply by supplying the appropriate constitutive re- equilibrium interbubble distance, controls the shape of the
lation. Thus, IPFEM requires no modi cation of the mesh- potential for any interbubble distance r r c . We normalize
ing, visualization, and analysis tools, whereas the general this distance by r c valued unity, so that r 0 is a dimensionless
quasicontinuum method requires signi cant effort to code quantity less than unity. In order to convert from these nor-
and implement. malized values to physical units in actual experiments, one
3 Dedication to an important phenomenon: The special multiplies normalized quantities by physical measures of
set of problems that IPFEM can study, e.g., the site and char- and r c e.g., critical load of value P, multiplied by / r c
acter of homogeneous nucleation occurring in nanoindenta- 1.55 J/0.003 m, to obtain force in N . We express the
tion, is suf ciently important and underdeveloped that a quantities that modulate the potential, and, as
simple method that achieves this goal with minimal effort
should be promoted. As one cannot truly separate a method
4 2
1 r0, 2 1 r0 . 7
from the problems it can solve, IPFEM is, by design, linked
to application of the criterion and the general phenomenon
of homogeneous defect nucleation. We observe more brittle behavior at larger r 0 . In all the
In this section, we discuss prediction of elastic instability following results, r 0 is taken to be 0.85, and plots of poten-
via the criterion for both MD and IPFEM simulations of tial and depth are shown in these reduced units. At equilib-
nanoindentation in two- and three-dimensional crystals. rium, the potential energy per bubble is 3, the volume per
3 r 2 /2. At equilibrium, the raft is an isotro-
bubble is 0 0
A. MD simulation of indentation on a 2D half plane pic elastic medium with elastic constant matrix,
In order to verify the assumptions and implementation of
the criterion, we conducted 2D MD simulations of the 2 0
bubble raft experiments outlined in Gouldstone et al.16 As
2 0
C, 8
we are not interested in the mechanical response of an actual
soap bubble raft, but instead regard it as a generic model of 0 0
close-packed 2D crystals, we arbitrarily choose the following
potential with a smooth, short-range cutoff for ease of imple- whereby the fourth order C i jkl tensor is reduced to a 3 3
mentation, matrix due to symmetry in an isotropic elastic system in
which Voigt indices 1,2,3 denote xx, y y, xy stress/strain com-
4
r 1 2,
r1 r1 ponents, respectively. The strain energy per bubble at uni-
Vr 6
form Lagrangian strain is
0, r1
104105-6
QUANTIFYING THE EARLY STAGES OF PLASTICITY . . . PHYSICAL REVIEW B 67, 104***-****
order of 400. Therefore max, the maximal oscillation fre-
T
E 3 V r0 1,0 1 2 1,0
quency, is 3 . We use the six-value Gear predictor-corrector
algorithm33 to integrate the trajectory of atoms moving on
T
V r0 1/2, 3 /2 1 2 1/2, 3 /2
the interatomic potential surface, with an integration time
step of 0.003, corresponding to about 100 force evaluations
T
V r0 1/2, 3 /2 1 2 1/2, 3 /2 . 9
per period for the hardest oscillation, ensuring numerical sta-
bility. The thermostat scheme34 of Berendsen et al. is used,
Using the fact that V ( r 0 ) 0, we may expand the above as
with the temperature-coupling time constant chosen to be 3,
2
V r2 3 3 making the system strongly dissipative such that vibrations
0 xx xy yy
2
E xx
2 4 2 4 due to kinetic energy are quickly restored to the level desig-
nated by the speci ed temperature; this minimizes re ection
2
3 3 of oscillations from boundary conditions. At the beginning of
xx xy yy
4 2 4 a simulation, the indenter is placed such that its perimeter
exactly touches the rst row of bubbles, the indenter force is
V r2 9 2 2 2
9 3
3 initially zero but it increases immediately upon penetration.
0 xx xx yy yy xy, 10
2 8 4 8 2 The bottom layer of bubbles is xed, constraining the raft to
be deformed by the indenter. PBCs are used in the x direction
from which we can derive that
in-plane . In the MD simulation, the indenter proceeds in
displacement control, moving toward the substrate by a
9V r2 3V r2 3V r2
0 0 0
small, xed increment every time step. In the conjugate gra-
C 11 C 22, C 12, C 33,
8 8 8
0 0 0 dient minimization, the indenter advances by 0.02 every step,
11 after which the entire raft of bubbles is relaxed except for the
or if we substitute bottom layer, so the net force on each bubble is zero. In both
simulations, the force V ext( r ) exerted on the bubbles is
8 summed and is correlated with the indenter displacement. We
2
V r0 12 r0 1 2, 12 found it more advantageous computationally to explore the
2
1 r0
general phenomenon using molecular dynamics with smaller
we obtain bubble rafts of thousands of bubbles. However, in order to
verify the criterion with maximum accuracy, we utilized
23 conjugate gradient minimization on larger bubble rafts com-, 13 prising 100 000 bubbles.
2
1 r0
Figure 4 shows the change in far- eld indentation load
and the Poisson s ratio is therefore 1/3. When we de ne the and min as a function of imposed indenter displacement for
equilibrium separation distance r 0 0.85, we obtain a 2D MD simulation comprising 112 400 bubbles. The in-
153.96; this value of under zero far- eld applied load- denter is cylindrical with radius R of 160, d of 0.8, and the
ing is the value of min . Any shear wave with k w mini- interatomic potential is given by Eq. 6 . The critical inden-
mizes ( k, w) equally well, with no preference along crys- tation depth for the rst homogeneous defect nucleation
tallographic directions, prior to the application of far- eld event is 14.84, in reduced units; conversion to physical units
loading. requires multiplication of this value by a real cutoff length
The indenter is implemented as an external potential on r c . The temporal correspondence of discontinuous load re-
the bubbles, laxations and the vanishing of min is excellent for the rst
load discontinuity, and the same trend is observed for ve
Rr d subsequent discontinuities. Through visualization of the
rR
exp,
atomic coordination number, we have veri ed that each load
d Rr
V ext r 14
relaxation corresponds exactly with the nucleation of one
0, r R,
dislocation dipole. This result indicates that the criterion
where r is the distance between any bubble to the origin of can predict ac
Contact this candidate