The effect of impurity niobium on diffusion of self-interstitial atom and self-diffusion by interstitial mechanism in zirconium: an atomistic simulation

The effect of niobium on diffusion of self-interstitial atom and self-diffusion by interstitial mechanism has been simulated by the molecular dynamics method in zirconium at a temperature up to 1000 K. It has been shown that even a small amount of niobium impurity in the HCP Zr matrix entails significant qualitative and quantitative changes in the diffusion coefficient of a self-interstitial atom and self-diffusion by the interstitial mechanism.


Introduction
For the last decades, zirconium-niobium alloys (1%-2.5%Nb) have been successfully exploited in nuclear engineering [1] as a material of fuel cladding and fuel assembly elements in pressurized water reactors (PWRs). When operating, the alloys suffer from strong neutron irradiation leading to the change in their properties. The diffusion mobility of point defects produced in atomic displacement cascades is an important characteristic that determines the process of radiation damage of zirconium and its alloys.
There are two methods of atomistic simulation employed to study the migration of point defects in a solid: molecular dynamics (MD) simulation and Kinetic Monte Carlo (KMC) method. Also, two combined MD-KMC methods have been proposed in Refs. [2,3]. Numerous studies (see [2][3][4][5][6][7][8][9][10]) highlight an anisotropic behavior of migration of the point defects in α-Zr. Such anisotropy is shown to play a key role in the radiation growth of zirconium [4,5,11,12]. A self-interstitial atom (SIA) migrates mainly in the basal plane (0001). This feature agrees with the last ab initio studies [10,13,14], where the SIA position in the basal plane is found as the most energetically preferable. Also, it has been shown that the axial diffusion coefficient is well described by the Arrhenius law whereas, in the basal plane, the diffusion coefficient sometimes does not obey the Arrhenius law. Woo et al [8] and Tikhonchev et al [3] have shown that some interatomic potentials predict a small decrease in basal diffusion coefficient with an increase in temperature. The influence of doping impurities on the migration of point defects in zirconium has not been investigated so far. In [15] reporting on atomistic simulation of radiation damage of Zr-Nb alloys, it has been shown that the presence of niobium impurity in the α-Zr can significantly change the SIA diffusion properties. This phenomenon is explained by the high positive binding energy of the Nb atom with the SIA configuration in the HCP Zr matrix. However, only room temperature (T=300 K) is considered in [15].
Here, we report on the atomistic simulation of the SIA-configuration diffusivity and self-diffusion by the interstitial mechanism in the Zr-0.5%Nb binary alloy at the temperature up to 1000 K. The molecular dynamics (MD) method with a semi-empirical n-body interatomic potential is used for simulation. The ZrNb system is considered as a random substitutional solid solution. The considered niobium concentration (0.5 wt.%) is close to the Nb fraction in α-Zr phase of E110 ® and M5 ® alloys (Zr-1%Nb). The results are compared with the data obtained for the pure HCP Zr, some of them have already been reported in paper [3]. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Methods
The semi-empirical n-body interatomic potential for a Zr-Nb system proposed by Lin et al [16] is used for simulation in Zr-Nb system. It was developed on a base of the embedded atom method [17] and fitted to the results of experiments and first-principles calculations for pure elements and some binary compositions. A slightly modified version of this potential has been employed earlier to simulate the displacement atomic cascades in ZrNb alloy [15]. In our study, we also apply this modified version. In [15], this potential is shown to predict the energetically preferable positions of SIA for pure HCP Zr in the basal plane, which agrees qualitatively with the ab initio calculations [14]. However, this potential slightly overestimates the values of the SIA formation energy and predicts the basal crowdion (BC) as the most preferable interstitial position. This prediction is in disagreement with the ab initio results [14], where the basal octahedral position is the most stable. Since the SIA diffusion is studied in the presence of Nb impurities, the energetic properties of Nb atoms in the zirconium matrix containing the SIAs should be properly reproduced. The used potential provides substitution energy of a single Nb atom in Zr bulk of 0.47 eV and binding energy of 0.57 eV between the Nb atom and the SIA configuration in Zr. The corresponding ab initio evaluations are 0.61 and 0.68 eV, respectively [18]. Therefore, the chosen potential provides a reasonable agreement with the main ab initio results. Noteworthy, this potential does not reproduce the α↔β phase transformation of zirconium: HCP-phase is stable up to the melting temperature. Nevertheless, this feature is not critical because the temperature rage considered in the present study lies below the transition temperature of 1136 K [19].
The HCP crystallites of pure Zr and Zr-0.5Nb alloy were considered. The crystallite had a size of ∼77×78×57 Å and contained ∼15000 atoms. Periodical boundary conditions were used. For alloy, each atom of a crystallite was set as the Nb atom with probability of 0.005 and as the Zr atom with probability of 0.995. Since the crystallite is sufficiently large, this simple way allows us to construct a random alloy appropriate for the considered problem. Therefore, neither coherent potential approximation (CPA) [20] nor Special Quasirandom Structures (SQS) [21] approach has been employed. The MD simulation was performed using the velocity Verlet numerical scheme [22] with a time step of 1-2 fs. The crystallite temperatures of 300, 400, 500, 600, 700, 800, 900, and 1000 K were considered. First, the lattice parameters were evaluated for each temperature by simulating an isothermal-isobaric (NPT) ensemble with a zero pressure followed by relaxation of the crystallite with an additional atom placed in the BC interstitial position for 10 ps at the given temperature and lattice parameters. The Berendsen thermostat was used for the temperature control. Then, the crystallite behaviour was simulated as a microcanonical (NVE) ensemble. No temperature control was performed in this stage. The temperature fluctuations observed during the simulation did not exceed 1%. The simulations were carried out in several statistically independent parallel flows. The total simulated time of the SIA migration was within the range of 400-800 ns depending on temperature.
During simulation, the SIA migration and self-diffusion by interstitial mechanism were studied. Change of the SIA position was traced by the method of Wigner-Seitz cells. The cell containing two atoms was interpreted as a place of SIA. The coordinates of the lattice site corresponding to this cell were considered as the current coordinates of SIA. The SIA diffusion coefficient was determined by the Einstein formula The diffusion was two-and one-dimensional, respectively. Also, the sums of the squared displacements of all atoms in the basal plane and along the [0001] axis were calculated: where N is the number of atoms, t is the time, Δt is the time interval. It is known (see, e.g., [23] and references therein) that the average values of R xy 2 and R z 2 are described by the following equation.

Results
The obtained evaluations of the SIA-diffusion coefficient and self-diffusion coefficient for pure Zr are given in figures 1 and 2, respectively, in the Arrhenius coordinates. Noteworthy, figure 1 summarizes our results obtained earlier in [3] to be compared with the new results. A slight difference is observed for the radial diffusion coefficient Da SIA in the temperature range of 300-600 K. At the temperatures above 600 K, the value  of Da SIA is well described by the Arrhenius ratio with the activation energy of Eact=0.03 eV. The evaluations of Dc SIA are well described by the Arrhenius equation at Eact=0.09 eV for all temperature ranges under consideration. The values of self-diffusion coefficients Da are well described by the Arrhenius equation with a the parameters changing at the same temperature of 600 K. Activation energy of migration in the basal plane is evaluated as 0.025 and 0.04 eV for the temperature intervals of 300-600 and 600-1000 K, respectively. The calculated values of Dc are well described by the single Arrhenius equation at Eact=0.07 eV. Both SIA migration and self-diffusion run mainly in the basal plane.
The following feature of the SIA migration in Zr-0.5%Nb alloy has been demonstrated. If the initial SIAconfiguration consists of the Zr atoms only, it migrates relatively fast until it is trapped by the Nb atom. Mixed Zr-Nb SIA-configuration migrates distinctly slower and is quite stable due to high positive binding energy. Then, one or more niobium atoms joint to SIA-configuration. In that case, an interstitial Nb atom and one or more other Nb atoms are observed in vicinity. Since the mobility of such SIA-configuration is even lower, the diffusion coefficients are evaluated only for the temperature of 600 K and higher. The obtained values of   1 and 3 allows us to conclude that mobility of SIA in alloy is significantly lower than in pure Zr. This phenomenon becomes more significant with the temperature decrease. At 600 K, the difference in the diffusion coefficient reaches three orders of magnitude. It is worth mentioning, that there is an important qualitative difference in anisotropy of SIA migration. Namely, Dc<Da for pure Zr and Dc>Da for the case of alloy.
The calculated values of the self-diffusion coefficient for the Zr and Nb atoms are plotted in figures 4 and 5, respectively. The results show that self-diffusion coefficient of the Zr atoms in alloy is 2-3 orders of magnitude lower than for pure Zr and almost isotropic. The migration of the Nb atoms is anisotropic at Dc>Da. The temperature dependence of coefficients Da and Dc is reasonably approximated by the Arrhenius equation at the activation energy of ∼0.6 eV. At the temperature of 600 K, the values of Da and Dc for niobium atoms are close in order of magnitude to Dc and Da values for the Zr atoms in pure Zr, respectively. The difference between these values increases with the temperature and gets 30-fold increase at 1000 K.

Conclusion
MD simulation of the single SIA migration in pure HCP Zr and Zr-0.5%Nb alloy has been performed. Evaluations of the SIA diffusion coefficient and self-diffusion coefficient by the interstitial mechanism have been obtained. The significant qualitative and quantitative difference is demonstrated for the obtained values in the cases of pure Zr and alloy. The SIA mobility in the alloy is significantly lower than for pure Zr due to high positive binding energy between SIA-configuration and Nb atom. Thus, SIA has been fast trapped by the Nb atoms. In pure Zr, the SIA migrates mainly in a basal plane, whereas in alloy it migrates mainly along the [0001] axis perpendicular to the basal plane. In the alloy, self-diffusion of the Zr atoms is almost isotropic and significantly lower than in pure Zr. The Nb atoms, in contrast, have demonstrated high mobility at the temperature above 600 K. They migrate more preferably along [0001] axis than in (0001) plane. It should be noted that this property agrees qualitatively with the experimental results on diffusivity of interstitial Cu, Cr, Fe and Ni in α-Zr [24].
The obtained results are believed to be useful for investigation of microstructure evolution of ZrNb alloys under operation in nuclear plants.