Deterministic Phonon Transport Predictions of Thermal Conductivity in Uranium Dioxide With Xenon Impurities

One of the fundamental quantities of interest in the safe and efficient operation of nuclear reactors is thermal conductivity $(κ)$ of the fuel, which greatly influences heat transfer throughout the structure of the nuclear core and into the coolant [1,2]. In addition to heat transfer, thermal conductivity is coupled to many other processes in the reactor core. Shifting thermal gradients have a strong influence on the macroscopic cross sections of interaction for neutrons [1]. These cross sections are of high importance since they dictate rates of nuclear fission, neutron absorption, and scattering within the fuel. As temperatures increase, so do the effects of Doppler broadening, which alter neutron scattering and absorption. While this is not a new phenomenon, reactor operators must be keenly aware of the effect temperature has on the absorption and scattering behavior of neutrons. The focus of our work is to develop a predictive computational tool which simulates thermal transport in heterogeneous nuclear fuel in service and under irradiation, with fission product defects.

Thermal conductivity in nuclear fuel is currently obtained through empirical relationships which have been experimentally determined from measurements made during the past 60–70 years [3]. Thermal resistance measurements are performed on nuclear fuel with operating histories, i.e., irradiated fuel and values are obtained at specific temperatures and isotope concentrations. This approach does not consider the constantly changing concentration of isotopic byproducts in the fuel, nor has it historically provided appropriate values across a wide range of reactor operating conditions without significant interpolation error [3,4]. Reactor designers and operators, however, rely on interpolation to fill in the gaps which can be a significant source of uncertainty in predictions of reactor performance. Nuclear reactors are constructed conscious to this attribute and thus have a significant conservatism and safety margin [2].

A better predictive approach to the computation of thermal conductivity could reduce these margins, creating improved performance and economics without compromising safety. A predictive simulation tool could also reduce the reliance on experiment for the development of new fuels for advanced reactors.

Development of nuclear reactors is an ongoing process. The current generation of power reactors is light water moderated and use uranium dioxide $(UO2)$ fuel. In the future, generation IV nuclear reactors will continue to use solid fuel. Uranium-based tristructural-isotropic particles are used in prismatic block high temperature gas reactors, while other fuels such as uranium-molybdenum and uranium-carbide are in development for a plate or pellet-based application [5,6]. Experimental measurements of thermal conductivity will likely be performed on these new fuels, and could result in more empirical correlations used to predict their performance under irradiation. The need for a predictive computational tool to reduce the reliance on destructive thermal conductivity measurements is steadily increasing.

We are developing a deterministic computational framework for simulating phonon transport. When supplied with appropriate information (temperature, isotopic concentration of fission products, and material properties), this framework would predict thermal conductivity in heterogeneous nuclear fuel with an operating history. The fission product significantly hindering thermal transport in UO₂ is widely accepted as xenon [7], a noble gas which accumulates in UO₂ over its operational life-cycle. The bulk of thermal conductivity characterization is done with molecular dynamics (MD) methodology, which models energy flow explicitly through atomic motion. MD is effective in predicting thermal conductivity, but is only able to model small systems of atoms due to the significant computational cost of the method.

To this end we leverage the code Rattlesnake, which solves a second-order form of the Boltzmann transport equation (BTE) using the self-adjoint angular flux (SAAF) formulation [8] using continuous finite element (CFEM) or discontinuous finite element spatial discretization [9]. Rattlesnake is a member of the multiphysics object-oriented simulation environment (MOOSE) [10] architecture developed and maintained by Idaho National Laboratory. We have shown Rattlesnake may be adapted to simulate phonon transport, and demonstrates promise in connecting transport phenomena at the nanoscale to properties which may be used at the macroscale [11]. We motivate the use of the SAAF formulation of the transport equation and discuss the advantages and disadvantages of the numerical solution of this equation in comparison with solvers for the traditional first-order integro-differential form of the equation.

Our numerical solution technique involves traditional source iteration (SI, a Richardson iteration) methods combined with a robust linear algebraic solver to solve the systems of discretized equations generated by the second-order BTE. With the application of preconditioning, we are able to rapidly solve these systems with tremendous savings in computational cost relative to traditional methods. Additionally, we have the capability to apply nonlinear diffusion acceleration as these simulations become very acoustically thick to yield an even greater convergence acceleration. As such, this approach opens the door for BTE simulation to model heat transport in relatively large systems that contain realistic and statistically representative material microstructures. This approach can potentially become a truly practical and predictive tool to nuclear engineers and material scientists.

For brevity, we have suppressed the independent variables in many of the terms of the equations. The phonon phase space density is $fω=f(r,Ω,ω,p),$ where r is the spatial coordinate $r≡(x,y,z), Ω$ is the unit vector denoting the direction of travel $Ω=(ϕ,θ)$⁠, and ω is angular frequency. p is polarization; the geometrical orientation of phonon travel is transverse in two directions (T) or longitudinal (L). Group velocity v_g is related to the propagation speed of phonons, which can have either acoustic (A) or optic (O) modes. However, for this work, we assume a single phonon speed averaged over the acoustic modes and polarizations at varying temperatures, an assumption for the transport of gray phonons which is addressed in Sec. 2.1.

The scattering kernel $f˙ω|scatt$ contains nonlinear operators and includes contributions from processes such as anharmonic phonon interactions or material defect scattering. Other contributions to Eq. (2) can include thermal boundary resistance (TBR) or defect scattering. We use a weak formulation of the phonon transport equation, in which we include upwinding terms to describe the interface condition of TBR, which we develop later.

where v_g has been brought to the right-hand side to obtain Λ, the phonon mean free path (the product of group velocity and relaxation time).

and can be applied to a phonon frequency distribution to yield phonon radiant intensity $I(r,Ω̂)$ which has units of $W·m−2·sr−1$⁠. The phonon frequency distribution is multiplied by $vgℏωD(ω)$⁠, summed over all phonon branches and polarizations and integrated over all possible frequencies (limited by the vibrational frequency of the medium). Here, $ℏ$ is the reduced Planck's constant, and $D(ω)$ is the phonon density of states.

Equation (6) is the equation of phonon radiative transfer (EPRT) [13]. The allure of the gray phonon EPRT is that the consequences of ballistic transport between geometric scattering features are determined explicitly but the entirety of the local intrinsic and extrinsic scattering physics is lumped into a single parameter that may be obtained empirically, from first principles, or some hybrid mixture of both. The next level of approximation is to model transport explicitly including the spectrum of participating phonon frequencies.

The central limitation of the gray approach is its inability to model transport in anisotropic materials or across strongly frequency selective boundaries. Further refinement includes modeling frequency-dependent transport across the phonon spectrum, but to still treat collision terms through the single relaxation approximation. From the point of view of our efficient SAAF solution method, the two approaches are numerically identical with independent equations coupled only through a single integral term. For simplicity of demonstrating the SAAF approach, we limit the scope of this work to the gray phonon model. While we have chosen to demonstrate the effectiveness and efficiency of our numerical approach in the context of frequency-independent transport, the extension to the solution of the frequency-dependent BTE with the relaxation time approximation is trivial since frequency appears as a parameter in the equation, i.e., the solution in each frequency group is decoupled from all other groups.

Morel and McGhee [8] outline an algebraic technique for the derivation of the self-adjoint form of the neutron transport equation. From a computational perspective, the SAAF formulation is advantageous, since the full angular flux becomes the unknown. Using reflecting boundary conditions becomes easier with the availability of the full angular flux, as the incoming and outgoing directions are coupled in the same manner as the first-order form of the transport equation. Through the application of a CFEM spatial discretization, the matrices are symmetric positive definite, which allows for the use of solution techniques such as the preconditioned Krylov family of solvers [10]. The MOOSE framework uses CFEM as a spatial discretization and by default employs Jacobian-free Newton Krylov [14] with preconditioning as a nonlinear iterative solution method.

From Eq. (8), the change in the phonon intensity at a point has two contributions: a streaming term from the spatial variation in intensity, and a collision term due to the deviation of the radiance from the equilibrium phonon radiance $I0(r)$⁠. Due to the implication of the single mode relaxation time approximation, the phonon radiative equilibrium intensity will be defined with a condition of zero heat generation, $∇·q=0$⁠. This suggests that a phonon radiative equilibrium could exist at all possible frequencies and provides some justification for the gray media formulation [13].

where T_b is a constant driving boundary temperature. For reflecting conditions, outgoing angular phonon radiance is defined as the reflection of the incoming angular phonon radiance, merely undergoing a directional change, i.e., $I(r,Ω)=I(r,Ω′)$ where we map $Ω$ to $Ω′$⁠.

The angular variable $Ω$ is discretized via the discrete ordinates method, sometimes referred to as the “S_N method.” The transport equation becomes a set of N equations for the radiant intensity in each discrete angle. Solutions to the EPRT generate angular phonon radiative intensity, and Rattlesnake employs level-symmetric quadrature to numerically integrate this quantity over solid angle, to obtain “moments” of the radiance. The S_N method has advantages in heterogeneous media over the spherical harmonics method $(PN)$⁠, in that the $PN$ angular moments are tightly coupled, and their solution requires more computational resources [15]. The recent application of a hybrid $SN−PN$ scheme to discretize the angular variable in the frequency-dependent phonon BTE has been shown to exhibit slow convergence in homogeneous silicon [16] which suggests that the P_N angular discretization approach may not be optimal for this type of problem.

The discrete ordinates SAAF transport equation has the feature that in each quadrature direction, a linear system of equations arising from the spatial discretization of an elliptic operator is solved for the angular intensity. This means that software for solving the diffusion approximation to transport can be readily converted to solve the transport equation. The treatment of voids and certain boundary conditions require special care, however. In contrast, the solution of the first-order form of the transport equation involves transport “sweeps” [15], where incident intensities from the problem boundary, and interior sources along the way, are propagated through the spatial mesh along the direction of travel to the exiting mesh boundary. The ordering of this mesh sweep is angle and problem dependent, and in multiple dimensions cyclic graphs are possible. Developing sweep algorithms that scale to large numbers of processors is an active research area, whereas efficient parallel solvers for elliptic equations are much more mature.

Other researchers have taken different approaches to computing heat flux, incorporating quantities obtained through MD simulations (phonon group velocity, wave vectors, and angular frequencies [17–19]).

Though heat flux and temperature gradient are computed over the entire domain, the dimension of interest is one where a temperature gradient is applied. This methodology is used in MD simulations and is repeated here.

To evaluate the effectiveness of Rattlesnake as a deterministic transport engine, we compared our numerical solutions of temperature and thermal conductivity in room temperature, homogeneous silicon to the work of Yilbas and Bin Mansoor, who performed deterministic phonon transport simulations in silicon under equivalent conditions [20]. They used a forward and backward finite difference discretization scheme to solve the BTE for phonons.

Yilbas and Bin Mansoor modeled a thin film of silicon, a common configuration which is used in many phonon transport simulations (with both deterministic and Monte Carlo methodologies) as a benchmark problem. Material properties of room temperature silicon are well known, and transport behavior at the nanoscale has been studied [13]. Material properties for this model were obtained from the open literature [13,21], and are listed in Table 1.

We use the generalized minimum residual [23] method preconditioned with algebraic multigrid [24] to solve the linear system, with an iterative convergence criteria of $ϵ=10−6$⁠, and an S₈ angular quadrature. Computation times for coarse and fine mesh cases were approximately 8 and 57 s each. Simulations were performed on a single core 2.8 GHz Intel i7 CPU with 16 GB RAM. The coarse and fine mesh solutions are within 10⁻⁵ of each other; for simulations of homogeneous media, it may be appropriate to use a more coarse spatial mesh to decrease computation time. The behavior of the nondimensional temperature solution for coarse and fine mesh cases obtained with Rattlesnake agree well with those from Yilbas and Mansoor, shown in Fig. 1. The temperature profile has a slight curvature, which is influenced by spatial discretization. Improperly scaled finite elements do not produce the desired solution behavior. As the acoustic thickness increases, phonon scattering regimes shift from ballistic to diffuse. In an acoustically thin medium, where $L≈Λ$⁠, phonons leaving a colder boundary are in the ballistic scattering regime, and propagate far across the medium to reach the hotter boundary causing the material temperature to be smaller than that associated with the prescribed incident intensity. In an acoustically thick medium, where phonons are in the diffuse scattering regime, this effect is significantly diminished. These are boundary scattering effects and are well characterized in simulations of phonon transport [13,20].

Du et al. [7] investigated the effects of xenon presence on thermal conductivity of UO₂. MD methods were used to simulate the effect of various concentrations and geometric configurations of xenon in the UO₂ lattice for a range of temperatures. They concluded that randomly dispersed xenon in the fuel matrix has a more significant impact on thermal conductivity than quantized xenon bubbles. At higher temperatures, phonon–phonon scattering from normal and Umklapp processes becomes a main contributor to the suppression of thermal conductivity, and heat transport is locally disrupted at the xenon defect. The phonon mean free path $(Λ)$ in UO₂ becomes shorter at high temperature and diffuse scattering dominates. Xenon concentration was limited to about 1% of volume in the simulations.

We model a selected problem from Du et al., computing temperature, heat flux, and thermal conductivity in a cell of UO₂ with a bubble of xenon in the center, in the absence of grain boundaries. The behavior of xenon in a UO₂ lattice at high temperature and pressures has been reported to vary widely, and available experimental data are minimal. Computational studies have been performed to investigate expected temperatures and pressures of xenon using various methodologies which drew varying conclusions [25–27]. We performed MD simulations to determine the properties of xenon at specific pressure and temperatures; these data are used to compute Λ in the bubble. Bates [28] performed thermal conductivity measurements on stoichiometric, unirradiated UO₂ for a large array of temperatures. We extract Λ from measured thermal conductivity in the Bates study, and perform simulations using the same bubble geometry to determine the impact of xenon on κ using values of Λ from Bates. We compare these results to the bulk κ of pristine UO₂ measured by Bates and to κ computed using parameters from Du et al.

The role of thermal boundary resistance at the UO₂–Xe interface is investigated, as incoming delocalized waves may reflect diffusely or specularly off the xenon bubble and have a significant effect on the local thermal conductivity. We introduce the diffuse mismatch model (DMM) and present a method to characterize thermal boundary resistance at heterogeneous defects in three-dimensional (3D).

We use Rattlesnake to simulate phonon transport in a cell of UO₂ with a xenon bubble in the cell center and no grain boundaries. The spatial domain $D$ is a rectangular cell, 25 nm along the z-axis with a cross section of 3.8 nm × 3.8 nm, consistent with the geometry used by Du et al. The xenon bubble has a radius of 1 nm and accounts for approximately 1% of the total volume. The finite element geometry is constructed with CUBIT, an unstructured mesh consisting of 100,689 tetrahedral elements (Fig. 2). The linear system is solved with the algebraic multigrid-preconditioned generalized minimum residual method with convergence criteria of $ϵ=10−6$⁠. We performed simulations with an angular quadrature order of S₂₄, a large amount of ordinates helps to mitigate ray effects [29], which can occur in deterministic phonon transport simulations.

Du et al. reported κ in UO₂ with Xe at 300 K, 800 K, and 1500 K. Where possible, we replicate their simulation conditions and report dimensionless temperature, heat flux, and thermal conductivity. We extract Λ from values of thermal conductivity for unirradiated UO₂ as documented by Du et al. In each simulation, a 1 K temperature difference is applied along the z-axis. We use the same mesh to perform additional simulations for different values of temperature using Λ extracted from experimental values of thermal conductivity measured by Bates [28]. Simulations using the mean free path from Bates are performed independently of the Du et al. simulations, in order to gain insight into the effect Λ has on overall heat flux, temperature gradient, and thermal conductivity.

The only material property entering into the gray BTE in Eq. (8) is the phonon mean free path, Λ, and so we need to determine these for UO₂ and Xe. However, it is also necessary to determine each material's phonon radiance, $I0(T)$⁠, as a function of temperature in order to impose the correct scalar flux at the external boundaries, and to set the transmission coefficients at the internal boundary.

For the UO₂, v_g and C_v of acoustic modes are computed from first principles calculations, and then Λ is chosen so as to reproduce the experimentally measured values of κ in unirradiated UO₂ [28]. The calculation of v_g and C_v in UO₂ is described in detail in Sec. 3.2.1. For Xe, the κ, v_g, and C_v are computed from classical molecular dynamics simulations, and similarly used to infer Λ in the Xe.

The effective group velocity and effective mean free path of the gray phonons were assumed to be isotropic in both the UO₂ and the Xe. In order to understand the degree to which transport is ballistic in either region we consider the acoustic thickness, ζ, in each domain. This is the domain size scaled by the material's effective phonon mean free path. Acoustic thickness of UO₂ is the ratio $ζUO2=DUO2,z/ΛUO2$⁠, where $DUO2,z$ is the distance between the hot and cold sides of the UO₂ cell. Similarly, the acoustic thickness of the Xe, $ζXe=DXe/ΛXe$⁠, describes the diameter of the bubble relative to the effective (gray) phonon mean free path in Xe; acoustic data and mean free path for simulations using values from Du et al. are contained in Table 2.

The transport properties of Xe are strongly tied to the Xe pressure, which in turn is set by the surface tension of the UO₂/Xe interface and the bubble size. In a 2 nm diameter bubble of Xe in UO₂, the pressure is estimated to be between 2 and 5 GPa, and so in this work, the Xe pressure was assumed to be 3 GPa at all temperatures studied. At this pressure, Xe is either solid or liquid across the temperature range that we study here, and so the use of a gray phonon model of transport is justified in the Xe. Across the range of temperatures we study, we hold the size of the bubble fixed (with radius of 1 nm), meaning that as the temperature is changed the number of Xe atoms in the bubble is not constant. This means that our sweep of simulations does not represent the change in thermal conductivity due to heating UO₂ containing Xe bubbles from 300 to 1500 K. However, it does provide us insight into the ballistic versus diffusive contributions to thermal resistance from 1 nm bubbles at different temperatures.

For UO₂, the effective group velocity, volumetric specific heat, and radiance of gray phonons was computed by averaging the properties of the three acoustic branches of the phonon dispersion over the entire Brillouin zone. The phonon dispersion was computed on a 50 × 50 × 50 q-point grid using Phonopy [30] based on the structure symmetry of UO₂ (Fd-3m) with interatomic force constants computed from first principles. The UO₂ simulations were executed with the plane-wave basis projector augmented wave method within the density functional theory framework as implemented in the Vienna ab initio simulation package [31–33]. A plane-wave energy cutoff of 600 eV was employed in the local density approximation [34], with a 6 × 6 × 6 Monkhorst-Pack k-point grid. A $2×2×2$ super-cell of the UO₂ unit cell (with 4 O and 8 U atoms) was used for all calculations, including calculation of the force constants. The perfect super-cell was found to be relaxed to a $1×10−3 eV/Å$ ionic tolerance and a $1×10−5$ eV electronic tolerance. UO₂ is antiferromagnetic, but there exist ferromagnetic solutions to the Kohn–Sham equation, and Vienna ab initio simulation package can get trapped into a ferromagnetic state. To prevent this from happening, the magnetic moment tag was selected to ensure that alternating uranium atoms in the structure had opposing spins, and the spin of the oxygen atoms was set to zero. We further used a Hubbard parameter, U, of 4.50 eV, and a Hund's exchange parameter, J, of 0.50 eV. The resulting electronic density of states (Fig. 3) agrees with that of Wang et al. [35]. The phonon dispersion (Fig. 4), also in good agreement with that obtained in the same reference [35].

This velocity and the heat capacity are only very weakly temperature dependent over the temperature range of interest, and thus, we approximate it by their average values of $1764 m·s−1$ and $1.007×106 J·m−3·K−1$⁠, respectively. Note that this effective velocity of the gray phonon bath is approximately 0.45 of the mean speed of sound obtained from the same UO₂ calculations. It is in effect the average group velocity of acoustic phonon modes weighted by the thermal energy in the mode. Similarly, the effective specific heat is not the true specific heat of UO₂, but only the contribution to its specific heat from the acoustic modes and the computed value for this is is close to the high temperature limit for the acoustic modes of $12kB/a3$⁠.

The thermal conductivity and transport properties of xenon under high pressure were computed from classical molecular dynamics simulations performed using the large-scale atomic/molecular massively parallel simulator package [36]. Interatomic forces were modeled using a Lennard–Jones potential with large-scale atomic/molecular massively parallel simulator parameters $ϵ=0.00425 eV$ and $σ=4.29 Å$⁠, with all interactions truncated after $20 Å$⁠.

Simulated systems of 10,000 Xe atoms under 3 GPa were prepared at a series of temperatures from 300 to 1700 K. This was achieved by equilibrating the system over a 500 ps simulation in the constant number of atoms, pressure, temperature ensemble before turning of the thermostat and barostat and simulating for a further 50 ps in the microcanonical ensemble. In these simulations, it was found that the Xe was solid at temperatures below $∼600$ K, and above that, remains liquid up to $1700 K$⁠. Once the systems were prepared, the system was simulated in for a further 1 ns in the microcanonical (NVE) ensemble, during which the thermal conductivity was computed using the Green–Kubo method [37]. The Green–Kubo method is founded on the fluctuation dissipation theorem to determine the thermal conductivity of a system from the lifetime of its natural thermal fluctuations during a simulation of the system at equilibrium. A minimum of six simulations were performed using different random starting configurations and the results averaged to obtain each thermal conductivity datum.

For v_g of Xe, we use the speed of sound computed at each temperature from the Xe's density and isentropic compressibility. The isentropic compressibility was computed by simulating adiabatic expansion. At each temperature, the system was first cycled in an NVE ensemble after which system dimensions were slightly increased, followed by another NVE cycling step. The compressibility was calculated from the differences in system volume and pressure before and after expansion. This set of simulations also served as the source of density data. The Xe density was also used to compute Xe's volumetric specific heat capacity at each temperature.

The computational approach for determining both thermal conductivity and speed of sound were validated by computing the pressures at slightly lower pressures for which there exists experimental data [38] and finding the properties to be in reasonable agreement. The computed density, thermal conductivity, v_g, and $ΛXe$ are plotted in Fig. 5, and clearly show a transition in properties between the solid and liquid Xe.

We must consider the resistive effect at material interfaces: the phenomenon of TBR, the ratio of a temperature discontinuity at an interface to the heat flux across that interface, due to a material difference at the junction. TBR is an extraordinarily subtle phenomenon and has some consideration in other deterministic phonon transport studies [39,40]. However, it is very important to consider in simulating phonon transport, and has been characterized in a number of other MD and Monte Carlo studies [41–43].

The physics of TBR are important to phonon transport because of how phonons behave when they encounter a physical interface between two adjacent materials. At this junction, phonons become subject to a phenomenon which manifests as a transmissive and resistive effect for phonons penetrating an interface into another material. This physical effect occurs as the intrinsic properties of material change. Phonons define the internal energy of a material; when they cross a boundary from one material into the next, the change in their contribution to internal energy as well as change in their velocity must be considered.

We develop the DMM [44] in a deterministic framework for 3D general geometries. We assume all phonons are diffusely scattered at the interface, with outgoing radiance emitted isotropically, and that scattering destroys the correlation between the wavevectors of incident and outgoing phonons; the probability that a phonon will scatter into a given side of the interface is independent of the phonon origin. Enforcing these conditions makes the probability of scattering into a given side proportional to the phonon density of states on that side, additionally constrained by the principle of detailed balance.

In Fig. 6, we denote “−” as the upwind radiant flux and “+” as the downwind radiance. The scalar radiance is identified by $Iα;β±(rint)$⁠. We must solve for the upwind and downwind scalar radiance at the interface $rint$ to characterize the effects of TBR in our transport simulation. In the transport equation, we solve for the angular radiance, which is integrated over solid angle to solve for scalar radiance.

At the location $rint$ with its unit vector normal to the interface which points from material α to material β denoted as $nα→β$⁠, we can write conservation equations which define the flow of phonons immediately at both sides of the interface. On the β side of the interface, phonons which flow away from the interface into material β come from two sources: they are transmitted through the interface from material α, and reflected from those incident on the surface from the β side.

A similar expression describes the upwind diffuse radiance of phonons flowing from material β into material α. These new expressions for the scalar radiance are now the isotropic emission sources on either side of the interface and are distributed at each direction of outgoing angular radiance at the interface $Iα;β(rint±,Ω)$⁠.

The effective implementation of this model allows for the description of localized heat flux and temperature around defects in heterogeneous structures. The DMM is a relatively crude descriptor of TBR, an improved model of interface physics would necessitate the inclusion of anharmonic effects, and would also need to be adaptive to phonon frequency selection over the interface itself [45]. We computed transmission and reflection coefficients for UO₂ and Xe from the material properties determined in Secs. 3.2.1 and 3.2.2; these are shown in Fig. 7. At all temperatures, the phonon radiance is approximately two orders of magnitude larger in the UO₂ than in the Xe; the boundary is highly resistive to phonons flowing from UO₂ into the Xe bubble with approximately 40% transmitted at 300 K, decreasing sharply as temperature increases.

We report simulated heat flux and thermal conductivity for an array of temperatures with two sets of values for $ΛUO2$⁠; one generated from the MD results of Du et al. (denoted $ΛDu$⁠), the other extracted from experimentally measured values of κ in unirradiated samples of UO₂ from Bates (denoted $ΛBates$⁠). In both cases, we use material properties for Xe based on MD simulations we performed, as well as the same spatial mesh. We observe the effects of thermal boundary resistance at the Xe bubble, which play a role in the amount of overall thermal resistance the bubble provides. In all cases, the presence of xenon lowers the thermal conductivity in the UO₂.

Figure 8 shows all values of simulated κ and the measured pristine κ. In the upper half of Fig. 8, we compare to Bates and follow a similar trend showing decreased thermal conductivity with increasing temperature. This is further affected by the presence of the xenon bubble; κ is reduced by approximately 30–55% over the temperature range with the sharpest difference occurring at lower temperature.

In the lower half of Fig. 8, we compare our κ to that of Du et al.; while we follow a loose trend in the shape of the curve, we experience a large discrepancy in simulated values. We under-predict κ in simulations with and without a Xe bubble, which may have multiple causes. Numerical results of κ from this study are compared to those of Du et al. in Table 3. Our computed group velocity in UO₂ is lower by approximately a factor of 2 compared to Du et al.; it is no surprise that our κ with Xe is also lower by approximately the same factor. We justify this exclusively for a gray approach in that we assume all phonons do not travel at the speed of sound in UO₂; indeed from the calculations in Sec. 3.2.1 it is clear that the group velocity is averaged, but has contributions from phonons with varying Λ lumped into a single term. Some phonon modes may be highly sensitive to frequency, and this detail is washed out in the gray approach. This assumption may explain some of the underestimation in the ballistic effects. In addition, we may not be capturing some intrinsic phonon scattering, and we currently do not have anharmonicity built into the relaxation time; there is no way to discern these effects. The relaxation time $τeff$ has many contributions: Umklapp processes, defect scattering, boundary scattering, and the resonant scattering of phonons. Classic MD can be used to characterize many of these processes, but MD is not able to capture certain quantum effects at low temperatures, such as specific heat [46], and this may be a reason for the discrepancy.

Dimensionless temperature for all simulated temperature is shown in Fig. 9. The influence of the xenon bubble is clear, and jumps at the interface are observed; these are more pronounced at lower temperatures when phonon scattering is highly ballistic. This effect results in localized negative temperature gradients, which was also observed by Yang in Si nanowires [43], and the gradient monotonically shifts toward zero with increasing temperature. Note that while the temperature gradient becomes negative, the net flux in this region is still positive (see Fig. 10)—heat is apparently flowing uphill. This result is counterintuitive when thinking of heat flow in the diffusive limit, but is entirely consistent with a ballistic picture of transport in which the energy of the phonon gas at any point contains nonlocal information.

The centerline heat flux $q(r)$ shown in Fig. 10 has been normalized to the 300 K value; $q(r)$ is inversely proportional to temperature and experiences a steady decline as temperatures increase. As temperature increases, $ΛUO2$ decreases and diffuse scattering becomes more prevalent, which contributes to the reduction in heat flux. As a result of TBR, large portions of the phonon radiance are reflected at the xenon bubble, decreasing local $q(r)$ by resisting the flow of phonons.

Heat flux in the UO₂ region remains approximately constant along the temperature gradient and changes drastically at the Xe bubble. We observe this effect in Figs. 10 and 11, where heat flux is severely depressed in the region local to the Xe bubble. The presence of a single xenon bubble does not significantly impact average $q(r)$ in the domain, but it does affect the behavior of the local heat flux. Du et al. established this by performing MD simulations which include multiple Xe bubbles and Xe randomly dispersed in the UO₂ matrix [7]. Ray effects are observed at lower temperatures when phonon scattering is highly ballistic (slight oscillations in $q(r)$⁠) but vanish as scattering becomes diffuse. With increasing temperature, the heat flux is gradually suppressed, and the effects of Xe on the local heat flux become less significant. $ΛUO2$ decreases by approximately a factor of 3 between the temperature extremes, and this decrease is more detrimental to heat flux than compared to the presence of a singular bubble of Xe. Additional Xe bubbles would have a greater negative effect on bulk thermal conductivity.

Table 4 contains the total number of source iterations required for convergence, and total acoustic thickness of the spatial domain over the range of temperatures. As ζ increases, required iterations decrease; this is counterintuitive as purely scattering thermal radiation and phonon transport simulations tend to require more iterations for convergence with increasing ζ. This effect is potentially related to the oscillations experienced in the ballistic scattering regime, where Λ is on the order of the entire spatial domain (Casimir limit).

We have presented the features of a 3D, generalized geometry radiation transport code modified to simulate phonon transport in a gray formulation. We have implemented the physics for thermal boundary resistance in 3D to describe phonon transport behavior at localized defects in the material. We have presented our deterministic transport results for a simulation of a 3D domain of UO₂ with a Xe impurity and compared them against MD results for similar geometry and simulation parameters [7]. Additionally, we use $ΛUO2$ extracted from experimentally measured pristine UO₂ for the same simulation setup, and mimic the shape of the κ curve in Ref. [28] but with lower values of κ due to the Xe presence.

The transport method we have presented is trivially extendable to simulate a multifrequency phonon spectrum using input variables derived from density functional theory (DFT) simulations, consistent with our discussion in this work (Secs. 3.2.1 and 3.2.2). The use of deterministic methods to simulate phonon transport is an underdeveloped aspect of the phonon transport community. Classical molecular dynamics simulations and DFT electronic structure calculations can provide detailed information about material properties, dispersion relations, and thermal conductivity but do so only at a local or nanometer scale. It is well understood that resistive processes also arise from the mesoscale structure of materials and these cannot be captured efficiently from atomistic calculations (a point that is reinforced by the results in this work).

By coupling Rattlesnake to MD and DFT methodologies, we show a way to bridge the gap between the atomistic and engineering scales. This multiscale method provides a framework for rapid prediction of the engineering-scale thermal conductivity in materials with evolving microstructures. Such a tool is particularly imperative for modeling nuclear fuels and their surrounding structural materials in which the thermal conductivity of the materials is a central component of the system performance.

We acknowledge Idaho National Laboratory for their support and thank Sebastian Schunert, Yaqi Wang, and Daniel Schwen for their valuable assistance. We also thank David Andersson of Los Alamos National Laboratory for providing MD data for comparison. For the DFT and MD simulations, this work used the Extreme Science and Engineering Discovery Environment (XSEDE).

• •

  National Science Foundation (Grant No. ACI-1548562).