Simulations of reactions between irradiation induced 〈a〉 -loops and mixed dislocation lines in zirconium

Journal of Nuclear Materials 462 (2015) 126–134

Contents lists available at ScienceDirect

Journal of Nuclear Materials journal homepage: www.elsevier.com/locate/jnucmat

Simulations of reactions between irradiation induced hai-loops and mixed dislocation lines in zirconium K. Ghavam, R. Gracie ⇑ Civil and Environmental Engineering, Department University of Waterloo, Waterloo, Canada

a r t i c l e

i n f o

Article history: Received 30 May 2014 Accepted 3 March 2015 Available online 19 March 2015

a b s t r a c t The interaction between mixed dislocations gliding in the basal plane with elliptical prism plane vacancy loops characteristic of irradiation induce loops is studied using Molecular Dynamics simulations. A methodology to create elliptical loops in the prism plane with major axes in the ½0 0 0 1 direction and with 21  0 direction is described. Dislocation reactions involving mixed dislocations, Burgers vectors in the ½1  21  0 with a Burgers vector of a=3½2  1 1 0 and gliding in the basal plane in the ½1 0 1  0 direcparallel to ½1  1 0 0Þ or tion are investigated. The reaction of the mixed dislocation and vacancy loops in prism planes ð1  0 1 0Þ and with Burgers vector a=3½1  21  0 leads to the creation of two jogs in the vacancy loop. These ð1 loops present a significant barrier to the motion of the mixed dislocation gliding in the basal plane.  0Þ prism plane with Burgers vector The reaction of the mixed dislocation and a vacancy loop in the ð1 0 1  1 1 0 leads to the cleaving of the vacancy loop and the absorption of part of the vacancy loop into a=3½2 the gliding dislocation. This loop presents a relatively weaker barrier to the motion of the mixed dislocation gliding in the basal plane. Strain rate is observed to have a large impact on the shear stress-displacement behavior of the reactions, but not on the reaction mechanisms. Crown Copyright Ó 2015 Published by Elsevier B.V. All rights reserved.

1. Introduction Damage caused by irradiation has long been known to have a significant effect on the plastic behavior of zirconium (Zr) and Zr-alloys used in the nuclear industry; however, to date there have been few studies at the atomistic scale, which elucidate the reaction mechanisms of irradiation induced dislocation loops with glissile dislocations. In this article, the results of Molecular Dynamics simulations are shown the illuminate the reaction mechanisms between glissile dislocations in the basal plane and irradiation induced vacancy loops in the prism planes. Non-irradiated zirconium deforms primarily by prism slip of hai-type dislocations and twinning [1]. At temperatures above 623 K, basal slip of hai-type dislocations is also present to a limited degree. Slip on pyramidal planes is less common, even at elevated temperatures. Prism slip occurs more favorable than basal slip because the stacking fault energy on the prism plane is lower than that on the basal plane [1,2]. During irradiation, high energy neutrons impact surface atoms causing cascading collisions and these cascades lead to the generation of vacancy and interstitial defects and defect clusters. These defects evolve at elevated operational temperatures to form ⇑ Corresponding author. E-mail address: [email protected] (R. Gracie). http://dx.doi.org/10.1016/j.jnucmat.2015.03.007 0022-3115/Crown Copyright Ó 2015 Published by Elsevier B.V. All rights reserved.

vacancy and interstitial dislocation loops [3]. Both irradiation induced vacancy and interstitial loops in the prism planes with  0i have been observed [3,4]. Based on Burgers vectors of a=3h1 1 2 the TEM observations, these loops are elliptical in shape with the major axis along the hci-axis and the minor axis residing in the basal plane [3]. The ratio of the lengths of the minor and major axes (b/a) for the interstitial loops is about 0.94 and is independent of the size of the loop. For vacancy loops, the ratio of b/a decreases from 0.8 for small loops (a = 10 nm) to a nearly constant value of 0.6 for loops with a > 100 nm. Kelly and Blake [5] reported a slightly different result, they observed dislocation loops with  0i, but they reported that the normal to Burgers vector a=3h1 1 2    1 direction, which the plane of the dislocation is close to the 2 0 2 is not reported in other experimental works. According to Kelly’s and Blake’s TEM observations, two-third of the loops were vacancy loops and the remainder were interstitial. They obtained these results for neutron irradiated zirconium after annealing at 725 K. The mechanical behavior of irradiated zirconium is significantly different than non-irradiated zirconium. The primary deformation mechanism is basal slip by dislocations of hai-type, rather than prism slip. Furthermore, while plastic deformation is relatively uniform in non-irradiated zirconium, localized deformation is observed in irradiated zirconium [6,7]. MD simulations have been extensively used to study the effects of irradiation on zirconium; however, to date MD simulations have

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

largely focused on the early effects of irradiation such as, displacement cascades, point defect diffusion, and defect cluster formation [8–14]. MD has been used to study the interaction of dislocations with small cascade induced defect clusters, consisting of less than 50 interstitials or vacancies [15–17]. These studies, conducted prior to 2007, used the 1995 potential developed by Ackland et al. [18]. In 2007, a more suitable atomistic potential for the simulation of Zr was developed by Mendelev and Ackland [19], the MA07 potential. It has been shown to more accurately estimate dislocation core structures, stacking fault energies, and the dynamic behavior of dislocations in the prism and basal planes than previous potentials [20]. Recently, Serra and Bacon [21] have used the MA07 potential to study the reaction of edge and screw dislocations gliding in prism planes with small irradiation induced interstitial loops, composed of up to 156 interstitials. MD has also been used to estimate stacking fault and dislocation core energies in other nuclear material, such a uranium oxide [22]. This article presents the first study of mixed dislocations gliding in the basal plane interacting with large elliptical vacancy loops in the prism planes, characteristic of irradiation induced loops. These dislocation reactions have not yet been studied in detail using MD. Section 2 describes the three cases of dislocation reactions between glissile dislocations in the basal plane with irradiation induced vacancy loops in the prism plane that will be studies. Section 3 describes a methodology to generate real-sized vacancy loops (a P 10 nm) in the prism planes, without the need to simulate displacement cascades nor defect diffusion processes. The MD results and discussion of the three cases are presented in Section 4 and Section 5, including a discussion on strain rate effects. Finally, in Section 6 the conclusions of this study are presented. 2. Model description Molecular Dynamics (MD) simulation is used to study the reaction mechanisms of glissile dislocations interacting with irradiation induced vacancy loops. In MD, since the forces acting on the atoms are derived from a potential, the selection of an appropriate potential and understanding its limitations is important. The Embedded Atom Method (EAM) potential proposed by Mendelev and Ackland (MA07) [19] is used. This potential is more accurate than previously developed potentials and accurately reproduces many properties of zirconium, including stacking fault energies and the structures of dislocation cores of dislocations in both the prism and basal planes. In this paper, the MD simulator LAMMPS has been used. The simulations were performed at 0.01 K using a micro-canonical NVE ensemble and a constant time step of 1 fs. The MA07 potential predicts that a perfect Zr HCP crystal will have lattice constants pffiffiffi a ¼ 3:2339 Å; b ¼ 3a, and c ¼ 1:5979a. We evaluated the stacking fault energy landscape for both the prism and basal planes using LAMMPS. For the basal plane, we found that the stacking fault energy is 12.4 meV/Å2 for a relative displacement of basal planes pffiffiffi    1 0 direction, which is the same as of a= 3 ¼ 1:85 Å in the 1 0 previously published [19,20]. For the prism plane, we found that the stacking fault energy was 8.4 meV/Å2 for a relative displacement of 0:5hai þ 0:14hci. This result is consistent with that reported in [20].1 Verification of these fundamental properties of Zr, was important for instilling confidence in the results reported later in this article. In this article, three different cases of the interaction of a dislocation gliding in a basal plane with an irradiation induced vacancy loop in the prism plane are studied. In this initial work, 1

In Ref. [19], the obtained value is 10.87 meV/Å2.

127

 21  0 direction with a only glissile dislocations parallel to the ½1  1 1 0 and traveling in the ½1  0 1 0 direction Burgers vector of a=3½2 are considered. The interaction of this glissile dislocation with  1 0 0g planes is considered. In each case, vacancy loops in the f1 the vacancy loops are elliptical with an aspect ratio (minor to major) of 0.8. The major axis of the loops has a length of 20 nm and is parallel to the hci-direction. Dislocation reactions are simulated in box-shaped domains, 21  0; ½1  0 1 0, and ½0 0 0 1 are such that the lattice directions ½1 orientated parallel to the x-, y-, and z-axes, respectively. The domain is 80  80  90 lattice units and is assumed to be periodic in the x- and y-directions. The bottom of the domain, z ¼ 0, is fixed in the z-direction. A shear velocity boundary condition, in the y-direction and with a magnitude of v y is applied to the top of the domain at z ¼ 90 lattice units. The initial velocity is applied to each atom, increasing linearly from zero at z ¼ 0 to a value of v y are z ¼ 90 lattice units. The following cases will be studied and are illustrated in Fig. 1. Case 1 The reaction of a vacancy  1 0 0Þ with Burgers vector ð1 dislocation. Case 2 The reaction of a vacancy  0 1 0Þ with Burgers vector ð1 dislocation. Case 3 The reaction of a vacancy  0Þ with Burgers vector ð1 0 1 dislocation.

loop in the prism plane  21  0 with the mixed a=3½1 loop in the prism plane  21  0 with the mixed a=3½1 loop in the prism plane  1 1 0 with the mixed a=3½2

In all cases, the domain size was selected to be a large as possible and in such a way as to limit the effect of the size of the simulation domain on the results. Simulations were carried out using parallel computations on 16–32 CPUs on SHARCNET.2

3. Model generation The generation of the initial equilibrium configuration of a simulation cell containing a vacancy loop characteristic of irradiation induced loops and a glissile dislocation in an HCP crystal is a non-trivial task. In nature, the vacancy loops are created by diffusion of vacancies generated by displacement cascades. The vacancies migrate to form voids, that once sufficiently large, collapse into dislocation loops. The simulation of the whole process from cascade, to the coalescence of vacancies into voids and the collapse of these voids into vacancy loops is not computationally feasibly, especially for large loops. While studies of dislocation interactions with irradiation induced loops in HCP crystals have been limited, there have been several studies concerning the interaction between the gliding dislocation lines and irradiation induced loops in BCC materials [23–25]. The methodology used in these earlier works has influenced the development of the methodologies used in this article to generate vacancy loops and mixed dislocation lines in HCP Zr. The methodology used here is somewhat different, in part, because the sizes of the irradiation-induced loops in this study are about ten times bigger and the dislocation lines are mixed. In this section, the methodologies used to avoid the explicit simulation of vacancy loop formation from the processes of displacement cascades, diffusion, and vacancy cluster formation and collapse are presented, specifically algorithms to (a) create a vacancy loop in the prism planes of a Zr crystal, (b) create a mixed dislocation in the basal plane of a Zr crystal, and (c) create both a 2

https://www.sharcnet.ca.

128

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

y [1 0 1 0]

y [1 0 1 0]

y [1 0 1 0] t2

b2

Vy

t2

a3

t2

b2

b2

t2

t2

a2 a1

Y0 Y0-b/6

b1

a/2

t1

t1

b1=a/3[2 1 1 0] Case 1

b1

x [1 2 1 0]

t1 Case 2

x [1 2 1 0]

Case 3

x [1 2 1 0]

 21  0 Fig. 1. A schematic of the initial and boundary conditions for cases 1–3 (left, middle and right). Illustration of the methology to create a mixed dislocation parallel to ½1  1 1 0 (left) and the division of a simulation box into two subdomains (red dashed boxes) one containing the mixed dislocation in the basal plane and with Burgers vector ½2 and the other containing the vacancy loop in the prism plane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

vacancy loop in the prism plane and a mixed dislocation in the basal plane of a Zr crystal. The algorithms were implemented into a MATLAB preprocessor, which creates LAMMPS input files with the necessary commands. 3.1. Methodology to generate vacancy loops This subsection describes the method used to generate an elliptical shaped vacancy loop in the prism planes, fag-planes, with a major axis in the hci-direction. The Burgers vectors of the generated loops are chosen to be in an hai-direction, such that the normal of the prism plane is oriented 30-degrees from the Burgers vector. The first step to generate a vacancy loop in an HCP lattice is to create a perfect HCP crystal, at least twice as large as the dimensions of the loop to be created. The second step is to remove all atoms from the simulation box that are contained in an elliptical shaped disk. The major axis of the disk should be parallel with the hci-axis and its length should be equal to the size of the loop that is to be created. The disk should be orientated so that its normal is parallel to the normal of the prism plane in which a vacancy loop is to be created. The aspect ratio of the disk should be the same as that of the vacancy loop to be created. The thickness of the disk is equal to the thickness of two adjacent layers of atoms in the prism plane. The third step is to collapse the void created in the previous step into an elliptical dislocation. This is accomplished by a constrained minimization of the energy of the remaining atoms in the crystal. The lattice is minimized under periodic boundary conditions and atoms are constrained to move only in the direction of Burgers vector. This will allow the atoms adjacent to the void to fill the void. Since the atoms are free to move in the basal plane, the atoms will 21  0 or the a=3½2  1 1 0 direction. tend to move in either the a=3½1 Thus, constraining the displacements of the atoms during the minimization of the lattice energy is extremely important. The last step of the methodology is to further minimize the energy of the lattice without any constraints. 3.2. Methodology to generate mixed glissile dislocations To study the reaction between glissile dislocations on the basal plane and vacancy loops typical of irradiation defects, a method was developed to create mixed glissile dislocations parallel to 21  0 and with Burgers vector equal to a=3½2  1 1 0. Edge disloca½1   tions parallel to ½1 2 1 0 and with a Burgers vector in the

pffiffiffi  0 1 0 direction are not considered, as these dislocations will 3a½1 disassociate into partial dislocations with Burgers vectors 1  2 0 and a=3½2  1 1 0. It is assumed that the methodology will a=3½1 be used to create mixed dislocations in lattices, where the lattice 21  0; ½1  0 1 0, and ½0 0 0 1 are parallel to the x-, y-, directions ½1 and z-axes, respectively.  21  0 with Burgers vector To create a dislocation parallel to ½1  a=3½2 1 1 0, first create a simulation box, where the lattice direc 21  0; ½1  0 1 0, and ½0 0 0 1 are parallel to the x-, y-, and tions ½1 z-axes, respectively, as illustrated in Fig. 1-left. Second, remove  0 1 0Þ that are in atoms from two adjacent planes parallel to ð1 the bottom half of the simulations box. The y-coordinates of these pffiffiffi atoms are y0 ¼ nb and y1 ¼ y0  b=6, where b ¼ 3a is the lattice unit cell length in the y-direction and n is a natural number (Fig. 1-left). Third, the atoms in the bottom half of the simulation box, with y < y1 are displaced in the x-direction a distance proportional to their y-coordinate; atoms at y ¼ y1 are displaced a distance of a=2, while atoms at y ¼ 0 are not displaced. This step is illustrated in Fig. 1-left. Fourth, the energy of the lattice is minimized under the constraint that displacements in the z-direction are zero and the boundaries are periodic in both the x- and y-directions. Following the above procedure, a mixed dislocation parallel to  21  0 and with Burgers vector a=3½2  1 1 0 is created. A dislocation ½1 with the opposite orientation, but with the same Burgers vector can be created by modifying the second step to remove atoms from the top half of the simulation box. Furthermore, dislocations 21  0 and with Burgers vector a=3½1 1  2 0 can also be parallel to ½1 created by appropriately modifying the third step. 3.3. Methodology to create a vacancy loop and a glissile mixed dislocation Additional model creation steps are required to create a simulation box with both a vacancy loop in the prism plane and a mixed dislocation in the basal plane. As the Burgers vectors of these dislocations are generally different, they have to be minimized under two different constraints. The domain is divided into two parts, in such a way that each subdomain will entirely contain either the vacancy loop or the glissile mixed dislocation; the division of the domain into two parts is illustrated in Fig. 1-left. To create a simulation box with both dislocations, the energy of the simulation box is minimized in three steps; first, the energy of the subdomain containing the mixed dislocation (red dashed box

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

on top in Fig. 1-left) is minimized following the methodology described in Section 3.2. Next, the energy of the subdomain containing the vacancy loop (red dashed box at bottom in Fig. 1-left) is minimized, following the methodology described in Section 3.1. Lastly, the energy of the whole system is minimized. 4. Results In this section, the dislocation reactions observed in the simulation of Cases 1–3 are discussed. In addition, the impact of the loading rate on the observed stress–strain curves is also discussed. In order to study the effects of loading rate, we considered five different values for the applied velocity v y : 0.56 Å/ps, 0.28 Å/ps, 0.056 Å/ps, 0.028 Å/ps, and 0.0056 Å/ps; the corresponding strain rates are 1:2  109 s1 , 6  108 s1 , 1:2  108 s1 , 6  107 s1 , and 1:2  107 s1 , respectively. The shear strain of the simulation box v is taken as yz ¼ Lzy t, where t is time and Lz is the height of the mobile part of the simulation box. Consequently, the strain rate v reported for each simulation is computed as _ ¼ Lzy . Here the stress of the simulation box is computed by dividing the sum of the symmetric per-atom virial stress tensor of each atom3 by the volume of the simulation box. It should be noted that the strain rates in the MD simulations are much higher than that used in regular mechanical testing which often prevents a quantitative comparison, though a qualitative comparison is still possible. 4.1. Reaction mechanisms of irradiation induced loops In this subsection, the reaction mechanism observed in the three cases are described. To simplify the discussion, let LN denote the mixed dislocation line gliding on the basal plane and let LP denote the vacancy loop. Due to the applied shear velocity on the upper surface of the domain, LN glides towards LP. The reaction mechanism is similar for different loading rates, and so we will discuss the mechanism with reference to simulation results obtained for an applied velocity of v y ¼ 0:028 Å=ps (strain rate of 6  107 s1 ).  1 0 0Þ prism plane 4.1.1. Case 1: vacancy loop in the ð1 Fig. 2 is a schematic of the reaction observed. Snap-shots of the actual MD simulations of the reaction are shown in Figs. 3 and 4. As LN nears LP, there is an observable spreading of the dislocation core of LP. There is also an observable spreading of the core of dislocation LN in the basal plane to enable the creation of the jog in LP, see Fig. 3c. During the reaction, LN cuts LP, leaving behind two jogs in LP in the basal plane both with magnitudes of a=3 in  1 1 0 directhe direction of the Burgers vector of LN, i.e., in the ½2 tion. Following the reaction, LP is no longer exclusively in a prism plane, since the reaction creates out of plane jogs. The nature of LN is unchanged due to the interaction, except for a small perturbation of LN from a straight line due to the drag imposed by LP. Following the reaction, LN leaves the simulation domain from one side and re-enters the simulation domain on the opposite side of the domain, because of the periodic boundary condition in the ydirection.  1 0Þ prism plane 4.1.2. Case 2: vacancy loop in the ð0 1 The reaction mechanism for Case 1 and 2 are quite similar. Snap-shots of the MD simulations of the reaction are shown in Fig. 5. From these images it can be observed that when LN reacts with LP, see Fig. 5b and c, there is an observable spreading of the core of dislocation LN in the basal plane to enable the creation of 3

See the compute stress/atom comment in the LAMMPS manual.

129

the jog in LP. As LN nears LP, there is also an observable spreading of the dislocation core of LP. During the reaction, LN cuts LP, leaving behind jogs in LP in the basal plane with magnitudes of a=3 in  1 1 0 directhe direction of the Burgers vector of LN, i.e., in the ½2 tion. In contrast to Case 1, the jogs created by the reaction are in the original plane of LP and so LP completely resides in the prism  1 0Þ following the reaction. The nature of LN is plane ð0 1 unchanged due to the interaction, except for a small perturbation of LN from a straight line due to the drag imposed by LP.  0Þ prism plane 4.1.3. Case 3: vacancy loop in the ð1 0 1 The reaction mechanism observed in Case 3 differs significantly from that of Cases 1 and 2. Snap-shots of the actual MD simulations of the reaction are shown in Figs. 6 and 7. In the two previous cases, LP remained stationary, while LN glided due to the shear stresses created by the applied load. In Case 3, a rotation of LP in the counter-clockwise direction about the hci-axis is observed so that some segments of LP move towards LN (in the positive y-direction), while LN moves towards LP (in the negative y-direction). Due to the rotation of LP, LN only interacts with one segment of LP at a time, instead of two segments simultaneously, as would be expect if LP were sessile. At the time of the reaction with LN, LP  1 0 0Þ plane. It is has rotated so that it is approximately in the ð1 noteworthy that the Burgers vectors of LN and LP are the same in this case. LN reacts with LP to create a single helical dislocation,  1 1 0. This reaction occurs denoted as LH, with Burgers vector a=3½2 between stages (b) and (c) in Fig. 6. The segments of LH, which had previously been LN, continue to be driven in the y-direction by the applied stress. After stage (c), LH reacts with itself. The result of the reaction is the creation of a new ‘D’-shaped dislocation loop, denoted as LD, which is smaller than the original vacancy loop, LP, but has the same Burgers vector as LP, as shown in Fig. 6d. Effectively, LD is comprised of a combination of the middle segments of LN and the top half of LP. The second dislocation created by the self-reaction of LH, denoted by LN2, is effectively comprised of the outer segments of LN and the lower part of LP. Following the reaction, LN2 continues to glide due to the shear stress induced by the boundary conditions, while LD remains fairly stationary. Figures containing MD results were are drawn using Ovito [26]. 4.2. Hardening due to Irradiation Induced Loops In this subsection, the hardening effect of irradiation induced loops is investigated at 0.01 K. The simulation cell behaves in a nearly plane strain manner do to the use of periodic boundary conditions and the constraint that xx ¼ 0; rxy  0; rxz  0, and rxx – 0. The shear stress, ryz , versus the simple shear strain, yz , is plotted in Fig. 8 for Cases 1–3 for a strain rate of _ ¼ 1:2  108 s1 . The shear stress–strain curves are plotted along with the corresponding curves for a defect-free cell (which has a purely elastic response) and for a simulation cell containing just the mixed dislocation line in the basal plane, but with no irradiation induced loop. The stress–strain response of Cases 1–3 at _ ¼ 1:2  108 s1 clearly fall between the two reference curves, as expected, with Cases 1 and 2 showing a generally greater amount of hardening than Case 3. The dislocation loop in Case 1 poses a greater barrier to a dislocation gliding on the basal plane than the dislocation loop in Case 2. This is likely due to the fact that the dislocation reaction in Case 1 leads to an out of plane jog in the vacancy loop, while in Case 2 the dislocation reaction leads to an in-plane jog in the vacancy loop. We note that the stress–strain behavior does not deviate as substantially as might be expected from the elastic defect-free case due to the high loading rate. The elastic shear modulus in the yz-plane estimated from Fig. 8 is 44.4 GPa, which is very similar to the value of 44 GPa reported in [19].

130

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

y [1 0 1 0] t1

z [0 0 0 1]

b1

Vy

t2 b2 t2

a2 a1 a3

y [1 0 1 0]

x [1 2 1 0] Fig. 2. A schematic of the dislocation reaction observed in Case 1.

Z[0 0 0 1]

x [1 2

1 0]

(a) y [1 01

(b)

(d)

(c)

(f )

(e)

0]

Fig. 3. Case 1 – Interaction between the glissile mixed dislocation line (LN) and sessile dislocation loop (LP) under a strain rate of 6  107 s1 . The plane of the loop is parallel to the ha3 i-direction. (a) t ¼ 0 ps – yz ¼ 0%; (b) t ¼ 140 ps – yz ¼ 0:84%; (c) t ¼ 210 ps – yz ¼ 1:26%; (d) t ¼ 220 ps – yz ¼ 1:32%; (e) t ¼ 240 ps – yz ¼ 1:44%; and (f) t ¼ 250 ps – yz ¼ 1:50%.

y[1010]

x[1210]

(a)

(b)

(d)

(c)

(f)

(e)

Fig. 4. Case 1 – Interaction between the glissile dislocation line and irradiation induced dislocation loop under shear displacement (top view). The plane of the loop is parallel to ha2 i direction. (a) t ¼ 0 ps – yz ¼ 0%; (b) t ¼ 140 ps – yz ¼ 0:84%; (c) t ¼ 210 ps – yz ¼ 1:26%; (d) t ¼ 220 ps – yz ¼ 1:32%; (e) t ¼ 240 ps – yz ¼ 1:44%; and (f) t ¼ 250 ps – yz ¼ 1:50%.

Z[0 0 0 1]

x [1 2

1 0]

y [1 0 1 0]

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 5. Case 2 – Interaction between the glissile mixed dislocation line (LN) and sessile dislocation loop (LP) under a strain rate of 6  107 s1 . The plane of the loop is parallel to the ha1 i-direction. (a) t ¼ 0 ps – yz ¼ 0%; (b) t ¼ 120 ps – yz ¼ 0:72%; (c) t ¼ 160 ps – yz ¼ 0:96%; (d) t ¼ 190 ps – yz ¼ 1:14%; (e) t ¼ 230 ps – yz ¼ 1:38%; and (f) t ¼ 250 ps – yz ¼ 1:50%.

Fig. 9a–c illustrates the stress–strain responses of each of the three cases at each of the five loading rates. It is clear that strain rate impacts the shear stress–strain response significantly. The shear stress–strain curves for the lowest strain rate shows a

significantly softer response than those at the highest strain rates. In all three cases, it can be seen that for _ P 1:2  109 s1 the response of the simulation cell is very similar to the elastic response. This occurs because the simulation cell is being

131

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

Z[0 0 0 1]

x [1 2

1 0]

y [1

010

]

(a)

(b)

(c)

(d)

(e)

Fig. 6. Case 3 – Interaction between the glissile mixed dislocation line (LN) and irradiation induced dislocation loop (LP) under a strain rate of 6  107 1/s. The plane of the loop is parallel to the ha2 i-direction. (a) t ¼ 138 ps – yz ¼ 0:83%; (b) t ¼ 144 ps – yz ¼ 0:86%; (c) t ¼ 148 ps – yz ¼ 0:89%; (d) t ¼ 152 ps – yz ¼ 0:91%; and (e) t ¼ 154 ps – yz ¼ 0:93%.

y[1010]

x[1210]

(a)

(b)

(c)

(d)

(e)

Fig. 7. Case 3 – Interaction between the glissile dislocation line and irradiation induced dislocation loop under shear displacement (top view). The plane of the loop is parallel to ha1 i direction. (a) t ¼ 138 ps – yz ¼ 0:83%; (b) t ¼ 144 ps – yz ¼ 0:86%; (c) t ¼ 148 ps – yz ¼ 0:89%; (d) t ¼ 152 ps – yz ¼ 0:91%; and (e) t ¼ 154 ps – yz ¼ 0:93%.

Fig. 10 plots the shear stress at which the dislocation lines start to glide in Case 1 for each loading rate simulated. The simulation results are fitted to a polynomial strain rate hardening law:

900 800

Elastic Case 1 Case 2 Case 3 Line

700

σyz (MPa)

600

rY ¼ rY0 þ a_ m where rY0 ¼ 70 MPa, a ¼ 0:4321, and m ¼ 0:3284. The quasi-static flow stress, rY0 ¼ 70 MPa estimated for the mixed dislocations simulated here is very close to the Peierls stress of 79 MPa for screw dislocations in the basal plane reported in [20]. While the simulated stress–strain responses presented are not representative of quasi-static behavior, key characteristics of these curves can be attributed to characteristics of the dislocation reaction and the nature of the barrier posed by irradiation loop. Fig. 11 illustrates the stress–strain curves for Case 1 for a strain

500 400 300 200 100 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

εyz (%) Fig. 8. Shear stress versus simple shear strain for Cases 1–3 plotted with the elastic response (Elastic) and the response of a single gliding dislocation (line) for a strain rate of 1:2  108 s1 .

deformed at such a high rate that the velocity of the gliding dislocation is insufficient to completely accommodate the imposed deformation. In contrast, there is a marked decrease in the shear stress as the applied velocity (or strain rate) is decreased. It is noteworthy, that the stress–strain response for the lowest strain rate are still too large for the results to be representative of quasi-static behavior. Simulation at lower strain rates require substantially more computational resources than are currently available to the authors.4

4 the results for a strain rate of 1:2  107 s1 for each of the 3 cases each took about 7 days to obtain using a 32 cpu cluster.

rates of 6  107 s1 and 1:2  107 s1 , where the points labeled b–f are associated with the stages of the dislocation reaction illustrated in Fig. 3b–f, respectively. Point (b0) denotes the stress where the dislocation line, LN, begins to glide. Point (b) corresponds to the beginning of the interaction on the dislocation line with the first segment of the loop, LP. Initially NL is pinned by LP, resulting in the increase in the applied stress from point b to point (c). At point (c), the stress–strain curve reaches a local maximum as LN overcomes the barrier posed by the first segment of the irradiation loop. The stress–strain response initially softens as LN glides between the two segments of LP, point (d). The stress then rises once again as LN interacts with the second segment of LP, point (e). Finally, the stress drops, moving from (e) to (f), as LN glides and overcomes the barrier posed by the second segment of LP. The stress–strain behavior of Case 2 closely follows that of Case 1 and so is not discussed in detail; Case 3, however, varies considerably. Fig. 12 illustrates the stress–strain curves for Case 3 at a strain rate of 6  107 s1 and 1:2  107 , where the points labeled (a)–(d) are associated with the stages of the dislocation reaction illustrated in Fig. 6a–d. The nonlinear segment of the stress–strain response before point (a) corresponds to the motion of both LN and LP.

132

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134 800

800 dε/dt=12E8 dε/dt=6E8 dε/dt=1.2E8 dε/dt=0.6E8 dε/dt=0.12E8

700

600

500

σyz (MPa)

σyz (MPa)

600

400 300

500 400 300

200

200

100

100

0 0

0.2

0.4

dε/dt=12E8 dε/dt=6E8 dε/dt=1.2E8 dε/dt=0.6E8 dε/dt=0.12E8

700

0.6

0.8

1

1.2

1.4

1.6

0 0

1.8

0.2

0.4

0.6

εyz (%)

0.8

1

1.2

1.4

1.6

1.8

εyz (%)

(a) Case 1

(b) Case 2

800 dε/dt=12E8 dε/dt=6E8 dε/dt=1.2E8 dε/dt=0.6E8 dε/dt=0.12E8

700 600

σyz (MPa)

500 400 300 200 100 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

εyz (%)

(c) Case 3 Fig. 9. Shear stress versus shear displacement curves for Cases 1–3 for different strain rates.

slightly following the reaction of LN and LP, creating a single helical loop LH. At point (d), the stress–strain curve reaches a local minimum following the self reaction of LH, creating a ‘‘D’’-shaped loop, LD, and another dislocation line LN2. Initially, as LN2 glides away from LD, there is little change in the stress–strain response, moving from (d) to (e).

500 450 Yield Stress Fit Curve

σyz (GPa)

400 350 300

5. Discussion

250

In MD simulations, the small size of the simulation box and the short simulated time lead to extremely high strain rates, which generally complicates the comparison of MD simulations directly with experiments. The effective strain rates for the MD simulations

200 150 100

1

2

3

4

5

6

7

8

Strain Rate (1/sec)

9

10

11

12 8

x 10

Fig. 10. The stress at which the line dislocation starts to glide as a function of strain rates.

Point (a) denotes the stress where LN begins to react with LP. Point (b) corresponds to the end of the interaction of the dislocation line with the first segment of LP. The stress–strain response softens

reported here vary from 1:2  107 s1 to 1:2  109 s1 . The reaction mechanisms are almost unaffected by strain rate, with the most noticeable difference being that bending of the gliding dislocation line LN during the reaction with loop LP is more pronounced at lower strain rates. While simulations at lower strain rates are desirable to better characterize the strength of the barrier that vacancy loops pose to gliding dislocations, they are beyond the computational resources currently available to the authors. Future studies should revisit these dislocation reactions at lower strain rates. The stress–strain curves reported in Figs. 8, 9, 11, and 12 are continually increasing, which is counter intuitive. This

133

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134 450

600

e

c

f

d

400

e

500 350 400

b Dislocation line starts to glide

b0

yz

250

σ (MPa)

σyz (MPa)

300

200 150

f c d

300

200

Dislocation line starts to glide b0

b

100 100 50 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

0 0

1.6

0.2

0.4

0.6

0.8

ε (%)

1

1.2

1.4

1.6

1.8

εyz (%)

yz

(a) ˙ = 6 × 108 s−1

(b) ˙ = 1.2 × 107 s−1

Fig. 11. Case 1 stress–strain curve at two strain rates. Points labeled b–f are associated with the images in Fig. 3b–f, respectively.

450

400

400

350

350

300

a

250

σyz (MPa)

σ (MPa) yz

300

e

bd

200 150

200

a

b e

150

d

100

100

50

50 0 0

250

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

ε (%) yz

(a) ˙ = 6 ×

108 s−1

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

ε (%) yz

(b) ˙ = 1.2 × 107 s−1

Fig. 12. Case 3 stress–strain curve at two strain rates. Points a to e are associated with images in Fig. 6a–e.

phenomenon occurs in part because the velocity of dislocation glide is too low for the dislocation motion alone to accommodate the deformation imposed by the boundary conditions. The imposed deformation, which is not accommodated by the motion of the gliding dislocation, must be accommodated by elastic deformation of the simulation cell. The presence of additional elastic deformations, which would not be present at quasi-static loading rates, is one reason why the stress–stress curves tend to continuously increase. For strain rates greater than 6  108 s1 , the stress–strain curves demonstrate a delayed softening behavior, which increases with strain rate. For Case 1 and strain rate 6  108 s1 , softening in the stress–strain response occurs at about 7.7% strain. For the lowest strain rate of 6  108 s1 softening following the dislocation reaction is observed for Cases 1 and 2 at about 1.5% and 1.4%, respectively. In Case 3 softening occurs at about 0.6% strain. The unexpected rise in the stress–strain curve for Case 3 following the initial reactions is attributed to LN2 interacting with LD for a second time, due to wrap-around effects caused by periodic

boundary conditions. This suggests that the results presented for Case 3 are influenced by the domain size. While the simulation of larger domains is desirable, it is beyond the current computational resources available to the authors. The increase in the applied shear stress observed in Cases 1–3 compared to the case of pure glide of the mixed dislocation without the loop are generally higher compared to those observed in the MD studies of Serra and Bacon [21]. This result is not unexpected. Here reactions with dislocation gliding on the basal plane are consider; whereas in this previous work, dislocation gliding on the prism planes were model. Secondly, the defects studies in this earlier article are significantly smaller than those studied here and so differences are to be expected. Lastly, Serra and Bacon [21] conducted their simulations at 300 K, while those presented here are at 0 K and so it would be expected that the critically resolved shear stress required to overcome obstacles in this previous work would be lower than those in the current work. Temperature effects were not investigated in the current study, but should be in the future.

134

K. Ghavam, R. Gracie / Journal of Nuclear Materials 462 (2015) 126–134

6. Conclusions Dislocation reactions between mixed dislocations gliding in  1 0 0g, characteristic of ð0 0 0 1Þ and elliptical vacancy loops in f1 irradiation induced defects in Zirconium, were investigated with Molecular Dynamics. Elliptical vacancy loops with their major axes parallel to ½0 0 0 1, lengths of 20 nm, and minor to major ratios of  1 0 0Þ or ð0 1  1 0Þ 0.8 were analyzed. When a vacancy loop in ð1   reacts with a mixed dislocation parallel to ½1 2 1 0 with a Burgers  1 1 0, the result is the creation of two jogs in the vector of a=3½2  1 1 0 direction. In these two cases, the vacancy loop in the ½2 vacancy loop is a relatively strong barrier to the motion of the  21  0 mixed dislocation. When a mixed dislocation parallel to ½1  1 1 0 reacts with a vacancy loop in with Burgers vector a=3½2  0Þ with Burgers vector a=3½2  1 1 0, part of the vacancy loop ð1 0 1 is absorbed into the mixed dislocation and part of the vacancy loop is left behind in the form of a smaller ‘‘D’’-shaped dislocation loop. In this third case, the vacancy loop is a relatively weak barrier to the motion of the mixed dislocation. In all cases, strain rate was demonstrated to significantly effect the stress–strain behaviour; however, the reaction mechanism were unaffected.

Acknowledgments Financial support for this research from Atomic Energy of Canada Limited (AECL) and a National Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant award to Dr. Gracie are acknowledged. Guidance provided by AECL scientists Malcolm Griffiths and Wenjing Li is greatfully acknowledged. The

authors’ would also like to thank Compute Canada and SHARCNET for access to high performance computing resources. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

A. Akhtar, A. Teghtsoonian, Acta Metall. 19 (7) (1971) 655–663. A. Akhtar, Acta Metall. 21 (1) (1973) 1–11. A. Jostsons, P.M. Kelly, R.G. Blake, J. Nucl. Mater. 66 (3) (1977) 236–256. M. Griffiths, J. Nucl. Mater. 159 (1988) 190–218. P. Kelly, R. Blake, Philos. Magaz. 28 (2) (1973) 415–426. F. Onimus, I. Monnet, J. Béchade, C. Prioul, P. Pilvin, J. Nucl. Mater. 328 (2) (2004) 165–179. H. Rosenbaum, G. Rieger, D. Lee, Metall. Trans. 5 (8) (1974) 1867–1874. C. Woo, X. Liu, J. Nucl. Mater. 393 (3) (2009) 513–517. D. Bacon, V. Vitek, Metall. Mater. Trans. A 33 (13) (2002) 721–733. S. Wooding, L. Howe, F. Gao, A. Calder, D. Bacon, J. Nucl. Mater. 254 (2) (1998) 191–204. M.S. Khan, M.S. Islam, D.R. Bates, J. Mater. Chem. 8 (10) (1998) 2299–2307. D. Bacon, F. Gao, Y.N. Osetsky, J. Nucl. Mater. 276 (1) (2000) 1–12. J.-P. Crocombette, D. Ghaleb, J. Nucl. Mater. 295 (2) (2001) 167–178. V. Kapinos, Y.N. Osetskii, P. Platonov, J. Nucl. Mater. 184 (2) (1991) 127–143. F. Gao, D. Bacon, L. Howe, C. So, J. Nucl. Mater. 294 (3) (2001) 288–298. R. Voskoboynikov, Y.N. Osetsky, D.J. Bacon, Mater. Sci. Eng.: A 400 (2005) 54– 58. R. Voskoboynikov, Y.N. Osetsky, D.J. Bacon, Mater. Sci. Eng.: A 400–401 (2005) 49–53. G. Ackland, S. Wooding, D. Bacon, Philos. Magaz. A 71 (3) (1995) 553–565. M.I. Mendelev, G.J. Ackland, Philos. Magaz. Lett. 87 (5) (2007) 349–359. H. Khater, D. Bacon, Acta Mater. 58 (8) (2010) 2978–2987. A. Serra, D. Bacon, Modell. Simul. Mater. Sci. Eng. 21 (4) (2013) 045007. P. Fossati, L. Van Brutzel, B. Devincre, J. Nucl. Mater. 443 (1) (2013) 359–365. Y.N. Osetsky, D.J. Bacon, Modell. Simul. Mater. Sci. Eng. 11 (2003) 427–446. D. Terentyev, P. Grammatikopoulos, D.J. Bacon, Y.N. Osetsky, Acta Mater. 56 (2008) 5034–5046. D.J. Bacon, Y.N. Osetsky, Math. Mech. Solids 14 (2009) 270–283. A. Stukowski, Visualization and analysis of atomistic simulation data with ovito–the open visualization tool, Modell. Simul. Mater. Sci. Eng. 18 (1) (2010) 015012.