**2. Computational Method**

In this work, MD simulations are performed using the LAMMPS (Large scale Atomic/ Molecular Massively Parallel Simulator) code (Sandia National Laboratories, Livermore, CA, USA) [15], which is a classical molecular dynamics (MD) code used to simulate the

atomic interaction of selected materials. The interatomic potential developed by Cooper, Rushton, and Grimes (CRG) is used in this work for U-U, U-O and O-O interactions [4], which has been proven to reliably predict the various thermo-mechanical properties within the temperature range of 300 K to 3000 K [8,16–18]. In order to accurately describe the properties of UO2, the Coulomb interaction is further included with the original pair. The computational box used in this work is a 10 × 10 × 12 extension of fluorite (CaF2) unit cells containing 14,400 atoms. The lattice parameter for the computational box is the equilibrated lattice parameter at the investigated temperature. The periodic boundary condition is applied in all directions.

Generally, the O/U ratio in all defects after a displacement cascade in UO2 is close to two [19], which is in agreement with the results presented by Devanathan et al. [20] and Van Brutzel et al. [21]. In order to create this structure of uranium dioxide with Frenkel defects and maintain the neutral charge of the system, uranium and oxygen atoms are removed from the system by keeping 1:2 ratio. The same amount of uranium and oxygen atoms are then randomly inserted at the octahedral interstitial positions of the face-centered cubic (fcc) cation sublattice [22]. Different from Frenkel defects, oxygen-antisites are created by substituting O atoms with U atoms. Similarly, uranium-antisites are created by substituting U atoms with O atoms. In order to maintain the stoichiometry of the defected system, the number of O-antisites is equal to that of U-antisites. In order to investigate the effect of defect concentration, UO2 structures with 1%, 2%, and 5% Frenkel defects and 1%, 2%, and 5% antisite defects are built for further simulations. It should be noted that in the present work the concentration is defined as the value before MD relaxation at given temperatures. The main reason is that the relaxations at different temperatures could result in different concentration values after full relaxation. In order to avoid this misunderstanding during the investigation of concentration effects in this work, the concentration value before MD relaxation is used. For each defect concentration the statistical results are made based on 3 samples by randomly creating Frenkel or antisite defects.

The simulation box is first relaxed by the conjugated gradient (CG) method at 0 K and further relaxed for 500 ps at the temperature of interest under NPT (constant number, pressure, and temperature) ensemble with zero external pressure. It should also be noted the system containing 5% defects requires longer equilibration time (600 ps) to reach the equilibrium state. The equilibration is checked by the dependence of total energy and volume of the system on simulation time, which both indicate that the system has reached the equilibrium state before further calculations of lattice constants, thermal expansion coefficient, and elastic modulus. Therefore, the results obtained in the present work are reliable and could provide useful information to understand the properties of UO2 after irradiation. The timestep of 1 fs is used for all simulation processes.

The lattice parameter as a function of temperature for different concentrations of Frenkel defects and antisite defects is first calculated. The thermal expansion coefficient (α) is then determined from the first derivate of the lattice parameter with respect to the temperature using the following Equation (1):

$$\mathbf{u}(\mathbf{T}) = \frac{1}{L} (\frac{\partial L}{\partial T})\_p \tag{1}$$

where *L* is lattice parameter and ( *<sup>∂</sup><sup>L</sup> <sup>∂</sup><sup>T</sup>* )*<sup>p</sup>* is the slope of the plot for the lattice parameter as a function of temperature at the given temperature [23]. It should be noted that, in the present work, the structure of UO2 is FCC and thus the thermal expansion coefficient is isotropic.

As the structure of uranium dioxide is cubic, three independent elastic constants (*C*11, *C*12, and *C*44) need to be calculated. These elastic constants can be calculated by applying elementary strain in six directions and measuring the changes in the six stress components. In this work, the strain to induce the deformation of the simulation box was set to be 10−5. Based on the dependence of the stress on the strain, these constants are calculated as described in [24]. Based on these three constants, the bulk modulus (*B*), shear modulus

(*G*), and Young's modulus (*Y*) can be calculated. The bulk modulus is calculated with the following equations [14].

$$B = (C\_{11} + 2C\_{12})/3\tag{2}$$

Furthermore, according to the Hill's suggestion [25], in order to determine the shear modulus, the shear modulus (*GV*) using the Voigt method and the shear modulus (*GR*) using the Reuss method need to be obtained, which can be determined using the following Equations (3) and (4), respectively.

$$G\_V = (C\_{11} - C\_{12} + \mathfrak{A}C\_{44})/5\tag{3}$$

$$\mathbf{G}\_R = (\mathbf{5(C\_{11} - C\_{12})C\_{44}})/(4\mathbf{C\_{44} + 3(C\_{11} - C\_{12})}) \tag{4}$$

Thus, the shear modulus (*G*) using Hill's method can be obtained as the arithmetic average of *GV* and *GR*.

$$G = (G\_V + G\_R)/2\tag{5}$$

Based on the calculated bulk modulus and shear modulus, Young's modulus is calculated with the following Equation (6).

$$
\mathcal{Y} = \mathcal{B}B\mathcal{G} / (\mathcal{B}B + \mathcal{G}) \tag{6}
$$

#### **3. Results and Discussion**

### *3.1. Effect of Defects on Lattice Parameter and Thermal Expansion Coefficient*

The lattice parameter (*L*) of pure UO2 as the function of temperature is plotted in Figure 1. The error bars in the figure correspond to the standard deviation calculated among the five different statistical lattice constants calculated at the given temperature. For comparison, the results from the VASP calculation by Wang et al. [26] and from the experimental measurement by Taylor et al. [27], Yamashita et al. [28], and Momin et al. [29] are also included in Figure 1. It is also clear that *L* increases linearly with increases in temperature from 0 K to 1500 K as investigated in this work. From these results, the results from the present MD agree better with the experimental value than those from VASP calculations. Thus, the MD method and the related empirical potential could be used for the present purpose for further simulations.

**Figure 1.** Dependence of lattice constant on temperature including the results from the present work, the experiment studies by Taylor et al. [27], Momin et al. [29], Yamashita et al. [28], and the VASP results by Wang et al. [26].

When defects are formed after irradiation, the effects of radiation defects are also simulated in this work. In Figure 2, the dependence of the lattice parameter of UO2 on defect concentration at different temperatures is provided for Frenkel defects (dash) and antisite defects (solid). From Figure 2, it is clear that the lattice parameter of the system increases with an increase in Frenkel defect concentration from 0 to 5% at all investigated temperatures in the present work, although the increases are limited around 1.0%.

**Figure 2.** Dependence of lattice constant of UO2 on defect concentration at different temperatures for Frenkel defects (dash) and antisite defects (solid).

The present results indicate that the formation of Frenkel pairs or antisite defects could increase the lattice constant quickly when the defect concentration is less than 2.0%, above which the lattice constant is almost a constant or increases or decreases slightly. Thus, the results are different from previous studies about the effect of porosity on the lattice parameter *L* of ThO2, which reported that the *L* of the system increased linearly with increases in porosity concentration. The reason for the above difference may be mainly from the property of anisotropic effects induced by antisite defects and interstitials in Frenkel pairs, which is different from the isotropic vacancy or porosity. It should also be noted here that in this work, only defect concentrations up to 5% are considered. If higher concentrations were considered, the lattice constant may change accordingly.

As stated previously, because of the high temperature within the fuel due to the fission reaction, thermal expansion coefficient is considered to be an important factor for modelling and predicting the nuclear fuel's behavior [13]. Figure 3 presents the effects of Frenkel defects and antisite defects on the thermal expansion coefficient of UO2 as a function of temperature. The uncertainty of the thermal expansion coefficient at different temperatures has also been calculated by changing the defect distribution but keeping the same concentration as stated in the computational method section. As shown in Figure 3, the uncertainty is limited for the three cases investigated in the present work. For comparison, the experimental results [27], MD derived values [30], and first principles data [31] for pure UO2 are included in the figure.

**Figure 3.** Thermal expansion coefficient for pure UO2 and for defective UO2 with different percentages of Frenkel defect (**a**), and antisites (**b**), as a function of temperature. For comparison, the experiment studies by Taylor et al. [27], the MD results by Cooper et al. [30] and the first-principles data by Yun et al. [31] are included in the figure.

Figure 3a shows that compared to the perfect system, the thermal expansion coefficient of systems containing 1% Frenkel defects has similar results to that of perfect systems at temperatures from 600 K to 1200 K, while at 1500 K the 1% Frenkel defects induce larger thermal expansion coefficients with an increase around 4.0%. When the system contains a higher concentration of Frenkel defects, e.g., 2.0% and 5.0%, the lower thermal expansion coefficients are observed. For example, when the concentration of Frenkel defects is 2%, the maximum reduction of thermal expansion coefficient up to 15% is observed at 1500 K. When Frenkel defect concentration is 5%, the maximum reduction is around 30% within 600–1500 K. Thus, it is found that the reduction of the thermal expansion coefficient increases with the increase in the concentration of Frenkel defects. Furthermore, it can also be observed from Figure 3a that when the defect concentration is low, e.g., 1%, the thermal expansion coefficient increases with an increase in temperature, while the higher concentration defects (2% and 5%) result in an increase in the thermal expansion coefficient for systems from 600 K to 1200 K, but decrease from 1200 K to 1500 K. In fact, according to Sun et al. [32], the thermal expansion coefficient can be described as a function of both temperature (T) and atomic restoring force of the system (F(r)), especially around the defect region. The F(r) is a function of atomic distance, which depends on the temperature and defect distribution in the computational box. Therefore, when the temperature changes, the change of F(r) together with temperature results in non-monotone temperature dependence of the thermal expansion coefficient.

Figure 3b presents the effect of antisite defects on the thermal expansion coefficient of UO2. Different from the effects of Frenkel defects, Figure 3b clearly indicates that for the concentrations of antisites investigated in the present work, the thermal expansion coefficient of the system increases with an increase in temperature. When the temperature is 600 K, the thermal expansion coefficient decreases with an increase in antisite defect concentration. When the temperature is 900 K or 1200 K, the thermal expansion coefficient has similar values for systems containing antisite defects less than 2%, but decreases around 25–30% when the antisite defect concentration is 5%. When temperature is 1500 K, the thermal expansion coefficient has the same value for systems containing 1% to 5% antisite defects, which is lower than that of the perfect system. Comparing the results shown in Figure 3a,b, it could be found that antisite defects have stronger effects than Frenkel defects on the thermal expansion coefficient of UO2.
