The application of an atomistic J-integral to a ductile crack Jonathan A. Zimmerman and Reese E. Jones Mechanics of Materials Department, Sandia National Laboratories, Livermore, CA 94551, USA Abstract. In this work we apply a Lagrangian kernel-based estimator of continuum fields to atomic data to estimate the J-integral for the emission dislocations from a crack tip. Face-centered cubic (fcc) gold and body-centered cubic (bcc) iron modeled with Embedded Atom Method (EAM) potentials are used as example systems. The results of a single crack with a K-loading compare well to an analytical solution from anisotropic linear elastic fracture mechanics. We also discovered that in the postemission of dislocations from the crack tip there is a loop size-dependent contribution to the J-integral. For a system with a finite width crack loaded in simple tension, the finite size effects for the systems that were feasible to compute prevented precise agreement with theory. However, our results indicate that there is a trend towards convergence.

keywords: Fracture mechanics, molecular simulation, J-integral, ductile fracture. 1. Introduction Among the most useful metrics developed for the analysis of defects in materials is the J-integral [1]. Based on the seminal work by Eshelby [2, 3], the J-integral is a path independent integral that estimates the energetic driving force acting to propagate an existing defect in a continuous medium. The J-integral is commonly used in numerical simulations of continuum mechanical deformation to indicate when a critical loading state has been achieved resulting in crack growth. The critical value of the J-integral corresponding to this onset of crack propagation, Jc , is referred to as the fracture toughness. The strength of the J-integral lies in its path-independence. This property enables the J-integral contour to be constructed in regions of a loaded body where only linear elastic deformation has occurred, yet still be accurate when non-linear, inelastic deformation is occurring in close proximity to the crack tip. Although many have used this capability to examine structures that contain large amounts of plastic deformation, it has also been used to examine the fundamental issue of dislocation nucleation/emission from an originally pristine crack tip. Rice [4, 5, 6, 7] pioneered this area by using material properties of Young’s modulus, Poisson’s ratio, and the theoretical unstable stacking fault energy (γus ) to predict the value of J necessary for dislocation nucleation (Jnucl ), and thus also whether such nucleation is energetically

The application of an atomistic J-integral to a ductile crack

2

favorable over brittle crack-tip extension for specific loading conditions and crystal orientations relative to the crack geometry. These original works have been expanded on by other researchers to address issues such as what effect elastic anisotropy has on the prediction of Jnucl [8, 9], whether the critical driving force needed to induce continued motion of an existing dislocation is less than its nucleation criterion [10], what the behavior of the J-integral for a system undergoing elastic-plastic deformation due to the presence of geometrically necessary dislocations [11], and the effect that pre-existing dislocations surrounding a crack tip can have on its driving force [12, 13]. Also of significance is the issue of crack-tip shape, and how it impacts the dislocation emission process and the estimate of Jnucl [14, 15, 9]. When considering atomistic simulations, this issue is complicated by the presence of surface stresses along the crack faces [15]. One of the issues that has not been examined is how dislocation emission affects the J-integral once the dislocation glides outside of the contour used for its estimation. For cases of gross plasticity, the expectation is that J loses its path-independence attribute due to unloading of the material in the region where it is being calculated. However, it is not clear that this expectation would hold for the emission of a single dislocation. Moreover, in face-centered cubic (fcc) materials it is often the case that a partial dislocation is nucleated, leaving behind a stacking fault trail. It is not obvious whether a stacking fault has a non-zero contribution to the J-integral, and to our knowledge such an estimation has never been done. In this paper, we use atomistic simulation at zero temperature to examine dislocation emission from a crack tip in both an fcc and a bcc material, and evaluate the J-integral both before and after emission occurs. Recently, we developed a method for calculation of the J-integral using a Lagrangian kernel-based estimator of the continuum fields of stress, energy density and displacement gradient from atomic data, and then integrating the resulting Eshelby stress over a contour surrounding a crack tip [16, 17]. For the case of a brittle material in which cleavage (crack extension) occurs, the Jintegral estimated using this method displays path independence, and good correlation with linear elastic fracture mechanics. Specifically, the trend in the calculated J-integral prior to propagation and the value at which propagation occurs both agree with theory and given confidence in the calculated fracture toughness. This method was originally developed for zero temperature analyses [16], but was later extended to treat systems at finite temperatures [17]. Here, we examine its applicability to ductile materials that preferentially exhibit dislocation nucleation over cleavage. We analyze two such geometries: a “semi-infinite” crack in a cylindrical system with boundary conditions consistent with application of a specified stress intensity factor, and a double-edged, finite width crack in a rectangular plate subject to periodic boundary conditions (thus emulating an infinite array of finite width cracks separated by the system width) and a remotely applied strain. Our results show that for the fcc cylindrical system, dislocation nucleation occurs at the theoretically expected value of J, and for the bcc system the trend agrees with theory but emission occurs at a lower loading, as in [18]. For the rectangular systems the value of J at nucleation is strongly influenced by both the

The application of an atomistic J-integral to a ductile crack

3

interplay of system size with boundary conditions, and crack-tip interactions. For both cases, we examine the behavior of the J-integral after dislocation emission occurs. 2. Construction of the J-integral In zero-temperature equilibrium, the Eshelby stress tensor S [2, 3] is given by S = W I − HT P ,

(1)

where W is the (internal) energy density, I is the identity tensor, H is the displacement gradient and P is the first Piola-Kirchhoff (PK) stress. The displacement gradient is a kinematic measure H = ∇X u defined in terms of the displacement u of material from a reference position X. It is well-known that the Eshelby stress acts as a material force conjugate to the change in the reference configuration due to the evolution of a defect. Given this fact, the J-integral has many applications to the mechanics of defects, specifically in characterizing the evolution and propagation of cracks and dislocations, see, e.g. [19, 20]. Rice’s directly related J-integral [1] is defined as a boundary integral of S, i.e. , Z Z J= S N dA = W N − HT PN dA , (2) ∂Ω

∂Ω

where the Eshelby stress acts on N, the outward normal to the surface ∂Ω enclosing the region Ω in the reference configuration. The J-integral is commonly interpreted as the resultant the forces causing the configurational change [3]. Note that, given two regions Ω1 ⊂ Ω2 , JΩ2 \Ω1 = JΩ2 − JΩ1 (3) holds. Given an atomistic representation of a material, as opposed to the classical continuum, we need estimates of stored energy density W , the displacement gradient H, and the first Piola-Kirchhoff stress P derived from atomic data in order to construct the Eshelby stress field. Hardy projection [21, 22, 23] has the advantage of being consistent with the continuum balance laws by extension of Irving and Kirkwood’s seminal work [21] to continuous kernels. These kernels ψ(x), also called “localization functions”, have R the basic properties of being: (a) positive ψ > 0, and (b) normalized Ω ψ dV = 1. As in our previous work [16, 17], we have chosen to use a Lagrangian description [23] since it is natural for the material frame orientation of Eshelbian mechanics. The stored energy density field W can be defined as a local, weighted average of the potential energy per atom φα X X
φα (t) − φX W (X, t) = ψ(X − X ) = φα (t) ψ(X − Xα ) − W (X) (4) α α α

α

where φX α = φα ({Xβ }) is the potential energy in a fixed reference configuration {Xβ }, P X and W (X) = α φα ψ(X − Xα ). (Here and in the following, Greek letters denote

The application of an atomistic J-integral to a ductile crack

4

the label or index of a particular atom in the system.) Given the choice of reference configuration, the field W (X) is equal to the potential energy density for the zero temperature, perfect lattice configuration Xα and is such that W (X, t) = 0 implies P = 0, see [23]. The potential energy per atom φα is simply the sum the bond energies split equally between all atoms involved in the bond, c.f. [24] In order to define the displacement gradient H, we must first define the displacement u. This is done in a mass-weighted fashion as P (xα (t) − Xα ) mα ψ(Xα − X) , (5) u(X, t) = α P α α mα ψ(X − X)

with mα being the mass of atom α, in order to be consistent with the momentum, which is given a axiomatic definition in [21, 22, 23]. To efficiently represent all the continuous P fields, we employ a partition-of-unity interpolation NI | I NI = 1. In particular, the displacement field is given by X X u(X, t) = uI (t)NI (X) = NI ψIα uα , (6) I

I,α

where uI = u(XI , t) are from samples of the displacement u, Eqn. (5), on a grid of points XI . Here, ψIa = ψ(Xα −XI ). A differentiation of this interpolation leads directly to the estimated displacement gradient X H = ∇X u = uI (t)∇X NI (X) . (7) I

Finally, the referential first Piola-Kirchhoff stress [23] for central force potentials can be defined in terms of the force fαβ between atoms α and β and the difference between their positions in the reference configuration Xαβ = Xα − Xβ as X P(X, t) = − fαβ (t) ⊗ Xαβ Bαβ (X) (8) α<β

R1 Here, the so-called “bond” function Bαβ = 0 ψ (λ(Xα − X) + (1 − λ)(Xβ − X)) dλ is required for consistency with the continuum balance of momentum, see [22, 23], instead of a simple average of the virial weighted by ψIα . For further exposition, see [16, 17]. 3. Results We consider: (a) a cylindrical system with a single crack and boundary displacements given by an anisotropic K-field, and (b) a rectangular system with a finite width crack loaded in simple tension normal to the crack plane. For both cases, the materials simulated are fcc gold (Au) and bcc iron (Fe) using the Embedded Atom Method (EAM) potential by Foiles et al [25] and Simonelli et al [26, 18], respectively. The gold potential exhibits the following properties: intrinsic stacking fault energy γsf = 4.73 mJ/m2 ,

The application of an atomistic J-integral to a ductile crack

5

unstable stacking fault energy γus = 95.76 mJ/m2 , and surface energy for {110} surfaces (110) γs = 976 mJ/m2 . We employ a0 = 4.08 ˚ A as the zero pressure, zero temperature lattice parameter for gold. The iron potential has a zero pressure, zero temperature lattice parameter of a0 = 2.866 ˚ A and the properties of unstable stacking fault energy {110} for h111i{110} slip γus = 738 mJ/m2 , unstable stacking fault energy for h111i{112} {112} (110) slip γus = 863 mJ/m2 , and surface energy for {110} surfaces γs = 1.436 J/m2 . The simulations were performed using the molecular simulation code LAMMPS [27] using a sequence of energy minimizations to emulate a zero-temperature loading process. For the gold systems, cracks are pointing in the ±[00¯1] directions, have crack plane normals in the [110] direction, and have fronts along the [1¯10] direction, corresponding to the x, y and z directions, respectively. For this orientation, theoretically the first dislocation emitted is either a partial dislocation with Burgers vector [¯1¯12] on the inclined (111) plane, or a partial dislocation with Burgers vector [112] on the inclined (¯1¯11) plane. The initial crack was created by omitting the interactions between atoms on opposite sides of the crack plane and behind the crack tip. As such, no spontaneous crack closure is possible. According to the theory by Rice [5], this configuration possesses an 2 2 energetic cost for dislocation emission of Jnucl = 8γus /((1 + cos θ)  sin θ) = 1.27 J/m , using the value of γus already given and a value of θ = arccos √26 ≈ ±35.264◦ for the {111} planes listed above. By comparison, the energy cost of surface cleavage is (110) Jc = 2γs = 1.95 J/m2 . Thus, Rice’s criterion [5] predicts that dislocation emission is energetically favorable over surface cleavage since Jnucl < Jc . For the iron systems, we use one of the orientations examined by Shastry and Farkas [18]: cracks pointing in the ±[100] direction, crack plane normal in the [01¯1] direction, and crack front along the [011] direction, which correspond to the x, y and z directions, respectively. For this orientation, Shastry and Farkas observed emission of h111i-type edge dislocations from crack tips for their dual-tip system on the inclined (21¯1) and (¯21¯1) slip planes [18]. Interestingly, for this orientation (θ = ±54.73◦ ) the theory by Rice predicts Jnucl = 6.56 J/m2 , a value far greater than Jc = 2.87 J/m2 . Thus, Rice’s theory predicts that surface cleavage is energetically favorable over dislocation emission. However, dislocation emission occurs for the dual-tip crack examined in [18], and also herein as we shall show. 3.1. Semi-infinite crack under K-field loading Using the cylindrical systems shown in Figs. 1 and 2, we apply the near-crack-tip displacement solution from anisotropic elasticity (see Appendix) to an annulus of boundary atoms. For this simulation, only mode I loading is applied so that the vector of stress intensity factors is k = {0, KI , 0} in Eq. (A-1) which was applied in increments of 0.005 for the fcc system and 0.01 for the bcc system. These values of KI are normalized with respect to twice the equivalent isotropic shear modulus, 2µ, where µ = 39.7095 GPa for gold and 47.7115 GPa for iron. The system is periodic in the out-of-plane direction, which is 3a0 thick. The elements of the mesh used to represent the continuum fields are

The application of an atomistic J-integral to a ductile crack

6

aligned with the unit cells of the lattice. Figs. 1 and 2 illustrate the arrangement of integration loops used to evaluate the J-integral. Note that loop 1 through 7 increase in size and loop 0 does not surround the crack.

(a)

(b)

Figure 1. Configuration of the K loaded fcc system (a) and close-up of the crack tip (b). Atoms are colored by their potential energy, which is in units of eV/atom. The mesh and loop 2 are shown for scale. Loop 1 is 8 × 8 elements on a side, loop 2 is 16 × 16 elements, loop 3 is 20 × 20 elements, loop 4 is 24 × 24 elements, loop 5 is 28 × 28 elements, loop 6 is 32 × 32 elements, and loop 7 is 36 × 36 elements. Loop 0 is 8 × 8 elements and is not around the crack.

Figures 3 and 4 show that the J-integral converges with system size despite the fact that the emitted dislocations become pinned at the boundary, as shown in Figs. 1 and 2, and exert long-range elastic effects on the crack region. For the fcc system, prior and up to emission, the calculated J-integrals have excellent agreement with theory, Eq. (A-1). From Fig. 3 it is clear that the first dislocation emits at J1 = 1.25 J/m2 and the jump is ∆J1 = 0.29 J/m2 using loop 2. The chosen loading increment is small enough to discern the two emission events. Interestingly, the fcc system seems to recover the unperturbed trend in J1 with KI2 after both emission events. Direct examination of the atomic configuration shows that only during the emission events is there a significant change in the tip conformation, see Fig. 1(b). For the bcc system, prior to emission the calculated J-integrals have excellent agreement with theory, Eq. (A-1), but the first dislocation emits at J1 = 3.42 J/m2 and the jump is ∆J1 = 0.20 J/m2 using loop 4. This behavior is consistent with the findings in [18] where emission occurs before crack propagation and before the theoretical limit for nucleation. Also, unlike the fcc system, the second emission occurs at much higher KI and J1 = 3.78 J/m2 values than the first emission and is particularly size-dependent. We ascribe this difference in behavior to the significant blunting at the atomic scale in the bcc system and the more drastic atomic rearrangement relative to the fcc system

The application of an atomistic J-integral to a ductile crack

(a)

7

(b)

Figure 2. Configuration of the K loaded bcc system (a) and close-up of the crack tip (b). Atoms are colored by their potential energy, which is in units of eV/atom. The mesh and loop 2 are shown for scale. Loop 1 is 8 × 8 elements on a side, loop 2 is 16 × 16 elements, loop 3 is 20 × 20 elements, loop 4 is 24 × 24 elements, loop 5 is 28 × 28 elements, loop 6 is 32 × 32 elements, and loop 7 is 36 × 36 elements. Loop 0 is 8 × 8 elements and is not around the crack.

1.6

Ro=28 Ro=56 Ro=84 theory

1.4

J / Jnucl

1.2 1 0.8 0.6 0.4 0.2 0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

2

(KI/2µ)

Figure 3. Fcc system size-dependence of the calculated J-integral for loop 2 for the K loaded system. Ro is the radius of the system in multiples of the lattice constant a. Also note that the theoretical response is only valid pre-nucleation.

The application of an atomistic J-integral to a ductile crack 1

Ro = 28 Ro = 56 Ro = 84 theory

0.8

J/Jnucl

8

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

2

(KI/2µ)

Figure 4. Bcc system size-dependence of the calculated J-integral for loop 2 for the K loaded system. Ro is the radius of the system in multiples of the lattice constant a. Also note that the theoretical response is only valid pre-nucleation.

which remains relatively sharp, compare Fig. 1b to 2b. Perhaps this shows a limitation of Rice’s observation [1] that the actual geometry of the crack tip does not affect J and that after emission its value is solely determined by K. Figures 5 and 6 show that the J-integral is path-independent prior to the emission of the first dislocation, but displays a path-dependency once a dislocation has been nucleated. Within this post-nucleation regime, the insets in these figures clearly show a monotonic trend with loop size in both the fcc system, which has stacking faults, and the bcc system, which does not. By referencing differences to loop 2 (loop 1 is clearly very close to the crack tip), we are able to calculate a per length contribution to J, shown in Figs. 7 and 8, that converges away from the crack tip, using Eq. (3). We believe this size dependent perturbation to the J1 values is due to a change in the appropriate reference configuration as the dislocations pass through the lattice, and not a contribution due to stacking faults (which are only present in the fcc system). Interestingly, the values in the fcc and bcc systems are comparable despite the different elastic properties and dislocation slip paths. Alternatively, this size effect could be due to nonlinear shielding between the crack tip and the emitted, trapped dislocations. Given the findings in [28], our smallest system (Ro = 28) may be displaying significant shielding, but for our largest system (Ro = 84), the shielding should be negligible. Another possibility to explain the size dependence is the back stress due to the dislocation/s which are pinned at different distances from the crack in the three systems. Examining the P22 stress field in the bcc system, in particular, before and after the two distinct emission events, we see in Fig. 9,

The application of an atomistic J-integral to a ductile crack 1.8

loop 0 loop 1 loop 2 loop 3 loop 4 loop 5 loop 6 loop 7 theory

1.6 1.4

J / Jnucl

1.2 1

9

1.2

0.8 1

0.6

0.8

0.4 0.2

0.6

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

2

(KI/2µ)

Figure 5. Path-dependence of the calculated J-integral for the largest K loaded fcc system (Ro = 84). Inset shows the two emission events.

1

loop 0 loop 1 loop 2 loop 3 loop 4 loop 5 loop 6 loop 7 theory

J/Jnucl

0.8

0.6

0.6

0.4

0.5 0.2 0.4 0.9

0 0

0.2

0.4

0.6

0.8 2 (KI/2µ)

1

1 1.2

1.1 1.4

Figure 6. Path-dependence of the calculated J-integral for the largest K loaded bcc system (Ro = 84). Inset shows the first emission event.

The application of an atomistic J-integral to a ductile crack 0.002

loop 3-2 loop 4-2 loop 5-2 loop 6-2 loop 7-2

0 -0.002 J / ( Jnucl l )

10

-0.004 -0.006 -0.008 -0.01 -0.012 -0.014 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

2

(KI/2µ)

Figure 7. Contribution to J per length ℓ = l/a, were l is the length of the dislocation path included in the loop, for the K loaded fcc system.

0.001 0

J / ( Jnucl l )

-0.001 -0.002 -0.003 -0.004 loop 3-2 loop 4-2 loop 5-2 loop 6-2 loop 7-2

-0.005 -0.006 -0.007 0

0.2

0.4

0.6

0.8 1 2 (KI/2µ)

1.2

1.4

1.6

1.8

Figure 8. Contribution to J per length ℓ = l/a, were l is the length of the dislocation path included in the loop, for the K loaded bcc system.

The application of an atomistic J-integral to a ductile crack

11

for the largest system, that: (a) there is relaxation after emission, (b) the stress field in between the crack tip and the dislocation/s is smooth and approximately equal to the background stress everywhere, and (c) there are discrete artifacts in the stress field that reflect the dislocation slip paths. Given these observations, it seems that the back

(a) before 1st emission

(b) after 1st emission

(c) before 2nd emission

(d) after 2nd emission

Figure 9. P22 stress field for the loaded bcc system (Ro = 84) directly before and after emission of the first and second dislocations. P22 values are in units of GPa.

stress, which should appear as a smooth, long-range perturbation to the crack-tip stress field, is minimal; whereas, the discontinuity in the appropriate reference configuration appears as a line-like perturbation to the expected stress field. Exact values of several stress components pre- and post-nucleation are given in Table 1 and show roughly the

The application of an atomistic J-integral to a ductile crack

12

same trend as the 22 component depicted in Fig. 9. This table also shows that the same conclusions hold for the middle sized system (Ro = 56). Ro = 56 before 1st dislocation emission after 1st dislocation emission before 2nd dislocation emission after 2nd dislocation emission Ro = 84 before 1st dislocation emission after 1st dislocation emission before 2nd dislocation emission after 2nd dislocation emission

P11 15.49 14.08 13.49 13.37

P12 -1.45 -3.03 -4.34 -2.67

P21 -0.72 3.29 -0.34 -1.75

P22 15.67 15.22 10.18 8.83

P33 7.40 7.67 4.64 6.52

15.54 12.92 13.33 13.01

-1.30 -3.92 -4.36 -2.65

-0.66 0.62 -0.35 -1.78

16.05 11.42 10.14 8.51

6.97 4.51 4.15 5.69

Table 1. Stress field values for the node closest to the crack-tip for the bcc system (Ro = 56, 84) before and after emission of the first two nucleated dislocations. Values given are in units of GPa.

3.2. Finite width crack in a rectangular plate under tension For this example, the atomic system is rectangular with a finite width crack at its center, refer to Fig. 10. System width in the x-direction is 2W , length in the y-direction is 2L, thickness in the z-direction is T . In all of our simulations, the systems are constructed of √ nx ×ny ×nz unit cells of dimensions ℓx ×ℓy ×ℓz , where ℓx = a0 , ℓy = ℓz = a0 2, and nx and √ A. ny are chosen such that L/W ≈ 0.99, nz = 3 and the thickness is T = 3a0 2 = 23.08 ˚ For one set of simulations, the crack’s width is 2ℓc where ℓc ∼ 4a0 = 16.32 ˚ A. Thus, the crack width to system width ratio (ℓc /W ) decreases with increasing system size. For a second set of simulations, the crack’s width is scaled such that ℓc /W = 0.2 for all systems. Periodic boundary conditions are used in the x- and z-directions, while fixed conditions are used in the y-direction to both apply a specified amount of strain and prevent dislocations from piercing the upper and lower boundaries. A strain increment of 0.02% is applied in the y-direction, with energy minimizations performed between application of successive increments. Figure 10 shows one such system, where nx = 40, ny = 28, and ℓc /a0 = 4, at various levels of strain and dislocation emission. Sub-figures (a)-(c) show the atoms colored according to a parameter calculated via common neighbor analysis (CNA), with green corresponding to atoms in fcc structured environments; red atoms corresponding to hexagonal close packed (hcp), i.e. atoms bordering a stacking fault; and white used for atoms at the free surfaces of the crack. These sub-figures show that partial dislocation emission occurs from both crack tips at a strain of 2.10%, leaving behind stacking fault trails. Then, at a higher strain of 2.33%, a second set of partial dislocations is emitted,

The application of an atomistic J-integral to a ductile crack

13

creating two more stacking faults. These qualitative features were present in all systems examined.

(a) 0% strain

(b) 2.10% strain

(c) 2.33% strain

(d) 2.08% strain

(e) 2.10% strain

(f) 2.33% strain

(g) 2.08% strain

(h) 2.10% strain

(i) 2.33% strain

√ Figure 10. Rectangular plate with a finite width crack (W/a0 = 40, L/a0 = 28 2, ℓc /a0 = 4). Figures (a)-(c) show per-atom common neighbor analysis (green:fcc, red:hcp, white:free surface. Figures (d)-(f) show per-atom potential energy (red: high, blue: low). Figures (g)-(i) show the 11 displacement gradient field (blue: high, red: low).

Sub-figures (d)-(f) show the system with atoms colored according to their value of potential energy, while sub-figures (g)-(i) show the overlapping FE mesh used for evaluating continuum fields of displacement gradient (H), 1st Piola-Kirchhoff stress

The application of an atomistic J-integral to a ductile crack

14

(P), referential energy density (W ) and Eshelby stress (S). Elements are of the same dimension as the unit cells, ℓx × ℓy × ℓz , and the same boundary conditions are used on the mesh as on the atomic system. Shown in Figure 10 are elements colored according to their interpolated value of the (1,1)-component of the displacement gradient tensor. These images clearly show discontinuities in H11 due to dislocation nucleation and slip having occurred. These continuum fields were integrated along rectangular loops surrounding the left (loops 1, 2 and 3) and right (loop 4) crack tips in order to quantify the J-integral at each loading step. We observe that J varies nearly linearly with the square of the applied strain, and that (as for the cylindrical crack example) J is path-independent prior to the first dislocation emission event, and even relatively path-independent prior to the second dislocation emission event. At both events, the value of J jumps discontinuously to higher values after the dislocations have been emitted and resulting stacking fault trails have been created. Figure 11 shows the J curves obtained by three different system sizes from the first set of simulations, i.e. ℓc /a0 = 4. Although there is some variation in the critical value of J at which the jump occurs, these values are fairly close to one another with no discernible pattern to their dependency on system size. It is the case, however, that the magnitude of the jump in J appears to be increasing as system size increases. 1

J / Jnucl

0.8

0.6

0.4

0.2

20x14 40x28 60x42

0 0

1

2

3

4 2

STRAIN x 10

5

6

7

8

4

Figure 11. Calculated J-integral for finite width systems with a fixed crack width (ℓc /a0 = 4).

By comparison, the second set of simulations (ℓc /W = 0.2) shown in Figure 12 displays different patterns. As system size increases, the critical value of J at which the first dislocation emission event occurs appears to increase, while the subsequent magnitude of the jump in J appears to be relatively constant.

The application of an atomistic J-integral to a ductile crack

15

1

0.8

J / Jnucl

0.6

0.4 20x14 40x28 60x42 80x56

0.2

0 0

1

2

3

4

5

2

STRAIN x 10

6

7

8

4

Figure 12. Calculated J-integral for finite width systems with similar crack widths (ℓc /W = 0.2).

The values of J just prior to, and after the emission of the first and second sets of dislocations are given in Table 2 for all fcc systems simulated. We note that as W/a0 20 40 60 80 20 40 60

L/a0 √ 14 2 √ 28 2 √ 42 2 √ 56 2 √ 14 2 √ 28 2 √ 42 2

L/W 0.99 0.99 0.99 0.99 0.99 0.99 0.99

ℓc /a0 4 8 12 16 4 4 4

ℓc /W 0.2 0.2 0.2 0.2 0.2 0.1 0.067

ǫ1nucl

st

Jpre-nucl

∆J

[%]

[J/m2 ]

2.39 1.67 1.39 1.23 2.39 2.10 2.06

0.890 0.928 0.960 1.009 0.890 0.855 0.868

nd

ǫ2nucl

Jpre-nucl

∆J

[J/m2 ]

[%]

[J/m2 ]

[J/m2 ]

0.067 0.075 0.064 0.076 0.067 0.191 0.249

2.72 1.86 1.51 1.33 2.72 2.33 2.29

1.296 1.277 1.262 1.309 1.296 1.336 1.422

0.149 0.092 0.078 0.070 0.149 0.222 0.302

Table 2. Values of J prior to emission of first and second dislocations for the various systems studied. By comparison, Rice’s theory predicts Jnucl = 1.27 J/m2 .

the system increases in its overall dimensions, the amount of strain needed to induce dislocation emission decreases. This reduction of strain required for dislocation emission with increasing system size is more dramatic as the crack width scales with system size (ℓc /W = 0.2) than for the case of a constant crack width (ℓc /a0 = 4), indicating that emission is strongly affected by both the periodic boundary condition and the crack tip interaction. This scaling of the atomic system had a much clearer trend towards convergence than the fixed crack width study shown in Fig. 11.

The application of an atomistic J-integral to a ductile crack

16

Regarding the nature of the crack-dislocation interactions for periodic systems, there can be nonlinear interactions but only for very small sized systems. At larger sizes, the interaction consists of a linear superposition of periodically-spaced crack and dislocation stress fields, and in the limit of infinite system width this interaction would go to zero. Our results given in Table 2 show these effects as the value of J at nucleation slowly approaches its theoretical value as system size is increased. We believe that for the smallest size studied (W = L = 20a0 ), our system is sufficiently large to be in the linear superposition regime; however, a more thorough future study is warranted to quantify the size at which nonlinear interactions prevail. The results for the bcc system were qualitatively similar and are omitted for brevity. 4. Discussion As in our previous work [16, 17], we see good agreement with elasticity theory down to the nanoscale given appropriate loading and/or a large enough system. In particular, for the fcc cylindrical system the calculated J value at emission was within 2% of the theoretical prediction and in the self-similar rectangular system, where the crack tips become increasingly separated, there is a clear trend in the calculated J towards the theoretical Jnucl . For the bcc cylindrical system the trend of J with loading follows theory but the emission of dislocations at J levels below the theoretical predictions for propagation or nucleations is consistent with [18]. Note that we tried to allow the dislocations to leave the explicitly modeled system and thus remove their long-range elastic effects on crack by fixing the interior core and allowing the boundary to relax, but this approach was unsuccessful. At present, we do not have a good explanation of the magnitude or sign of the jump in J at emission. As mentioned previously, there is a loss of path independence of J-integral we believe is due to change of the reference configuration that is appropriate for elastic response. This perturbation of the J values is loop size dependent and we showed that this per-length contribution converges to a well-defined value far enough away from the crack tip. For both corresponding fcc and bcc systems, the Burger’s vectors of corresponding dislocations in the cylindrical and rectangular systems are the same and yet the jump in the cylindrical system is negative and shows recovery to the perfect lattice trend between emission events and the jump in the rectangular system is positive and does not display recovery. We calculated the Peach-Kohler force on corresponding dislocations in both fcc systems and the x component is approximately 0.08 J/m2 , which is about a third of the jump in the cylindrical system and on the order of the jump in the rectangular system. So, clearly the jump in J cannot be simply ascribed to the J value associated with the dislocation that has left the integration loop. Finally, in the Introduction we posed the question as to whether a stacking fault provides a non-zero contribution to the J-integral. Our analysis of the cylindrical crack for the fcc and bcc systems suggests that it does not, at least not to any

The application of an atomistic J-integral to a ductile crack

17

significant degree. Our method for estimating the J-integral clearly shows a loss of path independence and non-zero, negative contributions to J after dislocation emission has occurred. However, this behavior is common to both the fcc (which produces stacking faults) and bcc (which does not) systems, with the magnitude of the contribution to J similar for both systems. This result suggests that it is the large, discontinuous displacements relative to the reference configuration characteristic of the slip process that are responsible for the non-zero values of J, rather than the presence of the stacking fault itself. More analysis in this area is clearly warranted, as our results are suggestive, but not conclusive. Acknowledgments We thank S. A. Lurie (Institute of Applied Mechanics, Russian Academy of Sciences) for suggesting this avenue of research and Prof. D. Farkas (Virginia Tech) for supplying the Fe EAM potential parameterization. This work was supported by the Engineering Science Research Foundation (ESRF) at Sandia National Laboratories. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. References [1] J. R. Rice. A path independent integral and approximate analysis of strain concentration by notches and cracks. J. Appl. Mech., 35(2):379–386, 1968. [2] J. D. Eshelby. The force on an elastic singularity. Phil. Trans. Roy. Soc. Lon. A, 244(877):87–112, 1951. [3] J.D. Eshelby. The elastic energy-momentum tensor. Journal of Elasticity, 5(3-4):321–35, 1975. [4] Y. Sun, J. R. Rice, and L. Truskinovsky. Dislocation nucleation versus cleavage in ni3 al and ni. In High-Temperature Ordered Intermetallic Alloys, volume 213 of Materials Research Society Symposium Proceedings, pages 243–248. Materials Research Society, 1991. [5] J. R. Rice. Dislocation nucleation from a crack tip: An analysis based on the Peierls concept. J. Mech. Phys. Solids, 40(2):239–271, 1992. [6] Y. Sun, G. E. Beltz, and J. R. Rice. Estimates from atomic models of tension-shear coupling in dislocation nucleation from a crack tip. Matls. Sci. and Eng. A, 170:67–85, 1993. [7] J. R. Rice and G. E. Beltz. The activation energy for dislocation nucleation at a crack tip. J. Mech. Phys. Solids, 42(2):333–360, 1994. [8] Y. Sun and G. E. Beltz. Dislocation nucleation from a crack tip: A formulation based on anisotropic elasticity. J. Mech. Phys. Solids, 42(12):1905–1932, 1994. [9] L. Hung and E. A. Carter. Ductile processes at aluminum crack tips: comparison of orbital-free density functional theory with classical potential predictions. Modelling Simul. Mater. Sci. Eng., 19(4):045002, 2011. [10] F. Cleri, D. Wolf, S. Yip, and S. R. Phillpot. Atomistic simulation of dislocation nucleation and motion from a crack tip. Acta Mater., 45(12):4993–5003, 1997. [11] M. X. Shi, Y. Huang, and H. Gao. The J-integral and geometrically necessary dislocations in nonuniform plastic deformation. Intl. J. Plasticity, 20:1739–1762, 2004.

The application of an atomistic J-integral to a ductile crack

18

[12] F. Cleri, S. Yip, D. Wolf, and S. R. Phillpot. Atomic-scale mechanism of crack-tip plasticity: Dislocation nucleation and crack-tip shielding. Phys. Rev. Letters, 79(7):1309–1312, 1997. [13] T. K. Bhandakkar, A. C. Chng, W. A. Curtin, and H. Gao. Dislocation shielding of a cohesive crack. J. Mech. Phys. Solids, 58:530–541, 2010. [14] G. E. Beltz, D. M. Lipkin, and L. L. Fischer. Role of crack blunting in ductile versus brittle response of crystalline materials. Phys. Rev. Letters, 82(22):4468–4471, 1999. [15] J. Schiotz and A. E. Carlsson. The influence of surface stress on dislocation emission from sharp and blunt crack in fcc metals. Phil. Mag. A, 80(1):69–82, 2000. [16] R E Jones and J A Zimmerman. The construction and application of an atomistic J-integral via Hardy estimates of continuum fields. J. Mech. Phys. Solids, 58(9):1318–1337, 2010. [17] R E Jones, J A Zimmerman, J Oswald, and T Belytschko. An atomistic J-integral at finite temperature based on Hardy estimates of continuum fields. J. Phys.: Cond. Matter, 23(1):015002, 2011. [18] V. Shastry and D. Farkas. Molecular statics simulation of fracture in α-iron. Modelling and Simulation in Materials Science and Engineering, 4:473–492, 1996. [19] G. A. Maugin. Material Inhomogeneities in Elasticity. CRC press, Boca Raton, FL, 1993. [20] M. E. Gurtin. Configurational forces as basic concepts of continuum physics. Springer, New York, NK, 2000. [21] J.H. Irving and J.G. Kirkwood. The statistical mechanical theory of transport processes. IV. the equations of hydrodynamics. J. Chem. Phys., 18:817–829, 1950. [22] R J Hardy. Formulas for determining local properties in molecular-dynamics simulations: Shock waves. Journal of Chemical Physics, 76(1):622–628, 1982. [23] J. A. Zimmerman, R. E. Jones, and J. A. Templeton. A material frame approach for evaluating continuum variables in atomistic simulations. J. Comp. Phys., 229(6):2364–2389, 2010. [24] E.B. Webb, III, J.A. Zimmerman, and S.C. Seel. Reconsideration of continuum thermomechanical quantities in atomic scale simulations. Mathematics and Mechanics of Solids, 13(3-4):221–66, 2008. [25] S M Foiles, M I Baskes, and M S Daw. Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Phys. Rev. B, 33:7983–7991, 1986. [26] G. Simonelli, R. Pasianot, and E. J. Savino. Embedded-atom-method interatomic potentials for BCC iron. In Materials Theory and Modeling, volume 291 of Materials Research Society Symposium Proceedings, pages 567–572. Materials Research Society, 1993. [27] LAMMPS : Large-scale Atom/Molecular Massively Parallel Simulator, Sandia National Laboratories: http://lammps.sandia.gov, 2009. [28] J. Song, W. A. Curtin, T. K. Bhandakkar, and H. J. Gao Dislocation shielding and crack tip decohesion at the atomic scale. Acta Mater., 58:5933-5940, 2010. [29] T. T. C. Ting. Anisotropic Elasticity: Theory and Applications. Oxford University Press, New York, New York, 1996. [30] G. Hattori, R. Rojas-D´ıaz, A. S´aez, N. Sukumar, and F. Garc´ıa-S´anchez. New anisotropic cracktip enrichment functions for the extended finite element method. Computational Mechanics, 2012.

Appendix A. Anisotropic linear elastic fracture solution The theoretical displacement solution [29, 9, 30] is given by r
2 Re ADB−1 k u= π

(A-1)

Herein, Re indicates the real part of the quantity in parentheses; and likewise Im the √  imaginary part. The matrix D = diag zα is a 3 × 3 diagonal matrix with elements

The application of an atomistic J-integral to a ductile crack 19 √ √ zα = x + pα y. Here, pα | Im pα > 0 are the three of the six complex conjugate roots of   det [M(p)] = 0 | [M(p)]ik = Ci1k1 + p (Ci1k2 + Ci2k1 ) + p2 Ci2k2 (A-2)

which is sixth-order in p and can be considered a generalized eigenvalue problem. Here, Cijkl are the components of the fourth-order elastic stiffness tensor of the material under consideration. Given the x, y and z orientations discussed in the beginning of Section 3, C has the following non-zero elements:

C2222 C1122 = C2211 = C1133 C2233 C1212 = C2121 = C1313 C2323

C1111 = c11 1 = C3333 = (c11 + c12 + 2c44 ) 2 = C3311 = c12 1 = C3322 = (c11 + c12 − 2c44 ) 2 = C3131 = c44 1 = C3232 = (c11 − c12 ) 2

Here, c11 , c12 and c44 are the elastic constants for a cubic material with the x, y and z axes in the h100i crystal directions. Using the eigenvalues pα of equation (A-2), three matrices can be formed Mα = M(pα ). Of the eigenvectors of each of Mα only one, aα , satisfies M(pα )aα = 0. The columns of the 3×3 matrix A is then constructed from each of these null-space eigenvectors as A = [a1 a2 a3 ]. Similarly, the matrix B = [b1 b2 b3 ] is composed of the complementary vectors bα = Zα aα where [Zα ]ik = Ci2k1 +pα Ci2k2 . In general, the vector k = {KII , KI , KIII } is composed of stress intensity factors applied far from the crack tip. Finally, the crack tip energetic driving force is given by
1 J1 = kT Re ıAB−1 k. 2

(A-3)