- Email: [email protected]
An atomistic modeling of the xenon bubble behavior in the UO2 matrix
Accepted Manuscript An atomistic modeling of the xenon bubble behavior in the UO2 matrix A. Jelea, R.J. -M. Pellenq, F. Ribeiro PII: DOI: Reference:
S0022-3115(13)01118-5 http://dx.doi.org/10.1016/j.jnucmat.2013.09.041 NUMA 47761
To appear in:
Journal of Nuclear Materials
Received Date: Accepted Date:
2 January 2013 20 September 2013
Please cite this article as: A. Jelea, R.J.. Pellenq, F. Ribeiro, An atomistic modeling of the xenon bubble behavior in the UO2 matrix, Journal of Nuclear Materials (2013), doi: http://dx.doi.org/10.1016/j.jnucmat.2013.09.041
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
An atomistic modeling of the xenon bubble behavior in the UO2 matrix A. Jeleaa,*, R. J. -M. Pellenqb,c,d and F. Ribeiroa a
Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN/SEMIA/LPTM, BP3, 13115 Saint-Paul-lez-Durance, France
b
Centre Interdisciplinaire des Nanosciences de Marseille, CNRS, Campus de Luminy, Marseille, 13288, France c
Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, 77 Massachusetts avenue, Cambridge, 02139, MA, US
d
2, The joint CNRS-MIT Laboratory, 77 Massachusetts avenue, Cambridge, 02139, MA, US
Abstract The behaviour of the xenon nanoinclusions/bubbles in the uranium dioxide (UO2) matrix and their influence on its swelling were investigated through atomistic simulation techniques. The pressure in bubbles of less than 2 nm in diameter, calculated using a virial equation that takes into account the xenon/matrix interactions, is larger than the pressure calculated in simulations of the equivalent density and temperature of super critical bulk xenon. The radial distribution function of confined xenon is characteristic of a dense ( > 4 g/cm3) glassy phase. The swelling of the UO2 induced by the intragranular bubbles is proportional to the Xe/U ratio but independent of the temperature.
Keywords: uranium dioxide, molecular dynamics, surface, porous material, nanoconfinement
*
Corresponding author. Tel.: +33 4 42 19 97 28 ; fax: +33 4 42 19 91 66
E-mail address : [email protected] (A. Jelea) 1
1. Introduction Uranium dioxide (UO2) is one of the standard nuclear fuels in fission reactors. The behaviour of this particular solid fuel and the impact of irradiation on its properties are the focus of numerous studies. The fission of uranium isotopes under irradiation produces noble gas species such as xenon (Xe) and krypton (Kr) at a rate of one noble gas atom for every three or four fission events, xenon being seven times more frequently produced than krypton [1]. These noble gases have a very low solubility in the UO2 matrix and consequently they precipitate into the fuel matrix as intragranular and intergranular nanoinclusions/bubbles, or are released outside [2]. The bubble formation into the fuel matrix alters its thermomechanical properties and can lead to an important swelling. Moreover, rapid and large variations of temperature, like those occurring during accidents, can lead to sudden expansion of the confined fluid [3] and consequently to fuel fragmentation. The fuel properties alteration, taking place during the normal functioning of the reactor, can be mostly viewed as a performance limiting factor. By contrast, in accident conditions they can have dramatic consequences for the safety of nuclear devices. This is the reason why numerous experimental [4-7] and theoretical studies [8-17] are devoted to understanding the impact of noble gas nanoinclusions on the fuel properties, in particular their thermomechanical properties. The experimental studies revealed the existence of small bubbles of 1-3 nm in diameter in fuels with a low burnup (when less than 2% of the uranium atoms underwent fission) [2]. These small bubbles can then coalesce at higher burnup (when more than 4% of the uranium atoms underwent fission) to lead to a bimodal bubbles population, with the largest bubbles of tens of nanometers in diameter [2]. Interestingly, xenon density in these bubbles (3.8 to 6 g/cm3) is characteristic to solid xenon [5]. Usually, the fission gas behaviour is modelled either by the van der Waals equation of state 2
or the ideal gas law, which have been found to be inadequate to describe high density xenon nanoinclusions. As a consequence new equations of state for Xe under high pressure and high temperature have been recently developed [18]. Some authors [10,13] evaluated the pressure in the bubble by comparing to the pressure obtained in a simulation of an equivalent density and temperature of the bulk noble gas, thus neglecting the gas/matrix interface effect (confinement effect). As underlined by some authors [19], the confinement effect can play a crucial role. The present work is an atomistic simulation study on the gas/matrix interaction when xenon nanoinclusions/bubbles are trapped in UO2 at temperatures characteristic to this nuclear fuel during the normal functioning of the reactor (1000 K) or during accidents (2000 K). This is achieved using classical molecular dynamics simulations. The paper is structured as follows: Section 2 gives some accounts of the simulation approach. In the first part of Section 3, the effect of the matrix on the confined xenon is studied. The amount of fuel swelling as a result of the intragranular gas bubbles is calculated and discussed in the second part of Section 3. All results are summarized in Section 4.
2. Models and methods 2.1. Models for the bubble/matrix systems All our calculations were carried out using 3D periodic boundary conditions. The pattern of the bubble/UO2 model is built in a very similar manner as that described in reference [20], starting with a 6×6×6 supercell of perfect UO2. An almost spherical cavity with the radius of about 0.45 nm is created in the middle of this supercell by removing a stoechiometric number of uranium and oxygen atoms. The resulting system has 854 uranium atoms and 1708 oxygen 3
atoms. A number of 15, 26, 35, 43, 50, 60 and 70 xenon atoms were introduced in the cavity using a dense body centred cubic (bcc) Xe arrangement, then the system was relaxed through molecular dynamics. The Xe/U ratios corresponding to these cavity fillings are: 0.018, 0.030, 0.041, 0.050, 0.059, 0.070 and 0.082 respectively. Usually, the gas bubbles distort the zone in the matrix that is situated close to the bubble/matrix interface. The supercell size was defined in order to avoid the interaction between the distorted zone surrounding the bubble and its corresponding images. Thus, there is an unperturbed matrix zone of at least 0.7 nm between the bubble and its images.
2.2. Equilibrium structures The properties of the systems in the present work are calculated at equilibrium. The equilibrium configurations are obtained through Newton-Raphson energy optimization and molecular dynamics methods [21]. Energy optimization allows finding the local minima closest to the starting configuration on the potential energy surface (PES). Molecular dynamics helps exploring larger parts of the PES by propagating Newton’s equations of motion through time using finite difference formalism. At 0 K, the quantity that is minimized during the energy optimization is the potential energy. To take into account the temperature impact on the equilibrium structure, a natural approach is to consider free energy minimization. The free energy of a structure at a given temperature is calculated through the determination of the vibrational partition function from the phonon density of states, within the quasi-harmonic approximation. This is usually assumed to be correct for temperatures lower than half of the melting point temperature of the considered material [22]. 4
For all our molecular dynamics simulations the timestep is 1 fs. The temperature and the pressure are kept at constant values using the Nose-Hoover thermostat and barostat algorithms [21]. The calculations are performed using the GULP 3.1 simulation package [22], which implements all the aforementioned methods.
2.3. Interatomic potentials Within our bubble/matrix systems there are three types of interactions: the in-matrix U/U, U/O and O/O interactions, the xenon/matrix Xe/U and Xe/O interactions and the Xe/Xe interactions. In order to describe the in-matrix interactions the potential proposed by Basak [23] was selected. Indeed, a systematic study carried out by Govers over 20 different potentials identified Basak’s potential [8] as one of the best, more particularly at high temperature, and Geng et al. [12] showed that this potential is transferable to structures far from the UO2 fluorite structure used to define its parameters. Moreover, a recent study [20] shows this potential to be able to correctly describe the thermomechanical properties of nanoporous UO2. The analytical expression of the Basak potential is given by the equation 1 : ⎛
rij ⎞ Cij zi z j e 2 ⎟− 6 + + Dij ⎧⎨ 1 − exp β ij rij* − rij ⎟ r ρ r ⎩ ij ⎠ ij ij ⎝
φij = Aij exp⎜⎜ −
[
( (
))]
2
− 1⎫⎬ ⎭
(1)
This potential is a sum of three terms: a Buckingham-Hill term describing the van der Waals interactions, a Coulomb term describing the ionic interactions and a Morse term for the
5
covalent bonds, as in the UO2 the U-O bonds have an ionic–covalent character. The parameters of the Basak potential are given in table 1. For the xenon/matrix interaction, we used the potential proposed by Grimes [9]. This potential has a Buckingham-Hill form. The CXeU and the CXeO parameters (the dispersion part parameters) were calculated by Grimes using the Kirkwood-Slater formula as described in reference [24], and the repulsive part parameters A and ρ were calculated using an electrongas method [25]. The Grimes potential was initially proposed in a shell-core form. In the present work, it is used in a rigid ion form since this one is far more effective for molecular dynamics simulations [8]. In the shell-core form smaller timesteps are needed and the number of interactions to compute increases by a factor of 4. Since the Buckingham-Hill potential has a non-physical zone of attraction at small distances and the present calculations are carried out for very dense systems and high temperatures, there is a risk to explore this non-physical zone. Therefore, the Buckingham-Hill form was replaced at short distance with a Born-Mayer repulsive potential. Thus, the potential has the following form: ⎛
rij ⎞ ⎟ BM ⎟ ρ ij ⎝ ⎠
φij = AijBM exp⎜⎜ −
if
rij < rc
(2)
and ⎛
rij ⎞ Cij ⎟− 6 ⎟ r ρ ij ⎠ ij ⎝
φij = Aij exp⎜⎜ −
if rij ≥ rc
Here rc is the inflexion point of the Buckingham-Hill function situated between r = 0 and the equilibrium distance r = rmin. The ABM and ρBM parameters of the Born-Mayer function were calculated from the conditions of continuity and derivability of the potential in rc. The Grimes potential choice for our study was made after it was assessed together with 6
three other largely used semiempirical potentials: the Geng potential [12], the Jackson-2 potential [8,26] and the Chartier potential [14]. We used some quantum DFT calculation results from literature [14-17] for comparison. The Xe incorporation energy in different sites of the UO2 matrix is one of the most appropriate assessment parameters. The considered sites can be easily identified on figure 1 representing the UO2 elementary cell : the octahedral interstitial site situated in the middle of the cube formed by the oxygen atoms, the uranium and the oxygen vacancies obtained by displacing an uranium or an oxygen far from their perfect crystal positions to form the Frenkel pairs, the UO divacancy obtained by displacing both an uranium and the nearest neighbour oxygen and the UO2 trivacancy (Schottky defect) obtained by removing an uranium atom and two nearest neighbour oxygen atoms. There are three trivacancy configurations : 145 (I), 124 (II), 123 (III). The incorporation energy was calculated at 0 K as the difference between the internal energy of the relaxed systems containing xenon and the internal energy of the xenon free systems. The patterns of the periodic models representing the xenon free systems and the systems containing xenon were built based on a 4x4x4 supercell of UO2 fluorite structure. The in-matrix interactions were calculated using the Basak potential. Table 2 gives the semiempirical potential as well as the DFT values for the incorporation energies. All the semiempirical potentials give the good order of the incorporation energies (as compared to DFT) with the interstitial site having the highest incorporation energy and the configurations II and III of the trivacancy having the lowest ones. The Grimes potential seems to give an overall better description than the other potential models. For the xenon/xenon interactions three semiempirical potentials were assessed: a Buckingham-Hill potential used in reference [27], a Lennard-Jones potential also from reference [27] and a potential proposed by Tang and Toennies [28]. In the Buckingham-Hill
7
potential the short distance part was replaced with a Born-Mayer repulsive potential in the same manner and for the same reason as for the Grimes potential. As a result, the potential has also the form given by equation 2. Some experimental results [29,30] were used for comparison. The first part of the assessment consisted in comparing the predictions of the density versus pressure ρ(P) curve at room temperature given by the three potentials with the experimental results (figure 2). The ρ(P) curve calculated with Buckingham-Hill is the only one to follow the experimental points for all the considered pressure range. The other two potentials produce curves that move away from the experimental points [29] at larger pressures. In the second part of the assessment, the bulk modulus evolution versus pressure K(P) is calculated at room temperature and compared with the experiment (figure 3). In the lower part of the pressure range (from 3 GPa to 6 GPa) the Lennard-Jones potential gives the best results, but for larger pressures it fails to account for the experimental points [30]. The Buckingham-Hill values seem again to best approach the experimental data in the considered pressure range. The ρ(P) curve was calculated through molecular dynamics simulations in NPT ensemble on a supercell of 500 xenon atoms. Within this supercell a solid/liquid interface was created in order to let the system the freedom to evolve toward the most stable phase for a given (T,P) set. 30 ps were used for the equilibration and another 30 ps were used to accumulate the results. At room temperature (298 K) and the considered pressure range from 5 to 50 GPa the system evolves toward the face centred cubic (fcc) structure. This result is in agreement with the phase diagram from reference [31]. The K(P) curve was calculated through free energy optimization on fcc xenon at 296 K and the pressure from 3 to 10 GPa. The free energy 8
optimization technique can be safely used since the considered temperature is less than half the melting point temperature in this pressure range [31]. All the parameters of the xenon/matrix and xenon/xenon potentials used in this study are presented in table 3.
3. Results and discussion 3.1. The characterization of the xenon in the bubble Estimating the state of the xenon confined in the form of bubbles in UO2 as well as the pressure in the bubble are a challenging task from both experimental [4-7] and theoretical [12,13] viewpoints. High resolution transmission electron microscopy (HRTEM) and energy dispersive X-ray spectroscopy (EDX) analyses, carried out in the outer region of high burnup UO2 pellets [4,5], showed the existence of almost spherical xenon inclusions/bubbles with diameters ranging from 1 nm at low burnups to 10 nm at high burnups. The xenon in the bubbles seems to be in a solid or in a near solid state with densities from 3.8 to 6.0 g/cm3 (for bubbles with diameters ranging from 4 to 10 nm). For these bubble radii and densities, at 700 K, the pressure in the bubble was estimated at 1.6 to 15 GPa. X-ray absorption spectroscopy (XAS) experiments [6,7] carried out on a set of xenon implanted UO2 samples indicated the existence of highly pressurized xenon inclusions/bubbles. The pressure in the bubble was estimated at 2.1 to 2.7 GPa at 5 K and at 3.0 to 3.7 GPa at 870 K for bubbles of approximately 2 nm in diameter. In these experimental studies the pressure in the bubble was estimated by equating to the pressure of the bulk noble gas at an equivalent density and temperature. The same method of estimating the pressure in the bubble was employed in theoretical works [10,13]. However this approach might not be very appropriate given the existence of the 9
xenon/matrix interactions. Moreover, the bubble size, on the order of a few nanometers, helps to increase these interactions. In the present work the pressure in the bubble is calculated “in situ”, taking into account not only the xenon/xenon interactions but also the xenon/matrix ones. This is achieved using the virial equation for pair-additive forcefields as recommended in reference [32] in the case of periodic boundary conditions. A local form of this equation can be obtained using as a starting point the local stress tensor. The σαβ component of the local stress tensor σ corresponding to a sub volume V inside of a material (an inclusion in a matrix for example) is calculated using the equation 3, where the summation runs over all the atoms i contained in the sub volume V [33] :
σ αβ =
1 1 mi viα viβ + ∑ V i∈V 2V
⎛
1 ∂U ij ∂rij
∑∑ ⎜⎜ − r i∈V j ≠ i
⎝
⎞α β ⎟rij rij ⎟ ⎠
(3)
Here mi is the mass of the atom i, α and β are the Cartesian xyz components, viα is the α component of the atom i velocity vector, rijα is the α component of the vector joining the atoms i and j and rij represents its norm. U is the interatomic pair potential. The summation in the second term on the right side runs over all the atoms j inside the range of the interatomic potential U. Thus, the atoms j may be contained in the sub volume V or may belong to the rest of the system surrounding this sub volume. The pressure in the bubble Pbubble is calculated as the average over time of one third of the trace of the local stress tensor σ[34]. Using this definition for the pressure and the equipartition theorem to relate the instantaneous system temperature T with the average over all the degrees of freedom of the kinetic energy, one finally obtains the virial equation that gives the pressure
10
in the bubble:
Pbubble =
1 N Xe kT + V 3V
⎛ ⎞ 1 ⎜ ∑ rij f ij + ∑∑ riL f iL ⎟ ⎜ ⎟ 2 i L ⎝ i< j ⎠
(4)
The first term on the right-hand side of the equation 4 is the kinetic (ideal) term whereas the second one is the potential (non-ideal) term. The non-ideal term was split into two parts: a xenon/xenon part and a xenon/matrix one. Hence rij and fij are the distance and the force acting between the xenon atoms i and j whereas riL and fiL are the distances and forces acting between the xenon atom i and the matrix atom L (oxygen or uranium). As the model is periodic, the sums run on all over the atoms in the cell as well as the image atoms from the neighbouring cells. T is the instantaneous temperature, k is the Boltzmann’s constant, NXe is the number of the xenon atoms in the bubble and V represents the instantaneous volume of the bubble. There are multiple ways for estimating the instantaneous volume of the bubble. The reference [35] presents a method, based on the tessellation of space into Voronoi and Delaunay polyhedra, that allows the exact calculation of volumes of cavities in sphere packings in three dimensions. A less computationally expensive way to estimate the bubble volume assumes that the bubble is spherical its radius being calculated [36] as the distance from the center of the bubble to the nearest uranium or oxygen atoms. In the present work we used the spherical shape hypothesis since high resolution transmission electron microscopy (HRTEM) images [4,5] show the intragranular bubbles to have almost spherical shapes. However, at the atomic scale there are irregularities in the bubble shape. For example some matrix atoms (oxygen and/or uranium) can be found closer to the bubble mass center than some xenon atoms. In order to take into account, to some extent, this fact we calculated the bubble radius as the arithmetic mean between a “cavity minimum radius” (Rca) and a “cluster
11
maximum radius” (Rcl). The “cavity minimum radius” is defined as the distance between the bubble mass center and the closest matrix atom (oxygen or uranium) and the “cluster maximum radius” as the distance between the bubble mass center and the furthest (xenon) atom in the xenon cluster. The figure 4 shows the two possible situations that can occur when calculating the bubble radius with this method: some matrix atoms (oxygen and/or uranium) are closer to the bubble mass center than some xenon atoms (left side of figure 4) and all the xenon atoms are closer to the bubble mass center than the matrix atoms (right side of figure 4). Given the relative sizes of the polyhedron closing the xenon cluster and the Rcl and Rca radius spheres and the fact that the atoms have finite volumes, in the former case our method might slightly underestimate the bubble instantaneous volume whereas in the latter it might slightly overestimate it. An instantaneous pressure was calculated at each configuration using the instantaneous temperature and the instantaneous bubble volume, then the pressure in the bubble was calculated as the average of these instantaneous pressures. We verified if a sufficient number of data points was used for the pressure average using the so called “block” average method [37,38]. Within this method one starts by calculating the standard error of the mean of all the data points, then two consecutive data points are grouped and the standard error of the mean of this new block of data is calculated. The new block of data contains half the number of data points of the initial set. The operation is repeated until there are not enough data points to compute a standard error. The original data set has highly correlated data points since the time between two data points (the timestep) is too short. As a consequence the standard error has a low value. By increasing the block size the data points are less and less correlated and the standard error of the mean increases with the number of block operations until so many data points are grouped that two consecutive points are uncorrelated. At this moment the standard
12
error reaches a plateau. Thus a plateau (figure 5) on the standard error of the mean versus number of block operations curve was used as an indicator that we have obtained enough statistics. Molecular dynamics simulations in NPT ensemble were carried out on systems with different numbers of xenon atoms, built as described in section 2.1, for temperatures of 1000 K and 2000 K. The pressure was set at 0 GPa because this facilitates the comparison with other theoretical results from the literature that were obtained using this pressure. Each system at 2000 K was relaxed during 80 ps. The systems were then “cooled” at 1000 K through two intermediate steps: one at 1700 K and the other at 1300 K. At each of these two temperatures the system was relaxed during 60 ps. At 1000 K the systems were relaxed again during 80 ps. This method was used because the systems relax faster at 2000 K than at 1000 K, thus a larger number of relaxation steps might be needed to obtain equilibrium structures at 1000 K. The pressure in the bubble as a function of the number of the xenon atoms in the bubble at 1000 K and 2000 K was then calculated by averaging over 40000 data points. Globally the block analysis showed this number of data points to be sufficient, but at 2000 K, for some cases, no plateau was found in the standard error versus number of block operations curve (see figure 5). However, the pressures calculated by averaging over a larger set of data points (140000 data points), for bubbles with NXe = 15, 43 and 70 (covering all the analyzed NXe range), were similar to those obtained using 40000 data points (see figure 6 bottom). We thus concluded that the pressure versus NXe curves were calculated with a sufficient accuracy. The pressure in the bubble was also calculated by equating to the pressure of the bulk xenon at an equivalent density and temperature. The results of the two methods are compared in figure 6. Both methods predict that the pressure in the bubble rapidly diminishes with increasing NXe until it reaches a minimum, then slowly increases. The pressures calculated “in 13
situ” are different (higher) from those corresponding to the bulk xenon. In reference [36] a similar result, but in terms of free energy, was obtained: the free energy of the xenon in a cavity is higher than the free energy of the bulk xenon at the same density and temperature for Xe densities of 3.8-6 g/cm3 and temperatures of 300-1050 K. The pressure in the bubble estimated with the ideal gas law is also represented in figure 6. It is much lower (< 1 GPa) than the pressures calculated with the two other methods that take into account the non-ideal terms. Moreover, it becomes clear from figure 6 that the difference between the “in situ” and the bulk xenon estimates of the pressure can be as important as the difference between the real gas and the ideal gas equation of state estimates. The difference between the ‘”in situ” and the bulk xenon pressure is a xenon/matrix interface effect. For large bubbles, where the surface interactions are far less important than the volume ones, one can expect the results of the two estimates to converge. As shown in figure 7 (top) the bubble radius is also a function of the number of xenon atoms in the bubble: the bubble radius increases with NXe but this increase slows down at high NXe. This is because, during the molecular dynamics relaxation, the surface that closes the
bubble evolves until the equilibrium is reached between the xenon pressure in the bubble on one side and the surface tension (attributed to the xenon/matrix interface) plus the external pressure on the other side. At the atomic scale, the xenon atoms push uranium and oxygen atoms into the matrix creating more space for themselves. This phenomenon was already observed by other authors. In reference [13] xenon atoms were introduced in UO2, one by one, starting from the perfect UO2 lattice. After adding each new xenon atom the resulting system was relaxed through molecular dynamics. At the end of each relaxation the noble gas atoms form bubbles of different sizes (function of the number of xenon atoms in the bubble) by displacing the oxygen and uranium atoms at the xenon/matrix interface. In reference [36]
14
the xenon atoms were introduced into an already built cavity into the UO2 (like in the present work) followed by a molecular dynamics relaxation of the system. The authors observed that, for sufficiently large numbers of xenon atoms in the bubble, the xenon atoms at the surface of the bubble created more space for themselves by forcing oxygen ions into the UO2 lattice. An equation relating the intragranular bubble radius (Rbubble) to the number of xenon atoms in the bubble (NXe) was proposed in reference [39]:
Rbubble = BN
1
1 3 Xe
where
⎛ 3 ⎞3 B = ⎜⎜ ⎟⎟ ⎝ 4πρ ⎠
(5)
Here ρ is the density of xenon atoms in the bubble in nm-3. As can be observed in figure 7 (bottom), the xenon density in the bubble, at both 1000 K and 2000 K, is a function of NXe, following the same evolution as the pressure: it decreases towards a minimum then increases at higher NXe. Considering that the ρ(NXe) curves can roughly be described by a third order polynomial and taking into account the equation 5, one obtains for Rbubble(NXe) the following equation: 1
Rbubble =
0.62035 N Xe 3
(aN
3 Xe
2 + bN Xe + cN Xe + d
)
1 3
(6)
where a = -0.00026, b = 0.0449, c = -2.3866 and d = 58.96. In figure 7 (top) we show that the equation 6 gives a good approximation of the molecular dynamics points. For almost all the cases we analyzed the density in the bubble is higher than 4 g/cm3, the
15
density of the solid face centred cubic (fcc) xenon [5]. This leads to the question of the state of the xenon confined in nanocavities at high temperature and pressure. The radial distribution function (RDF) shape analysis can give an insight into the state of the confined xenon. This method was already successfully used by other authors to characterize the pure xenon phase diagram [31] or even the xenon state into a bubble [13]. The RDF’s for the xenon in the bubble were calculated for systems with NXe = 15 up to 70 covering all the investigated range of conditions (figure 8). RDF’s (figure 9) were also built for bulk xenon at a density of about 6.1 g/cm3 (within the range of those found in the bubbles) and temperatures of 1000 K and 2000 K. Through molecular dynamics relaxation in NVT ensemble, the bulk xenon was found to be in a solid fcc phase at 1000 K and in a liquid phase at 2000 K, in agreement with the phase diagram from reference [31]. In order to facilitate the comparison with the RDF’s for the xenon in the bubble the RDF’s for bulk xenon were calculated only for a cluster of atoms found in a 0.6 nm radius sphere. The value of the sphere radius was taken to be representative for the bubble sizes analyzed in this study. The RDFs of the xenon in the bubble (figure 8) are very similar with the ‘liquid’ cluster RDF shown in the figure 9. Taking into account this observation and the fact that the xenon density in the bubble corresponds to the solid xenon, one can safely assume that the xenon in the nanocavities (bubbles), in the conditions analyzed here, is in a glassy phase, in spite of the fact that in some cases the (T,Pbubble) sets correspond to the solid fcc state on the bulk xenon phase diagram [31]. This is due to the fact that the xenon/matrix interaction alters the phase diagram of the xenon in the bubble through a confinement effect. In reference [36] a similar conclusion is reported : within the bubbles, the surface of the uranium dioxide restricts the rearrangement of the Xe atoms into the perfect fcc structure thus retarding the phase change to a solid fcc phase. 16
3.2. The calculation of the bubble induced swelling In this section we analyze the effect of the xenon/matrix interaction on the global behaviour of the nuclear fuel. A very important issue to be addressed is the nuclear fuel swelling during the normal functioning of the nuclear reactor or accidental situation. The swelling is the sum of three components: A swelling due to the UO2 matrix thermal expansion; A swelling due to the UO2 matrix fragmentation; A swelling due to the fission product accumulation inside the matrix, particularly the noble gases, in the form of bubbles. The last component, the bubble induced swelling, is supposed to become very important at high temperature. So it is crucial to quantify this swelling term and its thermal dependency. The bubble induced swelling was calculated using the equation 7 as given in the reference [13]:
S% =
V eq − V0eq × 100 V0eq
(7)
Here Veq is the equilibrium volume of the bubble/matrix supercell and V0eq is the equilibrium volume of the perfect UO2 supercell with the same number of UO2 groups as the bubble/matrix supercell. The equilibrium volumes were calculated through molecular
17
dynamics simulations in NPT ensemble for temperatures from 800 K to 2000 K and for a pressure of 0 GPa (no external stress). The swelling versus Xe/U ratio at 1600 K is plotted in figure 10 and compared with the results from reference [13]. The agreement between the results of the two studies is remarkable since in reference [13] a different set of interatomic potentials was used. At constant temperature the bubble induced swelling increases linearly with the Xe/U ratio. In order to verify if the system size has an effect on the calculated swelling we computed the swelling corresponding to a larger system. This system was built starting with a 7×7×7 supercell of perfect UO2. An almost spherical cavity with the radius of about 0.58 nm was created in the middle of this supercell by removing 16 UO2 units. The resulting system has 1356 uranium atoms and 2712 oxygen atoms. The cavity was filled with 35 and 70 xenon atoms as described in section 2.1. The corresponding Xe/U ratios are 0.0258 and 0.0516 respectively. The system was relaxed for 80 ps through molecular dynamics in NPT ensemble at T = 1600 K and P = 0 GPa. The swelling values calculated for these larger systems integrate well in the S%(Xe/U) curve from figure 10 showing that the system size has practically no effect on these results. The swelling versus temperature corresponding to three Xe/U ratios at P = 0 GPa is plotted in figure 11. It seems that the temperature has no effect on the bubble induced swelling. However, this conclusion should be restricted only to the intragranular bubbles since only these ones are analyzed here. Actually, the temperature dramatically increases the gaseous swelling but mainly due to the intergranular bubbles [40]. The intergranular bubbles grow faster and larger at high temperatures because the point defects and the fission gas atoms from the fuel grains can diffuse faster to the grain boundary and be trapped by the intergranular
18
bubbles.
4. Conclusions We have presented an atomistic simulation study on the behaviour of the intragranular xenon bubbles trapped in the UO2 matrix and their influence on the fuel swelling. The main conclusions that can be derived from this study are: 1) The xenon pressure in the bubble calculated “in situ” by taking into account the xenon/matrix interaction is systematically larger than the pressure calculated in simulations at an equivalent density and temperature of the bulk xenon. However, this discrepancy should vanish for large bubbles as a consequence of the surface (confinement) effect decrease with decreasing surface/volume ratio. Thus, the equations of state established for free xenon, whereas still appropriate for large bubbles, cannot be safely used to describe the xenon in bubbles of a couple of nanometres in diameter. 2) The bubble radius and the pressure in the bubble are functions of NXe, the number of xenon atoms confined in the bubble. The radius evolution is fast at low NXe but slows down at higher NXe. 3) The radial distribution function analysis is indicative of a glassy phase for the xenon in the bubble even for temperature/pressure points corresponding to solid fcc xenon on the bulk xenon phase diagram. This is also an effect of the xenon/matrix interface. However, a more detailed analysis would require the calculation of the diffusion coefficient of the confined xenon and a comparison with the corresponding bulk xenon ones. This will form the object of a future work dedicated to the state of the noble gases confined in UO2.
19
4) The bubble induced swelling of the fuel increases linearly with the quantity of xenon in the system (Xe/U ratio). However, for a given Xe/U ratio, the temperature doesn’t seem to have an important impact on the intragranular swelling in the 800K to 2000K temperature range.
Acknowledgements. The calculations were performed at the facilities of the CNRS/CINaM laboratory and IDRIS, the CNRS computing center. We would like to thank G. Tréglia from CNRS/CINaM for his useful suggestions.
20
References [1] J. Rest, J. Nucl. Mater. 402 (2010) 179-185. [2] P. Lösönen, J. Nucl. Mater. 280 (2000) 56-72. [3] K. Govers, C.L. Bishop, D.C. Parfitt, S.E. Lemehov, M. Verwerft, R.W. Grimes, J. Nucl. Mater. 420 (2012) 282-290. [4] K. Nogita, K. Une, J. Nucl. Mater. 250 (1997) 244-249. [5] K. Nogita, K. Une, Nucl. Instr. and Meth. in Phys. Res. B 141 (1998) 481-486. [6] P. Garcia, P. Martin, G. Carlot, E. Castelier, M. Ripert, C. Sabathier, C. Valot, F. D’Acapito, J.L. Hazeman, O. Proux, V. Nassif, J. Nucl. Mater. 352 (2006) 136-143. [7] P. Martin, P. Garcia, G. Carlot, C. Sabathier, C. Valot, V. Nassif, O. Proux, J.L. Hazemann, Nucl. Instr. and Meth. in Phys. Res. B 266 (2008) 2887-2891. [8] K. Govers, PhD Thesis, Université libre de Bruxelles, 2008. [9] R.W. Grimes, C.R.A. Catlow, Phil. Trans. R. Soc. Lond. A 335 (1991) 609-634. [10] D.C. Parfitt, R.W. Grimes, J. Nucl. Mater. 381 (2008) 216-222. [11] D.C. Parfitt, R.W. Grimes, J. Nucl. Mater. 392 (2009) 28-34. [12] H.Y. Geng, Y. Chen, Y. Kaneta, M. Kinoshita, J. Alloys Compd. 457 (2008) 465-471. [13] T.X. Feng, L.C. Sheng, Z.Z. He, G. Tao, Chin. Phys. B 19 (2010) 057102. [14] A. Chartier, L. Van Brutzel, M. Freyss, Phys. Rev. B 81 (2010) 174111. [15] G. Brillant, F. Gupta, A. Pasturel, J. Nucl. Mater. 412 (2011) 170-176.
21
[16] A.E. Thompson, C. Wolverton, Phys. Rev. B 84 (2011) 134111. [17] H.Y. Geng, Y. Chen, Y. Kaneta, M. Kinoshita, Q. Wu, Phys. Rev. B 82 (2010) 094106. [18] J.Y. Oh, Y.H. Koo, J.S. Cheon, B.H. Lee, D.S. Sohn, J. Nucl. Mater. 372 (2008) 89-93. [19] C. Alba-Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K.E. Gubbins, R. Radhakrishnan, M. Sliwinska-Bartkowiak, J. Phys.: Condens. Matter 18 (2006) R15-68. [20] A. Jelea, M. Colbert, F. Ribeiro, G. Treglia, R.J.-M. Pellenq, J. Nucl. Mater. 415 (2011) 210-216. [21] A.R. Leach, Molecular Modelling Principles and Applications, Pearson Education Limited, Harlow, 2001. [22] J.D. Gale, A.L. Rohl, Mol. Simul. 29 (2003) 291-341. [23] C.B. Basak, A.K. Sengupta, H.S. Kamath, J. Alloys Compd. 360 (2003) 210-216. [24] P.W. Fowler, P.J. Knowles, N.C. Pyper, Molec. Phys. 56 (1985) 83-95. [25] J.H. Harding, A.H. Harker, UKAEA Harwell Report, AERE-R-10425 (1982). [26] S. Nicoll, H. Matzke, C.R.A. Catlow, J. Nucl. Mater. 226 (1995) 51-57. [27] I.R. Brearley, D.A. MacInnes, J. Nucl. Mater. 95 (1980) 239-252. [28] K.T. Tang, J.P. Toennies, J. Chem. Phys. 118 (2003) 4976-4983. [29] A.N. Zisman, I.V. Aleksandrov, S.M. Stishov, Phys. Rev. B 32 (1985) 484-487. [30] S. Sasaki, N. Wada, T. Kume, H. Shimizu, J. Raman Spectrosc. 40 (2009) 121-127. [31] A.B. Belonoshko, Phys. Rev. B 78 (2008) 174109.
22
[32] M.J. Louwerse, E.J. Baerends, Chem. Phys. Lett. 421 (2006) 138-141. [33] P.S. Branicio, R.K. Kalia, A. Nakano, P. Vashishta, F. Shimojo, J.P. Rino, J. Mech. Phys. Solids 56 (2008) 1955-1988. [34] D.H. Tsai, J. Chem. Phys. 70 (1979) 1375-1382. [35] S. Sastry, D.S. Corti, P.G. Debenedetti, F.H. Stillinger, Phys. Rev. E 56 (1997) 55245532. [36] S.T. Murphy, A. Chartier, L. Van Brutzel, J.P. Crocombette, Phys. Rev. B 85 (2012) 144102. [37] H. Flyvbjerg, H.G. Petersen, J. Chem. Phys. 91 (1989) 461-466. [38] D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, 2001. [39] D.R. Olander, D. Wongsawaeng, J. Nucl. Mater. 354 (2006) 94-109. [40] J. Spino, J. Rest, W. Goll, C.T. Walker, J. Nucl. Mater. 346 (2005) 131-144.
23
Figure 1
24
Figure 2
25
Figure 3
26
Figure 4
27
Figure 5
28
Figure 6
29
Figure 7
30
Figure 8
31
Figure 9
32
Figure 10
33
Figure 11
34
Figure captions Fig. 1. The UO2 elementary cell. It contains 4 uranium atoms (large spheres) forming a face centred cubic like arrangement and 8 oxygen atoms (small spheres) forming a cube in the middle of the cell. The numbers help to define the trivacancy configurations in section 2.3. Fig. 2. Comparison between the density versus pressure ρ(P) curves calculated at 298 K with three xenon/xenon semiempirical potentials and experimental data. The semiempirical potentials are: a Buckingham-Hill (Br-BH) potential, a Lennard-Jones (Br-LJ) potential and a potential proposed by Tang and Toennies (TT). Fig. 3. Comparison between the bulk modulus versus pressure K(P) curves calculated at 296 K with three xenon/xenon semiempirical potentials and experimental data. The semiempirical potentials are: a Buckingham-Hill (Br-BH) potential, a Lennard-Jones (Br-LJ) potential and a potential proposed by Tang and Toennies (TT). Fig. 4. Calculation of the bubble radius as the arithmetic mean between Rca and Rcl. Rca is the distance between the bubble mass center and the closest matrix atom (oxygen or uranium) and Rcl is the distance between the bubble mass center and the furthest (xenon) atom in the bubble.
Two possible situations can occur when calculating the bubble radius with this method: some matrix atoms (oxygen and/or uranium) are closer to the bubble center than some xenon atoms (left) and all the xenon atoms are closer to the matrix center than the matrix atoms (right). Fig. 5. The standard error of the mean σ for the pressure in the bubble as a function of the number of block operations NB for 40000 time steps of simulation. A plateau can be observed at NB = 7 or 8 on the curves corresponding to systems at 1000 K. An example without a plateau is given for a system at 2000 K. Fig. 6. The pressure in the bubble as a function of NXe: top – systems at 1000 K, bottom – 35
systems at 2000 K. The error bars represent the sample standard deviation. The pressures calculated “in situ” taking into account the xenon/matrix interactions are compared with those calculated for bulk xenon at an equivalent density and temperature. For systems with NXe = 15, 43 and 70 atoms, at 2000 K, the pressure in the bubble was calculated by averaging over 40000 and 140000 (better statistics) data points. The pressures at an equivalent density and temperature given by the ideal gas law are also displayed. Fig. 7. Top – the bubble radius as a function of NXe. Bottom – the density of the xenon in the bubble as a function of NXe. The error bars represent the sample standard deviation. Fig. 8. Radial distribution functions (RDF) corresponding to the confined xenon in bubbles containing 15, 43 and 70 xenon atoms: top – at 1000 K; bottom – at 2000 K. Fig. 9. Radial distribution functions (RDF) corresponding to the bulk xenon in liquid and solid states. The RDF’s for bulk xenon were calculated for a cluster of atoms found in a 0.6 nm radius sphere. Fig. 10. The swelling as a function of the Xe/U ratio at 1600 K. The results of the present work are compared to those from the reference [13]. The “smaller system” is the system derived from a 6x6x6 supercell and the “larger system” is derived from a 7x7x7 supercell. Fig. 11. The swelling as a function of temperature for different Xe/U ratios.
36
Table 1 Parameters of the Basak potential used to describe the in-matrix interactions z
A(eV)
ρ(nm)
C(eV nm6)
D(eV)
r*(nm)
β(nm-1)
U
2.4
-
-
-
-
-
-
O
-1.2
-
-
-
-
-
-
OO
-
1633.390
0.0327
3.95 x 10-6
-
-
-
UU
-
294.745
0.0327
0
-
-
-
UO
-
693.870
0.0327
0
0.57745
0.2369
16.5
37
Table 2 The incorporation energies (in eV units) of xenon at different sites in the UO2 matrix as calculated using the semiempirical potentials of Grimes, Geng, Jackson (version 2) and Chartier for the xenon/matrix interactions and the Basak potential for the in-matrix interactions. Some DFT results from literature are given for comparison. Grimes
Geng
Jackson-2
Chartier
DFT[14-17]
Interstitial
17.1
22.8
20.3
12.8
9.7-12.0
U vacancy
5.8
5.8
7.0
7.5
2.0-5.8
O vacancy
13.4
19.6
17.0
10.2
7.5-9.1
UO divacancy
4.8
5.3
6.1
6.2
1.6
UO2 trivacancy (I)
4.9
5.0
5.9
6.5
5.4
UO2 trivacancy (II)
4.3
4.9
5.6
6.0
5.2
UO2 trivacancy (III)
3.7
4.8
5.2
5.1
0.2-5.2
38
Table 3 Parameters of the xenon/matrix and xenon/xenon potentials used within this study ABM(eV)
ρBM(nm)
A(eV)
ρ(nm)
C(eV nm6)
rc(nm)
XeU
2915.33
0.0410
6139.16
0.0340
7.184 x 10-5
0.11
XeO
136.70
0.0561
598.00
0.0426
1.0838 x 10-4
0.18
XeXe
1171.48
0.0466
4934.10
0.0356
2.9747 x 10-4
0.15
+
39
S0022-3115(13)01118-5 http://dx.doi.org/10.1016/j.jnucmat.2013.09.041 NUMA 47761
To appear in:
Journal of Nuclear Materials
Received Date: Accepted Date:
2 January 2013 20 September 2013
Please cite this article as: A. Jelea, R.J.. Pellenq, F. Ribeiro, An atomistic modeling of the xenon bubble behavior in the UO2 matrix, Journal of Nuclear Materials (2013), doi: http://dx.doi.org/10.1016/j.jnucmat.2013.09.041
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
An atomistic modeling of the xenon bubble behavior in the UO2 matrix A. Jeleaa,*, R. J. -M. Pellenqb,c,d and F. Ribeiroa a
Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN/SEMIA/LPTM, BP3, 13115 Saint-Paul-lez-Durance, France
b
Centre Interdisciplinaire des Nanosciences de Marseille, CNRS, Campus de Luminy, Marseille, 13288, France c
Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, 77 Massachusetts avenue, Cambridge, 02139, MA, US
d
Abstract The behaviour of the xenon nanoinclusions/bubbles in the uranium dioxide (UO2) matrix and their influence on its swelling were investigated through atomistic simulation techniques. The pressure in bubbles of less than 2 nm in diameter, calculated using a virial equation that takes into account the xenon/matrix interactions, is larger than the pressure calculated in simulations of the equivalent density and temperature of super critical bulk xenon. The radial distribution function of confined xenon is characteristic of a dense ( > 4 g/cm3) glassy phase. The swelling of the UO2 induced by the intragranular bubbles is proportional to the Xe/U ratio but independent of the temperature.
Keywords: uranium dioxide, molecular dynamics, surface, porous material, nanoconfinement
*
Corresponding author. Tel.: +33 4 42 19 97 28 ; fax: +33 4 42 19 91 66
E-mail address : [email protected] (A. Jelea) 1
1. Introduction Uranium dioxide (UO2) is one of the standard nuclear fuels in fission reactors. The behaviour of this particular solid fuel and the impact of irradiation on its properties are the focus of numerous studies. The fission of uranium isotopes under irradiation produces noble gas species such as xenon (Xe) and krypton (Kr) at a rate of one noble gas atom for every three or four fission events, xenon being seven times more frequently produced than krypton [1]. These noble gases have a very low solubility in the UO2 matrix and consequently they precipitate into the fuel matrix as intragranular and intergranular nanoinclusions/bubbles, or are released outside [2]. The bubble formation into the fuel matrix alters its thermomechanical properties and can lead to an important swelling. Moreover, rapid and large variations of temperature, like those occurring during accidents, can lead to sudden expansion of the confined fluid [3] and consequently to fuel fragmentation. The fuel properties alteration, taking place during the normal functioning of the reactor, can be mostly viewed as a performance limiting factor. By contrast, in accident conditions they can have dramatic consequences for the safety of nuclear devices. This is the reason why numerous experimental [4-7] and theoretical studies [8-17] are devoted to understanding the impact of noble gas nanoinclusions on the fuel properties, in particular their thermomechanical properties. The experimental studies revealed the existence of small bubbles of 1-3 nm in diameter in fuels with a low burnup (when less than 2% of the uranium atoms underwent fission) [2]. These small bubbles can then coalesce at higher burnup (when more than 4% of the uranium atoms underwent fission) to lead to a bimodal bubbles population, with the largest bubbles of tens of nanometers in diameter [2]. Interestingly, xenon density in these bubbles (3.8 to 6 g/cm3) is characteristic to solid xenon [5]. Usually, the fission gas behaviour is modelled either by the van der Waals equation of state 2
or the ideal gas law, which have been found to be inadequate to describe high density xenon nanoinclusions. As a consequence new equations of state for Xe under high pressure and high temperature have been recently developed [18]. Some authors [10,13] evaluated the pressure in the bubble by comparing to the pressure obtained in a simulation of an equivalent density and temperature of the bulk noble gas, thus neglecting the gas/matrix interface effect (confinement effect). As underlined by some authors [19], the confinement effect can play a crucial role. The present work is an atomistic simulation study on the gas/matrix interaction when xenon nanoinclusions/bubbles are trapped in UO2 at temperatures characteristic to this nuclear fuel during the normal functioning of the reactor (1000 K) or during accidents (2000 K). This is achieved using classical molecular dynamics simulations. The paper is structured as follows: Section 2 gives some accounts of the simulation approach. In the first part of Section 3, the effect of the matrix on the confined xenon is studied. The amount of fuel swelling as a result of the intragranular gas bubbles is calculated and discussed in the second part of Section 3. All results are summarized in Section 4.
2. Models and methods 2.1. Models for the bubble/matrix systems All our calculations were carried out using 3D periodic boundary conditions. The pattern of the bubble/UO2 model is built in a very similar manner as that described in reference [20], starting with a 6×6×6 supercell of perfect UO2. An almost spherical cavity with the radius of about 0.45 nm is created in the middle of this supercell by removing a stoechiometric number of uranium and oxygen atoms. The resulting system has 854 uranium atoms and 1708 oxygen 3
atoms. A number of 15, 26, 35, 43, 50, 60 and 70 xenon atoms were introduced in the cavity using a dense body centred cubic (bcc) Xe arrangement, then the system was relaxed through molecular dynamics. The Xe/U ratios corresponding to these cavity fillings are: 0.018, 0.030, 0.041, 0.050, 0.059, 0.070 and 0.082 respectively. Usually, the gas bubbles distort the zone in the matrix that is situated close to the bubble/matrix interface. The supercell size was defined in order to avoid the interaction between the distorted zone surrounding the bubble and its corresponding images. Thus, there is an unperturbed matrix zone of at least 0.7 nm between the bubble and its images.
2.2. Equilibrium structures The properties of the systems in the present work are calculated at equilibrium. The equilibrium configurations are obtained through Newton-Raphson energy optimization and molecular dynamics methods [21]. Energy optimization allows finding the local minima closest to the starting configuration on the potential energy surface (PES). Molecular dynamics helps exploring larger parts of the PES by propagating Newton’s equations of motion through time using finite difference formalism. At 0 K, the quantity that is minimized during the energy optimization is the potential energy. To take into account the temperature impact on the equilibrium structure, a natural approach is to consider free energy minimization. The free energy of a structure at a given temperature is calculated through the determination of the vibrational partition function from the phonon density of states, within the quasi-harmonic approximation. This is usually assumed to be correct for temperatures lower than half of the melting point temperature of the considered material [22]. 4
For all our molecular dynamics simulations the timestep is 1 fs. The temperature and the pressure are kept at constant values using the Nose-Hoover thermostat and barostat algorithms [21]. The calculations are performed using the GULP 3.1 simulation package [22], which implements all the aforementioned methods.
2.3. Interatomic potentials Within our bubble/matrix systems there are three types of interactions: the in-matrix U/U, U/O and O/O interactions, the xenon/matrix Xe/U and Xe/O interactions and the Xe/Xe interactions. In order to describe the in-matrix interactions the potential proposed by Basak [23] was selected. Indeed, a systematic study carried out by Govers over 20 different potentials identified Basak’s potential [8] as one of the best, more particularly at high temperature, and Geng et al. [12] showed that this potential is transferable to structures far from the UO2 fluorite structure used to define its parameters. Moreover, a recent study [20] shows this potential to be able to correctly describe the thermomechanical properties of nanoporous UO2. The analytical expression of the Basak potential is given by the equation 1 : ⎛
rij ⎞ Cij zi z j e 2 ⎟− 6 + + Dij ⎧⎨ 1 − exp β ij rij* − rij ⎟ r ρ r ⎩ ij ⎠ ij ij ⎝
φij = Aij exp⎜⎜ −
[
( (
))]
2
− 1⎫⎬ ⎭
(1)
This potential is a sum of three terms: a Buckingham-Hill term describing the van der Waals interactions, a Coulomb term describing the ionic interactions and a Morse term for the
5
covalent bonds, as in the UO2 the U-O bonds have an ionic–covalent character. The parameters of the Basak potential are given in table 1. For the xenon/matrix interaction, we used the potential proposed by Grimes [9]. This potential has a Buckingham-Hill form. The CXeU and the CXeO parameters (the dispersion part parameters) were calculated by Grimes using the Kirkwood-Slater formula as described in reference [24], and the repulsive part parameters A and ρ were calculated using an electrongas method [25]. The Grimes potential was initially proposed in a shell-core form. In the present work, it is used in a rigid ion form since this one is far more effective for molecular dynamics simulations [8]. In the shell-core form smaller timesteps are needed and the number of interactions to compute increases by a factor of 4. Since the Buckingham-Hill potential has a non-physical zone of attraction at small distances and the present calculations are carried out for very dense systems and high temperatures, there is a risk to explore this non-physical zone. Therefore, the Buckingham-Hill form was replaced at short distance with a Born-Mayer repulsive potential. Thus, the potential has the following form: ⎛
rij ⎞ ⎟ BM ⎟ ρ ij ⎝ ⎠
φij = AijBM exp⎜⎜ −
if
rij < rc
(2)
and ⎛
rij ⎞ Cij ⎟− 6 ⎟ r ρ ij ⎠ ij ⎝
φij = Aij exp⎜⎜ −
if rij ≥ rc
Here rc is the inflexion point of the Buckingham-Hill function situated between r = 0 and the equilibrium distance r = rmin. The ABM and ρBM parameters of the Born-Mayer function were calculated from the conditions of continuity and derivability of the potential in rc. The Grimes potential choice for our study was made after it was assessed together with 6
three other largely used semiempirical potentials: the Geng potential [12], the Jackson-2 potential [8,26] and the Chartier potential [14]. We used some quantum DFT calculation results from literature [14-17] for comparison. The Xe incorporation energy in different sites of the UO2 matrix is one of the most appropriate assessment parameters. The considered sites can be easily identified on figure 1 representing the UO2 elementary cell : the octahedral interstitial site situated in the middle of the cube formed by the oxygen atoms, the uranium and the oxygen vacancies obtained by displacing an uranium or an oxygen far from their perfect crystal positions to form the Frenkel pairs, the UO divacancy obtained by displacing both an uranium and the nearest neighbour oxygen and the UO2 trivacancy (Schottky defect) obtained by removing an uranium atom and two nearest neighbour oxygen atoms. There are three trivacancy configurations : 145 (I), 124 (II), 123 (III). The incorporation energy was calculated at 0 K as the difference between the internal energy of the relaxed systems containing xenon and the internal energy of the xenon free systems. The patterns of the periodic models representing the xenon free systems and the systems containing xenon were built based on a 4x4x4 supercell of UO2 fluorite structure. The in-matrix interactions were calculated using the Basak potential. Table 2 gives the semiempirical potential as well as the DFT values for the incorporation energies. All the semiempirical potentials give the good order of the incorporation energies (as compared to DFT) with the interstitial site having the highest incorporation energy and the configurations II and III of the trivacancy having the lowest ones. The Grimes potential seems to give an overall better description than the other potential models. For the xenon/xenon interactions three semiempirical potentials were assessed: a Buckingham-Hill potential used in reference [27], a Lennard-Jones potential also from reference [27] and a potential proposed by Tang and Toennies [28]. In the Buckingham-Hill
7
potential the short distance part was replaced with a Born-Mayer repulsive potential in the same manner and for the same reason as for the Grimes potential. As a result, the potential has also the form given by equation 2. Some experimental results [29,30] were used for comparison. The first part of the assessment consisted in comparing the predictions of the density versus pressure ρ(P) curve at room temperature given by the three potentials with the experimental results (figure 2). The ρ(P) curve calculated with Buckingham-Hill is the only one to follow the experimental points for all the considered pressure range. The other two potentials produce curves that move away from the experimental points [29] at larger pressures. In the second part of the assessment, the bulk modulus evolution versus pressure K(P) is calculated at room temperature and compared with the experiment (figure 3). In the lower part of the pressure range (from 3 GPa to 6 GPa) the Lennard-Jones potential gives the best results, but for larger pressures it fails to account for the experimental points [30]. The Buckingham-Hill values seem again to best approach the experimental data in the considered pressure range. The ρ(P) curve was calculated through molecular dynamics simulations in NPT ensemble on a supercell of 500 xenon atoms. Within this supercell a solid/liquid interface was created in order to let the system the freedom to evolve toward the most stable phase for a given (T,P) set. 30 ps were used for the equilibration and another 30 ps were used to accumulate the results. At room temperature (298 K) and the considered pressure range from 5 to 50 GPa the system evolves toward the face centred cubic (fcc) structure. This result is in agreement with the phase diagram from reference [31]. The K(P) curve was calculated through free energy optimization on fcc xenon at 296 K and the pressure from 3 to 10 GPa. The free energy 8
optimization technique can be safely used since the considered temperature is less than half the melting point temperature in this pressure range [31]. All the parameters of the xenon/matrix and xenon/xenon potentials used in this study are presented in table 3.
3. Results and discussion 3.1. The characterization of the xenon in the bubble Estimating the state of the xenon confined in the form of bubbles in UO2 as well as the pressure in the bubble are a challenging task from both experimental [4-7] and theoretical [12,13] viewpoints. High resolution transmission electron microscopy (HRTEM) and energy dispersive X-ray spectroscopy (EDX) analyses, carried out in the outer region of high burnup UO2 pellets [4,5], showed the existence of almost spherical xenon inclusions/bubbles with diameters ranging from 1 nm at low burnups to 10 nm at high burnups. The xenon in the bubbles seems to be in a solid or in a near solid state with densities from 3.8 to 6.0 g/cm3 (for bubbles with diameters ranging from 4 to 10 nm). For these bubble radii and densities, at 700 K, the pressure in the bubble was estimated at 1.6 to 15 GPa. X-ray absorption spectroscopy (XAS) experiments [6,7] carried out on a set of xenon implanted UO2 samples indicated the existence of highly pressurized xenon inclusions/bubbles. The pressure in the bubble was estimated at 2.1 to 2.7 GPa at 5 K and at 3.0 to 3.7 GPa at 870 K for bubbles of approximately 2 nm in diameter. In these experimental studies the pressure in the bubble was estimated by equating to the pressure of the bulk noble gas at an equivalent density and temperature. The same method of estimating the pressure in the bubble was employed in theoretical works [10,13]. However this approach might not be very appropriate given the existence of the 9
xenon/matrix interactions. Moreover, the bubble size, on the order of a few nanometers, helps to increase these interactions. In the present work the pressure in the bubble is calculated “in situ”, taking into account not only the xenon/xenon interactions but also the xenon/matrix ones. This is achieved using the virial equation for pair-additive forcefields as recommended in reference [32] in the case of periodic boundary conditions. A local form of this equation can be obtained using as a starting point the local stress tensor. The σαβ component of the local stress tensor σ corresponding to a sub volume V inside of a material (an inclusion in a matrix for example) is calculated using the equation 3, where the summation runs over all the atoms i contained in the sub volume V [33] :
σ αβ =
1 1 mi viα viβ + ∑ V i∈V 2V
⎛
1 ∂U ij ∂rij
∑∑ ⎜⎜ − r i∈V j ≠ i
⎝
⎞α β ⎟rij rij ⎟ ⎠
(3)
Here mi is the mass of the atom i, α and β are the Cartesian xyz components, viα is the α component of the atom i velocity vector, rijα is the α component of the vector joining the atoms i and j and rij represents its norm. U is the interatomic pair potential. The summation in the second term on the right side runs over all the atoms j inside the range of the interatomic potential U. Thus, the atoms j may be contained in the sub volume V or may belong to the rest of the system surrounding this sub volume. The pressure in the bubble Pbubble is calculated as the average over time of one third of the trace of the local stress tensor σ[34]. Using this definition for the pressure and the equipartition theorem to relate the instantaneous system temperature T with the average over all the degrees of freedom of the kinetic energy, one finally obtains the virial equation that gives the pressure
10
in the bubble:
Pbubble =
1 N Xe kT + V 3V
⎛ ⎞ 1 ⎜ ∑ rij f ij + ∑∑ riL f iL ⎟ ⎜ ⎟ 2 i L ⎝ i< j ⎠
(4)
The first term on the right-hand side of the equation 4 is the kinetic (ideal) term whereas the second one is the potential (non-ideal) term. The non-ideal term was split into two parts: a xenon/xenon part and a xenon/matrix one. Hence rij and fij are the distance and the force acting between the xenon atoms i and j whereas riL and fiL are the distances and forces acting between the xenon atom i and the matrix atom L (oxygen or uranium). As the model is periodic, the sums run on all over the atoms in the cell as well as the image atoms from the neighbouring cells. T is the instantaneous temperature, k is the Boltzmann’s constant, NXe is the number of the xenon atoms in the bubble and V represents the instantaneous volume of the bubble. There are multiple ways for estimating the instantaneous volume of the bubble. The reference [35] presents a method, based on the tessellation of space into Voronoi and Delaunay polyhedra, that allows the exact calculation of volumes of cavities in sphere packings in three dimensions. A less computationally expensive way to estimate the bubble volume assumes that the bubble is spherical its radius being calculated [36] as the distance from the center of the bubble to the nearest uranium or oxygen atoms. In the present work we used the spherical shape hypothesis since high resolution transmission electron microscopy (HRTEM) images [4,5] show the intragranular bubbles to have almost spherical shapes. However, at the atomic scale there are irregularities in the bubble shape. For example some matrix atoms (oxygen and/or uranium) can be found closer to the bubble mass center than some xenon atoms. In order to take into account, to some extent, this fact we calculated the bubble radius as the arithmetic mean between a “cavity minimum radius” (Rca) and a “cluster
11
maximum radius” (Rcl). The “cavity minimum radius” is defined as the distance between the bubble mass center and the closest matrix atom (oxygen or uranium) and the “cluster maximum radius” as the distance between the bubble mass center and the furthest (xenon) atom in the xenon cluster. The figure 4 shows the two possible situations that can occur when calculating the bubble radius with this method: some matrix atoms (oxygen and/or uranium) are closer to the bubble mass center than some xenon atoms (left side of figure 4) and all the xenon atoms are closer to the bubble mass center than the matrix atoms (right side of figure 4). Given the relative sizes of the polyhedron closing the xenon cluster and the Rcl and Rca radius spheres and the fact that the atoms have finite volumes, in the former case our method might slightly underestimate the bubble instantaneous volume whereas in the latter it might slightly overestimate it. An instantaneous pressure was calculated at each configuration using the instantaneous temperature and the instantaneous bubble volume, then the pressure in the bubble was calculated as the average of these instantaneous pressures. We verified if a sufficient number of data points was used for the pressure average using the so called “block” average method [37,38]. Within this method one starts by calculating the standard error of the mean of all the data points, then two consecutive data points are grouped and the standard error of the mean of this new block of data is calculated. The new block of data contains half the number of data points of the initial set. The operation is repeated until there are not enough data points to compute a standard error. The original data set has highly correlated data points since the time between two data points (the timestep) is too short. As a consequence the standard error has a low value. By increasing the block size the data points are less and less correlated and the standard error of the mean increases with the number of block operations until so many data points are grouped that two consecutive points are uncorrelated. At this moment the standard
12
error reaches a plateau. Thus a plateau (figure 5) on the standard error of the mean versus number of block operations curve was used as an indicator that we have obtained enough statistics. Molecular dynamics simulations in NPT ensemble were carried out on systems with different numbers of xenon atoms, built as described in section 2.1, for temperatures of 1000 K and 2000 K. The pressure was set at 0 GPa because this facilitates the comparison with other theoretical results from the literature that were obtained using this pressure. Each system at 2000 K was relaxed during 80 ps. The systems were then “cooled” at 1000 K through two intermediate steps: one at 1700 K and the other at 1300 K. At each of these two temperatures the system was relaxed during 60 ps. At 1000 K the systems were relaxed again during 80 ps. This method was used because the systems relax faster at 2000 K than at 1000 K, thus a larger number of relaxation steps might be needed to obtain equilibrium structures at 1000 K. The pressure in the bubble as a function of the number of the xenon atoms in the bubble at 1000 K and 2000 K was then calculated by averaging over 40000 data points. Globally the block analysis showed this number of data points to be sufficient, but at 2000 K, for some cases, no plateau was found in the standard error versus number of block operations curve (see figure 5). However, the pressures calculated by averaging over a larger set of data points (140000 data points), for bubbles with NXe = 15, 43 and 70 (covering all the analyzed NXe range), were similar to those obtained using 40000 data points (see figure 6 bottom). We thus concluded that the pressure versus NXe curves were calculated with a sufficient accuracy. The pressure in the bubble was also calculated by equating to the pressure of the bulk xenon at an equivalent density and temperature. The results of the two methods are compared in figure 6. Both methods predict that the pressure in the bubble rapidly diminishes with increasing NXe until it reaches a minimum, then slowly increases. The pressures calculated “in 13
situ” are different (higher) from those corresponding to the bulk xenon. In reference [36] a similar result, but in terms of free energy, was obtained: the free energy of the xenon in a cavity is higher than the free energy of the bulk xenon at the same density and temperature for Xe densities of 3.8-6 g/cm3 and temperatures of 300-1050 K. The pressure in the bubble estimated with the ideal gas law is also represented in figure 6. It is much lower (< 1 GPa) than the pressures calculated with the two other methods that take into account the non-ideal terms. Moreover, it becomes clear from figure 6 that the difference between the “in situ” and the bulk xenon estimates of the pressure can be as important as the difference between the real gas and the ideal gas equation of state estimates. The difference between the ‘”in situ” and the bulk xenon pressure is a xenon/matrix interface effect. For large bubbles, where the surface interactions are far less important than the volume ones, one can expect the results of the two estimates to converge. As shown in figure 7 (top) the bubble radius is also a function of the number of xenon atoms in the bubble: the bubble radius increases with NXe but this increase slows down at high NXe. This is because, during the molecular dynamics relaxation, the surface that closes the
bubble evolves until the equilibrium is reached between the xenon pressure in the bubble on one side and the surface tension (attributed to the xenon/matrix interface) plus the external pressure on the other side. At the atomic scale, the xenon atoms push uranium and oxygen atoms into the matrix creating more space for themselves. This phenomenon was already observed by other authors. In reference [13] xenon atoms were introduced in UO2, one by one, starting from the perfect UO2 lattice. After adding each new xenon atom the resulting system was relaxed through molecular dynamics. At the end of each relaxation the noble gas atoms form bubbles of different sizes (function of the number of xenon atoms in the bubble) by displacing the oxygen and uranium atoms at the xenon/matrix interface. In reference [36]
14
the xenon atoms were introduced into an already built cavity into the UO2 (like in the present work) followed by a molecular dynamics relaxation of the system. The authors observed that, for sufficiently large numbers of xenon atoms in the bubble, the xenon atoms at the surface of the bubble created more space for themselves by forcing oxygen ions into the UO2 lattice. An equation relating the intragranular bubble radius (Rbubble) to the number of xenon atoms in the bubble (NXe) was proposed in reference [39]:
Rbubble = BN
1
1 3 Xe
where
⎛ 3 ⎞3 B = ⎜⎜ ⎟⎟ ⎝ 4πρ ⎠
(5)
Here ρ is the density of xenon atoms in the bubble in nm-3. As can be observed in figure 7 (bottom), the xenon density in the bubble, at both 1000 K and 2000 K, is a function of NXe, following the same evolution as the pressure: it decreases towards a minimum then increases at higher NXe. Considering that the ρ(NXe) curves can roughly be described by a third order polynomial and taking into account the equation 5, one obtains for Rbubble(NXe) the following equation: 1
Rbubble =
0.62035 N Xe 3
(aN
3 Xe
2 + bN Xe + cN Xe + d
)
1 3
(6)
where a = -0.00026, b = 0.0449, c = -2.3866 and d = 58.96. In figure 7 (top) we show that the equation 6 gives a good approximation of the molecular dynamics points. For almost all the cases we analyzed the density in the bubble is higher than 4 g/cm3, the
15
density of the solid face centred cubic (fcc) xenon [5]. This leads to the question of the state of the xenon confined in nanocavities at high temperature and pressure. The radial distribution function (RDF) shape analysis can give an insight into the state of the confined xenon. This method was already successfully used by other authors to characterize the pure xenon phase diagram [31] or even the xenon state into a bubble [13]. The RDF’s for the xenon in the bubble were calculated for systems with NXe = 15 up to 70 covering all the investigated range of conditions (figure 8). RDF’s (figure 9) were also built for bulk xenon at a density of about 6.1 g/cm3 (within the range of those found in the bubbles) and temperatures of 1000 K and 2000 K. Through molecular dynamics relaxation in NVT ensemble, the bulk xenon was found to be in a solid fcc phase at 1000 K and in a liquid phase at 2000 K, in agreement with the phase diagram from reference [31]. In order to facilitate the comparison with the RDF’s for the xenon in the bubble the RDF’s for bulk xenon were calculated only for a cluster of atoms found in a 0.6 nm radius sphere. The value of the sphere radius was taken to be representative for the bubble sizes analyzed in this study. The RDFs of the xenon in the bubble (figure 8) are very similar with the ‘liquid’ cluster RDF shown in the figure 9. Taking into account this observation and the fact that the xenon density in the bubble corresponds to the solid xenon, one can safely assume that the xenon in the nanocavities (bubbles), in the conditions analyzed here, is in a glassy phase, in spite of the fact that in some cases the (T,Pbubble) sets correspond to the solid fcc state on the bulk xenon phase diagram [31]. This is due to the fact that the xenon/matrix interaction alters the phase diagram of the xenon in the bubble through a confinement effect. In reference [36] a similar conclusion is reported : within the bubbles, the surface of the uranium dioxide restricts the rearrangement of the Xe atoms into the perfect fcc structure thus retarding the phase change to a solid fcc phase. 16
3.2. The calculation of the bubble induced swelling In this section we analyze the effect of the xenon/matrix interaction on the global behaviour of the nuclear fuel. A very important issue to be addressed is the nuclear fuel swelling during the normal functioning of the nuclear reactor or accidental situation. The swelling is the sum of three components: A swelling due to the UO2 matrix thermal expansion; A swelling due to the UO2 matrix fragmentation; A swelling due to the fission product accumulation inside the matrix, particularly the noble gases, in the form of bubbles. The last component, the bubble induced swelling, is supposed to become very important at high temperature. So it is crucial to quantify this swelling term and its thermal dependency. The bubble induced swelling was calculated using the equation 7 as given in the reference [13]:
S% =
V eq − V0eq × 100 V0eq
(7)
Here Veq is the equilibrium volume of the bubble/matrix supercell and V0eq is the equilibrium volume of the perfect UO2 supercell with the same number of UO2 groups as the bubble/matrix supercell. The equilibrium volumes were calculated through molecular
17
dynamics simulations in NPT ensemble for temperatures from 800 K to 2000 K and for a pressure of 0 GPa (no external stress). The swelling versus Xe/U ratio at 1600 K is plotted in figure 10 and compared with the results from reference [13]. The agreement between the results of the two studies is remarkable since in reference [13] a different set of interatomic potentials was used. At constant temperature the bubble induced swelling increases linearly with the Xe/U ratio. In order to verify if the system size has an effect on the calculated swelling we computed the swelling corresponding to a larger system. This system was built starting with a 7×7×7 supercell of perfect UO2. An almost spherical cavity with the radius of about 0.58 nm was created in the middle of this supercell by removing 16 UO2 units. The resulting system has 1356 uranium atoms and 2712 oxygen atoms. The cavity was filled with 35 and 70 xenon atoms as described in section 2.1. The corresponding Xe/U ratios are 0.0258 and 0.0516 respectively. The system was relaxed for 80 ps through molecular dynamics in NPT ensemble at T = 1600 K and P = 0 GPa. The swelling values calculated for these larger systems integrate well in the S%(Xe/U) curve from figure 10 showing that the system size has practically no effect on these results. The swelling versus temperature corresponding to three Xe/U ratios at P = 0 GPa is plotted in figure 11. It seems that the temperature has no effect on the bubble induced swelling. However, this conclusion should be restricted only to the intragranular bubbles since only these ones are analyzed here. Actually, the temperature dramatically increases the gaseous swelling but mainly due to the intergranular bubbles [40]. The intergranular bubbles grow faster and larger at high temperatures because the point defects and the fission gas atoms from the fuel grains can diffuse faster to the grain boundary and be trapped by the intergranular
18
bubbles.
4. Conclusions We have presented an atomistic simulation study on the behaviour of the intragranular xenon bubbles trapped in the UO2 matrix and their influence on the fuel swelling. The main conclusions that can be derived from this study are: 1) The xenon pressure in the bubble calculated “in situ” by taking into account the xenon/matrix interaction is systematically larger than the pressure calculated in simulations at an equivalent density and temperature of the bulk xenon. However, this discrepancy should vanish for large bubbles as a consequence of the surface (confinement) effect decrease with decreasing surface/volume ratio. Thus, the equations of state established for free xenon, whereas still appropriate for large bubbles, cannot be safely used to describe the xenon in bubbles of a couple of nanometres in diameter. 2) The bubble radius and the pressure in the bubble are functions of NXe, the number of xenon atoms confined in the bubble. The radius evolution is fast at low NXe but slows down at higher NXe. 3) The radial distribution function analysis is indicative of a glassy phase for the xenon in the bubble even for temperature/pressure points corresponding to solid fcc xenon on the bulk xenon phase diagram. This is also an effect of the xenon/matrix interface. However, a more detailed analysis would require the calculation of the diffusion coefficient of the confined xenon and a comparison with the corresponding bulk xenon ones. This will form the object of a future work dedicated to the state of the noble gases confined in UO2.
19
4) The bubble induced swelling of the fuel increases linearly with the quantity of xenon in the system (Xe/U ratio). However, for a given Xe/U ratio, the temperature doesn’t seem to have an important impact on the intragranular swelling in the 800K to 2000K temperature range.
Acknowledgements. The calculations were performed at the facilities of the CNRS/CINaM laboratory and IDRIS, the CNRS computing center. We would like to thank G. Tréglia from CNRS/CINaM for his useful suggestions.
20
References [1] J. Rest, J. Nucl. Mater. 402 (2010) 179-185. [2] P. Lösönen, J. Nucl. Mater. 280 (2000) 56-72. [3] K. Govers, C.L. Bishop, D.C. Parfitt, S.E. Lemehov, M. Verwerft, R.W. Grimes, J. Nucl. Mater. 420 (2012) 282-290. [4] K. Nogita, K. Une, J. Nucl. Mater. 250 (1997) 244-249. [5] K. Nogita, K. Une, Nucl. Instr. and Meth. in Phys. Res. B 141 (1998) 481-486. [6] P. Garcia, P. Martin, G. Carlot, E. Castelier, M. Ripert, C. Sabathier, C. Valot, F. D’Acapito, J.L. Hazeman, O. Proux, V. Nassif, J. Nucl. Mater. 352 (2006) 136-143. [7] P. Martin, P. Garcia, G. Carlot, C. Sabathier, C. Valot, V. Nassif, O. Proux, J.L. Hazemann, Nucl. Instr. and Meth. in Phys. Res. B 266 (2008) 2887-2891. [8] K. Govers, PhD Thesis, Université libre de Bruxelles, 2008. [9] R.W. Grimes, C.R.A. Catlow, Phil. Trans. R. Soc. Lond. A 335 (1991) 609-634. [10] D.C. Parfitt, R.W. Grimes, J. Nucl. Mater. 381 (2008) 216-222. [11] D.C. Parfitt, R.W. Grimes, J. Nucl. Mater. 392 (2009) 28-34. [12] H.Y. Geng, Y. Chen, Y. Kaneta, M. Kinoshita, J. Alloys Compd. 457 (2008) 465-471. [13] T.X. Feng, L.C. Sheng, Z.Z. He, G. Tao, Chin. Phys. B 19 (2010) 057102. [14] A. Chartier, L. Van Brutzel, M. Freyss, Phys. Rev. B 81 (2010) 174111. [15] G. Brillant, F. Gupta, A. Pasturel, J. Nucl. Mater. 412 (2011) 170-176.
21
[16] A.E. Thompson, C. Wolverton, Phys. Rev. B 84 (2011) 134111. [17] H.Y. Geng, Y. Chen, Y. Kaneta, M. Kinoshita, Q. Wu, Phys. Rev. B 82 (2010) 094106. [18] J.Y. Oh, Y.H. Koo, J.S. Cheon, B.H. Lee, D.S. Sohn, J. Nucl. Mater. 372 (2008) 89-93. [19] C. Alba-Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K.E. Gubbins, R. Radhakrishnan, M. Sliwinska-Bartkowiak, J. Phys.: Condens. Matter 18 (2006) R15-68. [20] A. Jelea, M. Colbert, F. Ribeiro, G. Treglia, R.J.-M. Pellenq, J. Nucl. Mater. 415 (2011) 210-216. [21] A.R. Leach, Molecular Modelling Principles and Applications, Pearson Education Limited, Harlow, 2001. [22] J.D. Gale, A.L. Rohl, Mol. Simul. 29 (2003) 291-341. [23] C.B. Basak, A.K. Sengupta, H.S. Kamath, J. Alloys Compd. 360 (2003) 210-216. [24] P.W. Fowler, P.J. Knowles, N.C. Pyper, Molec. Phys. 56 (1985) 83-95. [25] J.H. Harding, A.H. Harker, UKAEA Harwell Report, AERE-R-10425 (1982). [26] S. Nicoll, H. Matzke, C.R.A. Catlow, J. Nucl. Mater. 226 (1995) 51-57. [27] I.R. Brearley, D.A. MacInnes, J. Nucl. Mater. 95 (1980) 239-252. [28] K.T. Tang, J.P. Toennies, J. Chem. Phys. 118 (2003) 4976-4983. [29] A.N. Zisman, I.V. Aleksandrov, S.M. Stishov, Phys. Rev. B 32 (1985) 484-487. [30] S. Sasaki, N. Wada, T. Kume, H. Shimizu, J. Raman Spectrosc. 40 (2009) 121-127. [31] A.B. Belonoshko, Phys. Rev. B 78 (2008) 174109.
22
[32] M.J. Louwerse, E.J. Baerends, Chem. Phys. Lett. 421 (2006) 138-141. [33] P.S. Branicio, R.K. Kalia, A. Nakano, P. Vashishta, F. Shimojo, J.P. Rino, J. Mech. Phys. Solids 56 (2008) 1955-1988. [34] D.H. Tsai, J. Chem. Phys. 70 (1979) 1375-1382. [35] S. Sastry, D.S. Corti, P.G. Debenedetti, F.H. Stillinger, Phys. Rev. E 56 (1997) 55245532. [36] S.T. Murphy, A. Chartier, L. Van Brutzel, J.P. Crocombette, Phys. Rev. B 85 (2012) 144102. [37] H. Flyvbjerg, H.G. Petersen, J. Chem. Phys. 91 (1989) 461-466. [38] D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, 2001. [39] D.R. Olander, D. Wongsawaeng, J. Nucl. Mater. 354 (2006) 94-109. [40] J. Spino, J. Rest, W. Goll, C.T. Walker, J. Nucl. Mater. 346 (2005) 131-144.
23
Figure 1
24
Figure 2
25
Figure 3
26
Figure 4
27
Figure 5
28
Figure 6
29
Figure 7
30
Figure 8
31
Figure 9
32
Figure 10
33
Figure 11
34
Figure captions Fig. 1. The UO2 elementary cell. It contains 4 uranium atoms (large spheres) forming a face centred cubic like arrangement and 8 oxygen atoms (small spheres) forming a cube in the middle of the cell. The numbers help to define the trivacancy configurations in section 2.3. Fig. 2. Comparison between the density versus pressure ρ(P) curves calculated at 298 K with three xenon/xenon semiempirical potentials and experimental data. The semiempirical potentials are: a Buckingham-Hill (Br-BH) potential, a Lennard-Jones (Br-LJ) potential and a potential proposed by Tang and Toennies (TT). Fig. 3. Comparison between the bulk modulus versus pressure K(P) curves calculated at 296 K with three xenon/xenon semiempirical potentials and experimental data. The semiempirical potentials are: a Buckingham-Hill (Br-BH) potential, a Lennard-Jones (Br-LJ) potential and a potential proposed by Tang and Toennies (TT). Fig. 4. Calculation of the bubble radius as the arithmetic mean between Rca and Rcl. Rca is the distance between the bubble mass center and the closest matrix atom (oxygen or uranium) and Rcl is the distance between the bubble mass center and the furthest (xenon) atom in the bubble.
Two possible situations can occur when calculating the bubble radius with this method: some matrix atoms (oxygen and/or uranium) are closer to the bubble center than some xenon atoms (left) and all the xenon atoms are closer to the matrix center than the matrix atoms (right). Fig. 5. The standard error of the mean σ for the pressure in the bubble as a function of the number of block operations NB for 40000 time steps of simulation. A plateau can be observed at NB = 7 or 8 on the curves corresponding to systems at 1000 K. An example without a plateau is given for a system at 2000 K. Fig. 6. The pressure in the bubble as a function of NXe: top – systems at 1000 K, bottom – 35
systems at 2000 K. The error bars represent the sample standard deviation. The pressures calculated “in situ” taking into account the xenon/matrix interactions are compared with those calculated for bulk xenon at an equivalent density and temperature. For systems with NXe = 15, 43 and 70 atoms, at 2000 K, the pressure in the bubble was calculated by averaging over 40000 and 140000 (better statistics) data points. The pressures at an equivalent density and temperature given by the ideal gas law are also displayed. Fig. 7. Top – the bubble radius as a function of NXe. Bottom – the density of the xenon in the bubble as a function of NXe. The error bars represent the sample standard deviation. Fig. 8. Radial distribution functions (RDF) corresponding to the confined xenon in bubbles containing 15, 43 and 70 xenon atoms: top – at 1000 K; bottom – at 2000 K. Fig. 9. Radial distribution functions (RDF) corresponding to the bulk xenon in liquid and solid states. The RDF’s for bulk xenon were calculated for a cluster of atoms found in a 0.6 nm radius sphere. Fig. 10. The swelling as a function of the Xe/U ratio at 1600 K. The results of the present work are compared to those from the reference [13]. The “smaller system” is the system derived from a 6x6x6 supercell and the “larger system” is derived from a 7x7x7 supercell. Fig. 11. The swelling as a function of temperature for different Xe/U ratios.
36
Table 1 Parameters of the Basak potential used to describe the in-matrix interactions z
A(eV)
ρ(nm)
C(eV nm6)
D(eV)
r*(nm)
β(nm-1)
U
2.4
-
-
-
-
-
-
O
-1.2
-
-
-
-
-
-
OO
-
1633.390
0.0327
3.95 x 10-6
-
-
-
UU
-
294.745
0.0327
0
-
-
-
UO
-
693.870
0.0327
0
0.57745
0.2369
16.5
37
Table 2 The incorporation energies (in eV units) of xenon at different sites in the UO2 matrix as calculated using the semiempirical potentials of Grimes, Geng, Jackson (version 2) and Chartier for the xenon/matrix interactions and the Basak potential for the in-matrix interactions. Some DFT results from literature are given for comparison. Grimes
Geng
Jackson-2
Chartier
DFT[14-17]
Interstitial
17.1
22.8
20.3
12.8
9.7-12.0
U vacancy
5.8
5.8
7.0
7.5
2.0-5.8
O vacancy
13.4
19.6
17.0
10.2
7.5-9.1
UO divacancy
4.8
5.3
6.1
6.2
1.6
UO2 trivacancy (I)
4.9
5.0
5.9
6.5
5.4
UO2 trivacancy (II)
4.3
4.9
5.6
6.0
5.2
UO2 trivacancy (III)
3.7
4.8
5.2
5.1
0.2-5.2
38
Table 3 Parameters of the xenon/matrix and xenon/xenon potentials used within this study ABM(eV)
ρBM(nm)
A(eV)
ρ(nm)
C(eV nm6)
rc(nm)
XeU
2915.33
0.0410
6139.16
0.0340
7.184 x 10-5
0.11
XeO
136.70
0.0561
598.00
0.0426
1.0838 x 10-4
0.18
XeXe
1171.48
0.0466
4934.10
0.0356
2.9747 x 10-4
0.15
+
39