Molecular Dynamics Simulations of Xe Behaviors at the Grain Boundary in UO₂

Abstract

In this study, we investigated the behavior of xenon (Xe) bubbles in uranium dioxide (UO₂) grain boundaries using molecular dynamics simulations and compared it to that in the UO₂ bulk. The results show that the formation energy of Xe clusters at the Σ5 grain boundaries (GBs) is much lower than in the bulk. The diffusion activation energy of a single interstitial Xe atom at the GBs was approximately 1 eV lower than that in the bulk. Furthermore, the nucleation and growth of Xe bubbles in the Σ5 GBs at 1000 and 2000 K were simulated. The volume and pressure of bubbles with different numbers of Xe atoms were simulated. The bubble pressure dropped with increasing temperature at low Xe concentrations, whereas the volume increased. The radial distribution function was computed to explore the configuration evolution of Xe bubbles. The bubble structures in the GB and bulk material at the same temperature were also compared. Xe atoms were more regular in the bulk, whereas multiple Xe atoms formed a planar structure at the GBs.

1. Introduction

Uranium dioxide (UO₂) has been widely used as fuel for nuclear reactors owing to its properties, such as there being no specific deformation when it is strongly irradiated, having an unchanged lattice structure at high temperatures, being non-volatile and being chemically unreactive with water [1]. During the operation of nuclear reactors, the nuclear fuel elements are subjected to a harsh working environment and numerous radioactive fission products are produced during the fission of the fuel assembly in the reactor core. With the development of burnup, solid and gas fission products are produced in the fuel elements and their volume is greater than that of the material before fission. The volume of the fuel element increases with the development of burnup, which is called irradiation swelling. The radiation swelling of nuclear fuel induces interactions between fuel pellets and cladding, resulting in radial deformation and transverse tension of the cladding tube, causing damage to the cladding tube, which seriously threatens the safe operation of the reactor [2]. Swelling caused by solid fission products is simple and increases linearly with burnup; the behavior of gas fission products is complex and this field has not been extensively studied.

Due to the numerous radiation and structural changes experienced during the life of UO₂ and storing, understanding and controlling the microstructural changes requires a comprehensive approach that considers all aspects of the material’s behavior, from basic radiation damage processes to longer-term changes in the material microstructure. To better understand the fission gas behavior, such as microstructural changes and swelling, many experimental and theoretical studies have been conducted [3,4,5,6,7,8]. It is difficult to analyze the behavior of xenon (Xe) atoms in UO₂ through experiments [9,10,11,12] because the growth mechanism of Xe bubbles and the entire physical process are not well understood. Therefore, intragranular rare bubbles under dynamic conditions, cannot be studied through experiments, making computational simulations the only choice. Niemiec et al. proposed a basic evolution equation describing the kinetics of the nucleation and growth phase transitions to study phase transformation or microstructure formation kinetics in physical systems originally composed of several grains [13]. Atomistic simulations are vital to a provide good insight into the atomic structure and damage mechanisms.

To date, all simulations have been conducted in monocrystals. Gadomski et al. [14] studied the kinetic anomalies occurring in nucleation and growth phenomena in complex systems, such as polycrystalline partly ordered alloys, quasicrystal line assemblies, and mesomorphs, to understand the kinetics of the evolution of the microstructure and the system during growth. UO₂ pellets are manufactured through traditional powder metallurgical processes; hence, they are polycrystalline materials composed of particles with a diameter of approximately 10 μm. The opening at the grain boundaries (GBs) is a crucial phenomenon since the energy of defects decreases near the GB and fission gas, which causes aggregation at the GB, and can be released to the outside of the fuel. As a result, the properties of UO₂ GBs at the atomic scale have been studied [15,16,17,18,19]. In this study, we compare the early behavior of Xe bubbles in the GBs and the entire block of UO₂. Previous numerical studies on this subject have mainly focused on metal systems, such as metallic tungsten [20,21,22]. Herein, we focus on the behavior of Xe atoms at the GBs and the bulk UO₂, including the migration of a single Xe interstitial atom, the nucleation and growth of Xe bubbles and the bubble pressure and expansion associated with bubble growth. This study serves as a good reference for higher length scale simulation models.

2. Simulation Method

2.1. Interatomic Potential

Molecular dynamics (MD) simulation is an effective means of studying the microworld, providing insight into the atomic structure and the mechanism of bubble growth. An MD simulation program, LAMMPS [23], was employed for all simulations in this study. However, the accuracy of MD simulations depends largely on the potential function of the atoms used to support the simulation.

Yang et al. studied the influence of different potential functions on the growth of high-pressure Xe bubbles [3]. With the continuous addition of Xe atoms, the pressure of the bubbles is divided into two stages. In stage I, where the bubble pressure is monotonically increased, the bubble characteristics and the microstructure evolution of UO₂ are relatively independent with the interatomic potential used; thus, the results of the five potential functions are approximately the same. In stage II, the bubble pressure is released and then fluctuates and the volume and pressure of the bubbles, as well as the evolution of the bubble configuration and the UO₂ matrix as a function of the number of Xe atoms, highly depend on the potentials used [3]. In addition, the formation energy of Xe at the same location in UO₂ varies significantly with the potential of UO₂ and Xe–UO₂ [3,24,25].

To ensure the accuracy of MD simulations herein, we evaluated different sets of potential functions and compared them with density functional theory (DFT) or DFT + U data to select a better comprehensive potential function for adding Xe to UO₂ [26,27,28,29,30].

The first and second sets of potentials have the same UO₂ potential and were developed by Basak et al. [31]. However, IPR [27] and Geng [32] potentials were used to interact with Xe–UO₂. In the third and fourth sets of potentials, the potential reported by Morelon [33] was used for the UO₂ matrix, whereas those reported by Chartier [29] and IPR [27] were used to describe the Xe–UO₂ interaction. In the fifth and sixth sets of potentials, the CRG potential [25] was used as a potential function of UO₂ and the interaction of Xe–UO₂ was described by Cooper [34] and IPR [27]. Therefore, the six sets of potentials are called Basak/IPR, Basak/Geng, Morelon/IPR, Morelon/Chartier, CRG/Cooper and CRG/IPR, respectively. Among them, three potential functions of UO₂, are paired with different potential functions of Xe–UO_2, because the interactions between Xe–U and Xe–O are determined in Xe–UO₂. Thus, even in the case of the same UO₂ potential function, we simulated it with different Xe–UO₂ potential functions and the results are different.

We can determine the accuracy of using interatomic potentials by comparing them with the energetics of relatively small defects calculated using MD [3]. We developed several physically relevant benchmarks to compare the six sets of potentials on an equal basis. We placed Xe at different sites in UO₂ (small, intermediate and large vacant sites). For small incorporation, we used the following defects: Xe in an interstitial (Int_Xe) and Xe on a uranium site with uranium in the neighboring interstitial (Sub_Xe-Int_U). The formation energies of the intermediate-sized incorporation sites are a Xe atom located at a U vacancy (Xe in Vu), an O vacancy (Xe in Vo), two nearest U and O vacancies (Xe in Vuo), or the three configurations of Schottky defect clusters (Xe in SD). Finally, we calculated the defect energetics of a single large incorporation site with Xe in a double Schottky defect cluster (Xe in 2SDs). To directly compare the MD values, a 2 × 2 × 2 supercell of the cubic fluorite unit cell (96 atoms for stoichiometric UO₂) with a lattice constant of 5.454 Å was used in the molecular statics calculations. The simulation results are compared with the data in

Table 1. Comparison of energy of xenon (Xe) point defects in uranium dioxide (UO₂) bulk calculated using various interatomic potentials at 0 K. Charge corrections for charge defects were not considered in the calculations herein.

	Basak/Geng	Morelon/Chartier	Basak/IPR	Morelon/IPR	CRG/IPR	CRG/Cooper	DFT
${\mathrm{E}}_{{\mathrm{Int}}_{\mathrm{Xe}}}^{\mathrm{f}}$	22.73	11.94	9.68	8.32	9.79	11.70	9.73 [27] 9.71 [28] 9.7–12.0 [26]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} \mathrm{Vu}}^{\mathrm{f}}$	5.76	5.31	1.84	0.70	1.72	3.60	2.0–5.8 [26] 1.95 [28]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} \mathrm{Vo}}^{\mathrm{f}}$	20.55	9.31	8.41	7.17	7.96	8.26	7.5–9.1 [26] 7.85 [28]
${\mathrm{E}}_{{\mathrm{Xe} \mathrm{in} \mathrm{V}}_{\mathrm{uo}}}^{\mathrm{f}}$	5.32	4.98	1.52	0.61	1.56	2.96	1.6 [26] 1.55 [28]
${\mathrm{E}}_{{\mathrm{Sub}}_{\mathrm{Xe}}-{\mathrm{Int}}_{\mathrm{U}}}^{\mathrm{f}}$	20.77	Unstable	11.89	10.76	10.29	12.27	11.33 [29]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} \mathrm{SD}}^{\mathrm{f}}$(1)	4.82	4.33	0.91	0.48	0.95	1.99	1.06 [27] 1.2 [29]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} \mathrm{SD} }^{\mathrm{f}}$(2)	4.92	4.77	1.53	0.67	1.68	3.02	1.83 [27] 1.8 [29]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} \mathrm{SD}}^{\mathrm{f}}$(3)	5.02	4.70	1.81	0.65	1.89	3.45	1.94 [27] 2.3 [29]
${\mathrm{E}}_{\mathrm{Xe} \mathrm{in} 2\mathrm{SDs}}^{\mathrm{f}}$	3.46	2.79	0.19	0.16	0.18	0.59	0.27 [27]
${\mathrm{E}}_{{\mathrm{Int}}_{\mathrm{Xe}}}^{\mathrm{m}}$	0.76	0.84	5.54	3.28	4.32	1.91	4.48 [30]

[3]. The difference in energy is small, the maximum error value is 0.28 eV and the difference is 0.1 eV. The energy difference may be attributed to the difference in the machines used. However, in general, the results show that the potentials used are adequate. Then, through the MD simulations, we used a larger supercell (12,000 atoms) for better convergence to calculate the energies of Xe atoms at different sites in the UO₂ lattice. In addition to the formation energies, we calculated the migration energy of interstitial Xe (${\mathrm{E}}_{{\mathrm{Int}}_{\mathrm{Xe}}}^{\mathrm{m}}$). The results are compared with selected DFT and DFT + U data in Table 1. The results obtained using the CRG/IPR potentials show the best consistency with the DFT results; thus, we chose CRG/IPR as the potential function for adding Xe to UO₂ for our subsequent simulations.

2.2. Model Building

The structure of GBs is essential for predicting the behavior of Xe bubbles. Up to now, most theoretical studies on the properties of GB have been based on molecular statics (MS) or DFT [35,36] without considering the influence of temperature. However, experimental observations have shown that GB microstructures and corresponding properties reflect the effect of temperature on GBs to a certain extent [37,38,39] because temperature changes may cause the relative sliding of grains at GBs [17,19]. Therefore, heat treatment of GBs is a necessary process to expand the understanding of GB. Hong et al. [40] reported that high temperature can cause complex ion transitions in some GBs. In this study, the symmetrically tilted GB, Σ5(310)/[100], in UO₂ was used as a GB model because the Σ5 symmetric tilt GB has been extensively studied both experimentally and theoretically [17,41,42,43]; thus, we could compare our results with the available data.

In this study, the symmetrically tilted GB Σ5(310)/[100] in UO₂ were formed by rotating perfect UO₂ crystals. Symmetrically tilted GB means that the crystals on both sides of the GB are tilted towards each other and the interface of the GB is symmetrical to the two grains. The specific process of constructing the GB is as follows. First, half of the perfect UO₂ crystal simulation box is rotated along the axis by half of the orientation difference angle. Then, the second crystal is constructed through mirror inversion, taking the first crystal as a reference. However, this method causes an overlap of the atoms at the crystal interface. Thus, we removed the overlapped interface atoms to maintain a neutral charge.

The annealed structures of the Σ5(310) GB based on MD are shown in Figure 1 and a common neighbor analysis was performed to identify the GB region and bulk region. The atomic configuration was viewed using OVITO software [44]. At 1000 K, the structure of the interface is conserved with the repetition of a triangular-like pattern (depicted by lines in Figure 1); however, at 2000 K, a structural change occurs with the boundary showing reduced diamond shapes arranged end-on-end. Previous simulation results show that, with an increase in temperature, the disorder degree of GB increases [45]; e.g., at 2000 K, the Σ5(310) GB has a more distorted triangular pattern [19].

2.3. MD Simulation Setup

For all MD simulations, the time step was set as 0.001 ps. For static relaxation, energy minimization was performed and the minimization algorithm was set as the conjugate gradient method (cg), with the stopping criteria for energy tolerance of 10⁻²⁵ s⁻¹ and a force tolerance of 10⁻²⁵ eV/Å.

We studied the characteristics of Xe atoms in UO₂ GB, including the formation energy of a single Xe atom in the GB, the migration of a single Xe atom in UO₂ GB and the nucleation of Xe bubbles, and compared them with those in the bulk. The specific simulation processes for various characteristics are different. Thus we constructed two simulation boxes of different sizes: one contains 38,290 atoms, with a size of 10 a₀ × 20 a₀ × 10 a₀ (a₀ is the lattice constant of UO₂) and was used to calculate the migration of a single Xe interstitial atom at the GB, and periodic boundary conditions were used in three directions of the box; the second is a 25 a₀ × 40 a₀ × 25 a₀ large system simulation box and was used to simulate the nucleation and growth of Xe bubbles and the energy of Xe atoms at the GB, and periodic boundary conditions were used in three directions of the box. Similarly, in the bulk, a cubic box of 10 a₀ × 10 a₀ × 10 a₀ containing 12,000 atoms was used with periodic boundary conditions to calculate the migration of a single Xe interstitial atom. The other system is 25 a₀ × 25 a₀ × 25 a₀ with 187,500 atoms and was used to simulate the nucleation and growth of Xe bubbles and the energy of Xe atoms. We performed the same simulation process at the GB and in the bulk to better compare the differences in the behavior of Xe bubbles. A cutoff distance of 1.1 nm was adopted and used for all potentials used in this study.

To determine the stable configuration of Xe atoms in the GB, the stable structures of Xe clusters with different sizes (the number of atoms less than six) were obtained by adding Xe atoms gradually and performing MS simulations for each configuration.

After determining the stable Xe site in the GB, the diffusion mechanism of a single Xe interstitial atom at the GB was studied. The migration of a single Xe atom at the GB was studied using MD and the nudged elastic band (NEB) methods [46,47,48]. Then, the mean square displacement (MSD) method was used to calculate the diffusion coefficient and diffusion energy barrier of a single Xe atom at the GB.

Finally, to study the evolution mechanism of Xe bubbles, the microcanonical NVT (constant volume and temperature) ensemble was used to simulate Xe bubbles in UO₂ GB. To describe the bubble characteristics during bubble evolution, the volume V and pressure P of the bubble after inserting each Xe atom were calculated by MD. The number of Xe atoms in the bubble is too small (only 50 Xe atoms) to obtain a good statistical result. V is the sum of the volumes of all Xe atoms, which were calculated using the Voronoi technology implemented in LAMMPS [49] and the pressure of the Xe bubble was calculated from the sum of the diagonal component of the atomic stress tensor and the bubble volumes as follows [8,50]:

$$P=-\frac{1}{3V}{\displaystyle \sum }_i^n\left[S_{11}\left(i\right)+S_{22}\left(i\right)+S_{33}\left(i\right)\right]$$

(1)

where n is the number of Xe atoms in the bubble, $S_{\alpha \alpha }\left(i\right)$ is the diagonal component of the stress tensor for atom i, V is the volume of all Xe atoms.

Before the first Xe atom was randomly introduced into the simulation box, energy minimization was performed and the minimization algorithm was set using the conjugate method (cg), with the stopping criteria for energy tolerance and force tolerance of 10⁻²⁵. Then, the first Xe atom was randomly added to the system. The initial system with Xe atoms was minimized to avoid long-distance movements of inadvertently overlapping atoms and the temperatures were set to 1000 K and 2000 K, respectively. The velocities of atoms were set and an additional 10 ps simulation was performed using microcanonical NVT (constant volume and temperature). Following this thermal equilibration (i.e., no change in the various properties of the system over time), Xe atoms were sequentially inserted one by one every 10 ps into the Xe bubble center of mass until 50 Xe atoms were contained in the bubble. After inserting each Xe atom, to avoid introducing artificial energy into the system resulting from atoms placed too close together, the system energy was minimized using the conjugate gradient (cg) and then simulated under NVT conditions for 10 ps at 1000 K and 2000 K.

3. Results

3.1. Formation Energy and Diffusion Behavior of a Single Xe Atom at the GB

The interaction between Xe atoms and the microstructure of UO₂ fuel is key in fission-gas release. To simulate the redistribution of fission gas atoms in the UO₂ microstructure, the interaction range between fission gas atoms and GB must be determined first, to determine the stable structure of migrated Xe atoms at the GB. We calculated the formation energy of a single Xe atom at different positions from the Σ5(310) GB.

There are two distinct regions of a Xe interstitial atom near the GB, as indicated by E in Figure 2. The first region is located at a distance of more than 1.5 a₀ for interstitial Xe, in which the energy is almost equal to the calculated formation energy of Xe interstitial atoms in the bulk (9.89 eV), indicating that the driving force for the segregation of Xe atoms to GB can be ignored. In the second region, where the distance is within 1.5 a₀ from the GB, E decreases, having a maximum value of approximately 6.88 eV, which indicates that interstitial Xe can be absorbed by GB. Thus, the interaction range between Xe interstitial atom and GB is approximately 1.5 a₀ from the plane of GB. Figure 2 also shows that the lowest energy point for interstitial Xe is located at the GB. Thus, it is energetically favourable for interstitial Xe to be absorbed by the GB region.

The behavior of fission gas in nuclear fuel is the main factor that determines the change in radiation swelling with fuel consumption. Therefore, it is necessary to study the migration of fission gas in UO₂ GB.

Herein, we calculated the migration of an interstitial Xe atom in UO₂ GB using NEB, i.e., a direct hopping mechanism from one octahedral site to another. NEB showed an energy barrier of 3.87 eV. Thus, whether migration is along or perpendicular to the GB direction, the calculated energy barrier is high. These results suggest that the migration of Xe between the two interstitial sites is not the main mechanism due to the higher energy barrier [51]. The relatively high energy barrier is mainly attributed to the lattice deformation caused by Xe movement between two adjacent interstitial sites, indicating that the octahedral interstitial position in UO₂ is highly strained and, therefore, energetically unfavourable [52,53].

In addition to NEB, we employed the MSD method to calculate the diffusion coefficients and diffusion energy barriers of a single Xe atom at the GB. A long simulation time and a short distance between the position were employed to obtain information on the atomic position. The results were obtained by dividing the tracks of the MSD and averaging them several times. As shown in Figure 3, we calculated the diffusion coefficient and diffusion energy barrier separately at the GB and in the bulk. The diffusion energy barrier of the Xe atoms at the GB was 1.40 eV and that in the bulk was 2.46 eV. The magnitude of the preexponential factor was used to describe the ease of diffusion of a Xe atom. We conclude that Xe atom diffuse more easily at the GB than in the bulk.

The diffusion barrier calculated using MSD and the migration barrier calculated using NEB are different. When NEB was used, the single Xe atom showed a high energy barrier in the interstitial site, which is consistent with that of the bulk. However, the calculation results with MSD are different. This is because Xe may not diffuse directionally in UO₂. The temperature accelerated dynamics simulations revealed that the dynamics of defect clusters strongly depends on their size and the diffusion direction is not one-dimensional [54].

3.2. Formation Energies of Small Xe Clusters at the GB

We calculated the formation energies of Xe atoms and atomic clusters at the GB and compared them with those in the bulk. The formation energies of Xe atoms at the GB are shown in

Table 2. Formation energies of Xe atoms at the GB and in the bulk UO₂.

N	1	2	3	4	5	6
$E_{\Sigma 5}^f$(eV)	1.3	3.8	6.34	10.91	13.41	15.59
$E_{Bulk}^f$(eV)	9.79	16.44	22.49	27.92	32.19	39.79

. The formation energy of the Xe bubble is defined as:

$$E^f\left(n\right)=E_{bubble}\left(n\right)-E_{bulk}-n{E}_{Xe}^{atom}$$

(2)

where $E_{bubble}\left(n\right)$ is the total energy of the simulation supercell containing n Xe atoms in the bubble, $E_{bulk}$ is the total energy of the supercell without bubbles and $E_{Xe}^{atom}$ is the isolated Xe atom energy, which was set to zero. As shown in Table 2, the energy of a single Xe interstitial atom at the GB was 1.3 eV, whereas that in the bulk it was 9.79 eV. The total energy difference (ΔE) was 8.49 eV. As the number of Xe atoms increases, ΔE increases. Whether at the GB or the bulk, the formation energy increases with an increase in the number of Xe atoms. At all stages, the formation energy of Xe clusters at the GB is smaller than that in the bulk, indicating that Xe bubbles are easier to form at the GB than in the bulk. This is because Xe atoms are at the GB plane, which improves the energy compared to that in the bulk. Similar to previous results [8], at smaller bubble sizes (1–5 Xe atoms), the Xe bubble nucleus at the GB has a much lower formation energy than that of bubbles in the bulk with a similar size.

3.3. Nucleation of Xe Bubbles at the GB

Early-stage nucleation and growth of Xe bubbles were simulated using the MD method. We simulated the nucleation and growth of Xe bubbles at 1000 and 2000 K. When the system reached a thermodynamic equilibrium after relaxing the system at a given temperature and volume, Xe atoms were regularly added to the bubble. For each Xe atom added, the system needed 10 ps to relax, to ensure there was no change in the properties of the system over time, i.e., the equilibrium state was reached. Thus, the Xe bubble growth was simulated until the bubbles could contain 50 Xe atoms.

Figure 4 lists the differences in the nucleation configurations of Xe bubbles at the Σ5(310) GB and the bulk. At 2000 K, Xe atoms were more regular in the bulk (Figure 4b). This is consistent with a previous report [3], which showed that Xe bubbles evolve into a glassy/amorphous state, a nearly face-centred cubic (FCC) solid structure and subsequently to a high density amorphous or glassy state. At the GB, Xe atoms are absorbed by the GB and multiple Xe atoms form a planar structure, which is similar to the case of rhenium (Re) atoms in metals [55]. The results are consistent with the figure showing the formation energies of an interstitial Xe at the GB. However, comparing the configurations at 1000 K, as shown in Figure 4a, we found that at a smaller bubble size (20 Xe atoms), Xe atoms at the GB remain at the GB plane. However, with the growth of Xe bubbles, most Xe atoms in the larger bubbles are located in the bulk-like region away from the GB plane. This is attributed to the increase in temperature, which causes a change in the lattice constant and the reconstruction of the GB. An increase in temperature decreases the interfacial energy and it is easy to form segregation at the GB, where the energy of Xe atoms in the bulk is higher than that of the atoms at the GB. Thus, Xe atoms spontaneously converge toward the GB.

Next, to analyze the characteristics of Xe bubbles, the volume (V) and pressure (P) of the bubbles after inserting each Xe atom were calculated. Figure 5 shows the relationship between the volume and pressure of Xe bubbles and the number of Xe atoms in the bubble. The inset of Figure 5a shows that the volume of Xe bubbles increases with an increase in temperature and, during bubble growth, the volume of Xe bubbles increases approximately linearly with an increase in the number of Xe atoms. However, at the same temperature, the volume of Xe bubbles at the GB is higher than that of the bubbles in the bulk. This is because GBs are surface defects in solid materials and there are many defects, such as vacancies, at the GBs. Thus, under the same conditions, the volume of Xe bubbles at the GB is larger than that of the bubbles in the bulk. Also, it reflects the influence of temperature on the volume of Xe bubbles, which is attributed to the increase in temperature, which changes the lattice constant and the reconstruction of GBs [19]. At 2000 K, the Σ5(310) GB consists of a more distorted triangular pattern, the middle gap is larger and the lattice constant increases. In Figure 5b, at the GB, the pressure of the Xe bubble is initially quite high and then drops with increasing Xe concentration until the number of Xe atoms reaches 10. When there are more Xe atoms in the bubbles, the pressure increases with an increase in the number of Xe atoms, which is consistent with that of Xe in the UO₂ bulk [4,8]. Then, with an increase in bubble growth, the bubble pressure generally decreases with some fluctuations involving several peaks and drops. This is because the bubbles may be far from the equilibrium, from which they grow because of rapid changes in temperature or insufficient local lattice vacancies [56]; thus, the smallest bubbles have a very high pressure, which results in a density comparable to that of solid Xe [57]. For larger bubbles, the density and pressure are relatively low. Experiments have shown that at the very early stages of bubble development, the bubble density is high with little evidence of deviation from circular bubbles. However, extensive bubble coalescence can greatly reduce the density of bubbles [58]. Under the same conditions, the pressure of Xe bubbles at the GB is less than that of the bubbles in the bulk. The internal pressure of Xe bubbles depends on not only the external stress on UO₂ but also the surface tension of the bubble voids [56].

In addition to volume and pressure, we evaluated the configuration of Xe bubbles at 1000 and 2000 K. Snapshots of the radial distribution function of Xe atoms at 1000 and 2000 K are shown in Figure 6. To avoid thermal fluctuations at 1000 and 2000 K, the positions of Xe atoms were obtained by taking an average total time of 50 ps before inserting the next Xe atom into the bubble. As shown in Figure 6, the peak is smaller and fuzzy at high temperatures because the atomic amplitude is large at high temperatures; thus, it is easy to deviate from the equilibrium position and approach the liquid state. With an increase in Xe atoms, Xe in the bubble evolves into a glassy/amorphous state, indicating that the behavior of Xe atoms in UO₂ is similar to that of a system with a hard spherical liquid in a solid. This observation is consistent with the results of Geng et al. [32]. For smaller bubbles (1–5 Xe atoms), the distribution of Xe atoms at the GB is more dispersed than that in the bulk, which also shows that Xe atoms at the GB diffuse more easily than those in the bulk, according to the MSD results. In the bulk, there is only one peak of Xe atoms at the beginning, indicating the uniform distribution of Xe atoms in the UO₂ matrix. With an increase in Xe atoms, the first peak gradually grows and the second peak appears, indicating the formation of Xe clusters and the emergence of the second layer of atoms in the clusters.

4. Conclusions

In this study, MD and static simulations were employed to investigate the energy and diffusion of a single Xe atom and small Xe clusters and the nucleation and growth of Xe bubbles at the UO₂ Σ5(310) GB. The results show that the formation energies of single Xe atoms and small Xe clusters at the GB are much lower than those in the bulk, indicating that impure Xe atoms are more stable at the GB and can precipitate into clusters. The diffusion activation energy of a single interstitial Xe atom at the GB is approximately 1 eV lower than that in the bulk. The energy barrier for interstitial Xe atoms to diffuse through the interstitial mechanism is high, which is mainly due to lattice deformation caused by the movement of Xe between two adjacent interstitial sites. The growth of Xe bubbles at the GB was simulated by sequentially inserting Xe atoms into a pre-existing Xe bubble. Then, the bubble volume and pressure were estimated for different numbers of Xe atoms. The pressure dropped with an increase in temperature, whereas the volume increased. At low Xe concentrations, the Xe bubble volume increased approximately linearly with an increase in the number of Xe atoms. We also compared the configurations of Xe atoms at the GB and in the bulk at the same temperature and found that Xe atoms are more regular in the bulk, whereas, at the GB, multiple Xe atoms form a planar structure.