First-Principle Studies on Local Lattice Distortions and Thermodynamic Properties in Non-Stoichiometric Thorium Monocarbide
1Engineering Research Center of Nuclear Technology Application, Ministry of Education (East China Institute of Technology), Nanchang 330013, China; qlwei11@ecut.edu.cn (Q.W.); wuyy20@ecut.edu.cn (Y.W.)
2School of Nuclear Science and Engineering, East China University of Technology, Nanchang 330013, China; 2022120631@ecut.edu.cn
3Engineering Technology Research Center of Nuclear Radiation Detection and Application, Nanchang 330013, China
4Institute of High Energy Physics, Chinese Academy of Sciences (CAS), Beijing 100049, China
5Spallation Neutron Source Science Center, Institute of High Energy Physics, Chinese Academy of Sciences (CAS), Dongguan 523803, China
6Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan 030006, China
*Correspondence: ybliu@ecut.edu.cn (Y.L.); wangbt@ihep.ac.cn (B.W.)Abstract
Thorium monocarbide (ThC) is interesting as an alternative fertile material to be used in nuclear breeder systems and thorium molten salt reactors because of its high thermal conductivity, good irradiation performance, and wide homogeneous composition range. Here, the influence of carbon vacancy site and concentration on lattice distortions in non-stoichiometric ThC1−x (x = 0, 0.03125, 0.0625, 0.125, 0.1875, 0.25, or 0.3125) is systematically investigated using first-principle calculations by the projector augmented wave (PAW) method. The energy, mechanical parameters, and thermodynamic properties of the ThC1-x system are calculated. The results show that vacancy disordering has little influence on the total energy of the system at a constant carbon vacancy concentration using the random substitution method. As the concentration of carbon vacancies increases, significant lattice distortion occurs, leading to poor structural stability in ThC1−x systems. The changes in lattice constant and volume indicate that ThC0.75 and ThC0.96875 represent the boundaries between two-phase and single-phase regions, which is consistent with our experiments. Furthermore, the structural phase of ThC1−x (x = 0.25–0.3125) transforms from a cubic to a tetragonal structure due to its ‘over-deficient’ composition. In addition, the elastic moduli, Poisson’s ratio, Zener anisotropic factor, and Debye temperature of ThC1-x approximately exhibit a linear downward trend as x increases. The thermal expansion coefficient of ThC1−x (x = 0–0.3125) exhibits an obvious ‘size effect’ and follows the same trend at high temperatures, except for x = 0.03125. Heat capacity and Helmholtz free energy were also calculated using the Debye model; the results showed the C vacancy defect has the greatest influence on non-stoichiometric ThC1−x. Our results can serve as a theoretical basis for studying the radiation damage behavior of ThC and other thorium-based nuclear fuels in reactors.
1. Introduction
With increasing demand for electricity and the depletion of uranium resources, the introduction of new nuclear fuels into the fuel cycle has become critical [1,2]. Thorium is a potential convertible nuclear energy resource which is which is approximately three to four times more abundant in the Earth’s crust than uranium [3,4]. In recent years, the development of a thorium fuel cycle has attracted considerable interest worldwide with the purpose of saving uranium reserves and further reducing the production of long-lived minor actinides [5,6]. Actinide carbides are considered to be one of the most promising nuclear fuel materials of generation IV reactors [7,8]. Recently, thorium-based carbides have attracted great attention because of their high melting points, corrosion resistivity, low thermal expansion coefficients, and high thermal conductivity [9]. Therefore, understanding the behavior and properties of thorium-based nuclear fuel is essential for exploring its potential application as nuclear reactor fuel material [10].
The Th-C system has two basic phases: thorium monocarbide (ThC) and thorium dicarbide (ThC2) [11]. Cubic (B1-type) ThC has a wide non-stoichiometric region in the carbon sublattice, ThC1−x (0 < x < 0.33) [12]. Fuel is irradiated in the reactor to produce non-stoichiometric ThC1−x, which may affect the thermodynamic performance of the fuel [13]. Experimentally, Satow et al. [14] concluded that the lattice parameter increases almost linearly with increasing carbon concentration between the compositions ThC0.68 and ThC0.95, while it remains a constant as carbon concentration is lower than ThC0.68 and greater than ThC0.95. Theoretically, ThC is metallic and structurally stable in the ground state [15]. The formation energy of carbon vacancies in ThC0.75 and ThC0.5 has shown that ThC can easily create carbon vacancies [16,17]. The relative stabilities of the fan-type and linear structures of gas-phase ThCn (n = 1–7) clusters were also investigated with DFT calculations by Yang et al. [18]. In addition, the high-pressure phase transition of ThC has been studied experimentally and theoretically. Yu et al. [19] experimentally revealed the phase transition of ThC from B1 to P4/nmm at ~58 GPa by synchronous X-ray diffraction. There is no phase transition in ThC under high pressures at 36 GPa [20] and 40–45 GPa [21], but a transitional P4/nmm phase is produced at 60–120 GPa [22], theoretically.
As an important nuclear energy material, it is well known that defects in ThC are unavoidable due to irradiation damage from high-energy neutrons. Therefore, it is necessary to study the structural stability of non-stoichiometric ThC. Most existing research focuses on stoichiometric ThC and its related phase transition at high pressures. There is less literature available on the lattice distortions and structural stability of non-stoichiometric ThC. In this study, considering the influence of the site and concentration of carbon vacancies on a non-stoichiometric ThC1−x (x = 0, 0.03125, 0.0625, 0.125, 0.1875, 0.25, or 0.3125) system, its lattice distortions, mechanical parameters, and thermodynamic properties were calculated.
2. Calculation Methods and Models
2.1. Calculation Method
The calculations were conducted using the VASP package [23] based on DFT [24], employing the projector augmented wave (PAW) method [25]. The exchange-correlation functional used to describe the interactions was the generalized gradient approximation described by Perdew, Burke, and Ernzerh (GGA-PBE) [26]. Twelve electrons (6s26p65f06d27s2) for Th and four electrons (2s22p2) for C were used as valence electrons in the ThC1−x system. Th contains only a small number of 5f states, and it is generally accepted that these states are itinerant: their nature does not need to be corrected with the Hubbard model [27]. Brillouin-zone integrations were carried out with Methfessel–Paxton [28] smearing with a width of 0.2 eV. Through convergence testing, the cutoff energy of atomic wave functions was set to 520 eV for all calculations. The Brillouin zone was sampled with a 9 × 9 × 9 k-point mesh for the 8-atom cell and a 5 × 5 × 5 k-point mesh for the 64-atom supercell using the Monkhorst and Pack (MP) scheme [29]; both meshes were proven to be sufficient for an energy convergence of less than 1.0 × 10−5 eV/atom and a force convergence of less than 0.02 eV/Å. The calculation details of p k-point mesh are shown in Appendix A.
2.2. Calculation Models
Under normal temperature and pressure conditions, ThC has the face-centered cubic structure of NaCl (B1), belonging to the crystal system. The atomic coordinates of Th are (0, 0, 0) and those of C are (0.5, 0.5, 0.5). Lattice parameters (a0 = 5.3510 Å, α = β = γ = 90°) were obtained from the optimized lattice structure. This is consistent with most theoretical values (5.335–5.388 Å) [7,13,17,22,30,31] and is close to the experimental values of 5.344 Å [14] and 5.430 Å [19]. ThC1−x (x = 0.03125, 0.0625, 0.125, 0.1875, 0.25, or 0.3125) with specific vacancy concentrations was created by random substitution method obeying Lowenstein’s rule [32], corresponding to the replacement of 1, 2, 4, 6, 8, or 10 carbon atoms with vacancies in a 64-atom supercell, respectively. The 8-atom unit cell structure and typical representatives of the 2 × 2 × 2 supercell structures of ThC1−x (x = 0–0.3125) are shown in Figure 1.
2.3. Crystal and Vacancy Formation Energies
ThC1−x can form from metal Th and the most stable graphite through the reaction. The formation energy per atom in the ThkCl supercell, Eform(ThC1−x), is expressed by Equation (1) [33]: where Etot(ThkCl) is the total energy of the ThkCl supercell and Etot(Th, Cg) is the energy per Th or C atom of each chemical species in its reference state. Here, the reference states are the ground state crystalline phases of Th and C, namely the thorium α phase and the carbon graphite phase. k and l are the numbers of Th and C atoms, respectively. According to this definition, a negative Eform means that the ThC1−x phase is thermodynamically stable, and the lower the formation energy is, the more stable the state is [34,35].
Vacancy formation energy (Evf) is obtained using Equation (2) [36]:
In Equation (2), positive values of Evf mean that the ThC1−x system is still stable as a result of the formation of carbon vacancies, i.e., stable non-stoichiometric phases are formed and vice versa but its stability will reduce. Certainly, such predictions are based only on thermodynamics and do not consider the kinetics of reactions.
2.4. Elastic Properties
ThC has cubic crystal system; the space group is 225, which has the highest symmetry degree among all crystal systems. Its independent stiffness matrix element number is only 3, that is c11, c12, and c44. In a cubic crystal system, three independent elastic constants satisfy the following relationship to maintain material stability, as prescribed by the Born–Huang criterion [37]:
Mechanical parameters, such as the bulk modulus (B), shear modulus (G), and Young’s modulus (E), are calculated to assess the influence of x on the structural stability of ThC1−x. These calculations are carried out by the Voigt–Reuss–Hill approximation [38] using the elastic constants (cij) of the crystal system, as shown in Equations (4)–(9) [39]. where BV and BR are the Voigt and Reuss bulk moduli, and GV and GR are the Voigt and Reuss shear moduli, respectively. B and G are arithmetic means of the Voigt and Reuss elastic moduli, expressed as:
From which Poisson’s ratio (ν) is given by: and the Zener anisotropy factor (A) [13] is given by:
The calculation method of the above parameters of other symmetrical structures can be referred to in Ref. [39].
2.5. Thermodynamic Properties
Thermodynamic properties are calculated on the basis of the Debye model. Debye temperature (ƟD) is an important fundamental parameter closely related to many physical properties, such as specific heat and melting temperature. At low temperatures, ƟD calculated from elastic constants is the same as that determined from specific heat measurements [13]. We calculated ƟD from the elastic constants using average wave velocity, vm, by the following common relation [13]: where vm is calculated by: and vl and vt are based on the elastic constant: where h is the Planck constant, k is the Boltzmann constant, NA is Avogadro’s number, ρ is the density of the crystal in g·cm−3, M is the molar mass of the crystal in g·mol−1, n is the number of atoms in a unit cell, and vl and vt are the longitudinal and transverse elastic wave velocities m·s−1, respectively.
The volumetric thermal expansion coefficient, αV (T), is then obtained from V (T) using: where V is the equilibrium volume at 0 K.
In addition, heat capacity and Helmholtz free energy are also calculated on the basis of the Debye model. At low temperatures, heat capacity is integrated to obtain [40]:
Helmholtz free energy A(T) is then obtained [40]:
E is the energy of ThC1x system, S is entropy, and T is absolute temperature, same as above.
3. Results and Discussion
3.1. Random Substitution
In order to find out the effect of carbon vacancy sites on the structural stability of non-stoichiometric ThC1−x, ten groups of models (represented by A, B, C, ..., I, and J) of each carbon vacancy concentration in ThC1−x (x = 0.03125, 0.0625, 0.125, 0.1875, 0.25, or 0.3125) were established using the random substitution method. The total energies of the optimized systems are shown in Figure 2.
As shown in Figure 2a, Etot of stoichiometric ThC is −552.069 eV, which is lower than the Etot of all non-stoichiometric ThC1−x. Etot of non-stoichiometric ThC1−x gradually increases as carbon vacancy concentration increases. These results indicate that ThC is the most stable structure. It is clear from Figure 2b that Etot is linearly related to carbon vacancy concentration: the smallest standard deviation (SD) reached 0.013% (group A), and the largest SD (group J) does not exceed 1.147%, with a coefficient of variation less than 0.032%. We conclude that vacancy-ordering effects [36] can be ignored in the non-stoichiometric ThC1−x system modeled using the random substitution method.
3.2. Structural Properties and Formation Energy
The lattice parameters, crystal formation energy, and carbon vacancy formation energy for different carbon vacancy concentrations of ThC1-x are shown in Table 1. −Δa/a0 is the change rate of lattice parameter a. As seen in Table 1, a equals to 5.3512 Å for ThC0.96875, 5.3427 Å for ThC0.875, and 5.3257 Å for ThC0.75, respectively, which are close to the existing experimental data of 5.3470 Å for ThC0.975, 5.3429 Å for ThC0.891 [41], and 5.31 Å for ThC0.70 [11]. When the range of x is from 0 to 0.3125, namely in the case of ThC → ThC0.6875, −Δa/a0 lies within 0.50%, while the values of α, β, and γ remain relatively stable, at 90.00° ± 0.25°, with distortion rates within ±0.28%. When x is less than 0.25, namely in the case of ThC → ThC0.75, the ratio of c/a is still equal to 1. It indicates that the crystal system can maintain a cubic structure, which is consistent with the result of Shein et al. [36]. When carbon vacancy concentration is increased, in the case of ThC0.75 → ThC0.6875, the ratio of c/a is equal to 0.997 and 0.994, respectively. These results show that significant structural distortions occur due to the ‘over-deficient’ composition of ThC1−x. At the same time, the change in lattice volume is calculated; when x = 0.125 and 0.3125, then the cell lattice volume decreases by 0.463% and 1.470%, respectively.
The relationship between the lattice parameter (a) and carbon concentrations (1 − x) is shown in Figure 3. The calculated lattice parameters are larger than the experimental values presented by Satow et al. [14], which may have been caused by the GGA algorithm, but they have the same trend. The lattice parameters we calculated decrease almost linearly with increasing carbon concentration between the compositions of ThC0.75 and ThC0.96875, while constant values were obtained for carbon concentrations lower than ThC0.75 and greater than ThC0.96875. The values of the boundary range are in good agreement with those of Satow et al. [14], which were obtained by three different experimental methods. This demonstrates the reliability of our calculation results. Furthermore, the occurrence of breaks at ThC0.75 and ThC0.96875 is considered to indicate the boundaries between the two-phase and single-phase regions. We can infer that ThC1-x belongs to the two-phase regions of Th + ThC and ThC + ThC2 when the C/Th ratio is lower than 0.75 and greater than 0.96875, respectively. Additionally, this is consistent with the experimental results of non-stoichiometric UC [42].
As shown in Table 1, all Eform of ThC1−x containing carbon vacancies are negative, and their values gradually increase as x increases. Eform of perfect ThC is −0.444 eV per atom, which was consistent with the calculation results of Shein et al. [12,36]. In addition, perfect ThC is the most stable compound. This result is consistent with the results calculated for uranium monocarbide (UC) [33]. Our calculated Eform of ThC0.75 was −0.358 eV per atom, which is also consistent with the result calculated by Shein et al. [36]. Calculation results for other non-stoichiometric ThC1−x (x = 0–0.3125) systems have not been reported in the literature.
Evf is positive, and its value increases as x rises, indicating that ThC can easily form carbon vacancies. However, compared with the stoichiometric ThC system, the stability of the non-stoichiometric system is reduced. Evf of ThC0.75 is equal to 0.172 eV, which is very consistent with the results obtained by Daraco et al. [17] using the GGA method with a 64-atom supercell, but lower than those obtained by Wang et al. [7] and Shein et al. [36] using 8-atom supercells, possibly due to the size effect [15].
3.3. Elastic Moduli
Elastic moduli are important parameters for characterizing the stability of materials [43]. We calculated the second-order elastic constants (cij) at the equilibrium lattice parameter by using the ‘stress-strain’ technique [44], as shown in Table 2. For stoichiometric ThC, c11, c12, and c44 are lower than the result of Aydin et al. [13], while our c11 and c12 are in good agreement with the theoretical analysis [15,45]. c44 is also consistent with [10] and within the range of [13,45].
Table 2 shows that ThC1−x (x < 0.25) is a cubic crystal. Each elastic constant matrix is determined by three variables (c11, c12, and c44), meeting the material stability condition of the Born–Huang criterion in Equation (3). It indicates ThC1−x (x = 0–0.25) is structurally stable. When x is greater than 0.25, there are six elastic constant variables (c11, c33, c44, c66, c12, and c13), and all of them also satisfy the stability conditions of a tetragonal crystal system: c11 > 0, c33 > 0, c44 > 0, c66 > 0, (c11 − c12) > 0, (c11 + c33 − 2c13) > 0, and [2(c11 + c12) + c33 + 4c13] > 0 [39]. It is indicated that ThC1−x remains stable when x is between 0.25 and 0.3125. However, c33 > c11 and c66 > c44, while c13 < c12 for ThC1−x (x = 0.25–0.3125). Those changes may affect the symmetry of the system. In conclusion, it can be seen that non-stoichiometric ThC1−x crystals can still maintain a stable structure despite significant lattice distortion for x = 0–0.3125.
All elastic constants decrease as vacancy concentration increases. The values of c11 are higher than those of c12 and c44. c11 represents elasticity in length, and longitudinal strain produces a change in c11. c12 and c44 are related to elasticity in shape, which is a shear constant, and transverse strain causes a change in shape [13]. As shown in Figure 4, c12 and c44 decrease more significantly than c11 as carbon vacancy concentration increases, and variation in c44 is perfectly linear from ThC0.96875 to ThC0.75. In contrast, shear constant c44 is important in NaCl structures because it is the modulus most sensitive to next-nearest neighbor, or atom-like, interactions [42]. Thus, c44 is expected to be the most sensitive to changes in carbon vacancy concentration in ThC1−x. These results also imply that ThC1−x, which exists as Th + ThC, can deviate from stoichiometry. This is consistent with the linear variation in lattice parameters with stoichiometry and can be explained by assuming that the structure of hypo-stoichiometric ThC1−x primarily contains free thorium, with some vacancies.
The bulk modulus (B), shear modulus (G) and Young’s modulus (E) were calculated using elastic constants (cij) and are shown in Figure 5. For stoichiometric ThC, the calculated bulk modulus is 131.61 GPa, which differs from the experimental value by 11.6 % (147 GPa for ThC0.95 at 300 K) [19], but the value agrees quite well (a difference of less than 1.0%) with the data calculated by Aydin et al. [13] (130.2 GPa) and Daraco et al. [45] (131.15 GPa). For non-stoichiometric ThC0.75, the calculated bulk modulus is 98.08 GPa, which differs from the experimental data by 10% (109 GPa for ThC0.76) [19,20].
As shown in Figure 5, B, G, and E all decrease as x increases. B changes by 10% for each 10% change in carbon atom concentration from the initial state to the final state. For comparison, Routbort et al. [42] studied the dependence of elastic moduli in UC on stoichiometry and found that the bulk modulus changes by 2% for each 10% change in carbon concentration. This indicates that ThC may be more prone to lattice distortion than UC when carbon vacancy defects are generated.
B/G ratios is also presented in Figure 6. According to the Pugh criterion [46], a material with a B/G ratio higher than 1.75 is considered ductile, while one with a B/G ratio lower than 1.75 is considered brittle [47]. We calculated the B/G ratio of ThC, which was found to be 1.82, thus indicating ductile behavior. The B/G of non-stoichiometric ThC1−x decreases as x increases, indicating that the ductility of non-stoichiometric ThC1−x decreases with an increase in carbon vacancies. When x is larger than 0.9375, non-stoichiometric ThC1−x would become brittle because B/G is less than 1.75.
In addition, Poisson’s ratio (ν) is a very important property for industrial applications because it provides more information about the characteristics of bonding forces rather than elastic constants [48]. As shown in Figure 6, the calculated ν value is equal to 0.27 for ThC at 0 GPa, and agrees with other theoretical values of 0.26 [13] and 0.28 [45]. It is concluded that interatomic forces are dominant in ThC. Moreover, according to the Poisson’s ratio criterion [49], in general, the ν value of a ductile material is approximately 1/3 and is less than 1/3 for a brittle material. As shown in Table 2, all of the values are less than 1/3 and decrease (from 0.27 to 0.24) as x increases. These are within the range (from 0.25 to 0.45) for typical metals, except for ThC0.6875. This indicates a reduction in ductility due to its ‘over-deficient’ composition.
In contrast, when the Zener anisotropy factor (A) is equal to 1.0, it indicates that a material is completely isotropic. There is no evident linear relationship between the Zener anisotropy factor and carbon vacancy concentration. As shown in Table 2, the calculated A value of ThC is 1.22, which is greater than the experimental results of 0.97 [45] and 0.98 [13], while the numerical result agrees well with the experimental value of 1.13 [15]. As shown in Figure 6, the values of non-stoichiometric ThC1−x (x = 0–0.3125) decrease as x increases, except for ThC0.75, and most of the values are close to 1, indicating that the anisotropy of ThC1−x (x = 0–0.3125) is small.
3.4. Debye Temperature and Thermal Expansion Coefficient
The relationships of the Debye temperature (ƟD), longitudinal wave velocity (νl), transverse elastic wave velocity (νt), and average wave velocity (νm) with (1 − x) for ThC1−x are shown in Figure 7. For stoichiometric ThC, the calculated vl, vt, and vm are 2618, 4648, and 2912 m·s−1, respectively. The results are in good agreement with those of the elastic constants calculated by Wang et al. (2657, 4709, and 2940 m·s−1) [15]. The resulting Debye temperature (ƟD) is 324.1 K, which is in good agreement with that calculated by Wang et al. [15] and Daroca et al. [45] (328 and 311 K, respectively) using the same method of elastic constants. This value is larger than the 298 and 280 K obtained by fitting the isochoric heat capacity curve at low temperatures [45,50] and the experimental value of 262 K [51] obtained on the basis of isobaric heat capacity measurements.
For non-stoichiometric ThC1−x, νl, νt and ƟD decrease as x increases. The change in νm is almost the same as that in νl. ƟD of non-stoichiometric ThC1−x is in the range of 320.0 to 269.6 K when x is from 0.03125 to 0.3125.
The thermal expansion coefficient (α) can be obtained from the temperature derivative of the lattice constant given in Equation (13). Variations in the thermal expansion coefficient with carbon vacancy concentration in the temperature range of 0–1000 K for ThC1−x are presented in Figure 8. It is noted that α rapidly increases with T at low temperatures and achieves saturation at approximately 350 K. In addition, α increases as x increases at a constant temperature, except for x = 0.03125. It is possible that carbon vacancy defects lead to the fracture of the covalent Th-C bond. The lattice volume reduction caused by a single carbon atom defect is not enough to offset the volume swelling caused by the bond fracture. The volume relationship in Table 1 can also explain this phenomenon: although the lattice constants are almost equal, the volume of ThC0.96875 is slightly larger than that of ThC. Accordingly, we can see that the α of ThC0.96875 is slightly less than that of ThC. As x increases further, it leads to increased lattice defects and non-uniformity, resulting in a further increase in the coefficient of thermal expansion.
As for non-stoichiometric ThC1−x, Ref. [52] reported that the value of the average linear thermal expansion coefficient for ThC0.96 was 8.5 × 10−6 K−1 (using αV = 3αl, we obtain αV = 2.55 × 10−5 K−1) between 974 K and 1174 K. Our average value of αV for the same range of temperature is 3.21 × 10−5 K−1 for ThC0.96875 (for a 64-atom supercell containing four carbon vacancies), and that for ThC0.76 was 6.6 × 10−6 K−1 (αV = 1.98 × 10−5 K−1) at 974–1104 K. Our average value of αV for the same temperature range is 3.91 × 10−5 K−1 for ThC0.75 (64-atom supercell containing four carbon vacancies). The theoretical values we calculated with the supercells were also larger than those obtained in the experiments. Here, asymmetry plays a more vital role than vibrational amplitude because the symmetry of a crystal also affects its thermal expansion. The worse the symmetry of the crystal, or the more defects in the crystal, the greater the coefficient of thermal expansion. For non-stoichiometric ThC1−x, the high concentration of carbon vacancies destroys its symmetry, which might be the cause for the overestimation of the thermal expansion coefficient in ThC1−x.
3.5. Heat Capacity and Helmholtz Free Energy
We also calculated the heat capacity Cv (Figure 9) and Helmholtz free energy A(T) (Figure 10) of ThC1−x using the Debye model. Compared with perfect ThC, the C vacancy defect had the greatest influence on the heat capacity and Helmholtz free energy of ThC0.96875. This can perhaps be explained with the volume relationship presented in Table 1: although the lattice constants are almost equal, the volume of ThC0.96875 is slightly larger than that of ThC. The effect of vacancy defects of other concentrations on heat capacity and Helmholtz free energy is irregular, which is not only related to the relative position of the C vacancy, but also affected by concentration. Although a numerical comparison cannot be made directly due to the inconsistency in the unit of the cell, the changes in the calculated results of free energy are consistent with those of ThC in [8].
4. Conclusions
In summary, the crystal energy, elastic parameters, and thermodynamic properties of ThC1−x (x = 0, 0.03125, 0.0625, 0.125, 0.1875, 0.25, or 0.3125) were studied using density functional theory (DFT) in conjunction with the random substitution method. The results showed that a vacancy-disordering effect is not evident, and it is feasible to adopt the random substitution method. We observed that ‘over-deficient’ carbon vacancies could affect the structural stability of ThC, even though the calculated elastic constants still satisfy traditional mechanical stability conditions. With an increase in carbon vacancy concentration, the lattice constant decreases and is distorted. When x is greater than 0.25, ThC1−x transforms from a cubic to a tetragonal structure owing to its ‘over-deficient’ composition. Moreover, ‘over-deficient’ carbon vacancies lead to a decline in the toughness and ductility, even leading to brittleness, of non-stoichiometric ThC1−x. Our calculated results can be used to analyze the stability of ThC fuel in the process of reactor combustion, and the method described in this paper can also be used for theoretical analysis of other thorium-based nuclear fuel in the future.
Acknowledgments
We acknowledge Li-yuan Dong for helpful discussions, and thank Saeed for revising the grammar in the manuscript.
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
The data presented in this study are available on request from the corresponding authors. The data are not publicly available due to ongoing research in the project.
Conflicts of Interest
The authors declare no conflict of interest.
Appendix Group
Appendix A
The convergence of the total energy (Etot) of the system is determined by two key computational parameters, namely the plane wave cutoff energy (Ecut) and the mesh of k-points, as computed with the MedeA-VASP package. To determine the proper value of Ecut and the number of k-points, Etot of the ThC system was tested with different Ecut (from 300 to 700 eV with a step of 50 eV) and k-points (3 × 3 × 3, 4 × 4 × 4, 5 × 5 × 5, 9 × 9 × 9, and 11 × 11 × 11). The results are shown in Figure A1. For all k-points, Etot can be stabilized when Ecut reaches 520 eV. Therefore, Ecut was set as 520 eV. For the same Ecut, Etot increases as the k-points increase. The difference in Etot using 9 × 9 × 9 and 11 × 11 × 11 k-point meshes is at most 0.005 eV. Therefore, the Brillouin zone is sampled by a 9 × 9 × 9 Monkhorst–Pack (MP) k-point mesh for 8-atom unit cells in the first step of structural optimization. Then, a 5 × 5 × 5 k-point mesh is selected to calculate the system energy and mechanical properties for 64-atom supercells.
| Phase | ThC | ThC0.96875 | ThC0.9375 | ThC0.875 | ThC0.8125 | ThC0.75 | ThC0.6875 |
|---|---|---|---|---|---|---|---|
| a/Å | 5.3510 5.352 [7] 5.341 [13] 5.344 [14] 5.351 [19] 5.388 [31] 5.3878 [36] | 5.3512 5.3470 a [41] | 5.3493 | 5.3427 5.3429 b [41] | 5.3351 5.325 c [14] | 5.3257 5.31 d [11] 5.312 [36] | 5.3248 5.292 e [14] |
| c/a | 1 | 1 | 1 | 1 | 1 | 0.997 | 0.994 |
| Δa/a0/% | 0 | 0 | −0.036 | −0.159 | −0.302 | −0.477 | −0.494 |
| V, unit cell, Å3 | 153.216 | 153.235 | 153.071 | 152.506 | 151.851 | 151.052 | 150.963 |
| Eform/eV | −0.444 −0.55 [12] −0.570 [36] | −0.438 | −0.425 | −0.407 | −0.383 | −0.358 | −0.337 |
| Evf/eV | 0.000 0.000 [36] | 0.013 | 0.038 | 0.075 | 0.123 | 0.172 0.15 [17] 0.29 f [7] 0.32 f [36] | 0.215 |
| Phase | c11 /GPa | c33 /GPa | c44 /GPa | c66 /GPa | c12 /GPa | c13 /GPa | B /GPa | G /GPa | B/G | E /GPa | ν | A |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ThC | 215.59 ± 0.60 | 79.72 ± 0.60 | 89.74 ± 0.43 | 131.61 | 72.49 | 1.82 | 183.84 | 0.27 | 1.22 | |||
| ThC [13] | 276.4 | 87.2 | 99.1 | 158.2 | 87.8 | 1.80 | 222.2 | 0.27 | 0.98 | |||
| ThC [15] | 222.49 | 80.41 | 92.03 | 135.52 | 74.34 | 1.82 | 188.54 | 0.27 | 1.13 | |||
| ThC [45] | 222.10 | 66.12 | 85.67 | 131.15 | 67.10 | 1.95 | 171.97 | 0.28 | 0.97 | |||
| ThC0.96875 | 214.73 ± 0.82 | 76.47 ± 0.58 | 84.15 ± 0.58 | 127.68 | 71.79 | 1.78 | 181.38 | 0.27 | 1.19 | |||
| ThC0.9375 | 213.53 ± 0.59 | 72.23 ± 0.41 | 78.24 ± 0.41 | 123.12 | 70.54 | 1.75 | 177.69 | 0.26 | 1.07 | |||
| ThC0.875 | 211.28 ± 0.83 | 66.14 ± 0.83 | 68.53 ± 0.59 | 116.47 | 68.63 | 1.70 | 172.08 | 0.25 | 0.90 | |||
| ThC0.8125 | 201.90 ± 0.78 | 59.08 ± 0.78 | 59.75 ± 0.55 | 107.20 | 63.93 | 1.68 | 159.99 | 0.25 | 0.83 | |||
| ThC0.75 | 180.51 ± 0.72 | 192.14 ± 0.72 | 54.62 ± 0.72 | 55.83 ± 0.72 | 60.90 ± 0.72 | 51.99 ± 0.72 | 98.08 | 58.40 | 1.68 | 146.18 | 0.25 | 0.91 |
| ThC0.6875 | 173.06 ± 0.63 | 180.81 ± 0.63 | 51.46 ± 0.45 | 51.95 ± 0.63 | 51.46 ± 0.45 | 46.95 ± 0.45 | 90.87 | 55.95 | 1.62 | 139.26 | 0.24 | 0.84 |