An interatomic potential for simulation of Zr-Nb system

Computational Materials Science 129 (2017) 259–272

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

An interatomic potential for simulation of Zr-Nb system D.E. Smirnova a,⇑, S.V. Starikov a,b a b

Joint Institute for High Temperatures, Russian Academy of Sciences, Moscow 125412, Russia Moscow Institute of Physics and Technology, Dolgoprudny 141700, Russia

a r t i c l e

i n f o

Article history: Received 1 October 2016 Received in revised form 7 December 2016 Accepted 8 December 2016

Keywords: Zirconium Niobium Zirconium alloys Molecular dynamics Point defects Diffusion Interatomic potentials

a b s t r a c t We report a new attempt to study properties of Zr-Nb structural alloys. For this purpose we constructed an angular-dependent many-body interatomic potential. The potential functions were fitted towards the ab initio data computed for a large set of reference structures. The fitting procedure is described, and its accuracy is discussed. We show that the structure and properties of all Nb and Zr phases existing in the Zr-Nb binary system are reproduced with good accuracy. The interatomic potential is appropriate for study of the high-pressure hexagonal x-phase of Zr. We also estimated characteristics of the point defects in a-Zr, b-Zr and Nb; results are proven to correlate with the existing experimental and theoretical data. In case of a-Zr the model reveals anisotropy of the vacancy diffusion, in agreement with previous calculations and experiments. The potential provides an opportunity for simulation of Zr-Nb alloys based on a-Zr and b-Zr. This conclusion is illustrated by the results obtained for the alloys with different niobium concentrations: up to 7% in case of hcp alloys and up to 50% for bcc alloys. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Zirconium alloys are widely chosen for structural applications in nuclear reactors. Originally, these alloys found their use as a material for fuel cladding tubes in water-cooled reactor. The choice was made due to low neutron absorption cross-section of zirconium, suitable corrosion resistance in water and appropriate mechanical characteristics [1]. Additionally, it was found that alloying with small amount of niobium (less than 5%) improves the corrosion resistance of the material without increase of the capacity for hydrogenation [2], the latter can cause formation of fragile hydride phase [3]. According to these requirements several alloy compositions are designed and manufactured. The most significant here are dilute Zr-Nb alloys, such as Zr-1%Nb (also marked as E110) and Zr-2.5%Nb (or E125). Structure of the given alloy is determined by its manufacturing procedure, which can include different heat-treatment and/or deformation [4,5]. Moreover, it has been observed that the Zr-Nb alloys can demonstrate microstructure changes under working conditions [6]. One of the reasons is that pure zirconium has three allotropic modifications existing at different pressures P and temperatures T. At ambient pressure and low temperature zirconium has a hexagonal close packed (hcp) structure, a-Zr. At T = 1135 K it transforms to body centred cubic (bcc) b-Zr. At T = 2128 K bcc ⇑ Corresponding author. E-mail address: [email protected] (D.E. Smirnova). http://dx.doi.org/10.1016/j.commatsci.2016.12.016 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

Zr melts. The third crystal modification, x-Zr, exists in the area of high pressures and has a hexagonal symmetry. According to the reported data, the pressure boundary of a $ x transition lies in a range between 2 and 6.5 GPa at room temperature [4]. One of the most discussed cases related to the structure changes in alloys concerns formation of the precipitates of different phases, for example, b-Nb or metastable b-Zr and x-Zr. At the same time, characteristics of the material (e.g. corrosion resistance) are strongly correlated with the alloys structure, so its deviations from the intended phase composition cause a change for the worse in durability. Understanding of these relations is required for extension of the materials lifetime, as well as for development of modified zirconium alloys with improved properties. To date many experimental and theoretical approaches have been used to investigate the correlation between structure changes and behaviour of the alloys. It has been stated that some of the precipitates have a positive impact on the alloys durability. For example, Ribeiro et al. [7] have mentioned that the configuration of Zr-1%Nb that consists of a-Zr matrix containing Nb in solid solution and fine b-Nb precipitates provides a combination of high creep resistance and small dimensional changes caused by irradiation. Experimental investigations of the role of precipitates in dilute Zr-Nb alloys have been also done in many works [6,8,9]. For example, the authors in [6] have established dependence of corrosion behaviour of Zr-Nb alloys on Nb content and cooling rate. It has been summarized (from the microstructural study and corrosion test) that the equilibrium Nb concentration below solubility limit (0.6 at.%Nb)

260

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

in a matrix plays more important role to enhance the corrosion resistance than the supersaturated Nb, the b-phase, and the precipitates. Thereby, many aspects can provide an impact on the resulting alloys characteristics. Nevertheless, one should keep in mind that experiments require material resources to determine nature of the exact precipitates in the given alloy. In such cases, it may be useful to involve theoretical approaches to support the experimental study. Previously several works have been aimed at theoretical investigations of structure of Zr-based alloys. Kharchenko and Kharchenko [2] have reported results of ab initio calculations that determined the structural properties of hcp and bcc Zr-Nb alloys with different Nb concentration. CALPHAD approach has been also applied for investigation of the Zr-Nb system. For example, Wang et al. [10] used CALPHAD model to trace concentration dependences of elastic constants of bcc solid solutions of Nb in Zr. The model parameters used in [10] have been based on available experimental and first-principles data. Molecular dynamics (MD) simulation has also found its application in investigation of Zr-Nb system. Previously MD has shown its advantages for study of structural characteristics and physical properties of matter at the atomistic level [11–14]. However, to apply this method to a specific material an appropriate interatomic potential is required [15]. In case of Zr-Nb system several interatomic potentials have been reported up to date. Most of them are aimed at study of pure Zr [16–18] and Nb [19]. Nonetheless, to our knowledge, two interatomic potential models have been developed for study of the binary alloys. Liang, Li and Liu have reported the potential for Zr-Nb system in a series of works [20,21]. In 2008 this research group has proposed a long-range empirical potential for Zr-Nb system, the model parameters have been summarized in [20]. The purpose of this potential was to investigate amorphization in Zr-Nb alloys and an interaction between Zr and Nb layers. The potential model adopted in [21] can reproduce melting temperatures of pure elements along with lattice and elastic constants, which have been used directly for fitting. The fitting base contains no data about point defects, results for defects energies were not presented. In 2013 Lin et al. developed an N-body potential for Zr-Nb system [22] on the basis of the embedded-atom method [23], the latter is also referred to as EAM. The test calculations show that the potential [22] is able to reproduce static properties of a-Zr, b-Zr and Nb, including characteristics of point defects (e.g. vacancies). However, the melting temperatures T m and thermal expansion coefficients of pure elements reported in [22] are underestimated in comparison with the experimental data. Thereby, it can be said that the potential models reported up to date for Zr-Nb system have different limitations. For example, they were not intended to be appropriate for study of high-pressure x-Zr. Also, previous MD works [20–22] did not pay significant attention to the behaviour of self-interstitial atoms (SIAs) in Zr and Zr-Nb alloys, while MD investigations of pure zirconium affirm the importance of both types of point defects [24]. A new interatomic potential that was developed in the present work is aimed on simulations of properties of different phases existing in Zr-Nb system, including binary alloys and x-Zr. For this purpose we use many-body interatomic potential model in an angular-dependent form. The paper is organized as follows. Section 2 contains information about the potential development technique and discussion about its accuracy. In Section 3 we show results obtained from simulations of zirconium phases with the constructed potential. Section 4 is devoted to study of pure niobium properties. In Section 5 the potential is validated by applying it for investigation of point defects characteristics and self-diffusion in Zr and Nb. In Section 6 we evaluate the potential application for simulation of the Zr-Nb binary alloys.

2. Development of a new interatomic potential for Zr-Nb In this work we develop the interatomic potential for Zr-Nb system based on the form of Angular-Dependent Potential (ADP) [25]. For the ADP the total potential energy U is given by the following formula:

U¼ þ

P i>j

uab ðrij Þ þ

P

iFa


i Þ þ 12 ðq

P

2 1X 1X 2 ðkkl Þ  m; i;k;l i i i 2 6

lki Þ2

i;k ð

ð1Þ

where

qi ¼

X

qb ðrij Þ; lki ¼

j–i

kkl i ¼

X wab ðrij Þr kij rlij ; j–i

X uab ðr ij Þr kij ; j–i

mi ¼

X kk ki :

ð2Þ

k

Here indices i and j enumerate atoms, while superscripts k; l ¼ 1, 2, 3 refer to the Cartesian components of vectors and tensors. Indices a and b denote the element types of atoms. The first term in Eq. (1) represents pair interactions between atoms via a pair potential u. The summation is over all j-th neighbours of i-th atom within the cut-off distance r cut ¼ 6:2 Å. The second term F is the embedding energy that is a function of the total electron density q. The two first terms in Eq. (1) give principal contribution to the system energy. Additional l and k terms introduce non-central interactions through the dipole vectors and quadrupole tensors. They are intended to penalize deviations of local environment from the cubic symmetry. The ADP potential form has been previously applied for study of various binary systems, such as, for example, Al-Cu [26], U-Mo [27] and U-N [28]. The ADP appears to be more sensitive, especially in case of non-cubic structures, compared with the original EAM potential. This conclusion is supported by the recent simulations that revealed local tetragonal distortions in the crystal structure of c-U and c-U-Mo alloys [29]. Thereby, we suppose that application of the angular-dependent model can improve accuracy of description for phases with hexagonal, tetragonal, orthorhombic, or more complex symmetry. In this work we constructed the potential for Zr-Nb system using the force-matching method [30]. This method provides a way to construct physically justified interatomic potentials from the fitting database, which does not contain experimental data. The idea is to adjust the interatomic potential functions to optimally reproduce per-atom forces (together with total energies and stresses) computed at the ab initio level for a fine-tuned set of reference structures. In this work optimization of the potential functions was performed with the help of the Potfit code [31,32]. The reference ab initio data were calculated using the DFT code VASP 5.2 [33]. Each of the reference structures (or so-called configurations) contains approximately 200 atoms in a simulation box with periodic boundary conditions (PBC). The exact number of atoms depends on the phase structure, density, and the number of defects included. The Brillouin zone was sampled with the 2  2  2 Monkhorst-Pack k-point mesh [34]. We tested the convergence of the energy in the studied structures using different k-points meshes, such as: 2  2  2; 3  3  3, and 3  4  4. According to these tests in our case the 2  2  2 mesh is suitable for bcc and hcp structures. The same setting is also extended to the case of hexagonal x-phase. The cut-off energy of a plane-wave basis set was equal to 520 eV. We use projector augmented wave pseudopotentials [35] included in the VASP package and the exchange-correlation functional within generalized-gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof. 3

Twelve electrons 4s2 4p6 5s1 4d for zirconium and eleven electrons

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272 4

4p6 5s1 4d for niobium were taken into account as valence electrons. Finally, for construction of the ADP we use 85 configurations containing 10,964 atoms altogether. These configurations represent hcp a-Zr, bcc b-Zr, hexagonal x-Zr, liquid Zr, bcc b-Nb, liquid Nb, bcc Zr-Nb alloys, and hcp Zr-Nb alloys. All these atomic configurations were obtained from classical atomistic simulations carried out at different temperatures and densities. Here, at the initial step of the work, we applied the potentials from [17,19]. MD and static calculations reported here and below were performed with the LAMMPS code [36]. Generally, the fitting contains the three steps, which reiterates to improve the description of the desired properties. (1) At the first step parameters of a new ADP potential are fitted to the reference ab initio database with the assistance of the Potfit code. As ‘‘parameters” here we mean coordinates of the spline nodes that are changed during fitting, starting from some trial version. According to the formalism of the ADP model we should find 13 potential functions to determine interatomic interactions in a binary system. The fitting procedure is directed by minimization of the target function Z given by the following sum:

Z ¼ Zf þ ZC ; Zf ¼

ZC ¼

Nf X X i¼1 a¼x;y;z 7NC X

ð3Þ ADP

DFT 2

W i ðf ia  f ia Þ ;

2

W j ðAADP  ADFT Þ : j j

ð4Þ

ð5Þ

j¼1

Z f and Z C denote two parts of the target function for forces and integral characteristics, respectively. In these equations the reference data are represented by per-atom forces f and integral characteristics A (energy and components of the stress tensor). N f is the total number of atoms for which the forces are taken into account during the minimization. NC is the total number of configurations. For each configuration we have one value of the energy and six components of the stress tensor. W i and W j are the weights attributed to different values taking part in the potential fitting. Index ‘‘DFT” denotes the reference value; ‘‘ADP” denotes the value computed with the fitted potential. (2) At the second step we perform primary verification of the potential through representation of the reference ab initio data (forces, energies and stresses). Generally, if we adopt the fitted potential version, we should proceed to step 3, otherwise, the fitting is repeated again. The potential developed in this work reproduces ab initio reference data with the average accuracy about 220.416 meV/Å for the forces, 17.442 meV for the energies and 0.1 GPa for the stresses. The numerical values of the spline nodes are summarized in the Appendix A.1 We also show the final potential functions in Fig. 1. To trace the precision of the potential fitting we determine the relative fitting errors for the forces averaged over each reference configuration:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 i N u i 1X tðF ADP  F DFT Þ ; DF ¼ 2 N i¼1 ðF iDFT Þ

ð6Þ

1 The potential in the LAMMPS setfl format can be downloaded from the NIST interatomic potential repository http://www.ctcms.nist.gov/potentials/.

261

here N is the number of atoms in a given configuration, F ADP is the force calculated with the developed potential, F DFT is the reference DFT force computed using VASP. Eq. (6) includes qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 absolute values of the forces: F ¼ ðf x þ f y þ f z Þ. The right panel in Fig. 2 shows the values of DF averaged on the configuration types. Quantitative distribution of the phases included in the fitting database is shown in left panel of Fig. 2. According to these figures we can say that the maximum DF corresponds to the high-pressure x-Zr. Nevertheless, an iterative potential fitting allows to improve representation of this phase: for the first examples of x-Zr included in the database DF was equal to 36%, while the fifth configuration (built with the iteratively refined test potential) was fitted with DF ¼ 18%. The following improvement of the x-Zr representation was not the purpose of this work because it deteriorates description of another zirconium phases, especially a-Zr. Thereby, the fitting results reflect the balance between optimal reproduction of different phases included in the database. We see that the accuracy of description lies within the range typical for the forcematching method, namely <20% (on average). (3) At the third step we apply the potential for classic MD simulations of the basic properties of Zr and Nb phases, together with Zr-Nb alloys. If the potential gives insufficient accuracy for description for the high-priority properties we edit initial set of configurations and repeat fitting. It is important to note that intermediate improved potential versions obtained in the process of fitting can be applied to prepare additional configurations during revision of the reference database. In this section we discuss only primary verification of the potential through quality of representation of the reference ab initio data. Information about the basic properties of Zr phases, Nb and Zr-Nb alloys calculated with a new potential for Zr-Nb system is summarized in the following sections in comparison with the existing experimental data and results of previous computations. 3. Basic characteristics of Zr phases The potential should be qualified by ability to reproduce known properties of the zirconium phases. In this section we compare our calculation results obtained with the existing experimental data or with the information computed from the first-principles. Also it is important to trace the difference between the new ADP and the potentials proposed for Zr and Zr-Nb previously by the other authors [17,22]. For this purpose comparison between different interatomic potentials is included in the discussion. In the following subsection we consider verification tests for hcp Zr. Next subsection deals with bcc phase, and after that we continue with discussion of high-pressure x-Zr. 3.1. Simulation of hcp a-Zr We evaluated how the ADP reproduce lattice parameters, elastic constants and cohesive energy of the hcp Zr. Table 1 contains calculated values for all of these properties, together with the experimental results (at T = 300 K). In our case none of the data gathered in Table 1 have been used directly in the fitting process. In spite of this fact, the accuracy of representation is about 0.5–1.2% for the lattice constants, about 2% for the cohesive energy and less than 22 % for C 11 ; C 13 , and C 33 . Significantly higher overestimation is seen only for C 12 . We also give the value of C 44 that can be compared to the DFT data reported previously. In Table 1 we summarize the characteristics obtained with another interatomic potentials – the one widely used for simulations of pure zirconium [17], and the potential aimed at study of Zr-Nb alloys [22].

262

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

Fig. 1. Potential functions of the created ADP. Types of the interatomic interactions related to each of the functions are indicated in the corresponding legend.

Fig. 2. Left panel: quantitative distribution of different configuration types used for the potential development. Right panel: relative force-fitting errors D averaged on various types of configurations. We use the following indication: a-Zr, b-Zr and x-Zr – hcp, bcc and hexagonal phases of pure Zr, respectively; L-Zr and L-Nb – liquid zirconium and niobium; D-Zr – bcc zirconium containing point defects; b-Nb – bcc Nb; a – hcp-based solid solutions of Nb in Zr; b – bcc-based solid solutions of Nb in Zr. Only the configurations with the attributed fitting forces are included in the figures.

Table 1 Properties of hcp Zr at zero pressure. Results calculated with ADP are presented in comparison with the existing experimental data and computational for the cohesive energy Ec [37], the equilibrium lattice parameters a and c [38], and elastic constants C ij [39].

Ec (eV) a (Å) c (Å) c=a C 11 (GPa) C 12 (GPa) C 13 (GPa) C 33 (GPa) C 44 (GPa)

Experiment

ADP (this work)

EAM [22]

EAM [17]

DFT [40]

6.32 3.232 5.148 1.6 143.4 72.8 65.3 164.8 –

6.44 3.216 5.21 1.62 150.5 120.5 65 129 8.5

6.25 3.2171 5.1621 1.6046 128.3 74.4 79.7 180.9 –

6.635 3.234 5.168 1.598 147 69 74 168 –

3.240 5.178 1.598 141.1 67.6 64.3 166.9 25.8

We also calculated the room-temperature isotherm of hcp zirconium. During this simulation the ratio c=a was not strongly fixed, so it could slightly change with the rise in pressure to satisfy the following condition: P xx ¼ P yy ¼ P zz ¼ P. Results of this computation are shown in Fig. 3, in comparison with the experimental

data presented by Zhao et al. [41]. Changes of the lattice parameters show the same behaviour as in the experiments. Experimental study [41] also contains information about volumetric thermal expansion aV of hcp Zr. According to the experiments aV is in the range between 1.9  105 K1 [42] and 2.5  105 K1 [41]. For

263

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

Fig. 3. Changes of the lattice parameters of hcp Zr depending on pressure. Left panel: volume of the lattice unit cell. Right panel: ratio of the lattice parameters c=a. 1 – calculations with the ADP potential at room temperature, 2 – experimental results from [41]. Dashed lines are plotted to guide the eye.

Table 2 High-temperature properties of bcc Zr at zero pressure. Results calculated with ADP are given in comparison with the existing experimental data for the lattice parameter a [41], elastic constants C ij [44] and melting temperature Tm.

a (Å) C 11 (GPa) C 12 (GPa) Tm (K)

Experiment

ADP (this work)

EAM [17]

EAM [22]

3.627 104 93 2128

3.627 92 87 1780

3.576 – – 2109

– – – 1580

Fig. 4. Isotherm of bcc Zr at T = 980 K. 1 – experimental data tabulated in [41], 2 – simulations with the ADP. Lines are plotted to guide the eye.

the model of Zr based on ADP the mean aV is found to be equal to 1.9  105 K1 (over the temperatures 300–1000 K), in agreement with experiments of Cleri et al. [42]. The value of aV computed in this work also agrees with the data measured from hightemperature investigations of linear thermal expansion of zirconium performed in [43]. According to these experiments mean (averaged over all directions) linear thermal expansion coefficient for a-Zr lies in the range between 5.8 and 6.8  106 K1. 3.2. Simulation of bcc b-Zr The potential developed in this work can be also used for simulation of high-temperature bcc Zr phase. We can see from Table 2 that the lattice parameter of the MD model is close to the experimental value (both are related to 970 K). To our knowledge, the

values of elastic constants C ij based of the experimental measurements for b-Zr have been reported only for three temperatures: 1188 K, 1483 K and 1883 K, see [44]. During discussion of this phase we can additionally rely on results of CALPHAD approach presented in [10]. Elastic constants obtained with the CALPHAD model are close to the experimental ones and equal to 99.7 GPa for C 11 and 87.9 GPa for C 12 at 1188 K. Elastic constants C 11 and C 12 computed in this work for 1188 K are within 6–11% from their experimental values reported in Table 2 for the same temperature. The prediction also agrees very well with the result computed with the CALPHAD. The table also contains information about the melting temperature T m . The value of melting temperature obtained by Lin et al. was obtained using the two-phase approach [11,45,46]. Such type of simulation means that one part of the 3D box (with PBC) is filled by the bcc lattice of zirconium atoms and another part is filled by liquid zirconium. Depending on the temperature and pressure the phase boundary can move or remain stable. In the latter case both phases coexist in equilibrium that gives the value of the melting temperature at the given pressure. In our work melting process was simulated through the heating of bcc zirconium model with an open surface [47]. Such model set up prevent overheating of the model and provides reliable value of the melting temperature at zero pressure. According to our simulation made the ADP predicts T m that is about 15% lower than the experimental value. Analyses of pressure changes in bcc Zr under compression at high temperature (970 K) shows that in general the MD model project the same trend as the existing experiments. Fig. 4 illustrates the comparison. For this simulation we use the model containing 2000 atoms of zirconium in a bcc structure, with PBC applied along x; y and z axes. We also estimated the value of linear thermal expansion coefficient alin for bcc Zr. According to the experiments [48,49] thermal expansion coefficient of bcc Zr is constant up to the melting point, and it is about 10  106 K1. The values measured for polycrystalline b-Zr in [43] are lower – the maximum here is about 7.8  106 K1. At the same time, experiments made in [50] show

264

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

Table 3 Properties of hexagonal x-Zr at zero temperature and zero pressure. Results calculated with ADP in this work are presented in comparison with the existing ab initio data for the lattice parameters a and c (both in Å) and elastic constants C ij (in GPa). The last string contains experimental data for lattice constants.

ADP [40] [51] [54] [58]

a

c

c/a

C 11

C 12

C 13

C 33

C 44

5.077 5.056 – 5.036 5.039

3.050 3.150 – 3.152 3.150

0.600 0.623 – 0.625 0.625

174 165.2 167 161.7 –

77 75.6 69.8 72.6 –

80 47.5 47.7 53.5 –

190 198.7 195.7 195.6 –

50.8 30.6 33.7 33.7 –

Table 4 Properties of bcc Nb at zero pressure. Results calculated with the ADP are presented in comparison with the existing experimental and calculation data for the lattice parameter a [60], cohesive energy Ec [37], elastic constants C ij and bulk modulus B [61].

a (Å) Ec (eV) C 11 (GPa) C 12 (GPa) C 44 (GPa) B (GPa) alin 106(K1)

Experiment

ADP (this work)

ab initio [19]

EAM [22]

EAM [19]

3.303 7.47 253 133 31 173 7 – 13

3.285 7.48 248 127 51.3 161 10.5

3.309 7.10 251 133 22 172 –

3.300 7.57 262.2 124.6 35.95 170.2 9.8

3.308 7.09 244 136 32 172 –

Fig. 5. Cold curve of bcc niobium: 1 – results of MD simulations of compression performed with the ADP, T = 5 K; 2 – results of ab initio computations performed in [19] at zero temperature; 3 and 4 – experiments [59] in multi-anvil and in DAC, respectively.

linear changes of alin with temperature: in the area of stability for this phase alin rises from 9  106 K1 to about 10.5  106 K1. Simulations with the ADP at T = 970–1800 K results in a constant alin that has an intermediate value equal to 8.1  106 K1. 3.3. Simulation of x-Zr Another important crystal phase of zirconium exists in highpressure region of the T-P phase diagram. Originally, x-Zr has hexagonal non-close packed structure. The data reported previously on experimental and calculated phase boundaries of x-Zr in pure zirconium allow to conclude that x-Zr appears in the range of pressures from 2.5 to 6 GPa (at room temperature). With following rise in pressure it becomes stable up to P  30 GPa, and than transforms to b-Zr [51]. For example, in recent experimental work reported in 2015 by Ono and Kikegawa [52] room-temperature x $ b transition in Zr was found at 33 GPa. Several attempts have been made to achieve theoretical description of x-Zr. Previously

different research groups [40,51,53,54] reported results of the ab inito calculations that treated structure and elastic moduli of x-Zr at zero temperature and pressure. To our knowledge, there is a lack of information about study of this phase with the help of interatomic potentials. At the same time, for titanium, whose phase diagram is similar to the zirconium one, the existing interatomic potential [55] helped to describe phase transformations between a-Ti, b-Ti and x-Ti, including the cases of shockinduced phase transitions [56,57]. The reason for the lack of attempts to create such model for zirconium may be attributed to the limited amount of data about x-Zr that are appropriate for the empirical potential fitting. At the same time, the forcematching method adopted in this work generally requires information only about the phase structure. That is one of the preferences of this potential development technique comparing to the empirical approach. Hence, the method may be applied for consideration of complex phases described by very few data, namely x-Zr. From the Table 3 we see that the results obtained by different authors for the elastic moduli of x-Zr are quite similar. According to the first-principles calculations made in [40], x-Zr is the most stable zirconium phase at T = 0 K. At T = 300 K situation changes and stability of a-Zr prevails. In our work also we created a model of x-Zr at zero temperature. Simulations were made with the hexagonal structure where positions of atoms satisfied the basis requirements for the x-phase. The model contains 1680 atoms in a simulation box with PBC. In our statical calculations with the ADP we see that x-Zr is stable at zero temperature, in agreement with previous theoretical results [40]. Elastic constants of the x-Zr simulated with the ADP are listed in Table 3, in comparison with the first-principles calculations. The combination of computed elastic moduli indicates mechanical stability of the phase. If we turn to the MD simulations at T = 300 K and zero pressure, it appears that x-structure becomes mechanically unstable at these conditions, in consistency with the experimental phase diagram. This fact suggests that one has to apply pressure to stabilize x-Zr at room temperature. The ADP potential predicts that xphase is stable in a range of pressures between 9 GPa and 12 GPa. This interval is determined from the pressure dependence of enthalpy computed in this work. During the computations the c=a ratio in x-Zr was kept constant, and such setting allows to avoid distortion of the pressure deviators. Thus, the examples presented above illustrate that the potential describes stability of x-Zr

265

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272 Table 5 The formation energy of different point defects in hcp Zr (in eV). All computational results are given for zero temperature. Defect type

ADP (this work)

DFT [67]

DFT [72]

EAM [24]

EAM [22]

Experiments [64]

Vacancy SIA

1.7 3.1

– 2.75

2.07 –

1.76 2.88

1.79 –

>1.5 –

at T = 0 K, as it was found in the first-principles computations. At the same time it allows to study this phase at higher temperatures (300 K) corresponding to an actual experimental observation. However, the pressure interval of stability determined for x-Zr in this work is rather narrow compared to another experimental and theoretical estimations. Summarizing the results obtained during investigation of Zr phases, we can say that the new ADP may provide background for further investigation of zirconium phase diagram. 4. Investigation of niobium properties We primarily tested the static properties of niobium at zero temperature and pressure: lattice parameters, cohesive energy, and elastic constants. All of these data are gathered in Table 4, together with information obtained from previous experiments and ab initio computations. For the sake of comparison results related to different interatomic potentials are also listed in the table. The ADP potential reproduces lattice parameter and cohesive energy with a high accuracy (see Table 4). Hierarchy of the elastic constants also indicates stability of bcc niobium phase at T = 0 K. Analysing the data gathered in the Table 4 one can say that the ADP developed in this work inferiors to the other models in description of ground properties of niobium. Nevertheless, this conclusion would not be generally correct because the ADP shows high accuracy in description of high-temperature and highpressure behaviour of niobium. For example, Fig. 5 illustrates pressure changes in zero-temperature niobium under compression. We also plotted the points obtained from the recent experimental compression [59] of bcc Nb at the room temperature. Two sets of the symbols denoted as 3 and 4 correspond to the Multi-anvil and DAC experiments, respectively. The cold curve calculated in this work agrees well with both sets of the experimental data. As an example of high-temperature property it would be important to discuss description of thermal expansion coefficient. According to the experiments [50,62], linear thermal expansion coefficient alin of niobium is found to have a temperaturedependent behaviour, the boundary values at T = 300 K and near the melting temperature are noted in Table 4. Simulations made with the ADP at different temperatures (T = 300 K–2000 K) result in a nearly constant alin that lies in the middle of the experimentally estimated interval (see Table 4). The value of alin obtained in this work is also consistent with the estimations based on results computed in [19] for pure Nb with the help of EAM potential. The melting temperature of bcc Nb computed at zero pressure is equal to 2500 K, while in experiments it is determined to be 2750 K. The difference between these values is less than 10 %. 5. Defects in pure zirconium and niobium phases Knowledge about driving forces leading to radiation damage is required since zirconium is used as a base for structural materials in nuclear engineering. Experimental and theoretical approaches have been practised to describe possible structure changes that can occur in real zirconium and its alloys due to formation of radiation defects. According to the experimental observations [63] the microstructure of Zr and its alloys during irradiation has proven to be increasingly complicated. One of the reasons for structure

Fig. 6. Diffusion coefficients computed along different directions in the hcp zirconium lattice. Results of the present work are shown by symbols. Diffusivities for a system with vacancy is shown for: 1 – diffusion in the basal plane and 2 – diffusion along c direction. Results for a system with a single SIA are also given for diffusion in the basal plane (3) and along the c-axis (4). Diffusivities along a and c directions computed in [74] for a system with vacancy are shown by lines 5 and 6, respectively.

changes is formation of the dislocation loops. It is worth mentioning that not only such difficult defect structures as a dislocation loop have been a subject of experimental studies. Characteristics of point defects (e.g. SIA and vacancy) are also important, as they can provide a key for understanding behaviour of more complex defect configurations. Migration and formation energies of a single f vacancy – Em v and Ev have been estimated for the hcp zirconium in the experimental works by Hood [64,65]. As for any experimental observation of single interstitial, to our knowledge, this information is not available for zirconium. Classical MD, together with the first-principles computations, have been also provided a help to understand the atomistic configurations of defects in zirconium phases. LeBacq et al. reported results of full-potential quantum calculations of vacancy formation energies for group IV metals, including Zr. Williame et al. [66], Samolyuk et al. [67,68] and Peng et al. [69] focused an attention on comparison between energetics of different locations of SIAs in hcp Zr. Formation energies of point defects have been also considered with application of different interatomic potentials [17,70]. On the one hand, such method has some limitations due to the fact

that the energies (for example, Evf ) are usually used directly for potential fitting, so they appears to be predetermined. Nevertheless, classical MD with an interatomic potential is able to provide an insight in defects diffusion. For example, the potential reported by Mendelev [17] has been adopted for detailed calculations of diffusion coefficients in hcp and bcc Zr [24]. Returning to the discussion of larger size of defects, one should mention a series of work dedicated to MD simulation of high-energy displacement cascades in a-zirconium [71]. We used an interatomic potential for estimation of energies related to formation and migration of point defects. For this purpose the well-known formalism is applied. We compute the energy of an ideal (defectless) system and compare it with the energy related to the system containing defect. Formation energies for SIA and vacancy are determined as follows:

266

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272 Table 6 The formation energy of different point defects in bcc Zr (in eV). Defect type

ADP this work

DFT [72]

EAM [24]

Vacancy SIA

2 2.1

2.3 –

1.88 2.05

work agrees with the other values predicted previously with the help of EAM potentials [22,24].

Fig. 7. Self-diffusion coefficients for the hcp zirconium lattice: 1 – results of this work obtained from MD with the ADP potential; 2 – data from the tracer diffusion experiments [75]; 3 – experimental results from Hood et al. [76] for self-diffusion in (4) and out () of the hcp basal plane; 4 – coefficients computed in [24] for selfdiffusion via SIA.

f

Ev ac ¼ EN1 

5.1.2. Defects diffusion We performed MD study of defect diffusion to estimate corresponding diffusion coefficients and migration energies. For this purpose a single defect was included in an hcp a-Zr lattice. For diffusion of SIA we have considered the temperature range from 600 K to 1100 K. The time length of the runs was equal to 2 ns. The systems with a single vacancy have been studied at temperatures from 800 K up to 1080 K, and requires more simulation time – from 3.5 up to 10 ns, depending on temperature. During simulations of diffusion we observe that migration of point defects appears to be anisotropic. Diffusivities depend on the given direction in the unit cell, so they are denoted as follows:

D? ¼

N1 EN ; N

ð7Þ

f ESIA ¼ ENþ1  Nþ1 EN : N

Here N is the number of atoms in the model of crystal lattice without any defects, EN is the energy of this model. EN1 corresponds to the model that contains a single vacancy and ENþ1 is for the model with a single SIA. In the present work we apply this calculation method at different conditions: at T = 0 K and at high temperatures. The type of calculation depends on the area of stability for each of the discussed phases. In the zero temperature case the energies included in Eq. (7) are obtained from statical computations, while for the high-temperature phases the energies were averaged over the MD run (see [24,27] for details). 5.1. Study of formation and migration of point defects in hcp a-Zr 5.1.1. Defects formation energies At first the defect formation energies were determined from statical relaxation at T = 0 K. The system contains 3072 atoms (in an ideal hcp structure) with PBC along all three axes. The values f of Evf ac and ESIA at zero pressure and T = 0 K are summarized in Table 5. We also traced how the formation energies change with rise in temperature from T = 620 K to 1180 K. For this case the energies from Eq. (7) were estimated over the MD run equal to f 1 ns. We see that for a-Zr ESIA seems to be constant (3.5 eV) in this

temperature range while Evf ac increases from 2.0 eV to 2.4 eV. We used the mean value equal to 2.2 eV to estimate self-diffusion coefficients in the following subsections. Statistical uncertainty of the resulted MD values at the given temperature is about 0.1 eV (averaged over results of several computation runs). According to the analyses [67,68] performed recently by means of electronic structure calculations within the framework of DFT (using VASP code) the minimum formation energy is related to the basal octahedral (BO) type of interstitial defect. From statical calculations at zero temperature with the developed potential we f see the same result. Moreover, ESIA of such defect is quite close to the value obtained for the same configuration by Samolyuk et al., see Table 5. The difference between these two values is about f

PN

2 i¼1 ðDxi

þ Dy2i Þ

4t

;

Dk ¼

PN

2 i¼1 ðDzi Þ

2t

:

ð8Þ

Here Dxi ; Dyi ; Dzi are displacements of i-th atom along x; y and z directions. In simulation x is parallel to the a-axis, and z is parallel to the c-axis (xy-plane coincides with the basal plane). N is the number of atoms in a simulated system with a single defect, and t is the calculation time. Diffusivities Dk and D? are equal up to a constant factor to the diffusion coefficients of point defects [73,27]. We see that vacancies show a tendency to diffuse more quickly in the basal plane, rather than along the c-axis. In other words it means that Em v ac is smaller for the migration in the basal plane than for the migration out of this plane. For SIA this effect is less pronounced in simulations, so D?  Dk . Difference in diffusivities is clearly seen in Fig. 6 where we plot D? and Dk for a systems containing vacancy and SIA. Anisotropy of diffusion in hcp Zr have been previously discussed in several theoretical works, for example in [68,74]. From Fig. 6 we see that our results for a vacancy diffusion are close to the computational results [74] obtained for the same defect type (based on thermodynamic integration). 5.1.3. Self-diffusion in a-Zr The calculation results discussed in the above subsection allow us to estimate self-diffusion coefficient for hcp Zr. To determine the prevailing diffusion mechanism we discuss self-diffusion guided by the different point-defects (vacancies and SIAs) separately. For this purpose we use the following equation:

Dself ¼ cv ac Dv ac þ cSIA DSIA :

ð9Þ

The two terms of the equation describe vacancy and SIA contributions. Here Dv ac=SIA are diffusivities of atoms in a system with a single vacancy or SIA. In this case we did not considered anisotropy of the self-diffusion, so the Dv ac and DSIA are averaged over all directions. cv ac=SIA are equilibrium concentrations of the given defects, which can be estimated from the defect formation energies:

cv ac=SIA ¼ exp 

F vf ac=SIA kB T

! ¼ exp

Svf ac=SIA kB

! exp 

Evf ac=SIA kB T

! ;

ð10Þ

where F vf ac=SIA and Svf ac=SIA are the variations in free energy and

10%. For a vacancy we see that Ev ac lies within the range deter-

entropy with the formation of the defect of the specified type, and

mined in the experiments. We also see that Evf ac computes in this

f kB is the Boltzmann constant. For Evf ac and ESIA we chose the values

267

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

Fig. 8. Left panel: diffusion of atoms in a bcc Zr via SIA (1) and vacancy (2). Right panel: self-diffusion coefficients of bcc Zr. 1 – experimental data from Lundy [78]; 2 – results obtained in this work for diffusion via SIA (from DSIA ); 2 – results obtained in this work for diffusion via vacancies (from Dv ac ).

Table 7 The formation energy of different point defects in bcc Nb (in eV). Defect type

ADP this work

DFT [19]

EAM [19]

EAM [22]

Experiments [82]

Vacancy SIA

2.6 4

2.72 3.95–4.89

3.10 3.83–4.50

2.75 –

2.6–3.1 –

computed in classical MD runs at non-zero temperatures. The entropy terms Svf ac=SIA were assumed to be the order of unity. In this f case we took expðSvf ac Þ equal to 3 and expðSSIA Þ equal to 1. The final self-diffusion coefficient is plotted in Fig. 7, in comparison with the experimental points [75,76] and previous MD results reported by Mendelev and Bokstein [24]. Comparison of the vacancy and interstitial contributions (see Eq. (9)) shows that in our work the final self-diffusion coefficient is mostly determined by vacancies. To support this conclusion we can say that the contribution from SIA is about 5 orders of magnitude lower in the whole studied temperature range. So the ADP potential predicts that vacancies are a leading defect for self-diffusion in hcp Zr. Fig. 7 shows that the resulted self-diffusion coefficients computed within this approach are close to the later experimental data from [76]. We also see that different interatomic potentials predict different nature of the self-diffusion in a-Zr. For example, the potential reported by Mendelev et al. [17,24] predicts that the interstitial atoms play more important role in self-diffusion than the vacancies. The line (4) plotted in Fig. 7 was obtained in [24] under this assumption. However, to our knowledge, there is no experimental information that could solve this contraversion and determine the actual character of self-diffusion in a-Zr.

5.2. Study of point defects and self-diffusion in b-Zr In case of bcc Zr experimental information about defects properties is even more limited than for hcp Zr. However, experiments show that the self-diffusion coefficient rises significantly after transition from a-Zr to b-Zr. Activation energy for self-diffusion in bcc Zr is also lower than for the normal bcc metals. These features are common for the group IV transition metals (Zr, Ti, Hf) and can be explained through the profile of the phonon spectra [77]. However, this is not the one specialty of self-diffusion in bcc Zr. The diffusion experiments of Lundy [78] performed for bZr also revealed unusual curvature of the Arrhenius plot appearing at T about 1500 K. According to the existing assumptions, this effect can be caused by different factors, such as influence of the phase transitions occurring in Zr, variation of the leading defect

Fig. 9. Lattice parameters of the hcp Zr-Nb alloys containing different amount of niobium (at 300 K). 1 – results of simulations with the ADP; 2 – experimental data from [83], 3 – experiments reported by Benites et al. in [84], 4 – experimental points from [85].

configuration or its formation energy, or interaction with the impurities [79,80]. Explanation of such anomalous behaviour of bcc Zr self-diffusion coefficient is a challenge that has not been

268

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

solved yet, so these fundamental effect deserves detailed investigations. 5.2.1. Point defects formation and diffusion in bcc Zr In this work we used the MD model based on ADP potential to estimate defect formation energies and self-diffusion coefficients in bcc Zr. As far as b-Zr is stable at T > 1135 K, it seems adequate to determine its defect characteristics at high temperatures. f However, if we try to obtain Evf ac and ESIA using Eq. (7) we see that at such high temperatures thermal fluctuations may obstruct averaging the energies. Based on this, we decided to estimate formation energies through the calculation of defects concentration (where it is possible). For this purpose we built a series of models representing bcc Zr with an open surface [24,27,81]. Each model contains 193,200 atoms total, including 39,200 atoms in a surface region and 154,000 atoms in a diffusion region (representing the bulk volume). Such model setting implies that the system does not contain preliminary created point defects, so they should be generated and regulated by the surface. Time-length of each calculation was equal to 4 ns. During the computation we traced the mean-squared deviations of atoms in the diffusion region and than converted it into self-diffusion coefficients, similarly to Eq. (8). We tested different temperatures in the area of stability of b-Zr: T was equal to 1450 K, 1510 K, 1630 K. Each of the models was instantly frozen, so type and number of defects can be traced and converted to the defect concentration. Than the defect formation energy was obtained from Eq. (10). The technique described above works well for SIA defects. However, to determine characteristics of vacancies we had to apply usual method described in the beginning of the section and including Eqs. (7). Resulted energies describing formation of point defects in b-Zr are summarized in Table 6. To our knowledge, in case of bcc Zr existing ab initio data are limited to the vacancy formation energy, while characteristics of SIAs were considered in the MD study [24]. Bokstein and Mendelev have shown that in case of Zr interstitial atoms can provide a contribution in self-diffusion. The result based on this assumption appears to be comparable with the existing experimental data on self-diffusion activation energy. From Table 6 one can see that f the values of Evf ac and ESIA obtained by Mendelev and Bokstein are very close. The same correlation is demonstrated in present work. As in case of hcp Zr, we have also studied defects migration in the bcc phase. Simulations of defects diffusion were carried out at temperatures T = 1000 – 1800 K and zero pressure. Each of the models contains 2000 atoms in a simulation box (with PBC), where a single defect was included. MD runs last from 3.4 up to 5 ns for a SIA and from 4.5 to 10 ns for a vacancy. Diffusion is isotropic, and coefficients were determined with the help of formulas similar to Eq. (8). Diffusivities of Zr atoms via vacancy and SIA are shown in Fig. 8, see the left panel. Using these results we estimated migration energies for point defects, dashed lines show the points involved in these estimations. The most high-temperature point was excluded due to proximity to the melting point. For bcc Zr we observe isotropic diffusion with Evmac = 0.35 eV, and

ESIA m = 0.17 eV. These results generally agrees with the migration v ac energies computed previously in [24], where Em = 0.46 eV, and ESIA m = 0.11 eV. According to the right panel of the Fig. 8, we can conclude that the ADP potential can describe temperature dependence of Dself in bcc Zr at 1500 K < T < 1800 K. The absolute values of Dself directed by the SIA diffusion are close to the experimental data. However, currently the potential does not reflect changes in curvature of the Arrhenius plot. We can analyse this result with respect to probable reasons of the anomalous diffusion, which were assumed previously for Zr [79]. Our interatomic potential did not reveal any

Fig. 10. Lattice parameters of the bcc Zr-Nb alloys. 1 – results of simulations with the ADP (T = 600 K); 2 and 3 – experimental data from [86,84], respectively (T = 300 K).

Fig. 11. Enthalpy of mixing estimated for the bcc Zr-Nb alloys. 1 – results of DFT computations performed in [87]; 2 – termodynamic estimations from [88]; 3 – DFT computations from [22]; 4 – simulation with EAM potential for Zr-Nb (from [22]); 5 – enthalpy computed with the ADP for Zr-Nb at 1000 K and zero pressure, this work. Symbols 6 and 7 show mixing enthalpy estimated in this work for the ordered alloys at 1000 K: star (7) corresponds to the fcc L12 Zr3Nb structure, and a triangle (6) shows result for the bcc B2 ZrNb.

influence of the phase transitions that could possibly lead to the change in diffusion coefficients. Also the ADP potential did not show variations of the defect formation energies, which could also be responsible for the anomalous character of the self-diffusion in Zr. Constant values of the formation energies were obtained previously in simulation with the other potential for pure Zr [24]. In this way we can suggest that the anomalous self-diffusion observed experimentally in bcc Zr requires further investigation that will be continued by the authors of this work and reported elsewhere. To conclude this subsection we should say that the developed ADP can be applied for simulation of point defects characteristics in bcc Zr. Herewith, results of the verification tests show that the ADP potential predicts SIA diffusion in bcc Zr and can partially describe temperature-dependence of the self-diffusion coefficient (at 1500 K < T < 1800 K).

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

Fig. 12. Zero-temperature energy estimated for the ordered Zr-Nb alloys with different structures. Symbols 1, 2 and 3 show results of DFT computations performed in [22] for L12 Zr3Nb, B2 ZrNb and L12 ZrNb3. Lines 10 , 20 , 30 are obtained for the same structures from simulations with the ADP potential (this work).

5.3. Formation of point defects in bcc Nb For niobium we obtained defect formation energies from statical computations at zero temperature. Previously it was reported that self-diffusion in Nb is determined only by vacancies [19]. Nevertheless, in this work we study formation of both types of point f ESIA

for niobium calculated defects. From Table 7 one can see that in present work agrees well with the minimum value computed by the first-principles techniques [19]. For a vacancy we observe agreement with the lower boundary of the available range of data gathered from experiments. In addition, the hierarchy of defects predicted for Nb with the help of ADP agrees with the theoretical data obtained previously by another authors. 6. Simulation of hcp and bcc Zr-Nb alloys Originally, structure of the given Zr-Nb alloy is determined by its manufacturing procedure, which includes heat treatment along with different types of deformation procuring desired properties. Full description of the alloy structure requires taking into account different aspects dealing, for example, with distribution of components and existence of precipitates. In the present work we considered Zr-Nb solid solutions based on bcc and hcp allotropic modifications of Zr. As it was discussed above, structural alloys aimed on application in nuclear field, usually contains less than 2.5% of Nb, and crystal structure of the alloys is generally based on hcp Zr. Here we apply the developed ADP potential for investigation of hcp alloys containing up to 7 at.% of niobium. Structure of such alloys is determined by the two lattice parameters (ahcp and chcp ), which show a tendency to become smaller with addition of niobium. According to the experiments made by different research groups, the ratio between these values is kept constant and close to the one for pure zirconium. Experimental data for ahcp ; chcp , and chcp =ahcp obtained at 300 K are summarized in Fig. 9. We create several models of random solid solutions containing various amount of niobium. Each model has an hcp structure with the unit lattice parameters that provide a zero pressure. The results of MD simulations with the ADP are also shown in the Fig. 9. The lattice parameters of the simulated alloys decrease with addition of niobium. Calculated values of ahcp are generally consistent with the experimental data, while chcp is overestimated. Nevertheless, in case of chcp the slope of the parameter changes in the studied concentra-

269

tion range is quite similar to the one estimated from the experiments. The lowest panel in Fig. 9 shows that the ratio chcp =ahcp obtained from MD simulations keeps its constant value with rise of niobium content. This fact also approves that the created ADP provides correct description of lattice parameters of different hcp Zr-Nb alloys at the room temperature. Another verification test carried out in this work deals with simulation of bcc Zr-Nb alloys. For this purpose we built a series of models corresponding to the random solid solutions of Nb in bcc Zr. In these simulations structure of the alloys has been observed over a different range of concentrations – up to 50% of Nb. The temperature was kept equal to 600 K, and pressure was set to 0 GPa in each calculation point. For verification of the calculation results we present a set of experimental data obtained in [86,84] for metastable bcc alloys at room temperature. Comparison between experimental and calculation results is made in Fig. 10. The ADP clearly provides a linear reduction in the lattice parameter abcc in response to increasing of niobium content in the alloy. This fact matches well with the experimental data, and is consistent with the Vegard’s law. Moreover, the slope of this calculated linear dependence is the same as the one traced towards the experimental points (see the dash-dotted line plotted in Fig. 10). We also estimated enthalpy of mixing (Hmix ) for the bcc Zr-Nb random alloys. For this purpose we leaned on the regular solutions model. Enthalpy of mixing Hmix was determined as follows:

Hmix ¼ HA1x Bx  ð1  xÞHA  xHB ;

ð11Þ

here HA and HB are the internal energies corresponding to the pure bcc Zr and bcc Nb at P = 0 and T = 1000 K; HA1x Bx – energy of the given random bcc alloy calculated at P = 0 and T = 1000 K; x – atomic fraction of niobium in the alloy. HA1x Bx was estimated for the three reference alloy compositions: Zr25Nb75, Zr50Nb50 and Zr75Nb25. For each of the cases we performed a number of computational runs with different distribution of the Zr and Nb atoms in a system. The HA1x Bx that corresponds to a system with the lowest internal energy was included in the Eq. (11). The resulted curve for Hmix of the disordered alloy is plotted in Fig. 11. The estimated Hmix is positive in the studied concentration range. It also depends on the alloy composition and has a maximum for Zr50Nb50 equal to about 5 kJ/mole. To our knowledge, there are no experimental data about the mixing enthalpy in Zr-Nb system. Nevertheless, we can compare our estimations with the previous theoretical results [22,87,88]. Barannikova et al. computed mixing enthalpies for hcp, bcc and hypothetical fcc Zr-Nb alloys using ab initio methods. According to their results the maximum value of Hmix for random bcc alloys is about 11.5 kJ/mole [87]. Thermodynamical assessment performed in [88] provided the value that is significantly lower – about 3 kJ/mole. Both of these curves are plotted in Fig. 11. We also can see that the estimation obtained in the present work with the help of the ADP potential lies between these two reference values. This figure also contains results obtained by Lin et al. with application of the EAM potential designed in their work [22]. We see from Fig. 11 that the Hmix described by this EAM potential lies significantly higher than all other estimations, including the DFT computations made by Lin et al. As an additional verification we tested high-temperature stability of hypothetical ordered phases. For this purpose we simulated several ordered phases with cubic structure: L12 Zr3Nb, B1 ZrNb, B2 ZrNb and L12 ZrNb3. All of them were tested in zero-pressure MD runs with the ADP Zr-Nb potential, temperature was set equal to 1000 K. According to the simulations results we can see that B1 ZrNb and L12 ZrNb3 are mechanically unstable under the given conditions while B2 ZrNb and L12 Zr3Nb seem stable over the tested time length (25 ps). This stability gives us an opportunity

270

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

to check the mixing enthalpy of these ordered phases at high temperature (calculated with respect to pure bcc Zr and bcc Nb). From Fig. 11 one can see that Hmix of bcc structure B2 is higher than Hmix for the random solution, while the mixing enthalpy for fcc L12 Zr3Nb is lower. Here we should mention that comparison of the Hmix is not enough to make a comprehensive analysis of the relative stability of ordered and disordered alloys at high temperatures. In high-temperature cases one should compare Gibbs energies G taking into account the entropy term (G = H-TS). Temperature dependence of the Gibbs energy is mostly determined by the entropy factor, which tends to be larger in case of disordered alloys. In this work we do not consider calculations of the free energies leaving it for further work. However, to enlarge verification via simulation of the hypothetical ordered structures we calculated energies of L12 Zr3Nb, B2 ZrNb and L12 ZrNb3 (at zero pressure and temperature) and compare them with the predictions obtained in [22] by using firstprinciples methods. Fig. 12 shows that the ADP potential reproduce exactly the same energy hierarchy as the one found in ab initio: EZr3 Nb > EZrNb > EZrNb3 .

7. Conclusion We presented and discussed a new interatomic potential aimed on study of structure and properties of the Zr-Nb alloys. The chosen fitting technique based on ab initio data allows to consider various phases existing in the system. Thus, the potential has a wide application area. The structure and general properties of three Zr crystal phases (a; b; x) existing in the Zr-Nb system are reproduced with a good accuracy. Especially it is worth noting that the potential is appropriate for study of high-pressure hexagonal x-phase. The ADP predicts that this zirconium phase is more stable than the other two in the pressure range between 9 and 12 GPa (at room temperature). Formation and migration energies of point defects are calculated for pure zirconium (a-Zr and b-Zr) as well as for bcc niobium. The results correlate with the existing experimental and theoretical data. It is important to note that the potential allows to observe anisotropic diffusion of vacancies in a-Zr, in accordance with the experiments. We also computed self-diffusion coefficients for bcc and hcp Zr. In case of a-Zr we see agreement with the existing experimental points. For b-Zr temperature dependence of the self-diffusion coefficient agrees with the experiments only partially, in the limited temperature interval. The potential gives an opportunity for simulation of Zr-Nb alloys with different structure and niobium content. In this report we approve adequate description of the lattice parameters for alloys based on a-Zr and b-Zr. Results of the potential verification summarized in this work allow to conclude that the developed ADP can be adopted for further investigation of zirconium phases, niobium and Zr-Nb alloys. Such vast variety of applications creates an opportunity for simulation of different complex structures, for example, including precipitates of the different phases.

Acknowledgements The calculations were carried out on the computer clusters MVS-100 K and MVS-10P (Joint Supercomputer Center of RAS) and ‘‘Lomonosov” (Moscow State University). The work was supported by the Program for Basic Research of the Presidium of the RAS 2 (coordinator is G. I. Kanel) and by the personal RF President Scholarship CP-199.2016.2 (D. E. S.) We are grateful to Dr. Alexey Kuksin for useful discussions and help with the development of the

interatomic potential. We also thank Prof. G. Norman for inspiration and support. Appendix A. The spline nodes of developed interatomic potential The interatomic potential developed in present work is specified by the spline nodes coordinates. In the tables below all potential functions are listed with the appropriate accuracy that allows to reproduce the potential for further application (see Tables A8–A13). Table A.8 Pair potential functions describing interactions of the following types: Zr-Zr, Zr-Nb and Nb-Nb, in eV. r (Å)

uZrZr ðrÞ

uZrNb ðrÞ

uNbNb ðrÞ

2.000000 2.247059 2.494118 2.741176 2.988235 3.235294 3.482353 3.729412 3.976471 4.223529 4.470588 4.717647 4.964706 5.211765 5.458824 5.705882 5.952941 6.200000

2.617885E+00 1.162526E+00 3.650032E01 1.951015E02 1.702512E01 2.120906E01 1.854105E01 1.311871E01 1.029225E01 6.063654E02 1.176603E02 1.632798E02 2.058145E02 8.765689E03 1.902092E03 4.835123E03 3.502197E03 0.000000

1.119329E+00 4.791570E01 4.470297E02 2.564771E01 3.104731E01 2.819147E01 2.144447E01 1.202948E01 1.783750E02 2.705878E02 2.087403E02 1.904803E02 1.610182E02 1.040849E02 4.474590E03 1.855767E03 3.503145E03 0.000000

– 1.607431E01 4.755347E01 4.819580E01 4.201582E01 3.296576E01 2.359412E01 1.131440E01 1.144510E02 6.600576E02 5.093681E02 3.808615E02 2.714059E02 1.883069E02 8.931215E03 3.550315E03 1.707270E03 0.000000

Table A.9 Electronic density functions q. r (Å)

qZr ðrÞ

qNb ðrÞ

2.000000 2.247059 2.494118 2.741176 2.988235 3.235294 3.482353 3.729412 3.976471 4.223529 4.470588 4.717647 4.964706 5.211765 5.458824 5.705882 5.952941 6.200000

1.083931E01 7.645973E02 6.107197E02 5.178194E02 4.204537E02 3.340608E02 2.531429E02 1.792409E02 1.371399E02 9.812501E03 6.458962E03 3.504674E03 4.574662E04 9.335764E04 2.349034E04 6.053589E04 4.136959E04 0.0000000

3.071360E01 1.565449E01 1.047126E01 6.781236E02 4.262434E02 3.294232E02 3.595219E02 3.155264E02 1.456222E02 3.332318E04 4.964953E04 1.296050E03 2.373624E03 2.036755E03 1.761418E03 9.907992E04 1.743727E04 0.0000000

Table A.10 embedded-atom functions F depending on the q, in eV. The F Nb ðqÞ is equal to 0 up to q = 0.205632.

q

F Zr ðqÞ

q

F Nb ðqÞ

0.0 0.072057 0.174898 0.277738 0.380578 0.483419 0.586259 0.689099 0.791939 0.894780

0.0 1.172797 3.822570 4.841703 5.134119 5.140360 4.969064 4.682909 4.277176 3.775535

– – 0.205632 0.308537 0.423781 0.539025 0.654268 0.769512 0.884756 1.000000

– – 0.000274 4.285508 5.181522 5.165299 4.949128 4.572652 4.131006 3.641580

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272 Table A.11 pffiffiffiffiffiffi ADP pair functions uðrÞ, in eV/Å. r (Å)

uZrZr ðrÞ

uZrNb ðrÞ

uNbNb ðrÞ

2.000000 2.466667 2.933333 3.400000 3.866667 4.333333 4.800000 5.266667 5.733333 6.200000

– 4.414937E02 4.182905E04 2.260217E02 9.401307E03 2.010107E02 5.177207E03 2.086232E03 1.180333E03 0.000000

8.294129E02 2.301147E02 2.571542E02 1.397535E02 5.448136E03 1.834118E05 1.168566E03 2.746916E03 7.201014E04 0.0000000

2.584908E02 4.100657E02 4.704426E02 2.637125E02 2.081906E02 8.916558E03 2.865232E03 4.732415E03 2.594833E03 0.00000

Table A.12 pffiffiffiffiffiffi ADP pair functions wðrÞ, in eV/Å2. r (Å)

wZrZr ðrÞ

wZrNb ðrÞ

wNbNb ðrÞ

2.000000 2.466667 2.933333 3.400000 3.866667 4.333333 4.800000 5.266667 5.733333 6.200000

– 6.352157E02 4.022481E02 1.894085E02 1.787584E02 9.381230E03 2.914564E03 7.687580E04 4.369823E05 0.0000000

1.369368E01 6.878460E02 4.030413E02 3.229168E02 1.722774E02 2.266340E03 2.919476E03 2.241204E03 4.868934E04 0.0000000

– 8.133130E02 4.903825E02 2.987643E02 1.483954E02 2.198034E03 3.790830E03 1.633282E03 4.423067E04 0.0000000

Table A.13 The derivatives of the potential functions necessary for setting up the boundary conditions. Here r0 (or q0 ) and rcut (or qcut ) are minimum and maximum coordinates 0 corresponding to the range where thepgiven function is p defined. ffiffiffiffiffiffi ffiffiffiffiffiffi 3 The values of u are 1 given in eV/Å; q0 are in Å ; u0 are in eV/Å2; w0 are in eV/Å ; F 0 is in eV.

u0ZrZr u0ZrNb u0NbNb q0Zr q0Nb

F 0Zr F 0Nb u0ZrZr u0ZrNb u0NbNb w0ZrZr w0ZrNb w0NbNb

r0 (Å)

rcut (Å)

q0

qcut

6.4975 2.6371 1.4902 0.1453 0.7125 – – 0.0326 0.3244 0.2495 0.0089 0.1654 0.0490

0 0 0 0 0 – – 0 0 0 0 0 0

– – – – – 11.9978 66.77 – – – – – –

– – – – – 5.0608 4.1304 – – – – – –

References [1] B. Cox, Some thoughts on the mechanisms of in-reactor corrosion of zirconium alloys, J. Nucl. Mater. 336 (2–3) (2005) 331–368. . [2] V.O. Kharchenko, D.O. Kharchenko, Ab-initio calculations for the structural properties of Zr-Nb alloys, Condens. Matter Phys. 16 (1) (2013) 13801, arXiv: arXiv:1206.7035v2 . [3] M. Veshchunov, V. Shestak, V. Ozrin, A new model of hydrogen redistribution in Zr alloy claddings during waterside corrosion in a temperature gradient, J. Nucl. Mater. 472 (2016) 65–75, http://dx.doi.org/10.1016/j.jnucmat.2016. 01.032. [4] A. Zhilyaev, I. Sabirov, G. Gonzalez-Doncel, J. Molina-Aldaregua, B. Srinivasarao, M. Prez-Prado, Effect of Nb additions on the microstructure, thermal stability and mechanical behavior of high pressure Zr phases under ambient conditions, Mater. Sci. Eng: A 528 (9) (2011) 3496–3505. . [5] B.B. Straumal, A.S. Gornakova, A.A. Mazilkin, O. Fabrichnaya, M. Kriegel, B. Baretzky, J.-Z. Jiang, S. Dobatkin, Phase transformations in the severely plastically deformed ZrNb alloys, Mater. Lett. 81 (2012) 225–228. . [6] Y.H. Jeong, K. Ok, H.G. Kim, Correlation between microstructure and corrosion behavior of Zr Nb binary alloy, J. Nucl. Mater. 302 (2002) 9–19.

271

[7] R.M. Ribeiro, C.B.A. Woyames, L.H. de Almeida, D.S. dos Santos, Effect of microstructure and addition of alloying elements on hydriding kinetics of ZrNb-based alloys, Int. J. Hydrogen Energy (2015) 1–10. . [8] H. Yang, Y. Matsukawa, S. Kano, Z. Duan, K. Murakami, H. Abe, Investigation on microstructural evolution and hardening mechanism in dilute Zr-Nb binary alloys, J. Nucl. Mater. 481 (2016) 117–124. [9] Q. Dong, H. Yu, Z. Yao, F. Long, L. Balogh, M. Daymond, Study of microstructure and precipitates of a Zr-2.5Nb-0.5Cu CANDU spacer material, J. Nucl. Mater. 481 (2016) 153–163, http://dx.doi.org/10.1016/j.jnucmat.2016.09.017. [10] X. Wang, L. Liu, M. Wang, X. Shi, G. Huang, L. Zhang, Computational modeling of elastic constants as a function of temperature and composition in Zr-Nb alloys, Calphad 48 (2015) 89–94. . [11] A.B. Belonoshko, Molecular dynamics of MgSiO3 perovskite at high pressures: Equation of state, structure and melting transition, Geochim. Cosmochim. Acta 58 (19) (1994) 4039–4047. [12] K.P. Zolnikov, A.V. Korchuganov, D.S. Kryzhevich, V.M. Chernov, S.G. Psakhie, Structural changes in elastically stressed crystallites under irradiation, Nucl. Instrum. Methods Phys. Res., Sect. B 352 (2015) 43–46, Proceedings of the 12th International Conference on Computer Simulation of Radiation Effects in Solids, Alacant, Spain, 8–13 June, 2014 . [13] D.K. Belashchenko, Molecular dynamics simulation of the structure and thermodynamic properties of liquid rubidium at pressures of up to 10 Gpa and temperatures of up to 3500 K, Russ. J. Phys. Chem. A 90 (9) (2016) 1707– 1716, http://dx.doi.org/10.1134/S0036024416090053. [14] O.L. Kuskov, D.K. Belashchenko, Thermodynamic properties of Fe-S alloys from molecular dynamics modeling: Implications for the lunar fluid core, Phys. Earth Planet. Inter. 258 (2016) 43–50. . [15] G.E. Norman, V.V. Stegailov, Stochastic theory of the classical molecular dynamics method, Math. Models Comput. Simul. 5 (2013) 305–333. [16] F. Willaime, C. Massobrio, Development of an N -body interatomic potential for hcp and bcc zirconium, Phys. Rev. B 43 (1991) 11653–11665, http://dx.doi.org/ 10.1103/PhysRevB.43.11653. [17] M.I. Mendelev, G.J. Ackland, Development of an interatomic potential for the simulation of phase transformations in zirconium, Phil. Mag. Lett. 87 (5) (2007) 349–359. . [18] M.J. Noordhoek, T. Liang, Z. Lu, T.-R. Shan, S.B. Sinnott, S.R. Phillpot, Charge-optimized many-body (COMB) potential for zirconium, J. Nucl. Mater. 441 (1-3) (2013) 274–279. . [19] M.R. Fellinger, H. Park, J.W. Wilkins, Force-matched embedded-atom method potential for niobium, Phys. Rev. B 81 (14) (2010), http://dx.doi.org/10.1103/ PhysRevB.81.144119. [20] S. Liang, J. Li, B. Liu, Solid-state amorphization of an immiscible Nb-Zr system simulated by molecular dynamics, Comput. Mater. Sci. 42 (4) (2008) 550–557. . [21] J. Li, Y. Dai, X. Dai, T. Wang, B. Liu, Development of n-body potentials for hcpbcc and fcc-bcc binary transition metal systems, Comput. Mater. Sci. 43 (4) (2008) 1207–1215. . [22] D.-Y. Lin, S.S. Wang, D.L. Peng, M. Li, X.D. Hui, An n-body potential for a Zr-Nb system based on the embedded-atom-method, J. Phys.: Condens. Matter 25 (20) (2013) 209501. . [23] M.S. Daw, M.I. Baskes, Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals, Phys. Rev. B. 29 (12) (1984) 6443–6453. [24] M.I. Mendelev, B.S. Bokstein, Molecular dynamics study of selfdiffusion in Zr, Phil. Mag. 90 (5) (2010) 637. [25] Y. Mishin, M. Mehl, D. Papaconstantopoulos, Phase stability in the Fe-Ni system: investigation by first-principles calculations and atomistic simulations, Acta Mater. 53 (15) (2005) 4029–4041. . [26] F. Apostol, Y. Mishin, Interatomic potential for the Al-Cu system, Phys. Rev. B 83 (2011) 054116, http://dx.doi.org/10.1103/PhysRevB.83.054116. [27] D.E. Smirnova, A.Y. Kuksin, S.V. Starikov, Investigation of point defects diffusion in bcc uranium and U-Mo alloys, J. Nucl. Mater. 458 (2015) 304–311. [28] A. Kuksin, S.V. Starikov, D.E. Smirnova, V.I. Tseplyaev, The diffusion of point defects in uranium mononitride: combination of DFT and atomistic simulation with novel potential, J. Alloy. Compd. 658 (2016) 385–394, http://dx.doi.org/ 10.1016/j.jallcom.2015.10.223. [29] S.V. Starikov, L.N. Kolotova, Features of cubic and tetragonal structures of UMo alloys: atomistic simulation, Scripta Mater. 113 (2016) 27–30. . [30] F. Ercolessi, J.B. Adams, Interatomic potentials from first-principles calculations: the force-matching method, Europhys. Lett. 26 (1994) 583. [31] P. Brommer, F. Gähler, Effective potentials for quasicrystals from ab-initio data, Phil. Mag. 86 (2006) 753–758. [32] P. Brommer, A. Kiselev, D. Schopf, P. Beck, J. Roth, H.-R. Trebin, Classical interaction potentials for diverse materials from ab initio data: a review of potfit, Modell. Simul. Mater. Sci. Eng. 23 (7) (2015) 074002. .

272

D.E. Smirnova, S.V. Starikov / Computational Materials Science 129 (2017) 259–272

[33] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B. 54 (16) (1996) 11169. [34] H.J. Monkhorst, J.D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13 (1976) 5188–5192, http://dx.doi.org/10.1103/PhysRevB.13.5188. [35] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projectoraugmented wave method, Phys. Rev. B 59 (3) (1999) 1758–1775. [36] S.J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. [37] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 2005. [38] D.R. Lide, CRC Handbook of Chemistry and Physics: A Ready-Reference Book of Chemical and Physical Data, CRC Press, Boca Raton, FL, 2009. [39] C. Fisher, E.S. Renken, Single-crystal elastic moduli and the hcp-bcc transformation in Ti, Zr, and Hf, Phys. Rev. 135 (2A) (1964) A482. [40] Y.-J. Hao, L. Zhang, X.-R. Chen, Y.-H. Li, H.-L. He, Phase transition and elastic constants of zirconium from first-principles calculations, J. Phys. Condens. Matter: Inst. Phys. J. 20 (23) (2008) 235230. . [41] Y. Zhao, J. Zhang, C. Pantea, J. Qian, L.L. Daemen, P.a. Rigg, R.S. Hixson, G.T. Gray, Y. Yang, L. Wang, Y. Wang, T. Uchida, Thermal equations of state of the a , b , and x phases of zirconium, Phys. Rev. B 71 (18) (2005) 184119, http://dx. doi.org/10.1103/PhysRevB.71.184119. [42] F. Cleri, F. Rosato, Tight-binding potentials for transition metals and alloys, Phys. Rev. B. 48 (1993) 22. [43] V. Petukhov, Thermal expansion of zirconium in the solid phase, High Temp. High Press. 35/36 (2003/2004) 15–23. [44] A. Heiming, W. Petry, J. Trampenau, M. Alba, C. Herzig, H.R. Schober, G. Vogl, Phonon dispersion of the bcc phase of group-IV metals. II. bcc zirconium, a model case of dynamical precursors of martensitic transitions, Phys. Rev. B 43 (1991) 10948–10962, http://dx.doi.org/10.1103/PhysRevB.43.10948. [45] J.R. Morris, C.Z. Wang, K.M. Ho, C.T. Chan, Melting line of aluminium from simulations of the coexisting phases, Phys. Rev. B. 49 (5) (1994) 3109–3115. [46] S. Starikov, V. Stegailov, Atomistic simulation of the premelting of iron and aluminium: Implications for high-pressure melting-curve measurements, Phys. Rev. B 80 (22) (2009) 22104. [47] A. Kuksin, G. Norman, V. Stegailov, A. Yanilkin, Surface melting of superheated crystals. atomistic simulation study, Comput. Phys. Commun. 177 (1–2) (2007) 34–37, http://dx.doi.org/10.1016/j.cpc.2007.02.073. [48] A. Heiming, W. Petry, W. Trampenau, W. Miekeley, J. Cockcroft, The temperature dependence of the lattice parameters of pure BCC Zr and BCC Zr-2 at.%Co, J. Phys. Cond. Mat. 4 (3) (1992) 727, http://dx.doi.org/10.1088/ 0953-8984/4/3/012. [49] A.F. Guillermet, Critical evaluation of the thermodynamic properties of zirconium, High Temp. High Press. 19 (2) (1987) 119–160. [50] G. Aurelio, A.F. Guillermet, G.J. Cuello, J. Campo, Structural properties and stability of metastable phases in the Zr-Nb system: Part II. Ageing of Bcc (b) alloys and assessment of b-Zr, Metall. Mater. Trans. A 34A (December) (2003) 2771–2779. [51] S. Zhang, X. Zhang, Y. Zhu, S. Zhang, L. Qi, R. Liu, First-principles investigations on elastic and thermodynamic properties of zirconium under pressure, Comput. Mater. Sci. 61 (2012) 42–49. . [52] S. Ono, T. Kikegawa, Determination of the phase boundary of the omega to beta transition in Zr using in situ high-pressure and high-temperature X-ray diffraction, J. Solid State Chem. 225 (2015) 110–113. . [53] C.W. Greeff, Phase changes and the equation of state of Zr, Modell. Simul. Mater. Sci. Eng. 13 (7) (2005) 1015–1027. . [54] B.-T. Wang, P. Zhang, H.-Y. Liu, W.-D. Li, P. Zhang, First-principles calculations of phase transition, elastic modulus, and superconductivity under pressure for zirconium, J. Appl. Phys. 109 (6) (2011) 063514. . [55] R.G. Hennig, T.J. Lenosky, D.R. Trinkle, S.P. Rudin, J.W. Wilkins, Classical potential describes martensitic phase transformations between the a; b, and x titanium phases, Phys. Rev. B 78 (2008) 054121. . [56] H. Zong, T. Lookman, X. Ding, C. Nisoli, D. Brown, S.R. Niezgoda, S. Jun, The kinetics of the x to a phase transformation in Zr, Ti: analysis of data from shock-recovered samples and atomistic simulation, Acta Mater. 77 (2014) 191–199. . [57] X. Zong, T. Lookman, X. Ding, S.-N. Luo, J. Sun, Anisotropic shock response of titanium: Reorientation and transformation mechanisms, Acta Mater. 65 (2014) 10–18. . [58] C.S. Barrett, T.B. Massalski, Structure of Metals, McGraw-Hill, New-york, 1996.

[59] Y. Zou, X. Qi, X. Wang, T. Chen, X. Li, D. Welch, B. Li, High-pressure behavior and thermoelastic properties of niobium studied by in situ x-ray diffraction, J. Appl. Phys. 116 (1) (2014) 013516, http://dx.doi.org/10.1063/1.4887436. [60] R. Roberge, Lattice parameter of niobium between 4.2 and 300 K, J. Less Common Met. 40 (1975) 161, http://dx.doi.org/10.1016/0022-5088(75)90193-9. [61] G. Simmons, H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook, second ed., The MIT Press, Cambridge, 1971. [62] K. Wang, R. Reeber, The role of defects on thermophysical properties: thermal expansion of V, Nb, Ta, Mo and W, Mater. Sci. Eng.: R: Rep. 23 (3) (1998) 101– 137, http://dx.doi.org/10.1016/S0927-796X(98)00011-4. [63] M. Griffiths, A review of microstructure evolution in zirconium alloys during irradiation, J. Nucl. Mater. 159 (1988) 190–218. [64] G.M. Hood, Diffusion and vacancy properties in a-Zr, J. Nucl. Mater. 139 (1986) 179–184, http://dx.doi.org/10.1016/0022-3115(86)90170-4. [65] G.M. Hood, Point defect diffusion in a-Zr, J. Nucl. Mater. 159 (1988) 149–175. [66] F. Willaime, Ab initio study of self-interstitials in hcp-Zr, J. Nucl. Mater. 323 (23) (2003) 205–212. . [67] G. Samolyuk, S. Golubov, Y. Osetsky, R. Stoller, Self-interstitial configurations in hcp Zr: a first principles analysis, Phil. Mag. Lett. 93 (2) (2013) 93–100. . [68] G.D. Samolyuk, A.V. Barashev, S.I. Golubov, Y.N. Osetsky, R.E. Stoller, Analysis of the anisotropy of point defect diffusion in hcp Zr, Acta Mater. 78 (2014) 173–180. [69] Q. Peng, W. Ji, H. Huang, S. De, Axial ratio dependence of the stability of selfinterstitials in HCP structures, J. Nucl. Mater. 437 (1-3) (2013) 293–296. . [70] C. Willaime, D. Massobrio, Development of an N-body interatomic potential for hcp and bcc zirconium, Phys. Rev. B 43 (14) (1991) 11653. [71] S.J. Wooding, L.M. Howe, F. Gao, A.F. Calder, D.J. Bacon, A molecular dynamics study of high-energy displacement cascades in a alpha-zirconium, J. Nucl. Mater. 254 (1998) 191–204. [72] O. Le Bacq, F. Willaime, Unrelaxed vacancy formation energies in group-IV elements calculated by the full-potential linear muffin-tin orbital method: invariance with crystal structure, Phys. Rev. B 59 (13) (1999) 8508–8515. [73] A.Y. Kuksin, D.E. Smirnova, Calculation of diffusion coefficients of defects and ions in UO2, Phys. Solid State 56 (6) (2014) 1214–1223, http://dx.doi.org/ 10.1134/S1063783414060201. [74] H. Wen, C. Woo, Temperature dependence and anisotropy of self- and monovacancy diffusion in a-Zr, J. Nucl. Mater. 420 (1-3) (2012) 362–369. . [75] J. Horváth, F. Dyment, H. Mehrer, Anomalous self-diffusion in a single crystal of a-zirconium, J. Nucl. Mater. 126 (3) (1984) 206–214. . [76] G.M. Hood, H. Zhou, R.J. Schultz, N. Matsuura, J.A. Roy, J.A. Jackman, Self- and Hf diffusion in a-Zr and in dilute, Fe-free Zr(Ti) and Zr(Nb) alloys, Defect Diffus. Forum 143–147 (1997) 49–54. [77] H. Mehrer, Diffusion in Solids: Fundamentals, Methods, Materials, DiffusionControlled Processes (Springer Series in Solid-State Sciences), Springer, 2007. [78] T. Lundy, Diffusion in body-centered cubic metals zirconium, vanadium, niobium, and tantalum (thesis), report ORNL-3617, Tech. rep., 1964. [79] J.M. Sanchez, D. de Fontaine, Model for anomalous self-diffusion in group-IVb transition metals, Phys. Rev. Lett. 35 (1975) 227–230, http://dx.doi.org/ 10.1103/PhysRevLett.35.227. [80] G. Neumann, V. Tölle, Self-diffusion in body-centred cubic metals: analysis of experimental data, Phil. Mag. A 61 (4) (1990) 563–578, http://dx.doi.org/ 10.1080/01418619008231935. [81] M.I. Mendelev, Y. Mishin, Molecular dynamics study of self-diffusion in bcc fe, Phys. Rev. B 80 (2009) 144111, http://dx.doi.org/10.1103/PhysRevB.80.144111. [82] H. Ullmaier, Atomic Defects in Metals, Springer, 1991. [83] B.D. Lichter, Precision lattice parameter determination of zirconium-oxygen solid solution, Trans. Met. Soc. AIME 218 (1960) 1015. [84] G.M. Benites, A. Fernandez, G.J. Cuello, J. Campo, Structural properties of metastable phases in Zr-Nb alloys, J. Alloy. Compd. 299 (2000) 183–188. [85] J.P. Guèrillon, Metaux, Corrosion, Industrie 557 (1972) 21. [86] G.B. Grad, J.J. Pieres, A.F. Guillermet, G.J. Cuello, J.R. Granada, R.E. Mayer, Systematics of lattice parameters and bonding distances of the omega phase in Zr-Nb alloys, Physica B 214 (1995) 433–435. [87] S.A. Barannikova, A.M. Zharmukhambetova, A.Y. Nikonov, A.V. Dmitriev, A.V. Ponomareva, I.A. Abrikosov, Influence of stresses on structure and properties of Ti and Zr-based alloys from first-principles simulations, IOP Conf. Ser.: Mater. Sci. Eng. 71 (1) (2015) 012078. [88] S.K. Menon, S. Banerjee, R. Krishnan, Application of free energy-composition diagrams in predicting the sequences of phase transformations in Zr-Nb alloys, Metall. Trans. a 9 (Issue 9) (1978) 1213–1220.