Phase-field simulation of the effect of grain boundary on fission gas migration in UO2 fuel

Grain boundaries are widely recognized as line defect sinks. During reactor operation, vacancies and fission gas atoms within UO2 fuel grains migrate to the grain boundaries, forming bubbles at these locations. In order to better understand the effect of grain boundaries on the migration of fission gas atoms, this study employed a phase-field model to simulate the nucleation and growth process of fission gas bubbles within the UO2 fuel grains and in the vicinity of grain boundaries. This study also investigated the degree to which grain boundaries affect fission gas atoms under different temperature conditions. The results indicated that in models containing grain boundaries, the nucleation of fission gas bubbles occurred earlier, as compared to models without grain boundaries. A noticeable bubble denuded zone was also observed adjacent to grain boundary interfaces. Furthermore, with increasing temperature, the bubble denuded zone becomes thicker. The effects of irradiation, vacancy diffusion, Ostwald ripening, as well as grain boundary trapping were discussed.


Introduction
During the operation of the reactor, various fission products such as inert gases Xe and Kr, are continuously produced as fission progresses.These two gases produce approximately 0.25-0.3atoms per fission [1], and the production rate of Xe is nearly 10 times that of Kr [2].Due to the low solubility of these two gases in UO2, most of Xe and Kr diffuse to form bubbles within or between crystals [3], [4].As bubbles increase and grow, intracrystalline bubbles can cause swelling of fuel pellets [3], and intergranular bubbles connect to form channels, causing fission gas to diffuse into free space along the channels, resulting in the release of fission gas [4].The generation and release of fission gas is one of the key constraints for improving the burnup level of reactors.The accumulation of fission gas in reactor fuel rods, especially when the reactor burnup reaches a certain level, can cause an increase in internal pressure and stress in the reactor fuel rods,ultimately resulting in deformation of the fuel rod cladding and the rupture of the cladding, causing the release of radioactive materials [5].Therefore, studying the release laws of fission gases can assist in improving fuel design, enhancing reactor safety, and increasing overall reactor economics [6].Studies on the release behavior of fission gases often employ the rate theory model, which uses the spatial mean-field method to improve computational efficiency but ignores the evolution of various microstructure features.Compared to the rate theory method, the phase-field method can not only investigate thermodynamic and dynamic factors, but also track the dynamic evolution process of the model.Based on Ginzburg-Landau theory, it establishes a mathematical model that can accurately describe the evolution of the system over time through partial differential equations [7], to obtain the instantaneous state of the system in time and space.The phase field method avoids the tracking problem of complex interfaces by incorporating different interface descriptions into its formulas, and is easily coupled with other physical fields (such as noise field, velocity field, composition field, etc.) [8].There have been many related studies using the phase-field method to simulate the release behavior of fission gases in UO2 fuel.Hu et al. [9] constructed a phase field model containing vacancy concentration field and interstitial atom concentration field using the Cahn Hilliard equation.The model simulated the evolution of voids under irradiation conditions and obtained a series of parameters such as critical radius and quantity density of voids.Millett et al. [10] constructed a phase field model including the Cahn-Hilliard equation and the Allen-Cahn equation.This equation system distinguishes different phase boundaries by controlling order parameters, allowing for different dynamic and thermodynamic descriptions of different phases.Through this model, three stages of void development were obtained.In addition, the model also demonstrated that grain boundaries act as absorption sinks for defects and have strong absorption effects on vacancies and interstitial atoms.Void denuded and peak zones were observed near the grain boundaries.However, the parameters selected in this model are not for specific materials, and fission gas atoms were not taken into account when considering the role of grain boundaries as absorption sinks, which cannot quantitatively analyze the release law of fission gases near UO2 fuel grain boundaries.This study also used a phase-field model and referred to the above research.The nucleation and growth process of fission gas bubbles inside and near UO2 fuel grains under irradiation conditions were simulated using a parabolic approximation free energy functional, specific diffusion coefficient values, and micron scale grain size.The influence of grain boundaries on the release of fission gas in UO2 fuel was studied, and the degree of influence of grain boundaries on different temperatures was investigated.

Definition of independent variables
This model investigates the concentration distribution of point defect vacancies and fission gas generated by irradiation effects in polycrystalline materials.Two conserved field variables are introduced: cv(r, t), cg(r, t) which represent the vacancy mole fraction and fission gas mole fraction.This mole fraction c can be determined by c = ρVa, where ρ is the number density of vacancies or fission gas atoms, Va is the volume of U atoms in UO2 crystal cells, with a value of 0.0409 nm [11].As mentioned in the introduction, fission gases in fuel matrix forms bubbles through diffusion, and there are two stable phases: the solid fuel matrix phase and the fission gas bubble phase, using the order parameter η distinguish between the two phases: when η = 0 is the matrix phase, η = 1 is the bubble phase.There is an equilibrium concentration of vacancies and gases in the matrix phase and bubble phase.For the matrix phase, the equilibrium concentration of vacancies and gases is represented by the following equation: where f v E and f g E are the formation energies of vacancies and fission gas atomic, f v E = 5.1eV [12].The dissolution energy of Xe atoms in vacancies s g E =5.21eV [13], choose . kB is the Boltzmann constant, and T is the temperature.We take the gas equilibrium concentration , b eq g c = 0.494 [11], when cg = 0.494, the van der Waals free energy of the bubble phase is the lowest.Due to the fluctuation of gas concentration cg during the simulation process, in order to maintain cv+cg ≤ 1 in the bubble phase, the vacancy equilibrium concentration is taken as , b eq v c = 0.9cg, , b eq v c will decreases as cg increases.
This phase-field model focuses on the release behavior of fission gases in polycrystalline fuels, using order parameters ϕi (i=1, 2,..., P) distinguishes grains with different orientations.Within the grains with the i-th orientation ϕi=1, in other oriented grain phases ϕi=0.Near the contact surface of grains with different orientations, ϕi is approximately equal to 0.5.Through investigating Φ = , it is possible to distinguish between the interior and boundaries of the grain, when Φ = 1 is within the grain, Φ ≈ 0.5 is on the grain boundary.

Free energy functional
The phase-field method uses the total free energy functional F composed of the free energies of each phase and interface, which is consistent with the free energy of non-uniform systems defined by Cahn Hilliard [14].The free energy functional F form used in this study is as follows: Where f solid is the matrix phase free energy density function, f bubble is the bubble phase free energy density function, f poly is the polycrystalline free energy density function, and f grad is the gradient contribution term.h(η) = (1-η) 2 j(η) = η 2 is the weighting function.When η = 0, the contribution of f solid to the total free energy reaches its maximum, and f bubble has no contribution.When η = 1 the opposite is true.
To simplify the numerical calculation of the control equation, parabolic approximations were applied to the matrix phase free energy density function and the bubble phase free energy density function, which are similar to the free energy density function used by Li et al. [15].And make a simple correction to ensure that the free energy at the grain boundary is less than the free energy within the grain.The approximate expression of the matrix phase free energy density function f solid is as follows: , , , , The correction term is , 4 1 5 5 m fix f Φ , assuming that the free energy density of vacancies and gases at the grain boundary of the matrix phase is about nine tenths of that within the grain, this leads to a trend of diffusion of vacancies and fission gases towards the grain boundary.The curvature of the parabola is controlled by m v k and m g k in the expression of the free energy density function, as follows: .8408 × 10 12 J/m³, calculated at T=1276K.Similarly, the approximate expression for the free energy density function f bubble of the bubble phase is as follows: The This correction term is added based on the characteristic of the Helmholtz free energy of the bubble phase as the gas concentration increases.When the gas concentration cg approaches the equilibrium concentration of 0.494, the Helmholtz free energy will sharply increase [11].Assuming that the vacancy concentration cg approaches 1, the free energy also exhibits similar characteristics.Take the curvature of f solid , b g k = 2.25 × 10 9 J/m 3 .The variation of the free energy controlled by this curvature with concentration matches well with the Helmholtz free energy in the range of 0<cg<0.494.Another curvature b v k = 2 × 10 10 J/m³, this is assuming that the vacancy concentration cv in the matrix phase does not exceed the value set by 0.02 when bubbles nucleate.The variation of matrix phase free energy density and bubble phase free energy density with concentration is shown in figure 1.The expression for the free energy density function of polycrystals f poly is as follows: ( , ) 4 2 4 2 Where m is the free energy barrier coefficient, aGB and as are the diffusion interface coefficients, aGB = 1.2, as = 0.8.The expression for the gradient contribution term f poly is as follows: κv, κg, κη, κϕ are the gradient coefficients of vacancy, gas, and two order parameters.m = 3.0 × 10 7 J/m 3 κϕ = 3.38 × 10 -7 J/m 3 , κv = κg = κη = 1.69 × 10 -6 J/m 3

Control equation
The control equation of the phase-field model adopts the Cahn-Hilliard equation [16], [17] to control the conserved field variables cv and cg, and the Allen-Cahn equation [18], [19] to control the order parameters η, ϕ.The control equations for cv and cg are given by: ( , ) ( , ) ( , ) ( , ) Where Mv and Mg are the mobility of vacancies and gases, and their specific form are given by: Where Dv and Dg are the diffusion coefficients of vacancies and gases, assuming Dv = Dg.The fission gas diffusion coefficient is given by Turnbull's model [20][21]:  The control equations for order parameter η and ϕ are given by: F L t (19) ( 1 ) Where Lη and Lϕ is the mobility coefficient, Lη = Lϕ = 1.56 × 10 11 m 3 /(J s) [11].
Scaling the model without dimensionless processing, time transformation coefficient t* = 1 × 10 -5 , τ is the simulation time, corresponding to the actual time t = τ / t*, free energy density transformation coefficient f * = 0.1.The diffusion coefficient used in the simulation is D = D / t*, mobility is L = L / t*, and free energy density is f = f × f *.

Effect of grain boundaries on nucleation
This section discusses the effect of grain boundaries on the nucleation of fission gas bubbles under T = 1276K, 0 < τ < 550s.Construct two different phase-field models, one without grain boundaries (GB) and the other with a straight GB passing through the center.The conclusion is drawn by calculating the nucleation time, radius, and porosity of the simulated area of the bubble.The simulation area of the model is 20μm×10μm, and the model adopts periodic boundary conditions, and the initial concentration of vacancies and fission gases in the entire simulated area is , m eq v c , , m eq g c .Calculated by Eq. ( 18), The diffusion coefficient at T =1200k is D3 = 8.72 × 10 -21 m 2 /s.T = 1276K, 0 < τ < 550s, the evolution of bubble distribution over time is shown in figure 2. The red area in the figure represents grains, while the blue area represents bubbles.At τ =155s as shown in figure 2 (a)-(b), neither model forms bubbles.At this time, vacancies and fission gases are dispersed in the grain.20 seconds later, it can be observed that the model containing GB (figure 2 (d)), where GB absorb defects around them, and bubbles begin to nucleate.The distribution of defects in the grains becomes uneven, indicating the formation of nuclei.The model without GB on the left (figure 2 (c)) does not show this sign, and defects are still evenly dispersed in the grains.At τ =235s, it can be observed that both models have a large amount of nucleation.In the model without GB (figure 2 (e)), bubbles with closer distances fuse with each other to form larger bubbles, while another model (figure 2 (f)) has larger bubbles at the GB.The bubbles inside the grain have not yet fused but larger than those from the model without GB.In addition, obvious "bubble denuded zone" were observed near grain boundaries.When the simulation time τ reaches 500s, bubbles in both models have grown and fused.At this time, the bubbles distribution of the model without GB (figure 2 (g)) is similar to that in the grain of the model with GB (figure 2 (h)).At τ = 235s, the model without GB has an average bubble diameter of 0.08μm, and the model with GB has a much larger intragranular bubbles diameter, which is 0.18μm.At τ = 500s, bubble diameter of the model without GB grows to 0.33μm, and the model with GB has a similar but smaller intragranular bubble diameter at this time, which is 0.31μm.It seems that in early stage, GB will promote the nucleation and growth of intragranular bubbles.After a period of simulation, the bubble sizes of the two models did not differ significantly.However, due to the fact that GB act as defects sink to absorb defects within the grains, the bubble sizes in the model containing GB are slightly smaller than those in the model without GB.Intergranular bubbles of the model with GB always have larger diameter, at τ = 235s, it is 0.25μm, and 0.43μm at τ = 500s.And it is also because GB act as defects sink, visible "void denuded zones" can be observed near GB.Width is defined as the average distance between the grain boundary and all bubbles closest to the grain boundary in this study.At τ = 235s, the width of the zone is 0.65μm, at τ = 500s the width of the zone is reduced to 0.49μm due to the growth of bubbles.We further quantify the impact of GB on gas bubbles by using porosity.In our study, we defined porosity as the percentage of bubble area to the entire simulation area.In statistics, the grid points with η ≥ 0.8 are considered as bubble phase.When T = 1276K, the variation of porosity with simulation time is shown in Figure 3 (a), including models without GB, models with GB, and models with GB but excluding areas with a width of 2μm near GB.As shown in Figure 3, regardless of the model, the evolution of porosity over time can be divided into three stages: the first stage is maintained at the zero level; the second stage is rapid growth stage; the third stage is the stage of slow and stable growth.
Entering the second stage from the first stage signifies the beginning of nucleation, and entering the third stage indicates that no new bubbles are generated and existing bubbles continue to grow and fuse, the process of which is called Ostwald ripening [22].These three stages are similar to those described in ref. [23].
From the figure 3 (a), it can be seen that the fission gas in the model with GB began to nucleate at τ = 165s, earlier than models without GB, which began to nucleate at 215s.Again, it indicates that the presence of GB promotes nucleation within grains.During the period from 180s to 250s, the porosity of the model with GB was always 0.005 higher than that of the model without GB.At τ = 200s, the porosity growth of models with GB begins to slow down, and the porosity of the model without GB maintains rapid growth, the gap in porosity between the two models is gradually narrowing.When the simulation

Effect of temperatures on nucleation
The effect of temperature on the release of fission gas will be discussed, and the influence caused by GB will not be considered this section.Two higher temperature conditions have been chosen: T = 1476K and 1676K.The diffusion coefficients under these two temperature conditions calculated by Eq. ( 19) and ( 20) are D2 = 1.583 × 10 -19 m 2 /s, D1 = 5.476 × 10 -19 m 2 /s, respectively.
The distribution of bubbles is shown in Figure 4.By comparing figure 2 (c) with figure 4 (a)(b), it can be observed that at τ = 175s, the fission gas at 1276K and 1476K has not yet formed bubbles, while visible bubbles have been formed at 1676K, with an average diameter of 0.35μm.At τ =175s, the bubble diameter at 1476K is 0.23μm, and the average diameter of bubbles at 1676K is 0.43μm.At τ =500s, the bubble diameter at 1476K is 0.38μm, and the average diameter of bubbles at 1676K is 0.56μm.The average diameter of bubbles formed under high temperature conditions at the same time is significantly larger than that formed at low temperature.and the number density of bubbles is lower than that at low temperature.Same situation as described in Li's [24] simulation.In addition, under T = 1276K, the nucleation of fission gas occurs in 215s, and 195s under T = 1476K, 150s under T = 1676K.The nucleation of fission gas under high temperature conditions occurs earlier than under low temperature conditions.The variation of porosity with simulation time of the model without GB under T = 1276K, 1476K and 1676K is shown in Figure 5.Under high temperature conditions, the porosity of the model enters a growth state earlier, corresponding to earlier nucleation of bubbles.The increase in temperature leads to an increase in the diffusion coefficient of fission gases and vacancies.It is generally believed that a larger diffusion coefficient leads to the formation of larger bubbles in fission gas.When the simulation time reaches 300s, the difference in porosity of the model of three temperature conditions is small (<0.002).Considering that in our simulation, the source terms of vacancies and fission gases exist in the form of net generation terms, and the net generation term is a constant value.That is to say, the models under three temperature conditions have the same total concentration of fission gas and vacancies at the same time.So it is possible for models under three temperature conditions to have the same porosity at a certain time in our simulation.In addition, this is also related to our statistical method.When calculating the bubble phase, we consider all grid points with η ≥ 0.8 as bubble phases.If the grid points with η >0.7 or 0.9 are considered bubble phases, the statistical results will be significantly different.But no matter what statistical method is used, the porosity of the model under high temperature conditions always enters the growth stage earlier than that under low temperature conditions.

Effect of grain boundaries at different temperatures
In section 3.1, it was found that GB have a promoting effect on the nucleation of fission gas bubbles.In section 3.2, it was found that high temperature conditions also promote the nucleation of bubbles, while promoting the formation of larger diameter bubbles.So, the impact of GB on the release of fission gases under high temperature conditions will be discussed in this section.The variation of porosity with simulation time under T = 1476K and 1676K is shown in figure 3 (b)(c).It can be observed that the trend of porosity over time at 1476K is similar to that at 1276K, where models with GB nucleate earlier and porosity enters the growth state earlier.At T = 1476K, the fission gas in the model with GB began to nucleate at τ = 150s, not only earlier than the model without GB at T=1476K, but it also earlier than the model with GB at T=1276K, conforming to the conclusions in Sections 3.1 and 3.2.At T = 1476K, when the simulation time reaches a certain time, the difference in porosity between the two models is also very small which is similar to the situation at T=1276K.But the situation is different at T=1676K.From Figure 3(c), it can be observed that the time difference between the two models when the porosity enters the growth stage is less than 5s, much smaller than at T=1276K and T=1476K.In the model with GB, nucleation occurs 5 seconds earlier at 1676K than at 1476K, the difference is not significant either.Considering that in this model, the diffusion coefficient of fission gas is the main factor affecting bubble nucleation, while the influence of GB on bubbles is reflected in the total free energy of GB.The total free energy varies very little with temperature, while the diffusion coefficient varies greatly with temperature, The total free energy at GB at 1676K is less than twice that at 1276K, but the diffusion coefficient is nearly 100 times that at 1276K.So it is normal for the variation of porosity over time at 1676K to differ from that at 1276K.The distribution of bubbles is shown in Figure 6.With temperature increasing, the width of the "bubble denuded zone" increases.At τ = 235s, the width of the zone at T=1476K and 1676K are 1.03μm and 1.42μm, respectively.At τ = 500s, the width of the zone at T=1476K and 1676K are 0.89μm and 1.18μm, respectively.Different from Millett's work [23] , no peak zones of bubbles were observed in our simulation.Because in our model, to simplify calculations, vacancies are generated by simulating the net generation term of vacancies.The recombination term and interstitial atom concentration were not considered.So there is no difference in diffusion coefficient between vacancies and interstitial atoms and no peak zones formed.

Conclusions
In this study, the fission gas diffusion coefficient given in the Turnbull model is used for establish a phase field model for the evolution of fission gas bubbles.The model uses parabolic approximation of free energy to simplify calculations and improve convergency.By establishing the models with and without GB, the influence of GB on the formation of bubbles from fission gas and the degree of influence at different temperatures have been studied.
Research has shown that GB have a promoting effect on the nucleation of fission gas bubbles within grains, and the lower the temperature, the stronger the promoting effect.The promoting effect under high temperature conditions is not significant.In addition, the model also demonstrates the influence of temperature on bubble nucleation.The higher the temperature, the larger the diffusion coefficient, and the larger the bubble size and the smaller the number density at the same time.The "bubble denuded zone" also widens with increasing temperature.The model also represents the three stages of the formation of fission gas bubbles: incubation stage, nucleation stage, and growth stage.During the incubation stage, no bubbles are generated and the porosity does not increase; During the nucleation stage, a large number of bubbles form and preferentially form at grain boundaries, resulting in a rapid increase in porosity; During the growth stage, bubbles stop producing, existing bubbles continue to grow and fuse with adjacent bubbles, resulting in a stable increase in porosity.

Figure 1 .
Figure 1.The variation of free energy density with concentration (cv, cg) distribution, the abscissa of the lowest point of the curve is the equilibrium concentration.As cg increases, , b eq v c gradually decreases, and the image of f bubble,v (blue dashed line) will shift to the left.
diffusion of gases under three temperature conditions [29], D1 is the intrinsic diffusion under high temperature (T > 1650K), D2 is the radiation enhanced diffusion under intermediate temperature (1650K > T > 1350K), and D3 is the non-thermal radiation driven diffusion under low temperature (T < 1350K).In our simulation, different diffusion coefficients are used based on different temperatures.The control equation contains two source terms, Pv and Pg respectively represent the rate of vacancy and gas generation.The rate of fission gas Pg generation is calculated by Pg = a V FY nt the rate o FY , F F is the fission rate density and taken as 1.09 × 10 13 fissions/(cm ³ s)[11] based on the operating data of a typical pressurized water reactor.Y is the yield of fission gas produced by each fission, taken as 0.27.Calculated Pg = 1.2037 × 10 10 s -1 , assuming Pv = 20Pg = 2.4074 × 10 11 s -1 .ξv and ξg are simple random functions that represent the thermal fluctuations of vacancies and gases, in order to achieve non-uniform distribution of concentration.The specific form of ξv is ( , )v v v c R t r , Rv is a random function, ranging from -0.01 to 0.01, ξg is similar to ξv..