Application of DFT to nuclear fuel characterization

Fundamental properties: UO₂

In addition to its technical importance as the standard nuclear fuel material, UO₂ has attracted interest as a prototypical actinide compound to investigate the intriguing behavior of f-electrons in an open 5f shell. UO₂ crystallizes in the cubic fluorite CaF₂ structure and develops antiferromagnetic order below 30 K.Reference Laskowski, Madsen, Blaha and Schwarz^3, Reference Dorado, Freyss and Martin^19, Reference Dudarev, Nguyen Manh and Sutton²⁵ Various DFT studies have been conducted to investigate basic properties of UO_2, including its structural, chemical, thermodynamic, electronic, and magnetic properties.Reference Dudarev, Botton, Savrasov, Humphreys and Sutton^1–Reference Dorado, Amadon, Freyss and Bertolus⁸ In particular, advanced DFT-based approximations beyond the LSDA and GGA, such as LSDA/GGA+U, Reference Dudarev, Nguyen Manh and Sutton^25, Reference Baer and Schoenes²⁶ hybrid functionals,Reference Prodan, Scuseria and Martin^5, Reference Kudin, Scuseria and Martin^27, Reference Atta-Fynn and Ray²⁸ self-interaction correction (SIC),Reference Petit, Svane, Szotek and Temmerman^29, Reference Petit, Svane, Szotek, Temmerman and Stocks³⁰ and DFT+DMFT (dynamical mean field theory) approaches,Reference Yin and Savrasov²² have been developed to correctly capture the strong correlations of 5f electrons in UO₂ that are responsible for its insulating behavior with a bandgap energy of about 2 eV. Figure 1 shows calculated m-decomposed partial f-electron density of states (PDOS) of UO₂ obtained from GGA+U calculations.Reference Yun, Kim, Lim and Park⁶ The PDOS explains the role of the 5f electrons in the U–O bonding in conjunction with antiferromagnetic ordering in UO₂.

Figure 1.

(a) Spin up and (b) spin down m-orbital decomposed partial electron density of states (PDOS) of the uranium 5f electrons in UO₂. (c) Sketch of the seven 5f orbitals.Reference Petrucci, Harwood and Herring³¹ The first-principles calculations provide detailed information on the chemical bonding between uranium and oxygen.

The diagram at the bottom shows the seven 5f, m = –3 to 3 orbitals.Reference Petrucci, Harwood and Herring³¹ Among these orbitals, four spin-up orbitals (m = –3,–1,1,3) are partially occupied in UO₂ with lobes directed along either the x- or y-axis in the coordinate system and without a nodal plane perpendicular to the z-axis. The four 5f PDOS peaks are coincident with peaks in the O p-DOS (i.e., these four 5f orbitals hybridize with O 2p _x and 2p _y orbitals).Reference Yun, Kim, Lim and Park⁶ The electronic structure of UO₂ is still a challenging topic, and comprehensive DFT studies are being conducted that consider spin-orbit and crystal field effects in order to investigate the crystal field levels and non-collinear magnetic structure of UO₂.Reference Zhou and Ozoliņš³²

Materials properties

Defects

Defects are central to an understanding of atomic transport in nuclear fuels. In particular, it is known that point defects form a major diffusion channel for fission gas atoms in the fuel matrix. To model defect-containing materials, a supercell approach is required. A 2 × 2 × 2 supercell containing 96 atoms has been predominantly used in recent DFT calculations.Reference Yun, Kim, Kim and Park^12–Reference Dorado, Freyss and Martin¹⁹Figure 2a shows the oxygen sublattice in a 2 × 2 × 2 supercell of cubic UO₂, ThO₂, and PuO₂. The oxygen atoms are located at the corners of each cube, and actinide atoms are located at the center of alternative grey cubes.

Figure 2.

(a) Oxygen sublattice of a 2 × 2 × 2 supercell of actinide dioxides; (b) the migration of an oxygen vacancy V_O between two adjacent lattice sites; (c) the migration of an actinide vacancy V_M (M = U or Th); (d) the migration of an actinide vacancy V_M’ due to the effective movement of a joint vacancy, V_MO. OIS, octahedral interstitial site; a ₀, lattice parameter.

Table I reports formation energies of major intrinsic defects in UO₂. Here an oxygen Frenkel pair consists of an oxygen interstitial and vacancy, and a Schottky trio defect consists of one uranium vacancy and two oxygen vacancies formed in the UO₂ lattice. The five computed energies were obtained using DFT-based methods and a 2 × 2 × 2 supercell. Most of the calculated values are in good agreement with the experimental data (last column), which illustrates that parameter-free DFT calculations are reliable and accurate within an acceptable error range. Another important achievement of DFT studies is the improved understanding of atomistic diffusion. Table II gives DFT-computed migration energies of single vacancies in UO₂ and ThO₂. The migration energy was calculated from the energy difference between an initial configuration and a saddle point when a M (U or Th) or O vacancy, denoted V_M or V_O, respectively, moves from its initial lattice site to the nearest adjacent site, as shown in Figure 2b. The migration energy of V_U is calculated to be almost 50% higher than the experimental value, whereas that of V_O agrees very well. However, the overestimated result for V_U can be explained by considering the effective movement of a combined U and O vacancy, V_UO, as depicted in Figure 2d. The migration energy of V_U is reduced by about 1 eV due to the existence of an adjacent V_O. Hence, the calculations predict a detailed atomistic process in which an oxygen vacancy plays an important role to enhance the mobility of a uranium vacancy in UO₂. The migration energies of single vacancies in ThO₂ (see Table II) are predicted to be significantly higher than those in UO₂.Reference Yun, Oppeneer, Kim and Park¹⁵ The assisting influence of a V_O was not observed in the movement of a Th vacancy. DFT studies thus predict that point defects are less easily created as well as less mobile in ThO₂ than in UO₂.

Table I.

Calculated and experimental formation energies in eV of intrinsic defects in UO₂.

The calculated resultsReference Iwasawa, Chen, Kaneta, Ohnuma, Geng and Kinoshita^16–Reference Dorado, Freyss and Martin¹⁹ were obtained using the generalized gradient approximation with added Coulombic correlation (GGA+U) approach, except for the resultsReference Yun, Oppeneer, Kim and Park¹⁵ obtained using GGA.

Table II.

ComputedReference Yun, Kim, Kim and Park^12, Reference Yun, Oppeneer, Kim and Park¹⁵ and experimentalReference Matzke³³ migration energies E _m in eV.

V_O and V_M indicate oxygen and actinide atom vacancies (M = Th or U), respectively, and E _m(V_M)’ indicates the effective migration energy of an M-vacancy within the joint V_M and V_O movement.

The diffusion of fission gases in the fuel matrix is closely connected to lattice defects, as these defects are known to act as a major diffusion channel of fission gas atoms. Therefore, based on the defect energetics, it is expected that fission gas atoms are less mobile in ThO₂ than in UO₂. Further calculations were conducted to understand diffusion behavior of fission gas atoms in ThO₂ and UO₂, and it was reportedReference Yun, Oppeneer, Kim and Park¹⁵ that Xe, which is one of the highest fractional released fission gases in UO₂, is more likely to diffuse in UO₂ than in ThO₂. Figure 3a and 3b show the electron charge distribution of UO₂ and ThO₂, respectively, and clearly illustrate the different behavior of Xe in the two oxides. In Figure 3a, Xe with index (1) is initially at V_U, and subsequently is seen to move during atomic relaxation when a V_O (2) is created at the nearest lattice site. In contrast, in ThO₂, Xe remains at the initial position of V_Th during atomic relaxation, despite the creation of a new V_O nearby, as shown in Figure 3b. A detailed schematic diagram depicting the Xe behavior in the two dioxides is plotted on the right side in Figure 3.

Figure 3.

Electron charge distribution of Xe-V_MO complex obtained through atomic relaxation in (a) UO₂ and (b) ThO₂. The magnitude of the electron density is indicated by a color grading, where red is used for the highest electron density and blue for the lowest one. Xe moves from the initial position at V_U to the center of V_UO when a V_O is created at the nearest lattice site, while in (b), it remains at the initial position of V_Th, despite the creation of a new V_O nearby. These processes are schematically explained in the right-hand panel.

In UO₂, the barrier-less movement of Xe during atomic relaxation is due to the strain energy of Xe caused by its relatively large atomic radius (Xe = 2.15 Å, U⁴⁺ = 1.01 Å, Th⁴⁺ = 1.07 Å).Reference Matzke³³ However, the strain-energy effect of Xe seems not enough to cause Xe movement in ThO₂. This is mainly associated with the strong resistance of Th against variation of its valency in ThO₂. Th in ThO₂ is most stable with 4+ valency, and electron redistribution is energetically unfavorable, which affects the formation and mobility of defect atoms and the diffusion of Xe atoms. These DFT results are supported by experiments showing that the Xe release is lower in UO₂ mixed with ThO₂, (U,Th)O₂Reference Kim, Park, Yun, Kim, Ryu, Song, Choo and Hong³⁴ than in UO₂. It is very important to understand the diffusion mechanism of fission gases in nuclear fuels, because fission gas is one of the parameters that severely affects fuel safety conditions, such as fuel swelling, thermal conductivity, and fuel rod pressure after it is released into the gap between fuel pellet and cladding. Controlling fission gas release is therefore a crucial factor to increase nuclear fuel burnup, which is the degree to which the fuel material is consumed in the nuclear reactor so that the fuel efficiency is improved. First-principles studies investigating the trap sites and diffusion mechanism of fission gases reported that Xe atoms prefer to accommodate the charge neutral UO₂ trivacancy among existing interstitial sites and various vacancies due to its charge neutrality as well as its large size, which are suitable to accommodate the large atomic radius of Xe.Reference Matzke^33, Reference Matzke³⁵ The trivacancy was moreover computed to play a central role in the Xe diffusion in UO₂ through the successive process of dissociation and association of the Xe-trivacancy complex,Reference Yun, Kim, Kim and Park¹² as shown in Figure 4.

Figure 4.

(a) Structure of a Xe-trivacancy complex, (b) dissociation of V_O from a Xe-trivacancy with binding energy of about 1.4 eV, (c) combination of a new V_O with the migration energy of 0.63 eV, (d) migration of V_U with the help of the nearest V_O with a migration energy of 2.19 eV, (e) formation of a new Xe-trivacancy complex by the spontaneous movement of Xe due to strain energy in the UO₂ lattice.Reference Anderson, Lezama, Uberuaga, Deo and Conradson⁹

It should be mentioned that accurate atomic transport modeling is a challenge in DFT research, as it is necessary to take temperature effects into account to achieve close energy values. Supercell size is also a sensitive parameter that affects energy calculations required for atomic transport; this aspect will be discussed in more detail later.

Thermodynamic behavior

The heat transfer from the fuel to the coolant is the central process to generate electricity from nuclear power. Knowledge of the fuel’s thermodynamic properties is therefore essential to model the fuel’s behavior under operation conditions at high temperatures. First-principles methods are increasingly being used to determine the thermodynamic properties of nuclear fuels. Recently, a DFT+DMFT study of the lattice dynamical properties of UO₂ and PuO₂ (i.e., their phonon spectra) was reported.Reference Yin and Savrasov²² The calculations suggested that an improved thermal conductivity of the fuel could be achieved through dedicated materials design, building on the observation that dispersive longitudinal optic phonon modes do not contribute to the heat transfer.Reference Yin and Savrasov²² Thus far, however, only a few DFT-based calculations of the phonon spectrum of actinide materials have been performed, employing either the linear response theoryReference Yin and Savrasov^22, Reference Dai, Savrasov, Kotliar, Migliori, Ledbetter and Abrahams³⁶ or the direct approach.Reference Parlinski^37, Reference Piekarz, Parlinski, Jochym, Oles, Sanchez and Rebizant³⁸ In the former approach, a perturbation of the atomic positions with the periodicity of the original lattice is treated with standard perturbation theory, and first-order corrections are computed. In the direct approach, or so-called supercell method, a finite displacement of several atoms in a supercell is employed to compute atomic force constants and the phonon spectrum. In this direct approach, it is important to achieve convergence of phonon frequencies with respect to supercell size. Generation of such data would be a valuable step for characterizing the thermal properties of nuclear fuels. An alternative first-principles approach is to use DFT to compute interatomic potentials, which are subsequently used in molecular dynamics (MD) simulations. This approach has already been successfully applied to describe the thermal expansion of dioxide fuelsReference Yamada, Kurosaki, Uno and Yamanaka³⁹ and the next-generation fuel uranium nitride (UN).Reference Chen, Wang, Lai, Li, Ao and Long⁴⁰

DFT calculations on nuclear fuel: Computational costs

The computationally heaviest step in DFT calculations for nuclear fuel materials is to perform the self-consistent field (SCF) total energy minimization (i.e., the iterative calculation of single-particle wave functions and the effective potential that depends on these wave functions until the starting and evaluated potential match). The costs of the SCF calculation depend mainly on the number of atoms in the chosen simulation cell. To obtain the electronic structure, cohesive properties, and magnetic properties, a calculation for a single unit cell suffices, and the required computation time is usually less than a few hours for the convergence of the SCF total energy with a single CPU. Table III shows the imposed parameter values and the computation time to calculate the electronic structure of UO₂. A tetragonal unit cell is used to model antiferromagnetically ordered UO₂, which consists of two U atoms with antiparallel magnetic moments and four O atoms. Once the self-consistent electronic charge density and wave functions are computed, other elements of the electronic structure, such as the electronic density of states and energy band structure, can be computed quickly.

Table III.

Computational details for self-consistent field (SCF) calculations of antiferromagnetic UO₂ using the VASP (Vienna ab initio simulation package) code.⁴⁵

The VASP code is an ab initio quantum mechanical calculation code using pseudopotentials and a plane-wave basis set. First column: The total energy calculation performed with the projector augmented wave (PAW) method and with spin-orbit coupling (SOC) and an additional Coulomb correlation U included. Second column: The formation energy of a uranium vacancy V_U using the PAW-GGA (generalized gradient approximation) method. The calculation includes relaxation of all atomic positions in the supercell until all atomic forces are minimized.

The use of a large simulation cell is essential for defect-related calculations, such as formation and migration energy of defects. The computation time required increases approximately as ~N ³, with N being the number of atoms in the simulation cell. Moreover, to simulate a system with defects, all atomic positions in the system have to be relaxed to obtain the overall minimal energy configuration. This implies that all atomic positions must be optimized until the forces acting on each of the atoms in the supercell vanish. This is usually an extremely time-consuming task. Table III summarizes the computational details for the calculation of the formation energy of a uranium vacancy V_U in UO₂ using the 2 × 2 × 2 supercell. The number of used CPUs was determined by carrying out an optimization test for the parallel computing platform. From Table III, it can be seen that the complete SCF total energy convergence with structural optimization requires a considerable computation time of 57 hours. Hence, the computational costs to simulate defect diffusion in nuclear fuels are massive. The reduction of computational costs is therefore an important key to increase the impact of first-principles techniques for nuclear fuel modeling. Nonetheless, in spite of the considerable computational efforts required, the costs of performing experimental studies on nuclear fuel materials are unequivocally higher.

Limitations and improvements to DFT modeling

Several challenges are currently encountered in analyses of nuclear fuels using DFT. The first area that requires additional research is extensions and improvements of the DFT+U methodology, since 5f states in actinide oxides are strongly correlated. Here, the U indicates the effective interaction parameter that is added to correct the on-site Coulomb repulsion interaction of the 5f electrons. The DFT+U method has become an important approximation in the study of nuclear fuel materials. It is, from a numerical perspective, a relatively simple extension beyond the standard GGA and can be applied in parallel computations to treat supercells with a large number of atoms. However, the DFT+U functional has the weakness that, depending on the initial density matrix in the self-consistent cycle, one can converge to a metastable state instead of the ground state.Reference Amadon, Jollet and Torrent⁷ A solution to this problem could be the recently proposed density-matrix controlling scheme,Reference Dorado, Amadon, Freyss and Bertolus⁸ which was successfully applied to reach the ground state of UO₂. The density-matrix controlling scheme is a definite improvement to correctly determine the ground state of actinide oxides within DFT+U calculations and could have a significant and positive impact on fuel materials research. However, it has not been applied yet for lattice defect calculations in a supercell. In addition to improvements of the DFT+U functional, it is also necessary that other approximations to the Hamiltonian, such as hybrid functionals,Reference Kudin, Scuseria and Martin^27, Reference Atta-Fynn and Ray²⁸ SIC,Reference Petit, Svane, Szotek and Temmerman^29, Reference Petit, Svane, Szotek, Temmerman and Stocks³⁰ and DFT+DMFT, Reference Yin and Savrasov²² be tested and developed further to calculate nuclear fuel materials, in combination with the development of parallelized versions of these approaches appropriate to the analyses associated with supercell domains.

A second demand on DFT-based techniques is to provide a better characterization of temperature-dependent properties in nuclear fuels. For many fuel materials properties, it is essential to take temperature-dependent effects into account. For example, fission gas diffusion occurs dominantly at temperatures over 1000°C.Reference Matzke³³ Nonetheless, standard DFT calculations can provide acceptable results for fission gas behavior.Reference Freyss, Vergnet and Petit^11, Reference Yun, Kim, Kim and Park^12, Reference Yun, Eriksson and Oppeneer^14, Reference Yun, Oppeneer, Kim and Park^15, Reference Dorado, Freyss and Martin¹⁹ For instance, the migration energy of Xe between two interstitial sites was calculated as 4.69 eV,Reference Yun, Eriksson and Oppeneer¹⁴ and this agrees with experimental data of 3.90 eV within 20%.Reference Matzke³³ However, when expansion of the lattice was accounted for, it was found that the migration energy was reduced to 3.78 eV,Reference Yun, Eriksson and Oppeneer¹⁴ which calls into question the validity of the nominal result. MerminReference Mermin⁴¹ generalized the Hohenberg-Kohn theorem (that a many-electron problem can be exactly mapped to a single-electron problem) to finite temperatures. Temperature-dependent DFT has been successfully applied to describe materials behavior over a wide range in temperature, up to almost 1000ºC.Reference Freyss, Vergnet and Petit^11, Reference Yun, Kim, Kim and Park^12, Reference Yun, Eriksson and Oppeneer^14, Reference Yun, Oppeneer, Kim and Park^15, Reference Dorado, Freyss and Martin¹⁹ Using temperature-dependent DFT in combination with MD simulations, Mattsson et al.Reference Mattsson, Sandberg, Armiento and Mattsson⁴² computed the anomalous temperature evolution of the vacancy formation energy in Mo up to about 2700ºC, and Root et al.Reference Root, Magyar, Carpenter, Hanson and Mattsson⁴³ studied xenon’s behavior at extreme conditions of high pressure and high temperature. Although T-dependent (temperature dependent) DFT can thus be employed to high-temperature dynamic and thermodynamic properties, it has not yet been applied to nuclear materials. Standard DFT-based MD simulations have already been performed.Reference Mattsson, Sandberg, Armiento and Mattsson^42–Reference Govers, Lemehov, Hou and Verwerft⁴⁴ These are, for example, well-suited to investigate the important role of fission gas diffusionReference Freyss, Vergnet and Petit^11, Reference Yun, Kim, Kim and Park^12, Reference Yun, Eriksson and Oppeneer^14, Reference Yun, Oppeneer, Kim and Park^15, Reference Dorado, Freyss and Martin¹⁹ in the nuclear fuel. DFT-based MD simulationsReference Govers, Lemehov, Hou and Verwerft⁴⁴ of the temperature-dependent diffusion coefficient and pathway were reported for UO₂. Furthermore, the defect jump rate and absolute concentration can be calculated as a function of temperature using the DFT-based MD.Reference Mattsson, Sandberg, Armiento and Mattsson⁴² MD is implemented in several DFT simulation packages.^{45, 46}

Current limitations on supercell size are a third area deserving of research attention. The scale of DFT-based modeling is, in general, limited to cells of the order of a few hundred atoms, and this limits the types of defect structures and interactions that may be considered. For example, to calculate the migration energy of fission gas atoms in a dilute concentration, a reasonably large supercell, consisting approximately of N times the unit cell for a defect concentration of 1/N, is needed to model a system containing several defects well separated from each other; otherwise, the defect-defect interaction between adjacent cells will be too large to obtain accurate energy values. Suitable supercell sizes can be determined from the geometry optimization of the supercell containing defects by carrying out total energy convergence tests, including atomic relaxation as a function of supercell size (i.e., the lattice coordinates of all the atoms in the supercell are modified during the total energy calculation until they converge to their equilibrium positions that are defined by the criterion that the forces acting on each atom are smaller than a lower value, usually 0.01 eV/Å). An efficient approximation has been developed recently that allows larger system sizes to be considered, allowing for an analysis of dislocation and dislocation-solute physics.Reference Trinkle and Woodward⁴⁷ A 540-atom supercell was used with periodic boundary conditions to the dislocation line direction. To improve computational efficiency for such a huge number of atoms, the cell was divided into three regions, and forces in the separate regions were relaxed, subject to a consideration of the long-range elastic response of the material.

Along with the development of algorithms focused on defects and actinides, it is also important to improve computational performance capacity through the development of efficient parallel computing algorithms. While not strictly focused on nuclear fuels, such generic improvements to DFT methodology may prove to have the largest impact of all on the ability to consider large numbers of atoms.

Summary and future role of DFT techniques in nuclear fuel research

First-principles modeling is now widely accepted as a tool not only to develop advanced fuels with improved efficiency in current reactors but also to design nuclear fuel materials of the future.Reference Grimes, Konings and Edwards⁴⁸ As already noted, for instance, density functional theory (DFT) studies predict that ThO₂ possesses qualities superior to UO₂ with regard to chemical stability and defect dynamics. Nitrides such as (U,Pu,MA)N and (Pu,MA,Zr)N, where MA is a minor actinide, are now being investigated using DFT as promising advanced nuclear fuel because of their high thermal conductivity and high melting point.Reference Muta, Kurosaki, Uno and Yamanaka⁴⁹ In a recent DFT study, point defects in uranium mononitride were investigated, calculating the defect formation energy as well as defect-induced electronic density redistribution.Reference Kotomin, Grimes, Mastrikov and Ashley^50, Reference Lu, Wang, Li, Shi and Zhang⁵¹ Uranium monocarbide and ternary alloys, including U-Pu-C, are now considered as possible candidates for next-generation reactor fuels because of the combination of a high melting temperature and high actinide density,Reference Vigier, Den Auwer, Fillaux, Maslennikov, Noël, Roques, Shuh, Simoni, Tyliszczak and Moisy^20, Reference Freyss²¹ and first-principles calculations have been performed to investigate the thermodynamic properties and phase diagram of these alloys.Reference Freyss²¹ Advanced DFT modeling can already provide valuable knowledge regarding the existence of potential candidate materials before experimental data, such as structure and basic properties.

During the last few decades, first-principles calculations performed within the framework of DFT have had a tremendous and positive impact on nuclear fuels research. It has become evident that DFT can significantly contribute to the nuclear energy industry by providing both insights and quantitative property data associated with fuel performance. Parameter-free DFT calculations have not only improved the reliability and accuracy of computed materials properties, but they also contributed to a better understanding of the microstructural character of nuclear fuels. However, the current DFT techniques need further improvement, with the primary need being associated with the ability to scale up to larger numbers of atoms and to develop finite temperature predictions of material behavior.