Accident tolerant composite nuclear fuels

Investigated accident tolerant nuclear fuels are fuels with enhanced thermal conductivity, which can withstand the loss of coolant for a longer time by allowing faster dissipation of heat, thus lowering the centerline temperature and preventing the melting of the fuel. Traditional nuclear fuels have a very low thermal conductivity and can be significantly enhanced if transformed into a composite with a very high thermal conductivity components. In this study, we analyze the thermal properties of various composites of mixed oxides and thoria fuels to improve thermal conductivity for the next generation safer nuclear reactors.


Introduction
The novel, composite fuels with enhanced thermal conductivity can withstand the loss of coolant for an extended time by allowing faster dissipation of heat, thus lowering the centerline fuel temperature.This consequently prevents the melting of a reactor core, minimizes oxidation of cladding and related hydrogen production and therefore prevents hydrogen explosion that occurred during Fukushima accident.The composite fuels allow reduction of the thermal gradients and thermal stresses that contribute to cracking, thereby increasing the longevity of the fuel.
Although many factors need to be investigated before nuclear fuel composites are adapted to the service in the harsh environment in the nuclear reactor the suitable fuels must have a high thermal conductivity.Traditionally metal oxides were used as fuel since they have a high melting point [1].Nitrides and Borides (with the appropriate isotopic composition that excludes strong neutron absorber 10 B) are also attracting attention due to their large thermal conductivity from the electronic contribution that increases with increasing temperature [2,3].In this work, we evaluate the thermal conductivity of ceramic composites where the heat carriers are primarily phonons.In our previous work [4] we have found that the phonon contribution to the thermal conductivity of thoria, (an alternative fuel with a high melting point) can be well reproduced up to very high temperatures by a simplified Slack model [5] with the input parameters calculated using first principles molecular dynamics.Silicon Carbide (SiC) and Beryllium Oxide (BeO) are materials proposed to be used in accident tolerant fuel.Within the lower temperature range where quasi-harmonic approximation (QHA) applies we have shown [6] using the first principle code Quantum ESPRESSO (QE) [7] that for SiC and BeO both Slack model and the ShengBTE (BTE) a solver for phonon thermal conductivity [8] predict similar thermal conductivities in agreement with the experiment except for cubic SiC where Slack model over-predicts thermal conductivity by a factor of two.Here, we investigate it further for other candidate components and evaluate the effect of high thermal conductivity of the composite fuel on the thermal performance.

Methodology
A straightforward and efficient interface (QE_nipy_advanced), developed in ipython by our group [9], is used to operate Quantum ESPRESSO (QE), an open source code for materials simulation [7].QE is a first principles code using density functional theory, plane waves and pseudopotentials; and it can predict ground state material properties in agreement with experiment provided the parameters like plane wave cutoff energies, k point optimal mesh and equilibrium structure are optimized accurately as described earlier [10].We also found previously [9,10] that new generalized gradient approximation (GGA) functional for solids, PBEsol used with the pseudopotentials generated for PBE or PBEsol functional lead to equilibrium lattice constants that agree very well with experiment and therefore it is used here in our new calculations.Equally good results are obtained with local density approximation (LDA) and norm conserved pseudopotentials which we have been used predominantly in our previous calculations [10].
We calculate lattice governed thermal conductivity using QHA as implemented in QE code and two methods as discussed above: Slack model [5] and BTE code [8].
Once the thermal conductivity is calculated for the limited temperature range within QHA, we fit it into an analytical form developed for phonons contribution as a function of temperature as: (1) where "a" and "b" parameters are independent of the temperature (T), and are associated with the contribution of phonon defects (a) and phonon-phonon (b) scattering to the mean path.The calculations and the related inverse fits (Eq. 1) are for fully dense (100 % TD) nuclear materials and therefore a porosity factor must be added [4] when comparing with experimental results [1].We demonstrated that this equation describes phonon assisted thermal conductivity for up to a very high temperature [4].
The phonon contribution to thermal conductivity as described by Eq. 1 decreases with temperature, and therefore it is of interest for reactor safety analysis to estimate the minimum thermal conductivity (kmin).We evaluate it here, using Young's modulus and Clark's model [11], from the equation: where Mav is the average mass per atom (Mav = M/(n NA)), where M refers to the molar mass, n is the total number of atoms per formula, NA is the Avogadro's constant), ρ is the density, and kB is the Boltzmann's constant.
The composite fuel can be fabricated as fibers or as particles of high-temperature conductivity component dispersed in fuel matrix.A simple analysis of thermal conductivity for these various arrangements for assumed here small volume fractions (5%-15%) of additivities shows that the thermal conductivities are comparable.We assumed therefore that the thermal conductivities of the composites could be approximated as an average over volume concentrations (C) of the thermal conductivities of their constituents (e.g. for ThO2-BeO):

Thermal conductivity
We carried out calculations of the lattice thermal conductivities as described in Section 2. In Fig. 1 we compare the theoretical results for the selected materials (SiC, BeO, and BAs) with experimentally observed one order of magnitude higher thermal conductivity than Urania and Thoria.According to the previous studies [6], we observed in Fig. 1 that with the input parameters derived using QE espresso code Slack model predicts one order of magnitude higher lattice thermal conductivity than Urania [1]  The total energy convergence of the BeO and structure in QE code was obtained using a kinetic energy cut-off of 75 Ry and the Brillouin zone integration over a Monkhorst-Pack of 8 × 8 × 5 mesh.Since we used the ultra-soft pseudopotentials, we need an augmented charge around the ion core.Hence the kinetic energy cut-off of plane waves basis describing charge density was taken twelve times the energy cut-off of the used wave functions.The dynamical matrix for phonon density of states of BeO was calculated on a mesh of 6 × 6 × 4 qpoints in the irreducible Brillouin zone respectively.While the input to Slack model was obtained using DFPT QHA as implemented in QE for ShengBTE code we used a supercell of 4 × 4 × 4 to calculate the third order interatomic force constants (IFCs) for BeO by solving BTE.The force cut-off distance was set such that the interaction range is up to the five nearest neighbours for BeO.A mesh of 6 × 6 × 4 q-points was used to calculate the second order IFCs using QE and needed to compute the kL of BeO.
For the simulation of the thermal properties of BAs and BP, we used available now PAW pseudopotentials developed with the new PBEsol functional.The optimal parameters for BAs were found to be 100 Ry a cut-off energy and mesh of 7 × 7 × 7 k-points, and we used cut-

CMPSE2017
off of plane waves basis describing charge density four times the plane wave energy cut-off.However, in case of BP, the convergence of energy was achieved with a mesh of 8 × 8 × 8 k-points and a plane wave cut-off energy of 100 Ry.Also, we observed that setting the plane waves basis describing charge density at a high value (1000 Ry) has led to a non-physical linear thermal expansion and therefore we have maintained the value at four times larger than the energy cut-off.Furthermore, we used a 5×5×5 q-points mesh for the second order force constants required for the thermal expansion calculation of BAs and BP.
A mesh of 5 × 5 × 5 q-points was used for BAs and BP to calculate the second order interatomic force constants (IFCs) used in thermal expansion calculations for the Slack model.The same mesh was used for BTE calculation for BAs with 5 × 5 × 5 supercells for third order IFCs.
We note however that due to the usage by the Slack equation of the thermal expansion coefficient, which is derived from the second derivative over the temperature of volume, the derived thermal conductivity shows numerical accuracy dependence.Therefore, for example, the calculations of the thermal conductivity of BAs, using thermal expansion coefficients calculated considering the same 14 deformations, 100 Ry plane-wave cut-off energy and different strains from 0.99 to 1.04 or 1.02 (indicated as 4% and 2% in Fig. 1) differ by a factor of two.
Additionally, we performed a comprehensive testing on the dependence of the number of deformations, cut-off energy and k points mesh on the thermal conductivity of BP predicted by Slack model (Fig. 2.).We can see that the result (indicated by dotted line) for strain 0.99-1.02,20 deformations, cut-off energy 80 Ry and 8 × 8 × 8 k point mesh significantly overpredicts the thermal conductivity, and this may be due to limited numerical accuracy at this strain.Therefore, we increased maximum strain range to 0.99-1.04.In Figure 3, we present for BP the selected results as indicated for a larger range of maximum strains 0.99-1.04,which were consistent with each other.We note that the calculated by Slack thermal conductivities for maximum strains 0.99-1.04agree with experiment [12] for our present calculations for BAs (Figure 1) and as shown in Figure 3 BP results also agree with experiments [13,14] at temperature 300 K and above.1) using BTE predicted the value an order of magnitude higher.Though our BTE results for BAs (Fig. 1) predict lower thermal conductivity than published later for a larger temperature range by the same authors (4288 Wm -1 K -1 versus ~10000 Wm -1 K -1 at 100 K) [15].The calculated here thermal conductivity at 100 K for BP using Slack method is between 3000 Wm -1 K -1 and 4000 Wm -1 K -1 as shown in Fig. 3, while it was estimated previously [15] to be also much higher ~10000 Wm -1 K -1 .The BTE calculations for BAs were calculated using the same PAW-PBEsol pseudopotentials and optimal parameters as for Slack model as described above and the same PBEsol functional.In the previous BTE calculations, local density approximation (LDA) was used with norm conserved pseudopotentials and only 80 Ry cut-off energy and 6 × 6 × 6 k points mesh, which we found not sufficient for the used here PAW-PBEsol pseudopotentials.Our comprehensive investigation indicates that the predicted lattice thermal conductivity by both BTE and Slack model for some materials may be very sensitive to the tuning optimization of parameters used in DFT calculation, selected type of DFT functional and the respective pseudopotentials.
In Figure 4 we present calculated here phonon spectra for BP and BAs as indicated.

BTE results of a very high thermal conductivity of BAs, originated from the limited thermal resistance of three phonon scattering, umklapp (e.g. two acoustic phonon frequencies (a) must be equal to that of optical phonon: a+a=o) and which is largely restricted in BAs due to a large gap between acoustic (a) and optical (o)
phonons remains to be confirmed experimentally.Acoustic phonons bunching also prevents umklapp involving three acoustic phonons (a+a =a) due to absence of high frequency acoustic phonons in BAs.As presented in Figure 4 the a-o gap is smaller in BP than in BAs, therefore its thermal conductivity is expected to go down with temperature according to Eq. 1 and as shown in Figure 3 and in agreement with experiment [13,14] and BTE calculation [15].Γ-X-Γ-L) and the corresponding phonons density of states in a.u. vs. the same frequency on the right panel.

Tables 1a,b show the results for SiC, both cubic and hexagonal, wurtzite BeO, cubic thoria, ceria, BP and BAs.
The presented in Table 1a parameters: a and b, are for the analytical expression (Eq.1), and in the last column we present thermal conductivity calculated using Eq. 1 for 400 K.The presented properties in Table 1b were calculated at 400 K, except for the minimal thermal conductivity.The minimal thermal conductivity was evaluated using Eq. 2 and both Young's modulus and density were calculated at 400 K.Both SiC and BeO have about three times larger minimal thermal conductivity than ceria and toria and comparable to that of BP, while BAs has smaller minimal thermal conductivity as predicted by Eq. 2. Table 1b shows that in general the lattice origin thermal conductivity is larger for compounds with a larger Debye temperature (θ), larger Young's modulus (and related minimal thermal conductivity) and smaller thermal expansion coefficient (α).This observation may be used as a guideline in searching for higher thermal conductivity materials.In general the thermal conductivity is higher for compounds with lighter atoms (average atomic mass is listed in Table 1a, column two).However, the thermal conductivity of thoria and ceria, are almost the same, while the average mass of atom in thoria is significantly higher (88 u versus 57 u).Slack also noted the correlation between the lattice thermal conductivity value and the Debye temperature [5].We used here generalised Slack model, where we include contribution to the thermal conductivity from both acoustic and optical phonons and all inputs to the model are calculated as described previously [4,10] and above.Slack model predicts large thermal conductivity for both BP and BAs and BP may be considered for a future application provided it has a proper isotopic contribution that is not absorbing neutrons.

Table 1a. The properties of selected compounds with the average atomic mass (in u) per unit cell listed in column two.
The lattice thermal conductivity [Wm -1 K -1 ], calculated at 400 K (column 5) using Eq. 1, and the respective a and b parameters (in W -1 mK and W -1 m units respectively) are shown in columns 3-4.

CMPSE2017
The solution of the equation quantifying steady state heat conduction for cylindrical symmetry (radius of the fuel rod equal to 0.006 m is assumed) provides the temperature profiles (T) and temperature gradients (dT/dr) in an operating fuel rod as a function of radius (r).We investigate various nuclear fuels: thoria, MOX, and the composites with one order of magnitude higher thermal conductivity constituents like BeO and BP.
The temperature profiles for various composites in the shape of a cylindrical fuel rod were calculated using the thermal conductivity evaluated by Eq. 4 and by solving the steady state heat conduction equation for cylindrical symmetry: (4) where "r" is the radial coordinate of a cylinder and "H" is the volumetric heat generation rate from fission as described previously [16].Eq. 4 is subject to the boundary conditions: zero temperature gradient at the centerline of the pellet and the fixed temperature value at the surface of the pellet, assumed to be the same as used before [17]: 870 K.  1a).
We also assume that the used fuel enrichment, which leads to the same linear power (48 kWm -1 ) does not affect thermal conductivity significantly and therefore can be neglected.
In Figs.5a,b, we present the respective solutions by Maple code in a fuel pellet with thoria fuel and its composite with BeO operating for one day with linear power of 48 kWm -1 .
We assumed that the fuels were fully dense and the performance of pure thoria was compared with 5vol%, 10vol% and 15vol% BeO/ThO2 composites.Replacement of 5 vol% BeO in thoria by SiC gives very similar, about 200 K centreline temperature reduction and similar reduction of temperature gradient [18,19].
The presented in Figure 6a centreline temperature reduction by replacing 5 vol% of MOX by BP predicts much higher reduction (1000 K) due to both lower thermal conductivity of MOX, which gives much higher centreline temperature for pure MOX (2200 K versus 1600 K) and predicted higher thermal conductivity of BP than that of BeO.  1a).
As shown in Figs. 5 and 6 the composite fuels using high thermal conductivity materials allow a significant reduction of the central fuel line temperature and the temperature gradient within the fuel.Therefore composite fuels would not only be safer, but would have higher longevity.
We have not included it here in our analysis but BN may be even better component of composite than BP since N is lighter than P and therefore has higher thermal conductivity [15].

Conclusion
The application of the discussed software enabled us to derive at the lattice thermal conductivity of various nuclear fuels using the thermo-mechanical properties obtained from first principles calculations.Slack model and the ShengBTE (BTE) a solver for phonon thermal conductivity predict similar thermal conductivities in agreement with the experiment except for cubic SiC where Slack model over-predicts thermal conductivity by a factor of two and BAs where BTE predicts order of magnitude higher thermal conductivity than Slack model and available currently experimental result.The studied here Borides (with the appropriate isotopic composition that excludes strong neutron absorber 10 B) like BAs and BP are predicted to have high thermal conductivity and BN may be even better component of composite since N is lighter than As and P. Furthermore, we calculated temperature profiles and gradients using a simplified, exemplary scenario of an operating fuel rod.
Our studies also help to identify important factors that contribute to enhancing lattice thermal conductivity of nuclear fuel and its implications.The findings are helpful in developing accident tolerant nuclear fuels.The composite fuels using high thermal conductivity materials allow a significant reduction of the central fuel line temperature and the temperature gradient within the fuel.Therefore, composite fuels would not only be safer but also have higher longevity.
The presented here multidisciplinary simulation may be implemented for educational purposes to enhance reactor safety analysis.

Figure 1 .
Figure 1.Comparison of the calculated by Slack and BTE methods lattice thermal conductivities of SiC, BeO and BAs.The shown experimental result is for BAs [12].

Figure 2 .
Figure 2. Thermal conductivity versus temperature with experimental values at various cut-off energies and k-points, and deformation number as indicated.The used strains were between 0.99 and 1.04 and 0.99 and 1.02, as indicated.The purple dash-dot curve agrees with experimental data, however, at a very low cut-off energy of 50 Ry.The best, desired curve is the black solid curve at 15 deformations, 100 Ry cut-off energy and mesh of 8 × 8 × 8 k-points.

Figure 3 .
Figure 3.The comparison of thermal conductivity as a function of temperature with experimental values [13, 14] at various deformation numbers, cut-off energies, and k-points.The best, selected curve is the black-dashed plot at 15 deformations.The maximum strain has been fixed to 0.99-1.04.The Slack model predicts the thermal conductivity of BAs in agreement to the experimental study by Bing et al. [12].However, the most interesting fact is that the thermal conductivity presented by Lindsay et al. [14], and ours own simulations (indicated by dot dashed red line in Fig.1) using BTE predicted the value an order of magnitude higher.Though our BTE results for BAs (Fig.1) predict lower thermal conductivity than published later for a larger temperature range by the same authors (4288 Wm -1 K -1 versus ~10000 Wm -1 K -1 at 100 K)[15].The calculated here thermal conductivity at 100 K for BP using Slack method is between 3000 Wm -1 K -1 and 4000 Wm -1 K -1 as shown in Fig.3, while it was estimated previously[15] to be also much higher ~10000 Wm -1 K -1 .The BTE calculations for BAs were calculated using the same PAW-PBEsol pseudopotentials and optimal parameters as for Slack model as described above and the same PBEsol functional.In the previous BTE calculations, local density approximation (LDA) was used with norm conserved pseudopotentials and only 80 Ry cut-off energy and 6 × 6 × 6 k points mesh, which we found not sufficient for the used here PAW-PBEsol pseudopotentials.Our comprehensive investigation indicates that the predicted lattice thermal conductivity by both BTE and Slack model for some materials may be very sensitive to the tuning optimization of parameters used in DFT calculation, selected type of DFT functional and the respective pseudopotentials.In Figure4we present calculated here phonon spectra for BP and BAs as indicated.