Maugis

L Verma

L Noirot

P Maugis

Modelling intra-granular bubble movement and fission gas release during post-irradiation annealing

Keywords: Fission gas, Nuclear fuel, Modeling, Intra-granular bubble

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Meso-scale

Introduction

Fission gases generated during irradiation in a nuclear fuel cause macroscopic phenomena like Fission Gas Release (FGR) and swelling in the fuel. These phenomena can affect the proper functioning of the fuel rod. FGR from the fuel increases the pressure in the fuel rod plenum, subjecting the cladding to additional stress, and also reduces the thermal conductivity of the fuel-cladding gap causing the fuel operating temperature to increase. Due to their low solubility in UO 2 , fission gases also precipitate into highly pressurized bubbles causing the swelling of the fuel. Swelling contributes to the fuel-cladding interaction, again exposing the cladding to higher stress and temperature conditions and ultimately hampering its lifetime. So, understanding fission gas behaviour is very important for optimal utilization of fuel rod during nuclear reactor operation.

In order to understand fission gas behaviour, in-pile as well as out-of-pile measurements are carried out in the nuclear fuel. It is, however, difficult to carry out measurements in a fuel under irradiation (in-pile) due to uncontrollable environment variables. Post-irradiation annealing (out-of-pile) tests are carried out to obtain data on FGR under controlled and monitored environment. One of the interesting issues during post-irradiation annealing tests has been the transport of intra-granular gas to the grain surface, as it is found to be significant while effective diffusion theory [START_REF] Speight | A calculation on the migration of fission gas in material exhibiting precipitation and re-solution of gas atoms under irradiation[END_REF] would predict no intra-granular gas release at all. Indeed, the effective diffusion coefficient is equal to D ef f = b b+g D Xe , where b and g are the re-solution and trapping probabilities per second, respectively, and D Xe is the intrinsic diffusion coefficient of the gas. In the absence of fission and if Xe is considered as insoluble in UO 2 , the re-solution probability, b, is nil. In other words, atomic gas is immediately trapped by the intra-granular bubbles, which are present in abundance [START_REF] Kashibe | Formation and growth of intragranular fission gas bubbles in UO 2 fuels with burnup of 6-83 GWd/t[END_REF]. Several mechanisms for the transport of intra-granular gas atoms outside the grain have been proposed. These include thermal resolution of gas atoms [START_REF] Macewan | Migration of xenon through a UO 2 matrix containing trapping sites[END_REF][START_REF] Macinnes | A model for the release of fission gas from reactor fuel undergoing transient heating[END_REF] or the inhibition of gas precipitation into bubbles [START_REF] Ronchi | On diffusion and precipitation of gas-in-solid[END_REF]. Other scenarios have to consider bubble movement. Evans [START_REF] Evans | Bubble diffusion to grain boundaries in UO 2 and metals during annealing: A new approach[END_REF] suggested that the rapid transport of gas atoms could be due to the movement of intragranular bubbles containing the gas in a vacancy concentration gradient.The mechanism of bubble movement in a vacancy gradient was further modelled by Veshchunov and Shestak [START_REF] Veshchunov | Modelling of fission gas release from irradiated uo2 fuel under high-temperature annealing conditions[END_REF] to analyze the fission gas release from irradiated UO 2 fuel under high-temperature annealing conditions. They also considered the role of dislocations in aiding the bubble movement.

In this paper, we focus on the mechanism proposed by Evans, and analyze the fission gas release due to the directed movement of gas bubbles in a vacancy concentration gradient. The experimental studies carried out by Kashibe et al. [START_REF] Kashibe | Formation and growth of intragranular fission gas bubbles in UO 2 fuels with burnup of 6-83 GWd/t[END_REF] showed a preferential coarsening of bubbles near the grain boundaries and showed that the bubble coarsening was associated with a decrease in the bubble number density. This coarsening was attributed to coalescence due to migration of bubbles. Thermal treatment of uranium oxide carried out by Zacharie et al. [START_REF] Zacharie | Thermal treatment of uranium oxide irradiated in pressurized water reactor: Swelling and release of fission gases[END_REF] also attributed the swelling and gas release to the coalescence of bubbles. Although there have been several studies focused on the behaviour of intra-granular bubbles and their post-irradiation treatment [START_REF] Lösönen | On the behaviour of intragranular fission gas in UO 2 fuel[END_REF][START_REF] Veshchunov | An advanced model for intragranular bubble diffusivity in irradiated UO 2 fuel[END_REF][START_REF] Li | Phase-field simulations of intragranular fission gas bubble evolution in UO 2 under post-irradiation thermal annealing[END_REF], we believe that a thorough analysis of FGR due to intragranular bubble movement incorporating bubble coalescence is missing. In a previous paper by the authors [START_REF] Verma | A new spatially resolved model for defects and fission gas bubbles interaction at the mesoscale[END_REF], a new model was presented for the interactions between fission gas bubbles and vacancies. The advantages and limitations of the model with respect to already existing methods like phasefield [START_REF] Li | Phase-field simulations of intragranular fission gas bubble evolution in UO 2 under post-irradiation thermal annealing[END_REF][START_REF] Tonks | An objectoriented finite element framework for multiphysics phase field simulations[END_REF] and cluster dynamics [START_REF] Barbu | Cluster dynamics modeling of materials: Advantages and limitations[END_REF] were also discussed. Simple analyses using the model showed that the model was capable of simulating the bubble movement and evolution in a vacancy concentration gradient. Here, we use the model to carry out a more extensive analysis of gas bubble movement and its impact on overall fission gas release.

The organization of the paper is as follows: In Section 2 are discussed the mechanisms of intra-granular bubble movement. The model used for the study, which was presented in detail in [START_REF] Verma | A new spatially resolved model for defects and fission gas bubbles interaction at the mesoscale[END_REF], is described in brief in Section 3. Section 4 presents the verification tests for the various aspects of the model such as the grid sensitivity, the diffusion calculation and the bubble movement. The results of the analyses for fission gas release are presented and discussed in Section 5 and the conclusion drawn from the analyses are presented in Section 6.

Mechanisms of intra-granular bubble movement

Directed movement in a vacancy concentration gradient

As stated before, Evans [START_REF] Evans | Bubble diffusion to grain boundaries in UO 2 and metals during annealing: A new approach[END_REF] proposed a mechanism for the transport of gas atoms by the movement of intra-granular bubbles containing the gas. According to him, during high temperature annealing, the influx of vacancies from the grain surface towards the grain would create a vacancy concentration gradient and this would impose a directed movement of intra-granular gas bubbles towards the boundary, eventually accelerating the gas release. He based this mechanism on the well-known mechanisms of movement of bubbles up vacancy gradients [START_REF] Nichols | Kinetics of diffusional motion of pores in solids: A review[END_REF] and efficiency of grain boundaries and surfaces as high temperature vacancy sources as recognized in the early works of Barnes et al. [START_REF] Barnes | The observation of vacancy sources in metals[END_REF][START_REF] Barnes | The generation of vacancies in metals[END_REF] and more recently in [START_REF] Kashibe | Formation and growth of intragranular fission gas bubbles in UO 2 fuels with burnup of 6-83 GWd/t[END_REF][START_REF] Chkuaseli | Modelling of fission gas spatial distribution in single grains of UO 2 fuel[END_REF]. The mechanism of directed bubble migration can be understood better with a pictorial representation shown in Fig. 1. From the figure, we see that the free surface -or grain boundary coated by large, already vented bubbles -generates vacancies which move towards the gas bubble and are trapped by it. This is due to the difference between the equilibrium concentration of vacancy at the vicinity of the free surface and the vicinity of the pressurized bubble [START_REF] Noirot | A method to calculate equilibrium concentrations of gas and defects in the vicinity of an over-pressured bubble in UO 2[END_REF]. The gas bubble may, in turn, emit vacancies towards other more pressurized bubbles. As a consequence, the bubble shifts towards the free surface/grain boundary. So, the bubble moves up the vacancy gradient and may reach the grain boundary and eventually cause the gas within the bubble to be released. Evans also carried out a quantitative assessment [START_REF] Evans | The role of directed bubble diffusion to grain boundaries in post-irradiation fission gas release from UO 2 : A quantitative assessment[END_REF] of the mechanism by using numerical calculations to simulate the phenomenon. He assumed a spherical grain with uniform distribution of over-pressurized gas bubbles. Annealing of the grain allowed an influx of vacancies from the grain surface inwards the grain, causing the bubbles to grow by migration and coalescence.

The model used by Evans postulates the existence of a front (defined by C v ≈ 0), where the bubbles undergo a local fractional swelling, ∆S, (a key parameter of his model) before the front progresses inwards the grain. In this model, the vacancies that arrive from the surface are supposed to be entirely trapped by the bubbles located on the front. This causes the vacancy gradient to last long enough to allow bubbles that are in the gradient to reach the grain boundary. However, the existence of this front has never been demonstrated. The FGR values are also found to be a function only of this fractional swelling in his model. Moreover, in this model, the bubbles which are further into the grain remain unaffected by the influx of vacancies while the bubbles at the front grow to equilibrium values. Though the bubbles closest to the grain boundary trap the vacancies first, it is not obvious that the bubbles a little further from the grain boundary do not grow at all before the closest bubbles reach equilibrium values. So, we require to develop more mechanistic models that would explicitly describe the bubbles and their migration in a vacancy gradient to access the impact of intra-granular bubble migration on the overall FGR.

Random (Brownian) movement of bubbles

A bubble can move in a solid without any driving force -in particular no temperature gradient -by the transfer of atoms around it . This transfer of atoms can be governed via volume or surface diffusion mechanisms. These mechanisms are random by nature.An advanced model for intra-granular bubbles diffusivity in irradiated UO 2 fuel via these mechanisms was proposed by Veshchunov and Shestak [START_REF] Veshchunov | An advanced model for intragranular bubble diffusivity in irradiated UO 2 fuel[END_REF] for the Van-der Waals gas in bubbles by using a non-linear adsorption law. We discuss the equations governing the diffusion of bubbles via random motion in the following sections.

Bubble diffusion by volume diffusion mechanism

The volume diffusion mechanism is the transfer of crystal atoms via vacancies in the bulk solid near the bubble and in the absence of any driving forces. The diffusion of bubbles is characterized by the diffusion coefficient, D vol b , which is defined according to the expression given by Olander [START_REF] Olander | Fundamental Aspects of Nuclear Reactor Fuel Elements[END_REF] as:

D vol b = 3 4π Ω r b 3 D v C eq v (bubble) ( 1
)
where Ω is the atomic volume, r b is the radius of the bubble, D v is the vacancy diffusion coefficient and C eq v (bubble) is the concentration of vacancy (at the vicinity of a bubble) in equilibrium with the bubble.

The relation that we use between D v and the uranium tracer coefficient, D U , is the following:

D U = D v * C eq v (P b = 0, κ = 0) * f (2)
where, f is the correlation factor that expresses that the U atoms do not exactly follow random-path theory because of the substitutional type of their diffusion. For this study, we have taken the value of f = 1. C eq v (P b = 0, κ = 0) is the concentration of vacancy (vacancy/site) at equilibrium in a solid limited by a flat surface (curvature, κ = 0) and no external pressure (P b = 0). This relation supposes that D U was measured when the concentration of vacancy was at equilibrium.

The volume self-diffusion coefficient,D U , can be calculated using Matzke's relation [START_REF] Matzke | Atomic transport properties in UO 2 and mixed oxides (U,Pu)O 2[END_REF]

D U = 6.5 * 10 -5 exp -5.6eV kT m 2 /s (3)
From Eq.11 in Section 3, we have the expression for C eq v (bubble) as:

C eq v (bubble) = C eq v (P b = 0, κ = 0) exp -Ω kT (P b -γ b κ) (4)
where, P b is the pressure in the bubble in Pa, γ b is the surface energy in J/m 2 and κ is the curvature of the bubble surface in m -1 . So, using Eq.2 and Eq.4 in Eq.1, we get the expression for D vol b as:

D vol b = 3Ω 4πr b 3 D U exp -Ω kT (P b -γ b κ) (5)
So, the random movement of bubbles via volume diffusion mechanism is characterized by the diffusion coefficient as given in Eq.5.

Bubble diffusion by surface diffusion mechanism

The surface diffusion mechanism is the transfer of crystal atoms on the surface of the bubbles. The coefficient of bubble diffusion via surface diffusion mechanism can be expressed from [START_REF] Olander | Fundamental Aspects of Nuclear Reactor Fuel Elements[END_REF] as:

D surf b = 3Ω 4/3 2πr b 4 D s (6)
where again, Ω is the atomic volume, r b is the radius of the bubble and D s is the surface diffusion coefficient.

Different values for the surface diffusion coefficient, D s , have been used by different authors in the past. The value for D s is usually obtained by either mass transport methods [START_REF] Reynolds | The surface self-diffusion of uranium dioxide[END_REF][START_REF] Henney | Surface-diffusion studies on UO 2 and MgO[END_REF][START_REF] Maiya | Surface diffusion, surface free energy, and grain-boundary free energy of uranium dioxide[END_REF] such as grain boundary grooving or scratch decay or by means of tracers [START_REF] Marlowe | Tracer study of the surface diffusivity of UO 2[END_REF][START_REF] Robertson | Surface diffusion of oxides (A review)[END_REF][START_REF] Olander | Interpretation of tracer surface diffusion experiments on UO 2 -roles of gas and solid transport processes[END_REF][START_REF] Zhou | Tracer surface diffusion on uranium dioxide[END_REF]. The available data on surface diffusion and matter transport driven by surface energy in ceramics has been reviewed by Matzke [START_REF] Matzke | Surface diffusion and surface energies of ceramics with application to the behavior of volatile fission products in ceramic nuclear fuels[END_REF]. The most accepted expression for D s recently [START_REF] Olander | Fundamental Aspects of Nuclear Reactor Fuel Elements[END_REF][START_REF] Bourgeois | Pore migration in UO 2 and grain growth kinetics[END_REF][START_REF] Veshchunov | An advanced model for intragranular bubble diffusivity in irradiated UO 2 fuel[END_REF][START_REF] Ahmed | Phase field simulation of grain growth in porous uranium dioxide[END_REF] is the one given by Maiya [START_REF] Maiya | Surface diffusion, surface free energy, and grain-boundary free energy of uranium dioxide[END_REF] or by Matzke [START_REF] Matzke | Surface diffusion and surface energies of ceramics with application to the behavior of volatile fission products in ceramic nuclear fuels[END_REF] in the temperature range of 1200-1800 0 C as given by:

D s = 50 exp -4.67eV kT m 2 /s (7)
Following results from Baker's data [START_REF] Baker | Materials for nuclear reactor core applications[END_REF], which showed an increase in the bubble diffusivity with increasing radius in the small radius regime, Mikhlin [START_REF] Mikhlin | The mobility of intragranular gas bubbles in uranium dioxide[END_REF] suggested a mechanism for the observed suppression of these small gas bubble diffusion mobility in solids. According to him, the interaction of ad-atoms with gas atoms, which are present in high density in small bubbles, causes a decrease of the ad-atom jump frequency on the bubble surface leading to the observed suppression of small gas bubble mobility. The expression for the bubble diffusion coefficient can be modified by incorporating a suppression factor, W b , which is equal to the probability that a region near an ad-atom is free of gas atoms and is expressed as:

W b = 1 - q V b N b (8)
where, q is the volume of the region with no gas atoms so that a jump can occur, N b is the number of gas atoms in the bubble of volume V b . We obtain the new values for the diffusion coefficient of bubbles by multiplying Mikhlin's suppression factor, W b , with the expression for D surf b . The expression for the bubble diffusion coefficient including Mikhlin's suppression factor becomes:

D surf b = 3 2π Ω 4/3 r b 4 D s 1 - q V b N b (9)
So, the bubble diffusion by surface diffusion mechanism use Eq.9 for the diffusion coefficient with Eq.7 for the value of D s .

Model description

To deal with the particular problem of bubble movement in a vacancy concentration gradient, the following assumptions and conditions have been adopted in the model [START_REF] Verma | A new spatially resolved model for defects and fission gas bubbles interaction at the mesoscale[END_REF]:

• We assume that the simulation starts after gas has been trapped in the bubbles. So, the Xe gas atoms are present inside the cavities and not in the solid. We also think that the irradiation defects largely annihilate with the temperature and the few that remain would be trapped in the extended defects. This is why we only consider thermal defects in the model.

• The UO 2 fuel is modeled as a mono-crystal containing spherical cavities.

• Each UO 2 is considered as an "atom" and Schottky defects as "vacancies" in the description of the model.

• Only vacancies are considered as point defects in the model for the present study. Indeed the auto-interstitial U or the anti-Schottky defect have large formation energies and can hardly be produced thermally.

• Bubbles are assumed to be spherical and remain so. This assumption is justified for the intra-granular bubbles at high temperature.

• No irradiation is considered in the model, i.e., it is applicable for annealing tests.

• The sources or sinks of vacancies are the surfaces in the model -surfaces of the bubbles and surface of the grain. In reality, dislocations could also be sources or sinks of vacancies, however, they have not been considered in the model presented here.

• The model is developed to function in both 2-D and 3-D, but for the results presented in this paper, only 2-D analysis has been done due to computation cost.

In the representation of the model, the domain is discretized on a regular mesh. Two types of independent variables are sufficient to represent the model and these are the spins and concentration fields. The first variables, spins, are integer numbers (0, 1 or 2) allocated to the cells (pixel in 2-D/voxel in 3-D) to categorize their position in the domain. The spin is 2 for the cells in the solid region, 0 for the cells in a cavity and 1 for the cells at the interface between solid and cavity. In the following sections, "2-cell", "1-cell" and "0cell" will refer to cells tagged with spin 2, 1 and 0, respectively. The 1-cells are partially solid and partially void.

The other variables, the concentration fields, represent the concentrations of different species in the model and are calculated by a system of ordinary differential equations. The field for a 2-cell is the concentration of vacancies, C v . For a 1-cell the field is the concentration of crystalline atoms, C a . All the concentrations are expressed in fraction of sites. C a in particular, is the number of U atoms in the cell divided by the number of U-sites in the cell (itself equal to the volume of the cell divided by Ω). C a can also be assimilated to the solid fraction (RS) in the interface cell. No fields are required for the 0-cells.

The limit condition imposed for the diffusion of vacancies in the solid is that the concentrations of vacancies at the interface cells will be the equilibrium concentrations of vacancies in the vicinity of the bubble that these cells are an interface of. These equilibrium concentrations of vacancies at the considered temperature in the proximity of a surface with curvature κ and exposed to a pressure P b are denoted as C eq v (P b , κ). The method to calculate the equilibrium concentrations of vacancies in the vicinity of over-pressurized bubbles has been presented by Noirot [START_REF] Noirot | A method to calculate equilibrium concentrations of gas and defects in the vicinity of an over-pressured bubble in UO 2[END_REF] and the expression for this is:

C eq v (bubble) = exp - ( v -s v T ) kT - Ω kT (P b -γ b κ) (10)
where, v is the formation energy of vacancies (more precisely, Schottky defects), s v is the excess entropy of vacancy formation, T is the annealing temperature in K, k is Boltzmann constant, Ω is the volume of one UO 2 -site and γ b is the surface tension at the solid-bubble interface (γ b depends on the temperature by the relation

γ b = 0.41 * (0.85 -1.4 * 10 -4 (T -273)) in J/m 3 [35]).
The excess entropy of vacancy formation, s v , is assumed purely vibrational and is taken as s v = 0 for this study. The expression for C eq v (bubble) in Eq.10 can be split in two where the first factor is the concentration of vacancy at equilibrium in a solid limited by a flat surface (curvature, κ = 0) and no external pressure (P b = 0):

C eq v (bubble) = C eq v (P b = 0, κ = 0)exp -Ω kT (P b -γ b κ) (11)
The pressure P b in the bubble is calculated using the Carnahan-Starling Equation of State [START_REF] Carnahan | Equation of state for non-attracting rigid spheres[END_REF] given as:

P b V N kT = (1 + η + η 2 -η 3 )/(1 -η) 3
where η is the packing fraction given by the expression:

η = π 6 N d 3 V
with d being the hard sphere diameter, calculated from the modified Buckingham potential as [START_REF] Brearley | An improved equation of state for inert gases at high pressures[END_REF]:

d = 4.45 * 10 -10 * 0.8542-0.03996 ln T 231.2
The bubbles in the model are characterized by the position of their centers, the total volume of the 0-cells included in the bubbles, the total volume of the bubbles, the gas content and the list of 1-cells forming the interface. In this study, we have considered the presence of gas in the bubbles, but we do not consider gas in the solid. We can define a bubble entirely within a cell. For such a bubble, the volume of the 0-cells would be 0 and the spin of the cell containing the bubble would be 1.

Methodology adopted in the model

The model simulates micro-structure evolution with time by taking into account the various aspects of interaction between defects and bubbles such as diffusion of defects and the update of bubble characteristics. A general algorithm adopted in the model is depicted in Fig. 2.

Note that we use an adaptive time step for the "macroscopic time", i.e., the one between two updates of the bubbles. The second condition for adapting the time step, namely the relative change in C veq for any bubble is the most demanding condition. Moreover, the diffusion calculation of vacancies is done using "microscopic time steps" during each "macroscopic time step". This microscopic time step must respect the Courant-Friedrichs-Lewy (CFL) condition for the numerical scheme that we adopted (discussed later).

(i) Initialization

The initial and boundary conditions as well as the physical properties are provided into the model. Table .1 describes the physical parameters considered in the model. The initial representation can be viewed using the visualization application PARAVIEW [START_REF] Ayachit | The ParaView Guide: A Parallel Visualization Application[END_REF]. In this paper, all the visualizations are represented in terms of RS, which is the solid ratio/fraction of a particular cell. RS = 1 (red region) means that the cell is completely solid and RS = 0 (black region) means the cell is completely void (cavity), with values of RS in between signifying an interface cell which is partly solid and partly cavity.

(ii) Diffusion of vacancies and crystal atoms

The diffusion of vacancies and crystal atoms is governed by the Fick's law. A balance equation is written for each 2-cell and 1-cell interacting

D U Self-diffusion coefficient [22] 5.82*10 -20 m 2 /s D v Vacancy diffusion coefficient 2.52*10 -13 m 2 /s † D U = C eq v (P b = 0, κ = 0)*D v *f T Isothermal annealing temperature 1600 • C kT Energy term 2.58*10 -20 J γ b
Surface tension at the solid-cavity 0.25666 J/m 2 interface D s Surface diffusion coefficient [START_REF] Matzke | Surface diffusion and surface energies of ceramics with application to the behavior of volatile fission products in ceramic nuclear fuels[END_REF] 1.39*10 -11 m 2 /s q Region without gas atoms 1.5*10 -27 m 3 in Mikhlin's suppression term Table 1: Physical parameters used in the model; † f, which is the correlation factor, is taken as f = 1 in this study. (Temperature dependent parameters are calculated at T = 1600 • C).

with its 2-cell neighbours as:

∂C j ∂t = - - → ∇. - → φ j (12)
Here -→ φ j represents the flux of the diffusing species, where j is either vacancy (when a 2-cell is considered) or crystal atom (when a 1-cell is considered). Fig. 3(a) represents the 0-cells (black), 1-cells (yellow) and 2-cells (red) in a domain. Eq.12 is solved for vacancies in the red domain with Dirichlet boundary conditions imposed on the yellow cells. It means that for this problem, there is a variable (C v ) for each 2-cell, while for the 1-cells, C v is imposed to C eq v (bubble). For the crystal atoms, Eq.12 is solved in each yellow cell with the Neumann boundary conditions. Since, we consider only vacancies as diffusing species, exchanging with a 2-cell, the flux of atoms is a consequence of the flux of vacancies and can be determined as:

- → φ a = - - → φ v (13)
On the red-yellow interface, the flux of atoms is imposed to the opposite of vacancy flux (Eq.13) and on the other interfaces, the flux of atoms

0-cells (empty) 1-cells (interface) 2-cells (solid)
𝐶 " #$ (Bubble 1) is nil (Fig. 3(b)). Eq.12 for the atoms of a 1-cell is just the integration with time of atoms that accumulated in the cell due to the emission of vacancies by the bubble. These equations are then discretized in space and time. In the present model, only the face neighbors have been taken into consideration using a 5-point stencil in 2-D for space discretization. Once the set of spacediscretized equations is written, a set of ordinary differential equations is obtained. An Euler explicit scheme is then used for time discretization.

𝐶 " #$ (Bubble 2) (a) (b) ∅ ' = -∅ " ∅ ' = 0

(iii) Bubble volume calculation

After the diffusion process, the total volume of the bubbles has to be calculated for each bubble. This is done using the volume of the 0-cells and the C a that has been calculated for each interface cell of the bubble. Then knowing the volume of each bubble, the temperature and the equation of state of the gas, we can determine the pressure, P b , as well as the new value of C eq v . The volume of the bubble and the relative change in the value of C eq v is checked and if the conditions are not satisfied (see Fig. 2), the time step is reduced by half and the diffusion calculation is repeated to meet the conditions for bubble volume and C eq v . (iv) Updating the new center of bubbles (target center)

If the flux of vacancies is not isotropic, the bubble is shifted from its position. The direction of bubble movement depends on the direction of vacancy absorption or removal from the bubble surface. The new center of the bubble is calculated as the barycenter of voids (i.e., the 0-cells in the bubble and the void fraction (1-RS=1-C a ) of the 1-cells). Since at each time step, the C a values are modified during the diffusion process, the new (target) centers of the bubbles are calculated accordingly. (v) Re-drawing the bubbles

Once the new characteristics of the bubbles are updated, the list of interface cells of each bubble also has to be updated. This means that the spins of some cells may be changed and the RS = C a (solid ratio in the interface cells) are re-calculated as if the bubble was re-drawn. The atom balance is checked and any error is accounted for during the volume update. If the target center of the bubble is too far from its previous center (meaning that their distance is larger than h), this process of re-drawing may be done in several steps, moving the bubble by shorter distances to mimic the simultaneous movement of bubbles.

Coalescence with other bubbles is also checked during the process. Once the re-drawing process is done for each bubble, a new time step can begin for the diffusion process.

(vi) Random movement of the bubbles

To incorporate the random movement of bubbles, the diffusion coefficient is chosen depending on the mechanism adopted (volume or surface mechanism). In 2-D, the most probable diffusion distance between the position of the bubble at time t = 0 and its position at ∆t is given by the expression:

d = 4D b ∆t (14)
where D b is the bubble diffusion coefficient and ∆t is the elapsed time.

We then choose a direction of movement by randomly generating an angle, φ [0, 2π[. However, we want the diffusion distance for a single jump to not exceed the length of one cell (defined as 'h' in the model). So, in order to restrict our jump distance to h, we impose another limit on the time step that must not exceed:

∆t = h 2 4D b (15)
The random displacement of the bubble is added to the directed displacement calculated in (iv) to determine the target center of the bubble. The bubble is then re-drawn in the way described earlier.

Verification

Verification for diffusion calculation

In order to verify the diffusion calculation in the model, a test was carried out in which the vacancies with an initial sinusoidal concentration profile were made to diffuse in the solid and the results obtained numerically were compared with the analytical results. The diffusion of vacancies in the solid follow the Fick's second law which is described by the equation:

∂C v (x, t) ∂t = D v ∂ 2 C v ∂ 2 x ( 16
)
An analytical solution of Eq.16 for a sinusoidal vacancy concentration profile can be obtained as:

C analytic v (x, t) = exp - 2πn a 2 D v t sin 2πn a x (17)
For the model, an initial concentration of vacancies was provided as a sinusoidal profile. A rectangular domain of size 128nm*64nm was taken with a grid size of 4 nm and ∆t = 0.99 ∆t CF L (to respect the CFL condition). Periodic boundary condition was imposed at all the boundaries. No bubbles were considered within the box for the test. One thing to note here is that this verification test is purely mathematical in the sense that C v (x, t) can take negative values (due to the sinusoidal profile).

The time evolution of vacancy concentration along the length of the domain was analyzed. As expected, with time the vacancies diffuse in the solid and the concentration profile becomes flatter and tends to zero with time. We defined the error as:

error = N C analytic v -C num v 2 N
where N is the number of cells in the domain. This error never exceeds 4.06*10 -4 in the calculation. Then we perform an analysis of the order of convergence in space and time. To determine the order of convergence in space, we choose a small ∆t and perform several calculations by varying 'h'. Each time the error defined above is calculated at the same time of C v evolution. Then we plot log(error) as a function of log(h). As expected for our numerical scheme, we obtain a straight line with a slope close to 2.0 (Fig. 4).

The dependence of error in function of h is found to be 10 -4.6864 * h 2.0176 . To determine the order of convergence in time, we perform several calculations with different ∆t and a fixed h. We plot log(error -10 -4.6864 * h 2.0176 ) in function of log(frac ∆t CF L ), with frac ∆t CF L = ∆t/∆t CF L . Fig. 5 shows the plot and as expected, we obtain aligned points with a slope close to 1.

The tests of order of convergence in space and time are found to behave as expected for the numerical scheme we adopt in this model.

Crystal atom balance

The process of re-drawing of bubbles is done using a recursive procedure to identify the new interface cells and their new solid ratio (RS=C a ). This can be done more or less precisely, but a high precision costs computation time. We used a procedure that keeps track of any error balance and avoid errors to cumulate with the time steps. Due to this procedure, the relative error in the crystal atom balance in the model is less than 4.5*10 -15 during the entire calculation.

Coalescence and space discretization

In the present formulation of the model, the problem of the coalescence of bubbles is not completely independent of the spatial discretization as the coalescence condition is "If two interface cells are neighbours and they belong to two different bubbles, then the two bubbles have to coalesce". Another approach that could be adopted to deal with the coalescence is to define a physical distance of coalescence, d coal (minimum acceptable distance between the bubble surfaces), and write the condition for coalescence differently. For instance, if d coal > 2h, then at least the second shell around each interface cell should be analyzed to check for coalescence. With such an algorithm, the solution, including the coalescence, should remain stable even if h is diminished. The development and tests of such an algorithm is beyond the scope of this work, but is feasible.

Our calculations, with h = 4 nm and the present algorithm corresponds to a d coal ∼8 nm which is reasonable since d -3 coal = 1.95 * 10 24 m 3 is of the order of magnitude of the maximum bubble densities observed.

Verification for bubble movement

As discussed earlier, a pressurized bubble will move up in a vacancy concentration gradient by trapping vacancies emitted by a free surface at lower pressure or higher curvature. In order to verify that the bubble movement simulated by the model is accurate, we carry out a test for the bubble movement in a vacancy concentration gradient and compare the velocity of the bubble to the theoretical value obtained from literature. For the test, we consider a domain of 128nm*32nm with a grid size of 4 nm. We consider two exterior regions, Region 1 and Region 2, with width 1 nm and 255 nm, respectively, on both ends of the domain box and a bubble of radius 7.15 nm within the box. In the model, the interfaces solid-Region 1 and solid-Region 2 are dealt with exactly as the solid-bubble interface. This means that a concentration C a of the crystal atoms is calculated at the interface cells and the surface can move accordingly. The only difference with the bubble surface is that this surface remains flat. We provide periodic boundary condition on the top and bottom and symmetric condition on the right and left boundaries.

We imposed constant pressures, P1 and P2, in the two regions to maintain a constant vacancy gradient in the solid throughout the duration of the test. 1 The properties of the different regions and the bubble for the test are presented in Table 2. A linear profile for initial vacancy concentration is provided. The evolution of the bubble with time is depicted in Fig. 6 in terms of RS which is the solid fraction/ratio in a cell. The figure shows that the bubble moves towards the region 2 by trapping vacancies which are emitted by the surface of Region 2, and by emitting vacancies towards the Region 1 which traps the vacancies arriving at its surface. The bubble also grows by trapping the vacancies. In the meanwhile, Region 2 reduces in size due to the loss of vacancies from its surface and the contrary for Region 1. The bubble can be seen to move up in the vacancy concentration gradient as expected. We compare the numerically calculated velocity of the bubble to the theoretical value. The velocity of a bubble moving in a vacancy concentration gradient is given as:

- → v b = 2D v - → ∇C v (18)
where, D v is the vacancy diffusion coefficient and -→ ∇C v is the vacancy concentration gradient (see Appendix A). Fig. 7 shows the theoretical as well as the calculated (numerically) values of velocity as a function of time. The theoretical value of velocity is a constant since we assume/impose a constant gradient of vacancy concentration. It can be seen from the figure that the numerically calculated value of bubble velocity is in quite good agreement with the theoretical value considering that no adjustment was made to obtain the numerical value. Moreover, the numerical velocity tends to the theoretical one when the grid is refined (from h = 4 nm to h = 2 nm).

Results and discussion

2-D analyses were carried out to determine the fission gas release due to intra-granular gas bubble movement to the grain boundary. First, the case with bubble movement in a vacancy concentration gradient was analyzed to assess Evans' mechanism. Later, the influence of the random movement of bubbles due to volume as well as surface diffusion mechanisms was analyzed without the diffusion of vacancies. Moreover, the overall fission gas release was analyzed by incorporating the random movement of the bubbles along with the directed motion in the vacancy concentration gradient. Finally, the influence of Mikhlin's suppression term on the overall FGR was assessed.

Fission gas release values from experiments for comparison

The analyses carried out are presented in this section with the primary focus being on the fractional FGR due to intra-granular bubble migration. We compare the values of FGR obtained from our analyses to the results from the experiments presented in the thesis of Valin [START_REF] Valin | Etude des méchanismes microstructuraux liés au relâchement des gaz de fission du dioxyde d'uranium irradié[END_REF]. She used a pellet of UO 2 fuel (enriched to 4.5% of U-235) from a PWR which had been irradiated for one cycle. The mean grain size was 11 µm. At the end of the irradiation, the burnup was 14.2 GWd/t U and there was nearly no FGR. The fragments of the whole pellet were annealled for 3 hours at 1600 • C. During the annealing, the specimen was swept by helium at a flow rate of 60 cm 3 /min. The Kr 85 release was continuously monitored by a spectrometer and the total quantity released was also determined by analyzing cold traps. The fraction of FGR was found to be ∼65%.

The value of FGR obtained from the experiments of Valin is, however, the overall FGR, from the intra as well as inter-granular bubbles. Earlier studies had been carried out using the MARGARET model [START_REF] Noirot | MARGARET: An advanced mechanistic model of fission gas behavior in nuclear fuel[END_REF] to discriminate the intra-granular and inter-granular retained gas. At the end of one cycle irradiation, we get approximately 10% of inter-granular gas. Supposing that the inter-granular gas has been completely released, leads to 55% of the total created gas coming from the grain in the FGR. This would mean that 0.55 0.9 = 61% of the intra-granular gas has been released during this annealing test. Moreover, for our analysis in this paper, we use a planar geometry rather than the actual spherical geometry of the grain. Going from spherical to planar geometry, we get an equivalent value for the FGR as ∼27% as compared to ∼61% FGR from intra-granular bubbles (Appendix B). So, our reference value for FGR from intra-granular bubbles during post-irradiation annealing of UO 2 for 3 hours at 1600 • C is 27% for comparison with the calculations.

Conditions for numerical analyses

For all the tests, we consider a planar domain of 5120nm*128nm with a grid size of 4 nm (i.e., the mesh is 1280*32 cells). We consider this domain size, as it represents the actual grain size radius (∼ 5µm). We consider a free surface and an exterior region at a low pressure and a large number of overpressurized bubbles within the grain. We try to respect a distance between the bubbles while generating the bubbles randomly in the grain that would be consistent with a bubble number density in the range 10 23 -10 24 bubbles/m 3 as observed in experiments [START_REF] Cornell | The growth of fission gas bubbles in irradiated uranium dioxide[END_REF][START_REF] Baker | The fission gas bubble distribution in uranium dioxide from high temperature irradiated SGHWR fuel pins[END_REF][START_REF] Kashibe | Formation and growth of intragranular fission gas bubbles in UO 2 fuels with burnup of 6-83 GWd/t[END_REF]. We reached an equivalent bubble density of ∼ 2.8*10 23 bubbles/m 3 . The bubble size is distributed randomly with the radius, r (0.4, 1) nm [START_REF] Kashibe | Formation and growth of intragranular fission gas bubbles in UO 2 fuels with burnup of 6-83 GWd/t[END_REF]44]. The volume per atom of gas within the nanobubbles was taken equal to Ω. This value is justified from the measurements of Nogita and Une [START_REF] Nogita | High resolution TEM observation and density estimation of Xe bubbles in high burnup UO 2 fuels[END_REF] and Thomas [START_REF] Thomas | Condensed-phase xenon and krypton in UO 2 spent fuel[END_REF], who measured Xe density inside intra-granular bubbles in UO 2 after base-irradiation in function of the bubble diameter. For 1-2 nm size bubbles, a value of 5g/cm 3 is consistent with these measurements. The mean mass of Xe produced by irradiation is 134.2g, which leads to a volume per atom close to Ω. We provide periodic boundary condition on the top and bottom boundaries and symmetric condition on the right and left boundaries. An isothermal annealing temperature of 1600 • C is used for the analysis. A domain with the provided conditions is depicted in Fig. 8 and the initial conditions used for the analyses are described in Table 3. Since the domain used for the analyses is quite large as compared to the bubbles, in order to have a better view of the distribution of the bubbles, we represent only an enlarged region of the domain (represented by dashed box in Fig. 9) towards the vicinity of the exterior surface in all the cases. We calculate the fractional FGR due to intra-granular bubble migration by the increase in the net content of gas in the exterior region before and after annealing, divided by the total gas contained initially in the bubbles.

5120

Enlarged region for analysis

Fission gas release due to directed bubble motion in a vacancy concentration gradient

Evans [START_REF] Evans | Bubble diffusion to grain boundaries in UO 2 and metals during annealing: A new approach[END_REF] had proposed that the directed movement of intra-granular bubbles in a vacancy concentration gradient could account for the large gas release observed during post-irradiation annealing tests. In order to test Evans' mechanism for the fission gas release, we consider a constant initial vacancy concentration of C ini v = 1.08*10 -12 in the solid, which is of the same order as the C eq v in the vicinity of the small bubbles2 . The simulation is carried out for a simulated time of 3 hours. We expect the bubbles in the vicinity of the exterior surface to grow by trapping the vacancies emitted from the surface and also to move towards the exterior surface the vacancy gradient. The distribution of the bubbles near the exterior surface after 3 hours is depicted in Fig. 10. Fig. 10 shows an enlarged view of the bubbles in the domain near the exterior surface at time t = 0 and at t = 3 h. It can be noticed that the bubbles grow by trapping the vacancies which are generated at the free surface. Some of the bubbles also move out of the grain into the exterior region. The number of bubbles in the domain decreased from 2647 bubbles at time t = 0 to 2634 bubbles at t = 3 h. The fractional FGR after 3 hours of annealing was found to be 0.112%, which is insignificant as compared to the values obtained from experiments if we compare it to our reference value (27%).

The analysis was continued for annealing up to 110 hours of simulated time in order to see any further evolution of the bubbles and the effect on the overall fission gas release. The bubbles grew further and up to a wider extent into the grain. The bubble evolution from time t = 0 to t = 110 h is depicted in Fig. 11. The number of bubbles further reduced to 2266 at time t = 110 h. However, the fractional fission gas release went up to only 1.204%, which is still insignificant as compared to the values obtained from experiments. The values for various parameters at time t = 0, 3 and 110 h for the case are presented in Table 4 and the mean radius of bubbles along the domain is depicted in Fig. 12. An interesting question is to evaluate the impact of D v on this result. In practice, we solve the equations using a non-dimensional time constructed as with l adim = 1 nm and the non-dimensional length as

t = 0 t = 3 h
t = 0 t = 3 h t = 110 h
t * = D v l 2 adim *
l * = l l adim
With this choice, the non-dimensional coefficient of vacancy diffusion, D * v = 1. This was done in order to accelerate the calculation of the diffusion (by avoiding multiplication by D v ). The interest in the question is also that in this case of directed movement in a vacancy gradient, where the kinetics is solely governed by D v , changing D v means that the t * is related differently to the real time. In this analysis, we showed that the vacancy gradient became almost flat before any significant gas release could take place. Taking another value for D v only changes the axis of time by a multiplying factor, however, the conclusion remains the same. So, having conducted the calculation up to a long duration is equivalent to having done a parametric study on D v , for the case of directed bubble movement and there is no impact on the results whatsoever.

In conclusion, it is evident from this analysis that the FGR from the grain cannot solely be due to the transport of gas via intra-granular bubble movement in a vacancy concentration gradient. Other mechanisms of gas release have to play a significant role in the overall fission gas release during post-irradiation annealing of UO 2 .

Influence of random bubble movement without vacancy diffusion

The influence of random movement of bubbles within the grain on the evolution of bubbles was tested. We present the analysis of random bubble movement via volume and surface diffusion mechanisms in the following section. In order to emphasize on the random movement only, we do not consider any driving force (diffusion of vacancies) in the solid. The simulation is carried out for a simulated time of 3 hours. To account for the variability due to randomness, all calculations are carried out 10 times and the mean of these values is used as the results.

Random movement of bubbles due to volume diffusion mechanism

The random movement of the bubbles by volume diffusion mechanism is governed by the bubble diffusion coefficient, D vol b , which is defined by the Eq.5. It was observed from the analysis that there was negligible movement of bubbles via volume diffusion mechanism. After 3 hours of annealing, the number of bubbles remained the same, 2647, as before. The mean value for the diffusion coefficient of bubbles, D vol b also remained unaltered at a very low value of 2.10*10 -27 m 2 /s. So, it could be concluded that there is no notable contribution of volume diffusion mechanism on the random movement of the bubbles and, thus, no additional impact on the fractional FGR.

Random movement of bubbles due to surface diffusion mechanism

The random movement of the bubbles by surface diffusion mechanism is governed by the bubble diffusion coefficient, D surf b , which is defined by the Eq.9. Using the value of q as 1.5*10 -29 m 3 , as proposed by Mikhlin for UO 2 , to evaluate the value of W b and multiplying the expression for D surf b with W b , we obtain the values for the diffusion coefficient of bubbles. The values of diffusion coefficient with the suppression factor are found to be very low and the mean value of D surf b at t = 0 is 3.79*10 -31 m 2 /s. There is no movement of bubbles because of such low values of the diffusion coefficients. After 3 hours of annealing, the number of bubbles remained the same, 2647, as before. The mean value for the diffusion coefficient of bubbles, D surf b also remained unaltered at a very low value of 3.79*10 -31 m 2 /s. So, from this analysis of the random movement of bubbles, we can conclude that neither the volume diffusion mechanism nor the surface diffusion mechanism can solely explain the observe FGR. However, since the Mikhlin's suppression term in the bubble diffusion coefficient suppresses the diffusion of small bubbles, an interesting scenario can further be considered if we provide the diffusion of vacancies into the grain and allow random motion of bubbles via surface diffusion including the Mikhlin suppression term. As the bubbles trap vacancies and grow, the suppression of diffusion of small bubbles would reduce. This would eventually enhance the diffusion of such bubbles and account for a higher fractional gas release.

Fission gas release due to combined random and directed bubble movement in a vacancy concentration gradient

We now try to analyze the overall fractional FGR via intra-granular gas bubble movement, taking into consideration the directed movement in the vacancy concentration gradient as well as the random movement of the bubbles due to surface diffusion mechanism including Mikhlin's suppression factor.

For incorporating the diffusion of bubbles due to surface mechanism, we determine the diffusion coefficient of bubbles by Eq.9. There exists a vacancy concentration gradient between the exterior surface (at low pressure) and the over-pressurized bubbles within the grain. The bubble distribution for this case at time t = 0, 3 and 110 h is depicted in Fig. 13.

As can be observed from Fig. 13, at time t = 3 h, the bubbles at the vicinity of the exterior surface have moved significantly as compared to time t = 0. At time t = 110 h, this bubble movement is extended further into the solid. The number of bubbles in the domain decreased from 2647 bubbles at time t = 0 to 2560 bubbles at t = 3 h and to 1879 bubbles at time t = 110 h. The mean diffusion coefficient of bubbles increased from its initial value of 3.76*10 -31 m 2 /s at t = 0 to 6.69*10 -23 m 2 /s at 3 hours and reached 3.91*10 -23 m 2 /s after 110 hours. The fractional FGR increased to 0.84% in 3 hours and to 2.32% in 110 hours. We notice that the bubble diffusion due to surface diffusion was enhanced as the vacancies were trapped by the bubbles and caused significant bubble movement. The reduction in number of bubbles was, however, attributed majorly to the coalescence of bubbles rather than the transfer of bubbles outside the grain and thus, the FGR is still too low. The values of various parameters at time t = 0, 3 and 110 h for the case are presented in Table 5 and the mean radius of bubbles along the domain is depicted in Fig. 14. The bubble diffusion coefficient by surface diffusion using Mikhlin's suppression term depends on two opposing factors depending on the radius of bubbles. The first term in Eq.9, reduces the bubble diffusion coefficient as the radius of the bubble increases. At the same time, increasing the radius causes the suppression factor term (W b ) to increase making the bubble diffusion coefficient to increase. Fig. 15 depicts the overall mean diffusion coefficients of bubbles at different distances within the domain. It can be observed that for time t = 3 h, the mean diffusion coefficient of bubbles in the vicinity of the exterior surface increases. This implies that the suppression factor term dominates the first term in the bubble diffusion coefficient expression. The extent of enhanced bubble diffusion coefficient is even further into the grain as time goes to t = 110 h. We note that the bubble diffusion coefficient enhancement due to increasing radius continues even after 110 hours of annealing time. This means that the contribution of the suppression factor term in aiding the bubble diffusion continues even at 110 hours of annealing, however, the fission gas release observed in 110 hours is merely 2.32%, which is not enough compared to the value from experiments.

Finally, we tried to continue the calculation until the concentration of vacancies in the solid reaches equilibrium value. This was done to check the extent to which the FGR values reach as long as the vacancy concentration gradient lasts. The evolution of vacancy concentration in the domain length with time is plotted in Fig. 16. It can be noted from the figure that the vacancy concentration values have not completely reached the equilibrium value. The simulation was discontinued after a simulation time of t = 3080 h because immediately after t = 3080 h, the two largest bubbles closest to the exterior surface coalesced with each other and formed a bubble with diameter larger than the width of the domain. The model did not consider any condition for such a scenario and thus, the simulation could not be continued further. At t = 3080 h, only 20 bubbles remained in the domain with a mean bubble diffusion coefficient of 1.39*10 -21 m 2 /s and a mean radius of 20.4 nm. At the end of the simulation, the vacancy concentration at the rightmost end of the domain was 2.31*10 -7 and at the leftmost end was 1.63*10 -7 . The FGR from the grain in 3080 hours was still found to be 2.158%, same as that at 110 hours. Note that this value is lower than the one presented at t = 110 h (2.32%) as the calculations up to t = 110 h were done 10 times to account for the variability due to randomness and the mean of these calculations was used as the FGR value, whereas, only one of these calculations was continued for the longer duration and it had a FGR of 2.158% at t = 110 h.

The bubbles grew by trapping the vacancies and coalesced among themselves but did not move out of the grain, thus, not contributing to the FGR. The bubble distribution after time t = 110 h and t = 3080 h is shown in Fig. 17. Even if we assume that the large bubble formed after coalescence moves out of the grain at some point, the FGR would reach 18.5%. So, even 0.00E+00 though the simulation could not be carried out till the equilibrium value for vacancy concentration was reached, it is still conclusive that the extent of FGR would still be less than the equivalent target value of ∼27% from experiments. Hence, from the analyses carried out in this paper, it can be concluded that, as far as crystal atoms diffusion at the surface and the Mikhlin's suppression term are reliable, the gas bubble migration in a vacancy gradient and its random movement cannot solely be responsible for the high value of fission gas release observed during post-irradiation annealing tests. This result does not depend on the diffusion coefficient of vacancies since we showed that even if the calculation is performed until vacancy equilibrium is reached, the FGR is too low as compared to the experiment.

Parametric investigation of Mikhlin's suppression term

In order to study its influence on the diffusion of bubbles, the Mikhlin's suppression term was excluded from the expression of bubble diffusion coefficient. This was done by using the value of q in Eq.8 as 0, making W b =1. The resulting expression which is given by Eq.6 is used to calculate the diffusion coefficient of bubbles. Here again, we do not consider the driving force of vacancy diffusion and consider only random movement by surface diffusion mechanism. The simulation was carried out for 3 hours of annealing. It was observed that the mean value of D surf b was very large. At time t = 0, the mean value of D surf b was found to be 6.04*10 -13 m 2 /s. Due to the high diffusion coefficient, the bubbles were observed to move quickly and coalesce with each other and the number of bubbles reduced from 2647 to just 5 bubbles at time t = 3 h. The D surf b values reduced gradually as the bubbles grew by coalescence and the mean D surf b at time t = 3 h was found to be 1.46*10 -18 m 2 /s. The mean radius of the remaining 5 bubbles was found to be 16.33 nm. The fractional FGR in this time was calculated to be 6.64%. The bubble distribution via surface diffusion mechanism at time t = 3 h in the entire domain is depicted in Fig. 18. Since the FGR values were comparatively higher for this case in 3 hours of annealing, this case was continued further for a longer duration. The simulation was carried out up to the time t = 5500 h. It was, however, realized that the remaining 5 bubbles also coalesced with each other leaving a single large bubble to move within the grain. This bubble can no longer coalesce with any other bubble, however, in principal, due to the periodicity conditions along the top and bottom boundaries, there exist bubbles similar to the one large bubble which it can interact with. So, in order to encounter the problem in a more realistic manner, we consider a domain 10 times larger in the vertical direction. For the new case, we consider the initial conditions to be the same as at time t = 3 h. The new domain is a 5120nm*1280nm box with grid size of 4 (i.e., a mesh of 1280*320 cells). Using the periodic conditions, the initial number of bubbles for the new case are 50. The mean radius of the 50 bubbles is 16.3 nm and the mean bubble diffusion coefficient is 1.46*10 -18 m 2 /s. The initial conditions of the domain with the bubbles (at t = 0, after the first 3 hours) are depicted in Fig. 19. The simulation was carried out and we observed that the bubbles moved randomly within the grain. The bubbles majorly coalesced with each other and at time t = 110 h, only 12 bubbles remained in the grain (Fig. 19). However, no bubbles moved out of the grain and thus, did not contribute to any further FGR. As the simulation continued, further coalescence occurred and also some bubbles moved out of the grain. By the end of the simulation time of 5500 hours, only 3 bubbles remained within the grain. The mean radius of these bubble was 57.96 nm and the mean bubble diffusion coefficient was 1.085*10 -20 m 2 /s. The overall FGR at t = 5500 h after the first 3 hours was calculated to be 28.67%. This FGR had been achieved at ∼3366 h and then the bubbles continued to move randomly within the grain but did not contribute to further FGR. The bubble distribution at time t = 5500 h is depicted in Fig. 19.

This investigation shows that even after neglecting the Mikhlin's suppression term, the random movement of bubbles do not explain the observed FGR obtained after 3 hours of annealing. The target values for FGR are reached in 3366 hours, which implies that in order to attain the target FGR in 3 hours, the value of D s required would have to be ∼ 10 3 times larger with no Mikhlin's suppression of surface diffusion. However, in our opinion, doing so would lead us too far from the literature recommendations.

Conclusions

Using a new model for fission gas bubble and vacancy interactions, the intra-granular gas bubble movement was analyzed for its impact on the overall gas release during post-irradiation annealing tests. This kind of mechanistic evaluation had not been done before, to our knowledge. The model was verified for various aspects like grid sensitivity, diffusion calculation and bubble movement and was found to simulate the various aspects with good precision. Intra-granular gas bubble migration and fission gas release was analyzed via the mechanism of directed movement in a vacancy concentration gradient as well as for the random movement without any driving force. A combined effect of directed and random movement due to surface mechanism was also analyzed. The conclusions for the different cases are as follows:

• The fission gas release due to directed movement in a vacancy concentration gradient was found to be 0.11% in 3 hours and 1.2% in 110 hours of annealing at 1600 0 C. This is far less than the values of fission gas release observed during our post-irradiation annealing reference experiment.

• The random movement of bubbles without the diffusion of vacancies was tested for volume and surface diffusion mechanisms. The bubble movement was negligible due to both, the volume diffusion mechanism and the surface diffusion mechanism with the Mikhlin's suppression term.

• The fission gas release due to combined directed movement and random movement by surface diffusion mechanism including Mikhlin's suppression term was found to be 0.84% in 3 hours and 2.32% in 110 hours of annealing time. Even the combined movement mechanism could not account for the observed high values of fission gas release during the post-irradiation annealing test. Even if the calculation was continued till the vacancy concentration gradient lasted, the value of FGR did not increase beyond 2.32%.

• Parametric study of the Mikhlin's suppression term to test its influence on the overall FGR was done. It was found that the FGR is sensitive to the Mikhlin's suppression term. However, even after neglecting the the Mikhlin's suppression term completely, the random movement of bubbles did not explain the observed FGR obtained after 3 hours of annealing.

The analyses show that based on bulk diffusion coefficient of vacancies and the surface diffusion of crystal atoms, and on the Mikhlin's suppression term, the directed movement of bubbles in a vacancy concentration gradient cannot account for the large FGR which is observed during the post-irradiation annealing test of UO 2 sample that we took as a reference, even if the directed movement is assisted by random movement of bubbles due to surface diffusion. Further work would be necessary to get a better knowledge about the suppression of surface diffusion, which seems to be a sensitive parameter.

Other mechanisms for gas atom migration within the grain such as influence of dislocations need to be addressed in detail to understand the fission gas release during post-irradiation annealing tests.

• C v ( r) is the vacancy concentration (fraction of sites) at a particular position r

• C v0 is the vacancy concentration (fraction of sites) at the origin of the position vectors

• ∇C v∞ is the vacancy gradient far from the bubble

• C eq v (bubble) is the equilibrium vacancy concentration at the vicinity of the bubble. It is the limit condition of the vacancy concentration field at the bubble surface.

• r B is the position of the bubble center • R is the bubble radius

• r = r -r B
We can re-write the above equation by adopting C v0 for C v0 + r B • ∇C v∞ , and taking the center of the bubble as the origin. This gives:

C v ( r) = C v0 + (C eq v (bubble) -C v0 ) R r + (1 - R 3 r 3 ) r • ∇C v∞
where the definitions that have changed are below :

• r is the position of a point, the center of the bubble being the origin

• C v0 is the vacancy concentration (fraction of sites) that would be at the origin (center of the bubble) if the bubble was not there Now, if we consider the spherical coordinates (r, θ, ϕ), and the z-axis placed in the direction of ∇C v∞ , as in The gradient of the vacancy concentration field is the sum of the gradients of the T 1 and T 2 fields which can be described as: The growth rate of the bubble is defined by :

∇(T 1) = -(C eq v (bubble) -C v0 ) R r 2
∂V ∂t = -Ω S ϕ v • dS = D v S ∇(T 1 + T 2) • dS
where V is the volume of the bubble, t is the time, Ω is the vacancy volume,

ϕ v = -1 Ω D v ∇(C v ( r))
is the vacancy flux (m -2 ), dS is the elementary surface vector oriented towards the exterior of the bubble, and D v is the vacancy diffusion coefficient.

dS = R 2 sinθdθdϕ e 1
It appears that the contribution to the growth rate of the bubble is only due to the T1 field, and this gives :

∂V ∂t = 4πRD v (C v0 -C eq v (bubble))
for each time of the simulation. So the numerical bubble velocity is also calculated in the fixed referential.

In several articles, the bubble velocity is said to be equal to 2D v ∇C v∞ . As Geguzin explains, this is the velocity of the bubble in a system of coordinates linked to the matter far from the bubble, in the case where the solid moves like in the test case, because of the vacancy flux that goes through it.

For a numerical simulation in 3D, the numerical velocity should be compared with 3D v ∇C v∞ . However, the simulations presented in this paper are 2D. That is why we need to find a theoretical expression in 2D.

In 2D, if we consider a bubble that has reached equilibrium with the vacancy field at the position where its center stands (C eq v (bubble) = C v0 ), the solution for the vacancy field in 2D would be:

C v ( r) = C v0 + (1 - R 2 r 2 )rcosθ ∇C v∞ T 2(r,θ)
In 2D, we verify that :

(C v ( r)) = 1 r ∂ ∂r (r ∂T 2

Figure 1 :

 1 Figure 1: Directed movement of a bubble in a vacancy concentration gradient.

Figure 2 :

 2 Figure 2: General algorithm adopted in the model.

Figure 3 :

 3 Figure 3: (a) Representation of 0, 1 and 2-cells in the domain; (b) Flux of atoms at various interfaces.

Figure 4 :

 4 Figure 4: Determining the order of convergence in space.

Figure 5 :

 5 Figure 5: Determining the order of convergence in time.

Figure 6 :

 6 Figure 6: Bubble movement and growth in a vacancy concentration gradient.

Figure 7 :

 7 Figure 7: Theoretical vs Numerical value of bubble velocity.

Figure 9 :

 9 Figure 9: Dashed box representing the enlarged region of the actual domain presented in the results.

Figure 10 :

 10 Figure 10: Directed bubble movement in a vacancy concentration gradient for 3 hours of simulated time (enlarged view near the vicinity of the external surface).

Figure 11 :

 11 Figure 11: Directed bubble movement in a vacancy concentration gradient for 110 hours of simulated time (enlarged view).

Figure 12 :

 12 Figure 12: Mean bubble radius along the domain for directed movement case.

Figure 13 :

 13 Figure 13: Bubble movement in a vacancy concentration gradient with random movement of bubble by surface diffusion for 110 hours of simulated time (enlarged view).

Figure 14 :

 14 Figure 14: Mean Bubble radius along the domain for combined movement case.

Figure 15 :

 15 Figure 15: Mean Bubble diffusion coefficient for bubbles along the domain.

Figure 16 :

 16 Figure 16: Vacancy concentration values along the grain at different times.

Figure 17 :

 17 Figure 17: Bubble distribution in the entire domain for the extended simulation.

Figure 18 :

 18 Figure 18: Bubble distribution due to random movement by surface diffusion mechanism at time = 3 h.

Figure 19 :

 19 Figure 19: Random movement of bubbles via surface diffusion in the new domain at time t = 0, t = 110 h and t = 5500 h, after the first 3 hours.

r 3 )Figure A. 20 :

 320 Figure A.20: Spherical coordinates in 3D.

2 = 3 r 3 )sinθ e 2

 2332 ∇C v∞ (1 + 2 R 3 r 3 )cosθ e 1 -∇C v∞ (1 -R

Figure A. 21 :

 21 Figure A.21: Polar coordinates in 2D.

Table 2 :

 2 Conditions used in the model.

	Type	Coordinate Initial width/	Imposed
		of center	Radius (nm) pressure (Pa)
	Region 1	(0,0)	1	9.19*10 9 (P1)
	Region 2	(512,0)	255	9.04*10 4 (P2)
	Bubble	(16,64)	7.15	-

Table 3 :

 3 Initial conditions used for the analyses.

		Periodic condition
	Symmetric condition	Free surface
				low P	128 nm
		nm	Over-pressurized bubbles
	Property	Description		Value
	T	Annealing temperature	1600 • C
	N b	Initial bubble number density	2.8*10 23 bubbles/m 3
	w	Exterior region width	255 nm
	r b	Small bubble Radius	(0.4 -1) nm
	Vol/at (ext) Volume per atom of gas in the Exterior region	286 nm 3
	Vol/at (bub)	Volume per atom of gas in the bubbles	0.0409 (= Ω) nm 3
	t	Simulation time	3/110 h

w Figure 8: Domain used for the analyses with various conditions.

Table 4 :

 4 t Evolution of various parameters with time for the directed movement case.

	time Number of Mean radius	Mean gas	FGR
	(hr)	bubbles	(nm)	content (at)	(%)
	0	2647	0.719	167.664	-
	3	2634	0.756	168.303	0.11 %
	110	2266	0.966	193.496	1.20 %

Table 5 :

 5 Evolution of various parameters with time for the combined directed and random movement case.

	time Mean D b Number of Mean radius	Mean gas	FGR
	(hr)	(m 2 /s)	bubbles	(nm)	content (at)	(%)
	0	3.76*10 -31	2647	0.719	167.66	-
	3	6.69*10 -23	2560	0.743	171.65	0.84
	110	3.91*10 -23	1879	0.899	231.09	2.32

Note that here, the solid moves due to the vacancy flux (Nabarro-Herring creep) only. However, as P1 and P2 are different, they should induce a resulting effect that would shift the solid much more quickly towards the right, if the solid is not somehow held in an other way. Here, we can imagine that the solid is a section of a very large and thin tube, with P1 inside the tube and P2 outside the tube.

The implicit hypothesis that we make is that the vacancies and auto-interstitials from the irradiation mainly annihilate with one another at the very beginning of the annealing

Acknowledgements

The authors acknowledge AMU, CEA, EDF and Framatome for the financial support of this research conducted in the framework of the PLEIADES project. The authors are also thankful to Gérald Jomard (CEA) for the valuable discussions and encouragement.

Appendix A. Demonstration for the bubble velocity in a vacancy gradient

Following Geguzin et al. [START_REF] Geguzin | Migration of Macroscopic Inclusions in Solids[END_REF], the vacancy field in the permanent state around a bubble placed in a vacancy gradient is the following :

where Now, to calculate the bubble velocity, let us define the center of the bubble at t as C, the center of the bubble at t + dt as C , a current point on the bubble surface as M . C may be defined as the center of voids of the bubble at t and the different new voids added to the bubble in each point M of the surface. With an arbitrary fixed origin point O, this gives :

This may be transformed, using OM = OC + CM = OC + r, into :

The bubble velocity v is obtained by :

Because of the spherical symmetry of the field T 1, it does not contribute to the bubble velocity. The contribution from T 2 can be calculated using the expressions of ∇(T 2) and dS to calculate ∇(T 2) • dS in function of r and θ first, and then using the expression of r :

to integrate the limits on ϕ from 0 to 2π, and on θ from 0 to π. The contributions along i and j are nil because of the cylindrical symmetry of the problem.

Finally, the expression for the bubble velocity is:

Note that this velocity has been calculated in the fixed referential. In the model, we calculate the new center of the bubble using the numerical equivalent of the relation A.1, and the numerical bubble velocity has been calculated by v = CC dt based on the output file that gives the bubble center position

We also have

The growth rate of the bubble is nil, which is consistent with the fact that the bubble is already at equilibrium.

The velocity of the bubble is calculated exactly as above

Finally, the expression of the bubble velocity in 2D and in a fixed referential is :

This latter expression is the suitable one for comparison with our numerical velocity in the paper.

Appendix B. Target value of FGR (From spherical to planar geometry)

Following the idea of Evans [START_REF] Evans | The role of directed bubble diffusion to grain boundaries in post-irradiation fission gas release from UO 2 : A quantitative assessment[END_REF], if the spherical grain is divided into shells and the FGR comes from the outermost shell of the grain, then a sphere of radius R ext would be depleted of gas bubbles up to the radius R 1 (Fig. The volume of FGR from the sphere can be represented as:

This volume represents 61% of total gas release, so Now, if we consider the same R 1 and R ext in a planar geometry, then the fraction of gas released would be given by 1 -R 1

Rext and this value is 0.2694. So, for a 61% FGR in a spherical geometry, we can roughly have an equivalent FGR of ∼ 27% in the planar geometry.