- Email: [email protected]

Contents lists available at SciVerse ScienceDirect

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

Charge-optimized many-body (COMB) potential for zirconium Mark J. Noordhoek, Tao Liang (梁涛), Zizhe Lu (芦子哲), Tzu-Ray Shan (單子睿) 1, Susan B. Sinnott, Simon R. Phillpot ⇑ Department of Materials Science and Engineering, University of Florida, Gainesville, FL 32611-6400, USA

a r t i c l e

i n f o

Article history: Received 8 March 2013 Accepted 1 June 2013 Available online 11 June 2013

a b s t r a c t An interatomic potential for zirconium is developed within the charge-optimized many-body (COMB) formalism. The potential correctly predicts the hexagonal close-packed (HCP) structure as the ground state with cohesive energy, lattice parameters, and elastic constants matching experiment well. The most stable interstitial position is the basal octahedral followed by basal split, in agreement with recent first principles calculations. Stacking fault energies within the prism and basal planes satisfactorily match first principles calculations. A tensile test using nanocrystalline zirconium exhibits both prismatic 0gh1 1 2 0i slip and pyramidal f1 1 2 2gh1 1 2 3i slip, showing the model is capable of reproducing f1 0 1 the mechanical deformation modes observed in experiments. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The hexagonal close-packed (HCP) structure presents unique challenges in describing interatomic interactions compared to other face-centered cubic (FCC) or body-centered cubic (BCC) materials. For example, the energetics of stacking faults in HCP metals is largely determined by the c0/a0 ratio. HCP metals with p a c0/a0 ratio close to or greater than the ideal ( (8/3) 1.633) have lower stacking fault energies on the basal {0 0 0 1} planes than 0g planes and are thus more susceptible to those on prism f1 0 1 basal slip. By contrast, HCP metals with a c0/a0 ratio much less than the ideal value have lower stacking-fault energies on the prism plane than on the basal plane, and thus undergo slip primarily on the prism planes. Zirconium is one such HCP metal and is the focus of this work. One important application of zirconium is in alloys commonly used for fuel claddings in nuclear reactors. These alloys were selected due to the low thermal neutron cross-section of Zr, their good corrosion resistance, and good mechanical properties [1]. Zircaloy 4 (typically Zr – 1.5% Sn – 0.2% Fe – 0.1% Cr) is the current industry standard, while ZIRLOTM (Zr – 1% Sn – 1% Nb) is an emerging Westinghouse-developed material. During the course of the fuel cycle, the fuel cladding undergoes a number of mechanical and chemical processes. Current challenges that limit cladding performance include grid-to-rod fretting (GTRF), pellet-cladding interactions (PCI), and corrosion [2]. Molecular dynamics (MD) simulations of zirconium and other metals commonly use the embedded atom method (EAM) or its modified derivatives. While ⇑ Corresponding author. 1

E-mail address: [email protected] (S.R. Phillpot). Present address: Sandia National Laboratories, Albuquerque, NM 87185, USA.

0022-3115/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnucmat.2013.06.004

these potentials provide a satisfactory characterization of many important properties in zirconium, their formalisms are incapable of describing nonmetallic environments and are thus not able to attack the issues of PCI or corrosion. Addressing such problems requires a simulation framework capable of handling each individual material system (metal, metal-oxide, metal-hydride), in addition to systems in which they co-exist and interact. The charge-optimized many-body (COMB) formalism is capable of providing a unified description of the interatomic interactions in systems with mixed bonding [3–5]. As a first step to developing a description of the Zr–ZrO2–ZrH2 system, here we develop a COMB description of Zr. For this work, charges for Zr are set to zero for the validation of the metal subsystem. Charged interactions will be included upon extension to the oxide and hydride systems. We show how the COMB potential for Zr is capable of providing a satisfactory description of the mechanical properties. Additionally, point defect phenomena are described and compared favorably with density functional theory (DFT) results and previously described interatomic potentials. 2. Potential form and fitting scheme In this section, we provide a summary of the pertinent chargeneutral atomic interactions within the third-generation COMB formalism; a detailed description including charged interactions may be found elsewhere [3]. The short-range interaction, UShort(rij), includes a summation over three attraction terms. The subscripts ij indicate interactions between atoms i and j.

U Short ðr ij Þ ¼

1 XX fF c ðr ij Þ½V R ðrij Þ bij V A ðr ij Þg 2 i j–i

ð1Þ

M.J. Noordhoek et al. / Journal of Nuclear Materials 441 (2013) 274–279

where:

3. Results

V R ðrij Þ ¼ Aij expðkij r ij Þ

ð2Þ

and

V A ðrij Þ ¼

3 X Bnij expðanij rij Þ

ð3Þ

n¼1

The bond order term, bij, is defined as:

"

!gi #21g N i X bij ¼ 1 þ fijk ðr ij ; rik ÞGðcos hijk Þ

ð4Þ

k–i;j

where: m

fijk ðr ij ; r ik Þ ¼ expðbij i ½r ij r ik mi Þ

ð5Þ

Gðcos hijk Þ ¼ b0 þ b1 cos hijk þ b2 cos2 hijk þ b3 cos3 hijk þ b4 cos4 hijk þ b5 cos5 hijk þ b6 cos6 hijk

ð6Þ

and

8 > >

9 1 rij 6 R > > = o r ij R 1 1 F c ðr ij Þ ¼ þ 2 cos p SR R < rij < S 2 > > > > : ; 0 r ij P S

ð7Þ

The cut-off function in Eq. (7) ensures a smooth, continuous function at the cut-off. The parameters developed in this work are shown in Table 1. This parameter set was determined by minimizing a penalty function, f(p), for each parameter set, p, with a least-squares method. The penalty function is calculated from weighted square absolute deviations of lattice parameters, elastic constants, and structural properties from target values obtained from experiments and electronic structure calculations. For this parameterization, nine structures consisting of five bulk structures, three stacking fault structures, and one surface structure were included in the fitting database.

Table 1 COMB potential parameters for Zr. COMB Aij (eV) B1ij (eV)

1122.16000 711.89600

B2ij (eV)

171.35700

B3ij (eV)

170.50400

k (Å1) a1ij (Å1)

1.67129 1.85834

a2ij (Å1) a3ij (Å1)

1.29329

b (Å1)

2.27174 0.444536 1.0 0.025674 0.025066 0.010837 0.002500 0.009815 0.002271 0.001545 4.1 4.2

gi mi b6 b5 b4 b3 b2 b1 b0 R (Å) S (Å)

275

1.22281

3.1. Bulk properties Table 2 compares the properties of zirconium determined from experiment, electronic-structure calculations, COMB, and various empirical potentials. The experimental 4 K a-axis lattice parameter is 3.23 Å [6]; COMB yields a 0 K value of 3.226 Å. The experimental c-axis lattice parameter is 5.15 Å; COMB yields 5.159 Å. For COMB, the bulk and shear moduli, as well as each individual elastic constant match experiment within 11%. Density functional theory (DFT) yields values of 227 J/m2 and 197 J/m2 for the energies of 0g stacking the basal (0 0 0 1) I2 stacking fault and the prism f1 0 1 fault, respectively [7]. COMB matches these values within 18%, with the ratio of the energy of the prism fault to the basal I2 fault of 0.81 compared to 0.87 for DFT. DFT results show that the surface energy is higher for the prism plane than for the basal plane. The COMB potential reverses this phase order with surface energies matching DFT values within 11%. COMB matches or improves upon the properties listed in Table 2 as compared to other empirical potentials. For all potentials, the lattice parameters and elastic constants satisfactorily match experimental results. Major differences arise when comparing stacking fault energies and surface energies. The analytic modified embedded atom method (AMEAM) [8] and long-range empirical potential (LREP) [9] severely underestimate these values, while COMB and EAM give values more in line with DFT. For stacking fault energies, COMB and EAM [10] closely match DFT, but COMB overestimates the values while EAM underestimates them.

3.2. Point defects It is known that irradiation-induced growth and creep plays an important role in material behavior when Zr-alloys are used as nuclear fuel cladding applications. Correctly describing the interstitial structure and formation energy is the first step for describing the microstructural evolution during irradiation. There are eight interstitial configurations commonly referred to in the literature, as shown in Fig. 1. The octahedral (O) and tetrahedral (T) sites lie be 0i tween basal planes and differ in their placement along the h1 1 2 direction. Their basal counterparts (BO and BT) lie directly below (above) within the (0 0 0 1) plane. The crowdion (C) interstitial is located between atoms in adjacent basal planes, while the basal crowdion (BC) lies between neighboring atoms in the h1 0 0 0i direction within the basal plane. Split (S) and basal split (BS) interstitials correspond to two atoms sharing one lattice position, with atoms being displaced along the h0 0 0 1i and h1 0 0 0i directions, respectively. In total, four interstitial sites lie within the basal plane while four are located between basal planes. Results of first-principles calculations showed a strong dependence on the conditions of the details of the computational approach [12–14]. In particular, early studies were limited to very small supercells, consisting of 4 4 3 repeat units containing 97 atoms [12–13]. In one study [12], the calculations were performed with constant volume and shape of the cell, allowing the internal relaxation of atoms. In others [13–14], the calculations were performed at constant volume per atom, while maintaining the shape of the supercell. The defect energies converge slowly with increasing system size [14]. In particular, it was found that the supercell should be at least 8 unit-cells along h1 0 0 0i for accurate results. Increasing the supercell size along h0 0 0 1i exhibited much smaller relative differences in the defect formation energy. Since no point defect information was included in the fitting database, all of the COMB values presented here are predictions of the potential. The COMB calculations for interstitial formation

276

M.J. Noordhoek et al. / Journal of Nuclear Materials 441 (2013) 274–279

Table 2 Properties of zirconium metal given by the COMB potentials compared with experiments, DFT calculations, and other empirical potentials. A * indicates that the select property was included in the fitting database. Properties

Experiment [6,11]

DFT [7]

EAM [10]

LREP [9]

AMEAM [8]

COMB

COMB/experiment [DFT] (% error)

a0 (Å)* c0 (Å)* c0/a0 Ec (eV/atom)* C11 (GPa)* C12 (GPa)* C13 (GPa)* C33 (GPa)* C44 (GPa)* B G DE (fcc–hcp) (eV/atom) DE (bcc-hcp) (eV/atom)* Basal I2 stacking fault (mJ/m2)* Prism stacking fault (mJ/m2)* Basal surface (mJ/m2) Prism surface (mJ/m2)*

3.23 5.15 1.5944 6.32 155 67 65 172 36 97 42

3.23 5.17 1.6006 N/A 156 61 62 166 26 94 39 0.040 0.078 227 197 1600 1660

3.234 5.168 1.598 6.635 147 69 74 168 44 100 42 0.054 0.103 198 175 1770 1930

3.23 5.15 1.5944 6.25 121.4 59.5 78.4 219.3 49.9 99.4 42.5 0.030 0.063 84 N/A 1097 1156

3.231 5.148 1.593 6.25 144 74 67 166 33.4 96.7 36.8 0.0049 0.0179 26 N/A 988 978

3.226 5.159 1.599 6.32 154 71 62 164 32 96 39 0.074 0.105 267 215 1770 1740

0 0 0 0 1 6 5 5 11 1 7 85 35 18 9 11 5

Octahedral (O)

Basal Octahedral (BO)

Tetrahedral (T)

Crowdion (C)

Basal Tetrahedral (BT)

Basal Crowdion (BC)

Split (S)

Basal Split (BS)

Fig. 1. Interstitial positions in an HCP lattice. The black atoms represent the interstitial, while gray atoms are normal HCP atoms. Sites in the top row lie out of the basal plane, while sites in the bottom row lie within the basal plane.

energies were performed using constant supercell volume while allowing for internal relaxation. Two supercell sizes of 10 3 4 and 30 20 20 repeat unit cells were used. Overall, the interstitial formation energies in COMB are slightly greater than those predicted by DFT, with the most stable sites matching within 0.8 eV. DFT calculations using large supercells showed the basal octahedral (BO) position to be the most stable, followed closely by the basal split (BS) [14]. COMB satisfactorily matches this trend, albeit with the energy differences between the two being smaller. The full results are shown in Table 3. Values listed as ‘‘Unstable” in Table 3 indicate that the interstitial spontaneously moves to another interstitial site or forms a complex structure. DFT calculations [14] suggest that interstitial sites within the basal plane have lower formation energies than their counterparts out of the basal plane. COMB matches this trend, with the exception being for the tetrahedral site, which has no DFT data for comparison. We have performed first-principles calculations for vacancy and divacancy formation energies using VASP [15]. The conditions for

the calculation are the same as those used in Ref. [12]. Vanderbilt ultrasoft (US) pseudopotentials within the generalized gradient approximation (GGA) were chosen [12,16–18]. The first Brillouin zone was sampled with a 3 3 4 grid; the energy cut-off was 225 eV. Here we used a larger supercell consisting of 5 6 3 repeat unit cells containing 180 atoms. The calculated lattice parameters for this configuration are 3.2314 Å in the a-direction and 5.1816 Å in the c-direction. Point defect formation energies were calculated using constant cell shape and volume, while internal relaxation of the atoms was allowed. This procedure resulted in a vacancy formation energy of 1.94 eV, which is similar to 1.86 eV obtained from a previous DFT study [12]. A schematic showing the only two possible structures of adjacent vacancies in an HCP lattice is given in Fig. 2. A distinction is made between vacancies next to each other within the basal (0 0 0 1) plane (in-plane) and those in adjacent basal planes (outof-plane). The DFT results show that divacancies within the same basal plane have positive binding energies, while those in adjacent

277

M.J. Noordhoek et al. / Journal of Nuclear Materials 441 (2013) 274–279 Table 3 Point defect formation energies. Units are in eV. Dashed lines indicate that no value was given.

Supercell dimensions Octahedral (O) Basal octahedral (BO) Crowdion (C) Basal crowdion (BC) Split (S) Basal split (BS) Tetrahedral (T) Basal tetrahedral (BT)

DFT (GGA) [12]

DFT (LDA) [13]

DFT (GGA) [14]

COMB

COMB

443 2.84 2.88 3.08 2.95 3.01 – Unstable 4.03

443 2.79 2.78 3.07 Unstable 2.80 2.90 – –

10 4 3 2.93 2.75 – – 3.02 2.87 – –

10 4 3 4.18 3.63 Unstable Unstable Unstable 3.66 3.89 4.31

30 20 20 4.12 3.53 Unstable Unstable 4.57 3.60 3.78 4.14

3.3. Stacking fault energy

Fig. 2. Schematic of (left) an in-plane divacancy and (right) an out-of-plane divacancy. The vacancies are denoted by open circles.

Table 4 Divacancy formation and binding energies. Units are in eV. DFT (this work) In-plane EFormation in2v Ebinding in2v Out-of-plane EFormation out2v Ebinding out2v

LREP [9]

AMEAM [8]

COMB

3.78 0.07

3.30 0.20

3.18 0.22

4.91 0.12

3.95 0.10

3.31 0.19

3.18 0.21

4.91 0.12

basal planes have negative binding energies. A negative binding energies means it is more energetically favorable for the vacancies to be well separated. For COMB and other empirical potentials, the divacancy binding energies are all positive and no significant differences for in-plane or out-of-plane divacancies are observed. The values predicted by COMB are approximately 1 eV larger than DFT, while LREP and AMEAM underestimate the DFT results by approximately 0.6 eV. The results are summarized in Table 4.

The stacking fault energy map for the (0 0 0 1) basal plane from COMB, DFT, and EAM [19] is shown in Fig. 3. Only the basal I2 stacking fault is included in the fitting database; the rest are predictions. The energies of the various translational configurations of the stacking faults were calculated by only allowing atomic relaxation in the direction perpendicular to the stacking fault plane 2 0i ¼ 1 h1 0 1 0i þ 1 h0 1 1 0i leading to [19]. The indirect path is 13 h1 1 3 3 saddle points in the center and corners of the maps in Fig. 3. The energy at these saddle points for COMB is 267 mJ/m2, which compares quite well with 213 mJ/m2 from DFT and 199 mJ/m2 from 0i partial EAM [19]. The unstable stacking fault energy for h1 1 2 dislocation is 346 mJ/m2, compared to about 260 mJ/m2 for DFT and 314 mJ/m2 from EAM [19]. 0g prism plane from The stacking fault energy map for the f1 0 1 COMB, DFT, and EAM [19] is shown in Fig. 4. Here, the direct path is 0i dislocation. The stable and unstable stacking the prismatic h1 1 2 faults are included in the fitting database. The energy for the stable prismatic stacking fault is 215 mJ/m2 from COMB, compared to 166 mJ/m2 from DFT and 145 mJ/m2 from EAM [19]. The unstable 0i dislocation is stacking fault energy for the prismatic h1 1 2 247 mJ/m2, compared to 194 mJ/m2 for DFT [19]. Thus, the difference between the stable and unstable stacking fault for prismatic 0i dislocation from COMB (32 mJ/m2) is very similar to the h1 1 2 DFT value (28 mJ/m2). In addition, the calculated unstable stacking 0i dislocation (247 mJ/m2) is only fault for the prismatic h1 1 2 slightly lower than the stable stacking fault energy in the basal plane (267 mJ/m2), which ensures that prismatic slip will be the predominantly observed slip system during deformation. The Legrand [20] criteria for the type of slip observed is based 66 cbasal on R ¼ CC44 cprism , where R < 1 leads to basal slip, R > 1 leads to prism slip, and R 1 leads to both types being activated. Our potential has a value of R = 1.61, compared with 1.44 for DFT [19], 1.20 for

Fig. 3. Basal (0 0 0 1) stacking fault maps for (a) DFT [19], (b) EAM [19] and (c) COMB.

278

M.J. Noordhoek et al. / Journal of Nuclear Materials 441 (2013) 274–279

0g stacking fault maps for (a) DFT [19], (b) EAM [19] and (c) COMB. Fig. 4. Prism f1 0 1

Fig. 5. (a) Common neighbor analysis (CNA) map [24] of 3D nano-crystalline Zr after tensile test with 8.8% strain. Red atoms are disordered atoms, blue atoms are normal HCP atoms. (b) Normal HCP atoms are removed from the strained structure. The orange atoms are atoms with coordination different from that of perfect HCP. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

an EAM potential [10], and 2.3 from Legrand [20]. Since the stacking fault energetics largely determine the mechanical response of a metal, the above analysis gives us confidence that the COMB potential for Zr will correctly reproduce at least the qualitative features of the mechanical response of Zr. We now test this explicitly. 3.4. Mechanical deformation In order to test the appropriateness of this potential for the study of the mechanisms of mechanical deformation, in this section we show preliminary results for a tensile test of a 3-D zirconium nanocrystal. The simulations were run in LAMMPS [21] using the COMB potential presented above. The nanocrystalline structure as shown in Fig. 5 contains 16 randomly oriented grains with average grain size of 17 nm, consisting of total of approximately 1.8 million atoms. This structure was annealed at various temperatures up to 900 K and then cooled to room temperature. This annealing process allows atoms within the grain boundary to diffuse to equilibrium positions. A tensile test was conducted at 300 K with a strain rate of 4.0 109 s1 using a timestep of 0.1 fs. An earlier study of Mg used a strain rate of 1.5 109 s1 [22], while an Al study used a strain rate of 2.0 108 s1 [23]. The high strain rates used in MD as compared to experiments

(1.0 103 s1) are due to the very short timescales accessible by MD. However, these high strain rates have been shown to describe the relevant physics in nanocrystalline metals [23]. In order to distinguish FCC like atoms from normal HCP atoms, common neighbor analysis (CNA) was used [24]. The structure was visualized with AtomEye [25]. Table 5 shows the twin and slip systems manifested experimen 0gh1 1 2 0i slip and the pyramitally in zirconium. Prismatic f1 0 1 dal f1 1 2 2gh1 1 2 3i slip are frequently observed in experiments [26], with prismatic slip being most easily activated. During the simulated tensile test, dislocation activity was observed in 15 of the 16 grains as shown in Fig. 5. Both prismatic and pyramidal dis 0gh1 1 2 0i is a full locations were also observed. The prismatic f1 0 1

Table 5 Slip and twinning modes of zirconium [28,29]. Properties

Zirconium

Slip systems

0gh1 1 2 0i Principle f1 0 1 2 0if1 1 2 2gh1 1 2 3i Secondary ð0 0 0 1Þh1 1

Twin systems

2gh1 0 1 1i 1st order f1 0 1 1gh1 1 2 6i 2nd order f1 1 2

M.J. Noordhoek et al. / Journal of Nuclear Materials 441 (2013) 274–279

2gh1 1 2 3i dislocation is a dislocation, while the pyramidal f1 1 2 partial. These simulation results are consistent with the above experimental observations. Twinning behaviors are also observed during the tensile test. It thus appears that the COMB potential can successfully describe the dislocation and twinning behaviors in zirconium. A full analysis of these deformation simulations will appear elsewhere [27]. 4. Conclusions

References [1] [2] [3] [4] [5] [6] [7] [8]

We have presented an empirical potential within the COMB formalism for zirconium that is capable of modeling a wide range of physical phenomena. Our model predicts that the basal octahedral interstitial is the most energetically favorable, in accordance with the most recent DFT study. Stacking fault energies and the relative values between the different types satisfactorily match DFT calculations reasonably well. With this ability to correctly model the stacking fault energies, we confirm that the potential is able to reproduce the various slip and twinning modes during mechanical deformation. Most importantly the COMB framework will allows us to extend this potential, so as to simulate materials structures of zirconium metal, with an oxide surface and hydride inclusions.

[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

Acknowledgements [23]

We are very happy to acknowledge very useful conversations with Donald Brenner and Stephen Foiles. We are also grateful to German Samolyuk, Stas Golubov, Yuri Ozetsky and Roger Stoller of ORNL for sharing the results of their DFT studies of defect energetics prior to publication. This research was supported by the Consortium for Advanced Simulation of Light Water Reactors (www.casl.gov), an Energy Innovation Hub (http://www.energy. gov/hubs) for Modeling and Simulation of Nuclear Reactors under US Department of Energy Contract No. DE-AC05-00OR22725.

279

[24] [25] [26] [27] [28] [29]

B. Cox, J. Nucl. Mater. 336 (2005) 331–368. K. Edsinger, C.R. Stanek, B.D. Wirth, JOM 63 (2011) 49–52. T. Liang, B. Devine, S.R. Phillpot, S.B. Sinnott, J. Phys. Chem. A 116 (2012) 7976. J. Yu, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 75 (2007) 085311. T.-R. Shan, B.D. Devine, T.W. Kemper, S.B. Sinnott, S.R. Phillpot, Phys. Rev. B 81 (2010) 125328. G. Simmons, H. Wang, Single Crystal Elastic Constants: A Handbook, MIT Press, Cambridge, MA, 1971. Y. Udagawa, Y. Masatake, H. Abe, N. Sekimura, T. Fuketa, Acta Mater. 58 (2010) 3927–3938. W. Hu, B. Zhang, F. Gao, D. Bacon, J. Phys.: Condens. Matter 13 (2001) 1193– 1213. Y. Dai, J.H. Li, B.X. Liu, J. Phys.: Condens. Matter 21 (2009) 385402. M.I. Mendelev, G.J. Ackland, Philos. Mag. Lett. 87 (2007) 349–359. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1986. C. Domain, A. Legris, Philos. Mag. 85 (2005) 569–575. F. Willaime, J. Nucl. Mater. 323 (2003) 205–212. G.D. Samolyuk, S.I. Golubov, Y.N. Osetsky, R.E. Stoller, Philos. Mag. Lett. (2012) 1–8.