Thorium dioxide is used industrially in high temperature applications, but more insight is needed into the behavior of the material as part of a mixed-oxide (MOX) nuclear fuel, incorporating uranium. We have developed a new interatomic potential model including polarizability via a shell model, and commensurate with a prominent existing UO2 potential, to conduct configurational analyses and to investigate the thermophysical properties of uranium-doped ThO2. Using the GULP and Site Occupancy Disorder (SOD) computational codes, we have analyzed the distribution of low concentrations of uranium in the bulk material, where we have not observed the formation of uranium clusters or the dominance of a single preferred configuration. We have calculated thermophysical properties of pure thorium dioxide and Th(1−x)UxO2 which generated values in very good agreement with experimental data.

Although the majority of nuclear power reactors operating today use uranium and/or plutonium-based fuels such as low-enriched uranium (LEU) or mixed oxide fuels (MOX), alternative nuclear fuels are increasingly considered as well.1 The most promising alternatives are fuels based on thorium, which is generally estimated to be between three and four times more abundant than uranium in the Earth’s crust, with large stores in India, the United States, and other areas around the world. Thorium compounds can be used to fuel nuclear reactors and they thus present an attractive alternative to the standard uranium and plutonium-MOX fuels.2 

Whereas thorium is a fertile material, it is not fissile and therefore a sustained conversion to the fissile 233U isotope can only occur in the presence of a neutron source. Although plutonium may serve as the source, there are several advantages to blending thoria with urania that have been confirmed experimentally, including changes in the decay heat, melting point, and thermal conductivity of the material, in addition to a decrease in the release rate of fission gases despite an actual increase in the production of fission gases.3 

Pure thorium dioxide (thoria) is a ceramic material with a melting point of 3651 K, the highest of any known oxide.4 Thoria has the fluorite structure, which is shared with the dioxides of, for example, cerium, plutonium, and uranium. Unlike uranium and plutonium oxides, ThO2 also represents the highest oxidation state of the material, which makes it exceptionally stable in the presence of oxygen or oxygenated water. This corrosion resistance has positive implications for the use of thoria in nuclear fuels.3 Since a neutron source is required for its application as a nuclear fuel, and as ThO2 and UO2 are isostructural, there is significant interest in thoria-urania solid solutions and both experimental and theoretical investigations have been carried out.5–8 As it is useful to study thoria in the context of uranium-doping, it is advantageous to have a reliable computational method to model these materials due to the safety hazards and high costs of conducting experiments on radioactive materials.

We are not aware of any robust theoretical investigations into the different possible configurations of uranium substitution in thoria, which might elucidate the distribution of the uranium atoms serving as a neutron source in a thorium dioxide nuclear fuel, and the thermodynamic properties of these configurations. It is important to understand how uranium behaves when inserted into thoria in order to understand the behavior of potential large-scale MOX fuel rods. Since thorium has no further oxidation states beyond Th4+, UO2 is more susceptible to oxidation than ThO2, and so any uranium clustering at surfaces or grain boundaries could lead to oxidation and thereby compromise the corrosion-resistance of the material. In this paper, we present our computational investigation into the numerous independent configurations obtainable at a variety of uranium concentrations in ThO2, using a new interatomic potential to describe the interaction between thorium and oxygen in this system.

In experimental test reactors using thorium fuels, typically only small concentrations, in the region of 5-15%, of uranium or plutonium are used. The role of U or Pu in these fuel blends is not as a fuel itself but rather as a neutron source for the production of 233U by 232Th. To model such small concentrations of uranium dopant in thoria, we need a sufficiently large system, which requires the use of affordable computational methods based on interatomic potentials. Early potentials, such as those used by Benson et al to study surface energies of thoria, were developed prior to the publication of experimental data on the bulk properties of thorium dioxide and therefore relied on estimating those properties from other ionic oxides.9 Later studies either did not publish the potential parameters or did not closely reproduce the available experimental data.10–12 While several good ThO2 potentials were published by Behera, et al in 2012, we have developed a thoria potential to be compatible with a leading UO2 potential13–15 and we have applied this model to the distribution of uranium in relatively large ThO2 systems.

The energies of different configurations of the mixed solid were evaluated using the General Lattice Utility Program (GULP), version 3.4.16,17 GULP implements the Born model of solids, which assumes that the ions in the crystal interact via short-range repulsive interactions, longer-range attractive interactions, and long-range Coulombic interactions.18 The short-range energy contribution due to interactions between each particle and all others is summed within a predetermined cut-off. The contribution of the long-range electrostatic forces is summed using the Ewald sum.19 The electronic polarizability of ions is included via the shell model of Dick and Overhauser in which each of the polarizable ions, in this case the oxygen atoms, are represented by a core and a massless shell, connected by a spring.20 The polarizability of the model ion is then determined by the spring constant and the charges of the core and shell. The advantage of using the shell model is that it better reproduces the elastic and dielectric properties of ionic solids than rigid ion models, though it increases the computational cost of the simulations.21 

The GULP code is widely used in solid state chemistry in simulations of materials using periodic boundary conditions. We have used GULP to derive the potential parameters and compute bulk properties for ThO2, with initial values taken from a suitable existing UO2 potential, and fitting to experimental crystallographic data. Details of the interatomic potentials used are given in section III A, equation (15). In order to make geometry predictions, the lattice energies are minimized with respect to the structural parameters, until the forces acting on the ions are all less than 0.001 eV Å−1. All structures reported are the result of constant pressure energy minimizations, where not only the ionic positions but also the cell parameters are allowed to vary to find the energy minimum. The external pressure was set to zero in all calculations.

The Th1−xUxO2 solid solutions were represented by a symmetry-adapted ensemble of configurations in a supercell, using the methodology implemented in the SOD (Site Occupancy Disorder) program,22,23 which has been employed previously in the simulation of a range of mineral solid solutions. Successful applications of SOD include a study of the thermophysical properties of carbonates, investigations into the mixing thermodynamics of oxides and a simulation of cation distribution and thermodynamics in sulfides.24–31 This program generates the complete configurational space for each composition in a supercell of the structure, before extracting the subspace of symmetrically inequivalent configurations, for which energies and other properties are evaluated. In this case a 2x2x2 thoria supercell (Figure 1) was used as the base simulation cell, and we have worked with uranium concentrations up to x = 0.16. The computational cost of investigating higher U concentrations in a system of this size was considerable and of limited practical concern, as the amount of uranium-doping required in thoria fuels is small. However, to investigate the full range of the solid solution and derive thermodynamic properties that can be compared with experimental values, a full configurational analysis was carried out on a 1x1x2 supercell.

FIG. 1.

a) 1x1x2 supercell, b) and c) 2x2x2 supercell face and side views of ThO2, where the Th atoms are represented by the blue spheres and the O atoms are represented by the red spheres.

FIG. 1.

a) 1x1x2 supercell, b) and c) 2x2x2 supercell face and side views of ThO2, where the Th atoms are represented by the blue spheres and the O atoms are represented by the red spheres.

Close modal

Once the configurational spectrum is obtained, it is possible to derive thermodynamic properties through further application of statistical mechanics. In the simplest formulation, the extent of occurrence of one particular configuration, n, in the disordered solid in configurational equilibrium at a temperature T can be described by a Boltzmann-like probability, Pn:

P n = 1 Z exp E n k B T
(1)

where En is the energy of configuration n and Z is the partition function given by the following equation:

Z = n = 1 N e x p E n k B T
(2)

The energy of the system is now calculated as a sum as follows:

E = n = 1 N P n E n
(3)

Then in the reduced configurational space, a new probability can be determined for an inequivalent configuration, m. This is calculated from the energy of that configuration, Em, and its degeneracy Ωm (the number of times that the configuration is repeated in the complete configurational space):

P m = Ω m Z exp E m R T = 1 Z exp E m ( red ) k b T
(4)

where m = 1, …, M (M is the number of inequivalent configurations), R is the gas constant, and Z is the partition function (Equation (5)), which guarantees that the sum of all the probabilities equals 1 and also gives access to the calculation of configurational free energies and entropies.

Z = m 1 M exp E m RT
(5)

Em(red) is the reduced energy and relates to the energy Em and the reduced, or degeneracy, entropy, Sm, via:

E m ( red ) = E m TS m
(6)
S m = k B T ln Ω m
(7)

From this relationship, it is possible to compare the inequivalent configurations as the degeneracy and the energies are considered. For two inequivalent configurations of the same energy, the one with the highest degeneracy will then have the higher degeneracy entropy and therefore a lower reduced energy and thus higher probability (equation (4)). It is possible in this way to obtain temperature-dependent configurational thermodynamics functions and equilibrium degrees of disorder.32–34 In this work we have focused on the distribution of uranium in a range of concentrations in the thoria lattice and the derivation of thermodynamic properties of the solid solution.

The configurational free energy G of the disordered solid can be obtained directly from the partition function:

G = k T ln Z
(8)

and any average observable, including the enthalpy H of the solution, can be estimated using configurational averaging:

H = m = 1 M P m H m
(9)

Having obtained the enthalpy and free energy of the solid solution, it is useful to evaluate the enthalpy of mixing:

Δ H mix = H Th 1 x U x O 2 1 x H bulk ThO 2 x H bulk UO 2
(10)

which can be compared with experimental calorimetric determinations, and the free energy of mixing:

Δ G mix = G Th 1 x U x O 2 1 x G bulk ThO 2 x G bulk UO 2
(11)

which can also be found experimentally from equilibrium composition measurements.

Technically important properties like the thermal linear expansion coefficient and heat capacity can be calculated:

α = 1 a 293 a T
(12)

where a is the lattice parameter of the cell at a given temperature, T, and a293 is the lattice parameter at 293 K. The values of a at seven temperatures were calculated directly in GULP followed by fitting of the values to a 3rd order polynomial (Table I) to enable the extrapolation to higher values. As noted, due to the computational cost of calculating the configurational thermodynamics for the full ThO2-UO2 solid solution using the large 2x2x2 supercell, a 1x1x2 cell was used instead for these calculations.

TABLE I.

Polynomial fit for ThO2 and UO2aT (Å)=c0 + c1T + c2T2 + c3T3.

Material c0 c1×105 c2×108 c3×1012 α293
ThO2  5.60793  2.79036  1.46253  -3.43684  5.608 
UO2  5.45337  2.54896  1.45446  -3.76733  5.462 
Th0.87U0.13O2  5.58811  2.93563  1.23247  -2.01114  5.598 
Th0.75U0.25O2  5.56882  2.71268  1.70802  -4.78603  5.578 
Material c0 c1×105 c2×108 c3×1012 α293
ThO2  5.60793  2.79036  1.46253  -3.43684  5.608 
UO2  5.45337  2.54896  1.45446  -3.76733  5.462 
Th0.87U0.13O2  5.58811  2.93563  1.23247  -2.01114  5.598 
Th0.75U0.25O2  5.56882  2.71268  1.70802  -4.78603  5.578 

GULP will also calculate certain thermodynamic properties of a system when a phonon calculation is performed. We have taken advantage of this fact to extract the specific heat capacities at constant volume of our system while computing phonons, where the constant volume heat capacity, Cv, is calculated from the vibrational partition function, Zvib:

C v = RT 2 ln Z vib T + T 2 ln Z vib T 2
(13)
Z vib = k points w k all modes ( 1 exp ( h v k T ) 1
(14)

We have used the thoria structure35,36 in our derivation of a Th-O Buckingham potential, VBuck(ri,j):

V Buck r ij = A ij e r ij ρ ij C ij r ij 6
(15)

where rij is the distance between ions i and j, and A, ρ, C, are variable parameters representing the pair-wise repulsion.

Thoria has a cubic unit cell with the space group Fm 3 ¯ m (N° 225), with four cation sites in the unit cell. The least-squares fitting against structure and other properties was done in GULP using experimental data from Idiri and Macedo and co-workers.35,36 While there has been extensive debate in the literature as to the correct value of the bulk modulus of thorium dioxide, in this work we have used the experimental value of 198 ± 2 GPa determined by Idiri et al using energy dispersive x-ray diffraction (EDXRD)35 Staun Olsen et al also studied the bulk modulus of thoria with high-pressure x-ray diffraction (XRD) and synchrotron radiation, finding a similar value of 195 ± 2 GPa.37 The oxygen-oxygen interaction in ThO2 is described by the potential published by Catlow, whereas initial values of A and ρ for the Th-O interaction were taken from the Catlow U-O Buckingham potential.15 Although U4+ (1.00 Å) has a smaller ionic radius than Th4+ (1.05 Å), due to the similarities between urania and thoria, these values provided a good starting point for the derivation of the Th-O potential, which is introduced in Table II.38 

TABLE II.

Parameters for the Buckingham potentials employed in this study.

Charge (e)
Species core shell Spring constant (eV Å−2)
Th  +4.000 
+4.000 
+2.400  −4.400  292.98 
Buckingham Potential  A (eV)  ρ (Å)  C (eV Å6
Th—Oshell  1281.775  0.3910  0.0 
U—Oshell  1217.800  0.3871  0.0 
Oshell—Oshell  22764.300  0.1490  112.2 
Charge (e)
Species core shell Spring constant (eV Å−2)
Th  +4.000 
+4.000 
+2.400  −4.400  292.98 
Buckingham Potential  A (eV)  ρ (Å)  C (eV Å6
Th—Oshell  1281.775  0.3910  0.0 
U—Oshell  1217.800  0.3871  0.0 
Oshell—Oshell  22764.300  0.1490  112.2 

Table III compares structural properties of ThO2 calculated using the newly developed potential with available experimental data. There is of course excellent agreement with the lattice parameter, which was used in the derivation, and good agreement with the elastic properties. By using a shell model potential, we are also able to correctly model the Cauchy violation observed in fluorite structures (C12≠C44). This potential will now be used to study the Th1−xUxO2 solid solutions.

TABLE III.

Bulk properties of thoria calculated with the new potential compared with the available experimental data and the percentage error between the two values.

ThO2 Experimental Value35,36 Calculated Value % Error in Calculated Value
Lattice Parameter (Å)  5.601  5.6001  0.00178 
Bulk Modulus (GPa)  198  206.4  4.24 
Shear Modulus (GPa)  98  94.4  −4.08 
Elastic Constant C11 (GPa)  367  389.5  6.14 
Elastic Constant C12 (GPa)  106  114.9  8.36 
Elastic Constant C44 (GPa)  79.7  65.8  −16.72 
ThO2 Experimental Value35,36 Calculated Value % Error in Calculated Value
Lattice Parameter (Å)  5.601  5.6001  0.00178 
Bulk Modulus (GPa)  198  206.4  4.24 
Shear Modulus (GPa)  98  94.4  −4.08 
Elastic Constant C11 (GPa)  367  389.5  6.14 
Elastic Constant C12 (GPa)  106  114.9  8.36 
Elastic Constant C44 (GPa)  79.7  65.8  −16.72 

It is important to have a thorough understanding of the temperature-dependent properties of any nuclear fuel in order to understand and predict its behavior. For example, the thermal expansion must be understood before the fuel rods and cladding are designed.39 If the fuel-cladding gap is either too small or too large, failure of the fuel cladding can occur. If the gap is too small, as the fuel expands, pressure is exerted on the cladding and deformation of the fuel can occur. If the gap is too large, the fuel may vibrate excessively and so damage the cladding.40 

We have calculated the value of the thermal linear expansion coefficient, α, at 6.4 × 10−6 K−1 at 300 K. This value is in line with the experimental values reported by Marples (7.3 × 10−6 K−1), Kempter and Elliot (8.2 × 10−6 K−1), and, at 293K, Yamashita et al (8.43 × 10−6 K−1)41–43 While our thermal expansion coefficient does underestimate the values found experimentally, this could be explained by the small difference in the initial lattice parameters used in the experimental systems and in our system. We have used a value for the lattice parameter, a, of ThO2 of 5.6 Å in our simulations, whereas Kempter, using XRD measurements, found an initial a value of 5.597 Å and Yamashita one of 5.58 Å, also using XRD. As ThO2 does not undergo much expansion with temperature, it is possible that by using a larger initial lattice parameter, our calculations failed to account for a small amount of expansion and so underestimated the coefficient of thermal expansion. It should also be noted that Kempter reported impurities in his original samples of ThO2 and these could have had a non-negligible effect on the experimental values reported.

A computational study based on Density Functional Theory (DFT) by Lu et al reports values at higher temperatures, i.e. α = 33–37 × 10−6 K−1, an order of magnitude above the values determined experimentally and that are estimated by our simulation.44 Martin et al have used a rigid ion potential in molecular dynamics simulations to calculate an α of 9.2 × 10−6 K−1 at 1500 K, which is in the same range as our value and the experimental values but this method does not account for the polarizability of the ions or the contribution of phonons in their simulation7 Similarly, using a Born-Mayer-Huggins potential with a partially ionic model, Ma et al have calculated an α at 1500 K of 11.9 × 10−6 K−1.45 These calculated values at 1500 K are lower than the experimental values reported at 1200 K of Yamashita et al and that at 1273 K of Kempter and Elliot. These findings indicate that interatomic potentials used to determine the coefficient of linear expansion tend to underestimate the values but are of the same order of magnitude as experiment and a distinct improvement on the values obtained from DFT for this material.

The equilibrium geometries and energies of the Th(1−x)UxO configurations were calculated both with and without considering vibrational effects. The SOD code identified a number of independent configurations for the uranium-doped thoria systems under investigation (Table IV). The number of inequivalent configurations increases significantly with each additional substitution until a U mole fraction of x = 0.5 is reached in the 1x1x2 supercell.

TABLE IV.

The total number of geometry configurations and the number of inequivalent configurations found with SOD using 1x1x2 and 2x2x2 supercells.

Supercell Th(1−x)UxO2 Mole fraction of Uranium, x Total Number of Configurations Total number of Inequivalent Configurations
  Th7U1O16  0.125 
  Th6U2O16  0.25  28 
  Th5U3O16  0.375  56 
1x1x2  Th4U4O16  0.5  70 
  Th3U5O16  0.625  56 
  Th2U6O16  0.75  28 
  Th1U7O16  0.875 
  Th31U1O64  0.031  32 
  Th30U2O64  0.063  496 
  Th29U3O64  0.094  4960  14 
2x2x2  Th28U4O64  0.125  35960  71 
  Th27U5O64  0.156  201376  223 
  Th26U6O64  0.188  906192  874 
  Th25U7O64  0.219  3365856  2706 
Supercell Th(1−x)UxO2 Mole fraction of Uranium, x Total Number of Configurations Total number of Inequivalent Configurations
  Th7U1O16  0.125 
  Th6U2O16  0.25  28 
  Th5U3O16  0.375  56 
1x1x2  Th4U4O16  0.5  70 
  Th3U5O16  0.625  56 
  Th2U6O16  0.75  28 
  Th1U7O16  0.875 
  Th31U1O64  0.031  32 
  Th30U2O64  0.063  496 
  Th29U3O64  0.094  4960  14 
2x2x2  Th28U4O64  0.125  35960  71 
  Th27U5O64  0.156  201376  223 
  Th26U6O64  0.188  906192  874 
  Th25U7O64  0.219  3365856  2706 

The variation of the lattice parameter, a, with the mole fraction, x, in the ThO2–UO2 solid solutions modeled follows Vegard’s law (Figure 3), in good agreement with the experimental measurements of Hubert et al and Banerjee et al5,46 We can further compare our models with these data by examining the cation-anion and cation-cation distances Th–O, U–O, Th–Th, and U–U (Figure 4 and Figure 5).

FIG. 3.

Comparison of the cell parameter vs. composition with Vegard’s law.

FIG. 3.

Comparison of the cell parameter vs. composition with Vegard’s law.

Close modal
FIG. 4.

Cation-anion distances decrease as the mole fraction, x, of the doped system increases. Black circles and triangles, this work; open triangles and circles, experimental data from Hubert et al.44 

FIG. 4.

Cation-anion distances decrease as the mole fraction, x, of the doped system increases. Black circles and triangles, this work; open triangles and circles, experimental data from Hubert et al.44 

Close modal
FIG. 5.

Cation-cation distances as a function of the mole fraction of uranium, x. Black circles and triangles, this work; open triangles and circles, experimental data from Hubert et al.46 

FIG. 5.

Cation-cation distances as a function of the mole fraction of uranium, x. Black circles and triangles, this work; open triangles and circles, experimental data from Hubert et al.46 

Close modal

Before we consider the distances in the full ThO2–UO2 solid solution, it is useful to examine the effect of a single uranium atom substitution in the thoria supercells. While the 1x1x2 supercells is of practical use in modeling the full solid solution, we observe that our values of the interatomic distances at x = 0.13 are higher than the experimental values at this mole fraction of uranium (Figure 4 and Figure 5). To investigate why this might be, we have investigated the effects of a single uranium substitution in both the 1x1x2 and 2x2x2 supercells (Figure 2). Of course a single substitution yields different mole fractions of uranium in the cells, where x = 0.13 in a 1x1x2 and 0.03 in a 2x2x2 supercell. However, we can see from the interatomic distances in Table V, that small changes are introduced to the lattice even at the lowest mole fraction of uranium studied, x = 0.03. While in this cell, the introduction of one U atom does not significantly affect the Th-Th or Th-O distances, the Th-U and U-O distances are smaller than the Th-Th and Th-O distances being replaced. Moving to the 1x1x2 cell, with one U substitution, all four of the possible cation-cation and anion-anion distances are affected by the presence of a uranium atom. Due to the shape of the 1x1x2 cell, a larger percentage of the oxygen atoms in the lattice are affected by the presence of uranium than in the 2x2x2 cell. These oxygen atoms are moved closer to uranium than they would be in the pure ThO2 lattice, thus affecting the other interatomic distances in the small supercell, resulting in an exaggerated effect when a single substitution is introduced into a small supercell. In a larger supercell with an equivalent mole fraction of uranium, there are many available configurations. We have calculated the interatomic distances in a 2x2x2 supercell for x = 0.13 as an average of all 71 inequivalent configurations (Table V). These distances reproduce more closely the experimental values of Hubert et al compared to the distances in the singly substituted 1x1x2 cell. Therefore, we consider that the overestimation of the interatomic distances at x = 0.13 in Figure 4 and Figure 5 is due to the effect of a single uranium substitution in a small simulation supercell. However, the 1x1x2 cell is able to provide a good representation of the solid solution at higher levels of uranium substitution, when the concentration is more balanced.

FIG. 2.

a) Th7U1O16 b) a Th31U1O64 supercells.

FIG. 2.

a) Th7U1O16 b) a Th31U1O64 supercells.

Close modal
TABLE V.

Interatomic distances (Å) for pure ThO2 and the 1x1x2 and 2x2x2 supercells with one uranium substitution (x = 0.13 and 0.03 respectively) and for the 2x2x2 supercell with 4 uranium substitutions (x = 0.13).

Th-Th Th-U Th-O U-O
ThO2  3.961  —  2.425  — 
Th7 U1 O16  3.978  3.970  2.437  2.422 
Th31 U1 O64  3.960  3.948  2.426  2.389 
Th28 U4 O64  3.952  3.940  2.423  2.386 
Th-Th Th-U Th-O U-O
ThO2  3.961  —  2.425  — 
Th7 U1 O16  3.978  3.970  2.437  2.422 
Th31 U1 O64  3.960  3.948  2.426  2.389 
Th28 U4 O64  3.952  3.940  2.423  2.386 

As noted, to obtain the interatomic distances, we have taken the average values over all the inequivalent configurations. Overall, our Th-Th distances are slightly overestimated compared to those distances found by Hubert et al at low concentrations of U, and slightly underestimated at high U concentrations. The U–U distances are more consistent but are underestimated across the full range of solid solutions. Hubert et al did not publish data on Th-U distances for us to use for comparison, but we have found in our calculations that the Th-U distances are greater than those of U-U but less than Th-Th for each value of x. In absolute terms all of the calculated cation-cation distances are extremely close to experiment and consistent with the fact that U is the smaller cation and that its incorporation into thoria creates stress in the lattice to maintain the four-fold coordination of the oxygen in the fluorite structure. The U–U separation is larger at low values of x. As typically only small amounts of uranium are added to a thorium-uranium MOX fuel, the larger U–U separation at these concentrations indicates that uranium does not cluster within the thoria lattice. The Th-O and U-O distances show similar behavior to the Th-Th distances, as expected from the smaller ionic radius of uranium compared to thorium.

To analyze the Boltzmann probability distributions, it is now useful to consider the larger 2x2x2 supercell. By using the larger cell, we can not only analyze the system at the low concentrations of uranium that are of particular interest in MOX fuels, but we can also more clearly observe any changes in the lattice that may occur on the addition of uranium. In Figure 6, with a mole fraction of uranium of x = 0.06, we show that configurations 4 and 5 together are clearly dominant, with probabilities of 38 and 41%, respectively. The corresponding structures may be seen in Figure 8. As the concentration of uranium and the total number of inequivalent configurations increases for Th29U3O64, there are no longer just one or two dominant configurations as illustrated in Figure 7.

FIG. 6.

Probability distribution of the independent configurations of Th30U2O64 at 500 K.

FIG. 6.

Probability distribution of the independent configurations of Th30U2O64 at 500 K.

Close modal
FIG. 8.

Th30U2O64 where Th is represented in blue, U in green, and O in red. a. Configuration 4 (38% probability); b. Configuration 5 (41% probability); c. Lowest probability structure.

FIG. 8.

Th30U2O64 where Th is represented in blue, U in green, and O in red. a. Configuration 4 (38% probability); b. Configuration 5 (41% probability); c. Lowest probability structure.

Close modal
FIG. 7.

Probability distribution of the independent configurations of Th29U3O64 at 500 K.

FIG. 7.

Probability distribution of the independent configurations of Th29U3O64 at 500 K.

Close modal

To further analyze the relationship between the structures of the inequivalent configurations and the probabilities that a given configuration will occur at equilibrium, it is useful to examine the energies of each symmetrically inequivalent configuration. In Table VI, the energies of the independent configurations obtained from GULP are listed relative to the lowest energy configuration of the Th30U2O64 (uranium mole fraction of x = 0.06) system. As shown in Table VI, the energies of the different configurations do not vary much, indicating that the U distribution can be expected to be fairly random. Following the simple arguments by McLean47 for metals and confirmed for ionic materials,48 that dopant incorporation is based on the reduction of elastic strain, the smaller U cation should be easily incorporated in the thoria lattice, which is confirmed by the lack of site preference shown by the U cations in the thoria lattice.

TABLE VI.

The number of times each inequivalent configuration occurs and the relative energies (Em) of each of these configurations, and the relative energies when the degeneracy of the system is considered (Em(red)) for the Th30U2O64 supercell.

Configuration Degeneracy (Ω) Em eV Em(red) eV
48  0.0147  0.0667 
48  0.0187  0.0707 
16  0.0993 
192  0.0101  0.0024 
192  0.0077 
Configuration Degeneracy (Ω) Em eV Em(red) eV
48  0.0147  0.0667 
48  0.0187  0.0707 
16  0.0993 
192  0.0101  0.0024 
192  0.0077 

However it is important to consider the degeneracy, Ω, of the system, which here is the number of symmetrically equivalent configurations that can be represented by one inequivalent configuration. The energy of the degenerate system, Em(red), considers the degeneracy of the system and therefore differs from the energy of the corresponding structure in the full configurational space (equations (6) and (7)). Therefore the probability of finding any one configuration must consider the energy of the degenerate system (Em(red)) and not simply the energy of the configuration alone.

For Th30U2O64, as noted in Table IV, there are 496 initial configurations which, in the reduced configurational space, yield five inequivalent configurations. The degeneracy and the relative energy of each configuration with and without considering the degeneracy (Em(red) and Em, respectively) are presented in Table VI. Configuration 3, which is the lowest probability configuration seen in Figure 6, is the structure with the lowest energy originally (Em = 0). However, as this configuration only occurs 16 times in total, compared to the Ω values of 48 for configurations 1 and 2 and 192 for configurations 4 and 5, configuration 3 has the highest Em(red) and therefore the lowest probability.

In the two dominant configurations at x = 0.063 and the lowest energy structure of x = 0.094 (Figure 9), there is a distortion of the lattice. This same distortion was observed in 2x2x2 supercells at the other concentrations considered, and can be explained by the changed cation-cation and cation-anion distances in the lower energy configurations as compared to the original lattice positions and the highest energy (lowest probability) structures. Due to the U-U and U-O distances being shorter than the corresponding Th-Th and Th-O distances, as uranium substitutions are made in the cell, oxygen atoms move closer to the U atom. This moves some of the nearby Th atoms away, as the distance between Th and oxygen atoms is preferably longer. These movements create local distortions in those regions of the cell, particularly where two uranium atoms are close together, as observed in Figure 8a and b and Figure 9. When the energies of the Th29U3O64 (U mole fraction of x = 0.094) are considered along with the degeneracy of each independent configuration (Table VII), Configuration 6 has the highest relative energy, both including the degeneracy of the configuration and without. The lowest energy configuration is Configuration 10, but the configurations with the lowest degenerate energy, Ed, are configurations numbered 13, 12, and 9. These three structures each occur 768 times in the full configurational space, the highest degeneracy of any of the 14 independent configurations. These configurations are therefore also the three highest probability configurations, with probabilities of 17%, 16%, and 14% respectively. However, as Figure 7. illustrates, no single configuration is dominant at a uranium mole fraction of x = 0.094.

FIG. 9.

Th29U3O64 a. Maximum probability; b. Lowest probability.

FIG. 9.

Th29U3O64 a. Maximum probability; b. Lowest probability.

Close modal
TABLE VII.

The number of times each inequivalent configuration occurs and the relative energies (Erel) of each of these configurations, and the relative energies when the degeneracy of the system is considered (Edrel) the Th29U3O64 supercell.

Configuration Degeneracy (Ω) Erel eV Edrel eV
96  0.0304  0.1122 
96  0.0158  0.0977 
384  0.0173  0.0394 
384  0.0123  0.0344 
384  0.0147  0.0368 
32  0.0388  0.168 
192  0.0212  0.0732 
192  0.0165  0.0685 
768  0.0188  0.0111 
10  384  0.0221 
11  256  0.0122  0.0517 
12  768  0.0102  0.0024 
13  768  0.0078 
14  256  0.0054  0.0449 
Configuration Degeneracy (Ω) Erel eV Edrel eV
96  0.0304  0.1122 
96  0.0158  0.0977 
384  0.0173  0.0394 
384  0.0123  0.0344 
384  0.0147  0.0368 
32  0.0388  0.168 
192  0.0212  0.0732 
192  0.0165  0.0685 
768  0.0188  0.0111 
10  384  0.0221 
11  256  0.0122  0.0517 
12  768  0.0102  0.0024 
13  768  0.0078 
14  256  0.0054  0.0449 

For all concentrations of uranium investigated within the 2x2x2 supercell, there is no single dominant configuration. Above x = 0.063, there are no longer even two dominant structures and instead several configurations share nearly equal probabilities and none of these probabilities exceed 17%. We do find that in the highest probability (lowest energy) configurations at all values of x, the U–U distances are smaller than these distances in the lowest probability (highest energy) configurations, but since there is a large number of possible independent configurations available in the system and no single one is dominant, including the lowest energy configuration, we predict that uranium atoms will be distributed throughout the cell without forming clusters.

From the average configurational enthalpies computed with the SOD code, using equation (5). we have calculated the enthalpies of mixing (ΔHmix) at several temperatures for a range of solid solutions of the 1x1x2 supercell system of Th(1−x)UxO2. By using a smaller supercell than that evaluated for the configurational energies and probability, the number of inequivalent configurations and hence computational cost was reduced and the full solid solution could be calculated. All values of ΔHmix are positive as shown in Figure 10, indicating that mixing of the two solids is endothermic. For an ideal solid solution, by definition ΔHmix = 0, whereas Th(1−x)UxO2 deviates from an ideal solid solution, with a maximum ΔHmix of approximately 2.5 kJmol−1 near the midpoint of the solid solution range. From the values of ΔHmix, we have calculated the values of ΔGmix (Figure 11), to see if mixing is spontaneous at temperatures above 298 K. To determine the vibrational contributions to the free energy of the system, the k-points were checked in GULP until the energy converged. We have used a k-points mesh of 2x4x4 to calculate the phonon frequencies of Th(1−x)UxO2. When vibrational effects are included at 600 K (Figure 12), ΔGmix is lowered, showing the importance of including vibrational contributions to ΔGmix to obtain more accurate values. On close examination of the mixing curves in Figure 12, a miscibility gaps is observed. However, the gap is only 0.1 kJmol−1 without consideration of vibrational effects, increasing to just 0.05 kJmol−1 when vibrational contributions are included. As this miscibility gap is so small in energetic terms, we do not expect to observe the formation of two separate domains at 600 K. Additionally, a slight shift towards the lowest values of x is observed in the ΔGmix curves in both Figure 11 and Figure 12 with the more thorium-rich Th(1−x)UxO2 systems having a lower ΔGmix than uranium-rich Th(1−x)UxO2.

FIG. 10.

ΔHmix calculated using a 1x1x2 supercell.

FIG. 10.

ΔHmix calculated using a 1x1x2 supercell.

Close modal
FIG. 11.

ΔGmix at different temperatures calculated for a 1x1x2 supercell, not considering vibrational effects.

FIG. 11.

ΔGmix at different temperatures calculated for a 1x1x2 supercell, not considering vibrational effects.

Close modal
FIG. 12.

ΔGmix at 600 K showing the configurational and the configurational and vibrational effects.

FIG. 12.

ΔGmix at 600 K showing the configurational and the configurational and vibrational effects.

Close modal

Understanding the heat capacity and thermal expansion is important in any nuclear fuel material in order to predict the temperature of the fuel in the reactor. We have calculated these properties for the full range of the ThO2-UO2 solid solution in the 1x1x2 supercell and compared our values to those found experimentally or using other computational methodologies where available. The results of the constant-volume heat capacity calculations for the pure materials and the solid solutions are presented in Figure 13 and agree closely experiment. As the Cv values of ThO2 we calculated are in line with the DFT values and have used a less compute-intensive methodology than the value found by Lu et al44 at all temperatures studied here, we expect that this methodology will be useful in determining Cv of even lower concentrations of uranium in the ThO2-UO2 solid solution, or other nuclear materials.

FIG. 13.

Constant volume heat capacity of Th(1−x)UxO2 solid solutions.

FIG. 13.

Constant volume heat capacity of Th(1−x)UxO2 solid solutions.

Close modal

The coefficient of thermal expansion, α, of pure thoria has been discussed in section III B but in MOX fuels it is important to understand the effect of both cations on the thermal expansion. Generation IV reactors are typically designed to withstand temperatures up to 1873 K, so we have calculated α from room temperature (298 K) to 1900 K to represent typical operating temperatures of this fuel type. In Figure 14, we present the results of Th(1−x)UxO2 where x = 0.125, 0.25, and 0.38 compared to ThO2. The values of α are in good agreement with the pure material at all temperatures calculated. There is particularly close agreement at low temperatures, with deviation from the pure material increasing above approximately 1000 K.

FIG. 14.

Coefficient of thermal expansion of Th(1−x)UxO2 solid solutions.

FIG. 14.

Coefficient of thermal expansion of Th(1−x)UxO2 solid solutions.

Close modal

We have developed a new interatomic shell-model potential that can be used to accurately model the structure and properties of pure ThO2. This potential can be combined with a leading UO2 potential to investigate uranium-doped thorium dioxide. In this study we have employed these potentials in a study of the distribution of uranium atoms in thoria supercells and we have closely reproduced available experimental data for the lattice parameters, coefficient of thermal expansion, and constant volume heat capacity.

We have shown that uranium atoms do not cluster in the bulk material but are instead distributed throughout the cell. In a Th(1−x)UxO2 solid solution the enthalpy of mixing is endothermic across the full range of the solution, but while the Gibbs free energy of mixing is favorable above 600 K, when vibrational effects are included in the calculations, the free energies of mixing show a very small miscibility gap. However, as the energy penalty is so small, we do not predict the formation of separate domains at high or low concentrations of uranium in the thoria lattice.

Via our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work made use of the facilities of HECToR and ARCHER, the UK’s national high-performance computing service, which is funded by the Office of Science and Technology through EPSRC’s High End Computing Programme. AS gratefully acknowledges funding from the Molecular Modeling and Materials Science Centre for Doctoral Training at UCL, and NHdL thanks the Royal Society for an Industry Fellowship and AWE for funding.

1.
D.C.
Crawford
,
D.L.
Porter
, and
S.L.
Hayes
,
J. Nucl. Mater.
371
,
202
(
2007
).
3.
J.S.
Herring
,
P.E.
MacDonald
,
K.D.
Weaver
, and
C.
Kullberg
,
Nucl. Eng. Des.
203
,
65
(
2001
).
4.
C.
Ronchi
and
J.P.
Hiernaut
,
J. Alloys Compd.
240
,
179
(
1996
).
5.
J.
Banerjee
,
S.C.
Parida
,
T.R.G.
Kutty
,
A.
Kumar
, and
S.
Banerjee
,
J. Nucl. Mater.
427
,
69
(
2012
).
6.
S.
Dash
,
S.C.
Parida
,
Z.
Singh
,
B.K.
Sen
, and
V.
Venugopal
,
J. Nucl. Mater.
393
,
267
(
2009
).
7.
P.
Martin
,
D.J.
Cooke
, and
R.
Cywinski
,
J. Appl. Phys.
112
,
073507
(
2012
).
8.
R.
Kandan
,
R.
Babu
,
P.
Manikandan
,
R.
Venkata Krishnan
, and
K.
Nagarajan
,
J. Nucl. Mater.
384
,
231
(
2009
).
9.
G.C.
Benson
,
P.J.
Freeman
, and
D. E.
,
J. Am. Ceram. Soc.
46
,
43
(
1963
).
10.
W.C.
Mackrodt
and
R.F.
Stewart
,
J. Phys. C Solid State Phys.
12
,
431
(
1979
).
11.
E.A.A.
Colbourn
and
W.C.C.
Mackrodt
,
J. Nucl. Mater.
118
,
50
(
1983
).
12.
M.
Nadeem
,
M.J.
Akhtar
,
R.
Shaheen
, and
M.N.
Haque
,
J. Mater. Sci. Technol.
17
,
638
(
2001
).
13.
M.
Osaka
,
J.
Adachi
,
K.
Kurosaki
,
M.
Uno
, and
S.
Yamanaka
,
J. Nucl. Sci. Technol.
44
,
1543
(
2007
).
14.
R.K.
Behera
and
C.S.
Deo
,
J. Phys. Condens. Matter
24
,
215405
(
2012
).
15.
C.R.A.
Catlow
,
Proc. R. Soc. A Math. Phys. Eng. Sci.
353
,
533
(
1977
).
16.
J.D.
Gale
,
J. Chem. Soc. Faraday Trans.
93
,
629
(
1997
).
17.
J.D.
Gale
and
A.L.
Rohl
,
Mol. Simul.
29
,
291
(
2003
).
18.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Oxford University Press
,
Oxford
,
1954
).
19.
20.
B.G.
Dick
and
A.W.
Overhauser
,
Phys. Rev.
112
,
90
(
1958
).
21.
D.G.K.
Upadhyay
,
Solid State Physics: Lattice Dynamics of Ionic Solids
(
Laxmi Publications
,
New Delhi
,
2008
).
22.
R.
Grau-Crespo
,
N.H.
de Leeuw
, and
C.R.A.
Catlow
,
J. Mater. Chem.
13
,
2848
(
2003
).
23.
R.
Grau-Crespo
,
N.H.
de Leeuw
, and
C.R.A.
Catlow
,
Chem. Mater.
16
,
1954
(
2004
).
24.
R.
Grau-Crespo
,
S.
Hamad
,
C.R.A.
Catlow
, and
N.H.
de Leeuw
,
J. Phys. Condens. Matter
19
,
256201
(
2007
).
25.
S.E.
Ruiz-Hernandez
,
R.
Grau-Crespo
,
A.R.
Ruiz-Salvador
, and
N.H.
De Leeuw
,
Geochim. Cosmochim. Acta
74
,
1320
(
2010
).
26.
Q.
Wang
,
R.
Grau-Crespo
, and
N.H.
de Leeuw
,
J. Phys. Chem. B
115
,
13854
(
2011
).
27.
S.
Benny
,
R.
Grau-Crespo
, and
N.H.
de Leeuw
,
Phys. Chem. Chem. Phys.
11
,
808
(
2009
).
28.
R.
Grau-Crespo
,
A.Y.
Al-Baitai
,
I.
Saadoune
, and
N.H.
De Leeuw
,
J. Phys. Condens. Matter
22
,
255401
(
2010
).
29.
J.
González-López
,
S.E.
Ruiz-Hernández
,
Á.
Fernández-González
,
A.
Jiménez
,
N.H.
de Leeuw
, and
R.
Grau-Crespo
,
Geochim. Cosmochim. Acta
142
,
205
(
2014
).
30.
S.
Haider
,
R.
Grau-Crespo
,
A.J.
Devey
, and
N.H.
de Leeuw
,
Geochim. Cosmochim. Acta
88
,
275
(
2012
).
31.
Y.
Seminovski
,
P.
Palacios
,
P.
Wahnón
, and
R.
Grau-Crespo
,
Appl. Phys. Lett.
100
,
102112
(
2012
).
32.
R.
Grau-Crespo
,
K.C.
Smith
,
T.S.
Fisher
,
N.H.
de Leeuw
, and
U. V.
Waghmare
,
Phys. Rev. B
80
,
174117
(
2009
).
33.
A.R.
Ruiz-Salvador
,
R.
Grau-Crespo
,
A.E.
Gray
, and
D.W.
Lewis
,
J. Solid State Chem.
198
,
330
(
2013
).
34.
K.C.
Smith
,
T.S.
Fisher
,
U. V.
Waghmare
, and
R.
Grau-Crespo
,
Phys. Rev. B
82
,
134109
(
2010
).
35.
M.
Idiri
,
T.
Le Bihan
,
S.
Heathman
, and
J.
Rebizant
,
Phys. Rev. B
70
,
1
(
2004
).
36.
P.M.
Macedo
,
W.
Capps
, and
J.B.O.
Wachtman
,
J. Am. Ceram. Soc.
47
,
651
(
1964
).
37.
J.
Staun Olsen
,
L.
Gerward
,
V.
Kanchana
, and
G.
Vaitheeswaran
,
J. Alloys Compd.
381
,
37
(
2004
).
38.
R.D.
Shannon
,
Acta Crystallogr. Sect. A
32
,
751
(
1976
).
39.
D.
Olander
,
Fundamental Aspects of Nuclear Reactor Fuel Elements
(
Energy Research and Development Administration
,
1976
).
40.
41.
J.A.C.
Marples
, in
Plutonium 1975 Other Actinides 5th Int. Conf. Plutonium Other Actinides 1975 Proc.
, edited by
H.
Blank
and
R.
Lindner
(
North Holland Publishing Company
,
Amsterdam
,
1976
).
42.
C.P.
Kempter
and
R.O.
Elliott
,
J. Chem. Phys.
30
,
1524
(
1959
).
43.
T.
Yamashita
,
N.
Nitani
,
T.
Tsuji
, and
H.
Inagaki
,
J. Nucl. Mater.
245
,
72
(
1997
).
44.
Y.
Lu
,
Y.
Yang
, and
P.
Zhang
,
J. Phys. Condens. Matter
24
,
225801
(
2012
).
45.
J.-J.
Ma
,
J.-G.
Du
,
M.-J.
Wan
, and
G.
Jiang
,
J. Alloys Compd.
627
,
476
(
2015
).
46.
S.
Hubert
,
J.
Purans
,
G.
Heisbourg
,
P.
Moisy
, and
N.
Dacheux
,
Inorg. Chem.
45
,
3887
(
2006
).
47.
D.
McLean
,
Grain boundaries in Metals
(
Clarendon Press
,
Oxford
,
1957
).
48.
N.H.
de Leeuw
and
S.C.
Parker
,
J. Chem. Phys.
112
,
4326
(
2000
).