Development of the mechanistic code MFPR for modelling fission-product release from irradiated UO2 fuel

Development of the mechanistic code MFPR for modelling fission-product release from irradiated UO2 fuel

Nuclear Engineering and Design 236 (2006) 179–200 Development of the mechanistic code MFPR for modelling fission-product release from irradiated UO2 ...

784KB Sizes 2 Downloads 80 Views

Nuclear Engineering and Design 236 (2006) 179–200

Development of the mechanistic code MFPR for modelling fission-product release from irradiated UO2 fuel M.S. Veshchunov a,∗ , V.D. Ozrin a , V.E. Shestak a , V.I. Tarasov a , R. Dubourg b,∗∗ , G. Nicaise b a

Nuclear Safety Institute of the Russian Academy of Sciences (IBRAE), Moscow, Russia b Institut de Radioprotection et de Sˆ uret´e Nucl´eaire, St. Paul lez-Durance, France

Received 25 February 2005; received in revised form 4 August 2005; accepted 4 August 2005

Abstract The models of the mechanistic code MFPR developed in collaboration between IBRAE (Moscow, Russia) and IRSN (Cadarache, France) are described. Exhaustive description of fission-gas behaviour in grain and out of grain is given in relation with individual validation results on analytical experiments under various conditions (steady irradiation, transient, post-irradiation annealing). It is shown that microscopic defects in the UO2 crystal structure can strongly influence fission-gas transport out of grains and release from fuel pellets. These defects include point defects such as vacancies, interstitials and fission atoms, and extended defects such as bubbles, pores and dislocations. A model for the dislocation generation and evolution in irradiated UO2 fuel has been developed and implemented in the mechanistic code MFPR along with a fuel densification model. Being combined with the set of equations describing evolution of point defects (vacancies and interstitials) and their interactions with fission-gas bubbles, a completely self-consistent consideration of the whole system of point and extended defects in irradiated UO2 fuel has been obtained. The mechanistic description of chemically active-elements behaviour is also presented. It is based on complex association of diffusion–vaporisation mechanism involving multi-phase and multi-component thermo-chemical equilibrium at the grain boundary with accurate calculation of fuel oxidation. Examples of application to the VERCORS 4 and 5 experiments show the possibilities of the code in the frame of severe-accident interpretation. © 2005 Elsevier B.V. All rights reserved.

1. Introduction

Corresponding author. Tel.: +7 095 955 2232. Corresponding author. Tel.: +33 4 421 99502. E-mail addresses: [email protected] (M.S. Veshchunov), [email protected] (R. Dubourg). ∗∗

0029-5493/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2005.08.006

Any model that attempts a realistic description of fission-gas release and fuel swelling as a function of fuel-fabrication variables and in a wide range of reactor operating conditions must treat them as coupled phenomena and must include various mechanisms


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

influencing fission-gas behaviour. The most advanced mechanistic codes FASTGRASS (Rest and Zawadzki, 1994) and VICTORIA (Heames et al., 1992) present various detailed models for evolution of intragranular bubbles formed from the solid solution of gas atoms in the UO2 matrix under irradiation conditions, however, they are essentially based on consideration of the equilibrium state of the bubbles. Such an approach radically simplifies the theory, since in this case the defect structure of the crystal (including point defects, such as vacancies and interstitials, and extended defects, such as dislocations) is almost completely excluded from consideration. Such simplification is well grounded only in steady-state irradiation conditions (with except of high burnups) which can be relatively well described by a simple parametric Booth model or other semiempirical correlation models. So, in these simple situations, the advantage of the mechanistic approach is not so visible. Under transient and/or annealing conditions the approximation of equilibrium bubbles is not anymore valid, and interactions of bubbles with point defects and dislocations become essential. Lacking mechanistic description of real defects (vacancies, interstitials and dislocations) and their interactions with bubbles, various artificial mechanisms were introduced in previous codes, in order to simulate more complicated regimes. These artificial mechanisms require “effective” parameters, which introduce strong uncertainties in the code predictions, and, as a result, the advantage of the mechanistic approach is partly lost. Trying to avoid any artificial tuning and introduction of artificial mechanisms were important reasons for development of the code module for fission-products release (MFPR). This code self-consistently describes evolution of various defects and their interactions with gas atoms and bubbles migrating out of grains, and new parameters characterizing the crystal defect structure also arise. However, being physically grounded, these microscopic parameters can be fixed from the analysis of available experimental data, and then used without any artificial tuning in further calculations. This is the main goal of this mechanistic approach performed in collaboration between the Institute for Nuclear Safety of Russian Academy of Sciences (IBRAE) and the French “Institut de Radioprotection et de Sˆuret´e Nucl´eaire” (IRSN), where the code is used to better understand the results of the experiments performed in

the context of radiological consequences of accidents. In particular, this approach will help in developing appropriate models for evaluating the risk of fissionproduct release in the event of a severe accident or design basis accident in a LWR. In addition, as shown by many observations and interpretations, the release of various chemically active fission products (FP) strongly depends on formation of separate phases, mainly at grain boundaries, and thermal behaviour of these phases in various gas atmospheres that determine the fuel oxygen potential and thus the stability of the different phases. The MFPR module includes modelling of chemistry effects on the behaviour of FPs within irradiated oxide fuel at high temperatures in the interval 500–3000 K. The ‘U–FP–O’ system is considered as a multi-phase system consisting of multi-component phases. Within the fuel matrix, the fission products migrate as atoms to the gas bubbles and solid phase precipitates at grain faces, and their mobility depends in particular on the extent of fuel oxidation characterized by the stoichiometry deviation. The release rate is thus proportional to partial pressures of FP vapours determined by chemical states of condensed FP species, and the major mechanism for release of the fission products is their vaporization at the gas/solid interface. In connection with that description, a model for fuel oxidation in steam/air mixtures was also developed. This paper will briefly describe, for an intact fuel geometry, the main models of the code both for fission-gases and chemically active elements and will give results of validation of individual models against separate-effects experiments and more global validation against the French VERCORS experiments (Ducros et al., 2001) performed for IRSN.

2. Physical processes and models Due to fission processes, a wide spectrum of fissionproduct elements is created in the UO2 matrix, the most important being combined into 15 element groups: Cs, Ce, I, Eu, Mo, Nd, Ru, Nb, Ba, Sb, Sr, Te, Zr, Xe, La. It is assumed that FP elements exist in the UO2 matrix in atomic form. All impurity atoms formed due to fission processes migrate to grain boundaries. In this way, fission-gas atoms can be captured by intragranular gas bubbles, which can also migrate to grain

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

boundaries. Part of the captured atoms can escape from bubbles by irradiation-induced and thermal re-solution processes. Chemical interactions between the FP elements and the dissolved oxygen result in formation of separate phases in solid precipitates on the grain boundaries and vaporization to the intergranular gas bubbles. These processes affect significantly the release rates of all FP elements including noble gases. Intergranular bubbles are represented by three groups: bubbles on grain faces (GF), on grain edges (GE) and in grain corners (GC). The GF bubble growth progresses up to the grain surface saturation, when interlinkage of the GF bubble and formation of grain face channels to the grain edges and corners occurs. Growth of the GE and GC bubbles leads to their interconnection by tunnels and formation of an open porosity structure. In the MFPR model, the process of FP release from solid fuel is divided into two consecutive stages: • intragranular FP transport from the bulk to the grain boundary accompanied by formation of solid precipitates and gaseous species in the intergranular bubbles, • accumulation of FP-bearing gases in the intergranular bubbles, and release to the open porosity through the system of bubbles on the grain boundaries and the network of channels and tunnels. The model currently ignores the transport delay in the FP releases into the gap related to the gas transfer through the channel network and through the open porosity. 2.1. Intragranular transport of FP elements In the transport problem a fuel grain is considered for simplicity as an isotropic sphere. The problem is formulated separately (and in different manners) for two subsystems: the noble gas (Xe, Kr) atoms dissolved in the matrix and filling the grain bubbles, and atoms of chemically active element and molecules of compounds formed by these elements. 2.1.1. Diffusion and release of chemically active FP elements In the MFPR chemistry model, irradiated oxide fuel including fission products and dissolved oxygen


is considered as a multi-phase system consisting of multi-component phases. According to experimental data (Kleykamp, 1985) and thermo-chemical calculations (Imoto, 1986; Cordfunke and Konings, 1988, 1993; Ball et al., 1989; Moriyama and Furuya, 1997), the major considered phases in the irradiated UO2 are the “fuel-FP” solid solution, the metal phase, the oxide phase of complex ternary compounds, the solid phase of CsI and the gas phase. Detailed description of the phase composition is given in Appendix 1. Zirconium, rare earth and alkaline earth elements are characterized by a high affinity for oxygen. Within the UO2 matrix, these elements exist mainly in the form of oxides. Zirconium and rare earth elements are partially or completely miscible with uranium dioxide to form a solid solution. The relatively low solubility of barium results in the majority of barium in the form of uranates, zirconates and molybdates precipitating as the multicomponent ‘grey’ phase, which also includes ternary compounds of strontium and caesium. Solubility of the metals Mo and Ru and their oxides in solid UO2 is extremely low. These metals (along with Tc, Rh, Pd) are present in the metallic (‘white’) inclusions of the oxide fuel. Within the fuel matrix the fission products migrate as atoms (Grimmes and Catlow, 1991) to the grain boundary (where they might be in gas bubbles or condensed precipitates forms), and their mobility depends in particular on the extent of fuel oxidation. In the current MFPR version it is assumed that separate solid phases (precipitates) are formed only on grain faces at the interface with gas face bubbles; solid precipitates are in equilibrium with the gas phase and with the solid solution. Partitioning of elements between chemical compounds and phases and the fuel stoichiometry are calculated self-consistently as they are determined by a balance between a series of oxidation/reduction reactions and solid-state diffusion. To formulate the transport problem, the complex system is separated into two subsystems. That is (1) solid solution (SS) of FP elements, their oxides and atomic oxygen dissolved in the UO2 matrix and (2) the subsystem ’solid precipitates–gas phase’ (SP/G) including the metallic “grey” phases and the phase of condensed CsI(c). Concentration profiles of FP elements in SS are described by the diffusion equations:   ∂ (1) −2 ∂ 2 ∂ (1) , (1) Y = Bi + r Di r Y ∂t i ∂r ∂r i

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

182 (1)

where Yi is the local volume concentration of the ith element in the subsystem 1—solid solution, Bi the rate of generation of the i-th element due to fission, and Di is the diffusion coefficient. Boundary conditions for these equations are determined by the thermo-chemical equilibrium within the SP/G subsystem and at the subsystem interface, which depends on fluxes of elements out of SS. The equilibrium composition of the phases is treated in terms of semi-ideal chemistry model in which phenomenological solid solubilities of FP elements are used, and the chemical potential of dissolved oxygen is described by the Lindemer–Besmann correlation (Lindemer and Besmann, 1985). More detailed description of the thermo-chemical approach used in MFPR is presented in Appendix 1. Solution of the above discussed complex problem yields partitioning of elements between various chemical states within the SP/G subsystem and between the two subsystems. The subsequent migration of the FP elements in the form of gaseous compounds to open porosity is related to the intergranular bubble growth kinetics. 2.1.2. Transport of fission gas The present version of the MFPR code includes two macro-models for intragranular transport of the fission product Xe: bi- and multi-modal models. The most detailed approach to describe the intragranular Xe gasbubble system is the multi-modal model, formulated in terms of the distribution function for the bubble sizes. In the bi-modal approach similar to that of the VICTORIA and FASTGRASS codes, the effects related to a finite width of the bubble size distribution are ignored. The basic space-time dependent variables are the volume concentrations of gas atoms Cg and bubbles Cb , the average number of gas atoms within a bubble Nb , and the average bubble volume Vb . It is convenient also to use dimensionless concentrations (number per uranium atom), which will be denoted by a small letter, e.g.: cg = Cg Ω,


where Ω is the atomic volume in UO2 . Additionally, MFPR includes self-consistent consideration of point defects (vacancies and interstitials) and dislocations, which mutually interact with each other and with gas bubbles and as-fabricated pores

during their evolution under irradiation or annealing conditions. These defects are described by three additional variables: vacancy concentration Cv , interstitial concentration Ci , dislocation density (length per unit volume) ρd , and concentration of atoms captured by dislocations Yd . Details of such an approach will be briefly presented in Sections 2.2–2.4. Transport equations for Cg and concentration of atoms-in-bubbles Yb = Nb Cb can be written as ∂Cg = Dg Cg − Fg→b − Fg→d + κG, ∂t


∂Yb = ∇(Db ∇Yb ) − ωbmg Yb + Fg→b − Fb→d ∂t


where Fg→b = kgg Fn Cg2 + kgb (Cg − Cgeq )Cb − bGYb ,


kxy = 4π(Dx + Dy )(Rx + Ry ),


x, y = g, b.

Different terms in the right-hand side of these equations correspond to the following physical effects, each being characterized by a particular kinetic parameter: • diffusion (Dg and Db are the diffusivities of gas atoms and bubbles, respectively); • gas atom generation by fission (G is the fission rate and κ = 0.241 is the probability of Xe formation by fission); • bubble biased migration (ωbmg ) along the temperature gradient; • bubble nucleation (with the nucleation rate Fn ); • gas atom capture by bubbles (Rg is the gas-atom radius and Rb is the bubble radius in the capture kernel kgb ); eq • thermal re-solution (Cg is the equilibrium gas concentration); • radiation induced re-solution (b is the re-solution probability factor). Modelling of these effects as well as sweeping of gas and bubbles by climbing dislocations (Fg→d and Fb→d ) is described in Sections 2.2–2.4 along with equations for the other MFPR variables (ρd , cv and ci ). Superposition of these effects (along with the van der Waals equation for the gas state) determines, in particular, the intragranular bubbles concentration Cb and mean radius Rb , as well as gas release to the grain faces.

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

In growing grains, the boundary conditions for Eqs. (4) and (5) are set at moving grain boundaries:   ∂Cg  = 0, Cg r=R (t) = Cδ ,  gr ∂r r=0  ∂Yb  = 0, Yb |r=Rgr (t) = 0. (7) ∂r r=0 Additionally, sweeping of intragranular bubbles by moving grain boundaries takes place. The MFPR model for the grain growth with the retarding effect of attached bubbles is presented in (Veshchunov and Tarasov, 2002).

Evolution of the vacancy and interstitial fields in grains is described in the mean field approximation (Brailsford and Bullough, 1981) in terms of dimensionless concentrations, cv and ci (number of vacancies and interstitials per uranium atom), by the following equations: c˙ v =

2 + kvgb )Dv cv

− αDi ci cv + Ke

+ Kb + Kp + (1 − ξ) K,

2 c˙ i = −(ki2 + kigb )Di ci − αDi ci cv − Kd + K,

(Brailsford and Bullough, 1981): 2 2 = ki,v kigb,vgb

ki,v Rgr coth(ki,v Rgr ) − 1 , 2 R2 /3] [1 − ki,v Rgr coth(ki,v Rgr ) + ki,v gr (10)

where Rgr is the grain radius. In accordance with (Brailsford and Bullough, 1981; Bullough et al., 1975; Dollins and Nichols, 1978), the total sink strength of vacancies and interstitials due to the extended defects are calculated as 2 ki,v = 4πRb Cb + 4πRp Cp

+ Zi,v (ρd + 2πRl Cl + 2πRvl Cvl ).


The rate of the vacancy thermal production has the form:

2.2. Evolution of point defects in irradiated UO2 fuel





where Dv and Di are the vacancy and interstitial diffusion coefficients, respectively, α is the recombination constant, kv2 and ki2 are the total sink strength of vacancies and interstitials into the extended defects (gas bubbles, pores, vacancy clusters and dislocations), 2 and k 2 are the grain-boundary sink respectively, kvgb igb strength for vacancies and interstitials, respectively, K is the Frenkel pair production rate (d.p.a. s−1 ) proportional to the fission rate G; ␰K is the rate at which vacancies are removed from solution to form vacancy loops, Ke is the rate of the vacancy thermal production, Kp is the rate of the vacancy irradiation-induced re-solution (knockout) from pores, Kb is the rate of the vacancy irradiation induced resolution from gas bubbles, and Kd is the rate of the interstitial absorption due to the interstitial loop generation. The grain-boundary sink strength for vacancies and interstitials is estimated from the relationship

Ke = Dv [4πRb Cb cvbs + 4πRp Cp cvps + Zv (ρd + 2πRl Cl + 2πRvl Cvl )cveq ],


eq vacancy concenwhere cv is the thermal-equilibrium ps tration, cvbs and cv are the boundary concentrations of vacancies at the bubble and pore surfaces, Cp and Rp are

the pore concentration and mean radius, respectively, Cvl and Rvl are the vacancy loops (clusters) concentration and mean radius, respectively, Cl and Rl are the interstitial-loop concentration and mean radius, respectively (see Section 2.3). The vacancy production rate from pores due to knockout by fission fragments passing through pores was described in (Stehle and Assmann, 1974) in the form: Kp = 8πR2p cp ληG,


where the number of vacancies knocked out of a pore per collision η = 100 and the length of the fission fragment path λ = 10−6 m, in accordance with estimations of (Dollins and Nichols, 1978). The vacancy production rate from bubbles due to knockout by fission fragments passing through bubbles, which was described in (Veshchunov, 2000) as the thermal annealing of bubbles in fission tracks, can be represented in the form: Kb = 4πRb Cb [Dv (cv − cveq ) − Di ci ] −

d(Cb Vb ) , dt (14)

i.e., as a difference between the rate balance of point defects sinking into bubbles and the real rate of bubbles


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

volume growth. In this case, the bubbles are considered at equilibrium (owing to the thermal annealing in fission tracks) and their growth is determined by absorption of gas atoms (as the rate-controlling step). The rate of the interstitial absorption due to the interstitial loop generation (in the bi-modal approximation for loops) is given by the equation: dCl Kd = πR2l B , dt


where B is the dislocation Burgers-vector length. The sink strength for interstitials is larger than that for vacancies due to the higher elastic interaction between dislocations and interstitials (Brailsford and Bullough, 1981): Zi = Zv (1 + 2ε),

0 < 2ε  1.


The effective uranium self-diffusion coefficient is evaluated as the sum of two (a-thermal and thermal) terms (Matzke, 1986):   EU (0) , (17) DU = AF + DU exp − T with A ≈ 1.2 × 10−39 m5 , DU = 2 × 10−4 m2 /s and EU = 64 200 K, T the temperature and F is the fission rate. On the other hand, DU can be represented as superposition of the vacancy and interstitial mechanisms with effective diffusivities Dv and Di : (0)

DU ≈ Dv cv + Di ci ,


which under steady irradiation conditions obey a relationship (Veshchunov, 2000):   kv2 ≈ D v cv , (19) Di ci ≈ Dv cv ki2 and therefore, DU ≈ 2Dv cv . 2.2.1. Effect of point defects on fission-gas bubbles The boundary concentrations of vacancies at bubble and pore surfaces are given by   ΩδPb bs eq cv = cv exp − , kT   ΩδPp cvps = cveq exp − (20) kT

with δPx = Px − Ph − Px =

2γ , Rx

Nx kT , (Vx − BXe Nx )

x = b, p,

where Ph is the external hydrostatic pressure, γ the effective surface tension of UO2 , Pb,p the gas pressure within a bubble or pore given by the van der Waals equation of state and BXe is the van der Waals constant for xenon. In the case of steady-state irradiation conditions the bubble volume, Vb , can be calculated using the equilibrium equation, δPb = 0, owing to thermal annealing of small intragranular bubbles in fission tracks (see Eq. (14)), and the corresponding boundary condition takes the form: cvbs = cveq . The nucleation factor, Fn , entering transport equations for the intragranular gas–bubble system, Eqs. (3) and (4), is usually interpreted as the parameter determining the probability that two gas atoms that have come together actually stick and form a bubble (Rest and Zawadzki, 1994). This parameter was introduced as a default value varying in a wide range from 10−4 to 10−7 . It was recently shown (Veshchunov, 2000) that in a more adequate interpretation, the factor Fn is proportional to the probability that a vacancy is located in a certain position (collision of two atoms), and therefore is equal to the vacancy concentration cv . Correspondingly, in the current version of the MFPR code with the vacancy field model activated, the nucleation factor is calculated in a self-consistence manner: Fn = cv .


2.3. Evolution of extended defects In the present model, generation of dislocation loops occurs mainly during the initial period of irradiation, concurrently with the fuel densification process, and therefore, may be significantly influenced by the kinetics of fuel densification (pore sintering). For this reason, comparative analysis of existing models for pore evolution under irradiation conditions was performed and the most adequate model was chosen for further

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

development and implementation in MFPR (Shestak et al., 2004). Consequently, in the new MFPR code version both processes of dislocation-loop generation and pore shrinkage are considered simultaneously and self-consistently with the point defects and gas-bubble evolution. 2.3.1. Pores evolution As demonstrated in (Dollins and Nichols, 1978), gas content Np and pressure Pp in large pores are essentially zero and, hence, are unimportant when calculating the change in the pore volume with respect to time: dVp = 4πRp [Dv (cv − cvps ) − Di ci ] − 8πR2p ληΩG, dt (22) where 



Ω exp kT

2γ + Ph Rp


2.3.2. Vacancy-cluster (loops) evolution A fraction ξ of vacancies created by fission are trapped in vacancy clusters (loops). The time variation of the vacancy cluster concentration is found from the equation (Dollins and Nichols, 1978): d 2πRvl Cvl = − [Zi Di ci − Zv Dv (cv − cveq )]Cvl dt Ωvl +

ξK , Ωvl

where Ωvl ≡

πR2vl B

(23) is the vacancy cluster volume.

In the bi-modal approximation the loop concentration obeys the following equation: αl Di ci2 dCl = , (24) dt Ω where αl is the dislocation-loop nucleation constant. It is assumed (Yoo, 1977) that dislocation loops are generated during the initial stage of irradiation until: Di ci − Dv (cv − cveq ) > 0.


Calculations of the loop size distribution (Hayns, 1975) confirm that the loop population achieves a steady-state concentration during some initial nucleation period. This conclusion is also in a fair agreement with direct observations in irradiated UO2 (Whapham and Sheldon, 1965; Whapham, 1966). The growth rate of the mean radius of dislocation loops obeys the following equations (Yoo, 1977): dRl 2π = [(1 + 2ε)Di ci − Dv (cv − cveq )] dt B  1   , if Rl < Ll ;  ln(8Rl /Rd ) × 1    , if Rl ≥ Ll , ln(Ld /Rd )


where Rd is the dislocation core radius estimated as Rd ∼ = 3B, the dislocation bias factor ε was introduced in Eq. (16), and Ld = (πρd )−1/2 ,

Ll = (3/4πCl )−1/3 ,

characterize the mean distance between dislocations and dislocation loops, respectively. The dislocation sink strength for vacancies is evaluated as (Yoo, 1977): Zv =

2π 4π2 Rl Cl + , ln(Ld /Rd ) ρd ln(8Rl /Rd )

ρd (t) = ρd (0) = Const, 2.3.3. Dislocation-loop evolution According to the present model, dislocations are generated under irradiation in the form of di-interstitial clusters (Hayns, 1975) and continuously grow by absorption/emission of point defects in two types, dislocation loops and dislocation network. A simple approach is developed by considering a bi-modal distribution of loops, their evolution and transformation into a dislocation network.


Zv =

if Rl < Ll ;

2π , ln(Ld /Rd )

ρd (t) = ρd (0) + 2πRl Cl ,

if Rl ≥ Ll ,


where ρd (0) is the initial network dislocation density. In this approach, it is assumed that when dislocation loops are large compared to the sphere of influence, i.e., Rl ≥ Ll , the dislocation loops may be better approximated as straight dislocations of the length equal to


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

irradiation effects at temperatures below ≈1500 ◦ C and thermal effects at temperatures above ≈1500 ◦ C. In order to improve the microscopic description of the fission-gas behaviour, new models were developed and implemented in the mechanistic code MFPR.

Fig. 1. Dislocation density as a function of burn-up. Comparison of MFPR predictions with measurements (Nogita and Une, 1994).

their circumferences by using the relationship that the corresponding dislocation line density is 2πRl Cl . In this approximation the growth rates of dislocation loops are given by the second part of Eq. (26) corresponding to the case Rl > Ll . A thorough validation of the defect model against various steady-irradiation experiments allows specification of the main microscopic parameters ε, ξ, zs , αl , etc. (Shestak et al., 2004), which can essentially differ from the corresponding values typical for metals, despite the basic physical mechanisms of the defect structure evolution being similar in different materials. Therefore, a similar formalism of the defect model (either for point defects (Brailsford and Bullough, 1981; Bullough et al., 1975), or for dislocations (Hayns, 1975; Yoo, 1977), originally developed for irradiated metals, with a new set of the microscopic parameters allows a mechanistic description of the defect system evolution in the ceramic material. In particular, results for the dislocation density evolution in irradiated UO2 obtained by the MFPR code are in a satisfactory agreement with measurements of (Nogita and Une, 1994) (see Fig. 1). In addition, these results provide initial conditions for simulation of annealing regimes, which are rather sensitive to the dislocation density attained during the pre-irradiation period of the annealing tests (see Section 2.4.4).

2.4.1. Intragranular bubbles under steady irradiation conditions at low temperatures ≤1500 ◦ C (irradiation effects) Theoretical analysis of intragranular bubblebehaviour shows that irradiation effects can produce strong limitations on the maximum bubble concentration attained under steady irradiation conditions at temperatures below 1500 ◦ C (Veshchunov, 2000). Owing to increase of the total sink of vacancies into the extended defects (bubbles and dislocations) in the course of their evolution, a strong reduction of the vacancy concentration at a late stage of irradiation takes place (in accordance with Eq. (8)). In its turn, according to Eq. (21), the bubble nucleation rate is directly proportional to the vacancy concentration and thus is strongly suppressed, leading to “saturation” of the bubble concentration (Veshchunov et al., 2000). After implementation of the full set of microscopic equations for evolution of point and extended defects and their interactions with gas bubbles (described in Sections 2.2 and 2.3), MFPR allows a satisfactory prediction of the intragranular bubble concentration saturation at a late stage of irradiation observed in recent tests with high burn-up fuel (Kashibe and Une, 1993) (Fig. 2).

2.4. Evolution of intragranular gas bubbles As demonstrated in (Veshchunov, 2000), currently existing models and codes generally underestimate

Fig. 2. Intragranular bubble diameter and concentration as a function of burn-up. Comparison of MFPR predictions with measurements (Kashibe and Une, 1993).

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

2.4.2. Intragranular bubbles under steady irradiation conditions at high temperatures ≥1500 ◦ C (thermal effects) The kinetic equation for the number Nb of gas atoms in a bubble with the radius Rb under irradiation conditions with consideration of the thermal re-solution of gas atoms from a bubble in the van der Waals approximation, takes the form: dNb = 4πDg Rb [cg − cgeq ], dt eq

(28) eq

where cg = 4πDg Rb [cg − cg ], and ϕ(P,T) = exp (BXe P/kT) is the function accounting for the gas-phase non-ideality (typical for small bubbles with Rb ≤ 5 nm) in the van der Waals approximation, P the pressure within the bubble, T the temperature, BXe the van der Waals parameter for Xe, and KS (T) is the Xe-atom solubility in UO2 .


The temperature dependence of KS is represented in the Arrhenius form:   ES ES − , (29) KS (T ) = KS (T0 ) exp kT0 kT where KS (T0 ) is the solubility at a fixed temperature T0 and ES is the activation energy of the solid solution. Numerical values of these two parameters were estimated in (Veshchunov et al., 2000) as KS (T0 ) ≈ 2 × 107 J−1 and ES = 3 eV on the base of available experimental data. Results of the MFPR code calculation (Veshchunov et al., 2000) confirmed that at high temperatures, the thermal re-solution process determines the onset of bubble formation under steady irradiation conditions. Namely, the thermal re-solution model predicts a complete suppression of the bubble nucleation process at T > 2000 ◦ C, in a qualitative agreement with observations of (Baker, 1977); this is shown by the results of MFPR calculations in Fig. 3a. Also, there is a signifi-

Fig. 3. Effect of thermal resolution on intragranular-bubble evolution (solid lines: model on, dashed lines: model off). Results of MFPR simulations of high-temperature tests (Baker, 1977) in (a), and low-temperature tests (Turnbull and Cornell, 1970) in (b).

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200


cant delay in the onset of the bubble formation at low temperatures T < 800 ◦ C, as observed in tests (Turnbull and Cornell, 1970) and reproduced by calculation in Fig. 3b. 2.4.3. Intragranular bubbles under transient conditions According to the model for intragranular bubbles (Nelson, 1969), irradiation-induced resolution of gas atoms occurs from a thin periphery layer (λ ≈ 1–1.5 nm) of bubbles, therefore, the resolution rate becomes inversely proportional to the bubble radius, Jres = bNb , where b ≈ b0 λ/(λ + Rb ), is the resolution probability, and Nb is the number of atoms in a bubble. However, the ejection of a gas atom into the surrounding matrix does not automatically result in its resolution, since those atoms would tend to return back to the bubble. In order to take this tendency into quantitative consideration, one can assume that the influx of the ejected atoms back to the bubble proceeds by diffusion within the re-solution layer δ ∼ 1 nm. Analytical consideration of this process results in the following formulation of the bubble growth equation (Veshchunov et al., 2000): dNb = 4πDcg Rb − b Nb , dt


with the renormalized expression for the resolution probability: b ≈ b0

λ δ . λ + Rb δ + Rb

Fig. 4. Size distribution function for intragranular bubbles calculated by MFPR with standard model for gas re-solution from bubbles.

2.4.4. Intragranular bubbles under annealing conditions Microstructure observations of irradiated fuel in the annealing tests confirm that gas release is accompanied by a rapid growth of the intragranular bubbles (up to hundreds of nm at high annealing temperatures) and a noticeable decrease (by several orders of magnitude) of the bubble concentration. This process of the bubble number decrease is usually associated with the Brownian motion of the bubbles leading to their coalescence


The improved model for gas re-solution from bubbles was implemented in the MFPR code. Numerical simulations of the bubble evolution behaviour under transient conditions of the tests (Ray et al., 1992) reveal a qualitative improvement of the calculation results in comparison with the standard approach (i.e. b ≈ b0 λ/(λ + Rb )), presented in Fig. 4. Indeed, the new model predicts formation of the “trimodal” bubble distribution function (i.e. single atoms and two populations of bubbles) with a much wider distribution of bubble sizes, Fig. 5, in accordance with the observations (Ray et al., 1992). The physical reasons for such a remarkable change in the bubble size distribution is connected with the onset of a new critical point in the nodal line Poincar´e formalism, as analyzed in (Veshchunov et al., 2000).

Fig. 5. Size distribution function for intragranular bubbles calculated by MFPR with advanced model for gas re-solution from bubbles.

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

into larger ones in the grain bulk and transport to the grain boundaries. During bubble growth and coalescence, extended defects such as dislocation loops uniformly distributed in the grain bulk, may act as the main sources of vacancies (necessary for the bubble equilibration) and afford the equilibrium concentration of the point defects in the crystal bulk. This explains dislocation creep and enhanced bubble growth by dislocation sweeping under annealing conditions, observed in (Whapham and Sheldon, 1965). In this situation a new mechanism for gas release due to dislocation creep emerges (Ozrin et al., 2002). This mechanism considers sweeping of bubbles and delivery them to the grain boundaries by climbing dislocation segments in the course of vacancy generation (necessary for equilibration of growing bubbles):   3ρd vd t dρd Fb→d 3vd ρd creep Fd→f = , dt , =− 2dgr t0 ρd dt 2dgr (32)

Fb→f = 2(Rcore + 2Rb )ρd vd Yb


The terms in the flux Fb → d corresponding to bubble sweeping by climbing dislocations includes the dislocation climb velocity, vd , which is determined by the amount of vacancies ejected by dislocations per unit time for equilibrating of intragranular gas bubbles:   cv (d) (d) vd = v0d ϕd (Rb , Nb ) 1 − eq , (34) cv


dislocations by bubbles has the form: (d)


(d) Nb


ρd . 2Yd


After an initial period determined by Eq. (37), a strong pinning of dislocations by swept bubbles saturates this source of point defects, and grain boundaries apparently become the dominant source of vacancies during the subsequent period of the annealing tests. In this situation a vacancy flux directed from the grain surface to its interior arises that induces bubble biased migration along the vacancy gradient in the opposite direction, as proposed recently by (Evans, 1997) in a qualitative approach. In order to handle this problem quantitatively, a simple analytical model was developed (Veshchunov et al., 2000). Neglecting bulk sources of point defects (i.e., dislocations), the model allowed it to be related quantitatively to the kinetics of bubble coalescence and swelling, in agreement with measurements (Baker and Kileen, 1987) (Fig. 6). Basing on this analysis, the MFPR transport equation for the gas-in-bubbles concentration Yb , Eq. (4), is modified by adding the term considering bubble biased migration with velocity vvac along the gradient of b the vacancy concentration cv (Geguzin and Krivoglaz, 1971): vvac b (r) = 2Dv

∂cv . ∂r


where the unperturbed velocity is given by v0d ≡

2πDu , B ln(Ld /Rd )


and the pinning by factor, ϕd , is determined by attached bubbles (Ozrin et al., 2002): (d)


ϕd (Rb , Nb , Yd , ρd ) =

1 (0) (d) (d) 1 + (B2 vd Yd /ρd Nb Db )

(36) (d)


where Rb and Nb are, respectively, the bubble radius and the mean number of atoms in bubbles captured by dislocations. The condition of complete coverage of

Fig. 6. Model prediction of gas release dependence on the relative mean bubble radius (Rb /R0 ) during annealing stage, in comparison with experimental observations (Baker and Kileen, 1987) at 1600 ◦ C at a variety of moments (50 s, 1.4 h, 24 h).


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

Fig. 7. (a) Dislocation velocity in a high temperature annealing regime (VERCORS 4) calculated by MFPR with different bubble diffusivities. (d) (b) Mean radii of intragranular bubbles Rb and bubbles on dislocations Rb in a high temperature annealing regime (VERCORS 4) calculated by MFPR with different bubble diffusivities.

Simulation of the annealing tests by MFPR shows that the dislocation sources acting in the initial period of the annealing stage significantly suppress gas release predicted by the Evans mechanism, and becomes dominant during the initial period. As a result, the dislocation mechanism provides so-called “burst release” observed in the initial stage of annealing, whereas the vacancy mechanism (considered by Evans) combined with the advanced model for bubble diffusivity (proposed by (Mikhlin, 1979)) provides a more gradual release in a late period of the annealing stage (when dislocations become immobile). This is illustrated by MFPR calculations of dislocation velocity and bubble size evolution during the high temperature annealing regime of the VERCORS 4 test (Ducros et al., 2001), represented in Fig. 7a and b for various models of bubble diffusivity (the standard and Mikhlin’s ones). Superposition of the mechanisms and models above described, allows an adequate description of various annealing tests (Zacharie et al., 1998), (Une and Kashibe, 1990), as shown in Figs. 8 and 9 (Ozrin et al., 2002). Calculations enable a good agreement with experimental data for the fractional gas release to be achieved, Fig. 8, and simultaneously for the value of intragranular bubble size, Fig. 9. Indeed, at the end of the annealing stage the measured mean bubble diameter was ≈55 nm (Une and Kashibe, 1990) whereas the calculated diameter is equal to 57 nm. A good agreement with the measured bubble concentration is also attained in these calculations. At high annealing temperatures (≥2000 ◦ C), after dislocation pinning, MFPR predicts an essential

decrease of the vacancy concentration in the central region of the grain. This leads to a strong overpressurisation of bubbles so that the thermal re-solution effect becomes significant (see Section 2.4.2). In particular, it reduces sinking of gas atoms into the bubbles and results in additional diffusion flux of gas atoms out of the grain. This high-temperature effect is illustrated by calculation of the annealing regime in the VI-3 test (Osborne and Lorenz, 1992) (Fig. 10). 2.5. Intergranular FP transport, release and swelling models Intergranular FP transport is described in terms of volume concentrations of the gas components con-

Fig. 8. Experimental data (Zacharie et al., 1998) and MFPR calculation results for Xe release as a function of annealing time at different temperatures.

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200


close to that developed by White and Tucker (1983). In the case of the grain face bubbles, the rate equation takes the form: d (i) (i) Y = Ff θf , dt f

θf ≡ θ(A∗ − Af ),


where Af = π(Rf sin θf )2 Xf ,

Fig. 9. MFPR calculation results for Xe atom and bubble concentrations and bubble radius as functions of annealing time at 1800 ◦ C. (i)

tained in the GF, GE and GC bubbles denoted by Yf , (i) (i) Ye and Yc , respectively. Total volume concentration of gas molecules composing bubbles is defined as a sum over all gas components:

Yx = Yx(i) , x = f, e, c. i

2.5.1. Model of FP release by bubble interlinkage The equations governing the FP transport from the system of the GF and GE bubbles to the open porosity are formulated based on a percolation type approach


is the projected area coverage of the grain face by the face bubbles, Rf is the radius of these bubbles, Xf is their surface concentration on the grain face (number of bubbles per unit surface area), θ(x) the Heaviside (i) step function, and Ff is the gas flux to the face bubbles (calculated in Section 2.5.2). The projected area coverage along with the face-bubble concentration and radius depends on time due to the coalescence effects (considered in Section 2.5.3). The critical-area coverage of the grain face: A∗ = 0.5,


determines the onset of formation of the grain-face interconnected channels to the grain edges. The rate equations for the concentration of the gas component in the GE and GC bubbles have the form: d (i) Y = Fe(i) θe , dt e (i)

d (i) Y = Fc(i) θe , dt c



where Fe and Fc are the gas fluxes to the grain edge and grain-corner bubbles, respectively; θ e ≡ θ(L* − le Nbpe ); L* is given by L* = Ledge − 2Rc ; Ledge the average edge length, Re and Rc are the radii of grain edge and grain-corner bubbles, respectively, le the effective length of edge bubbles and Nbpe is the number of edge bubbles per grain edge. Saturation of the grain edge and grain corner porosity takes place when these bubbles are just touching each other, that is Nbpe le = Ledge − 2Rc (White and Tucker, 1983)). This is the critical point for formation of the long-range interconnection of the grain-edge and grain-corner bubbles to the open porosity. Finally, the release rate of the i-th gas component is determined by Fig. 10. MFPR simulation of the VI-3 test: joint effect of the vacancy and thermal re-solution models after subsequent implementation of the models in the base version of the code.

d (i) (i) Y = [Fe(i) + Fc(i) + Ff (1 − θf )](1 − θe ). dt out


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200


2.5.2. Model for Xe grain-face diffusion transport The MFPR advanced model for the transport of Xe atoms self-consistently takes into account the effects of atom diffusion along the grain surface and irradiation-induced re-solution from intergranular bubbles (Veshchunov et al., 2000; Berdyshev and Veshchunov, 2002). It is shown that “circulation” of gas atoms collected by growing intergranular bubbles from the grain face and then returned back (by the re-solution process) into the grain matrix, makes intergranular bubbles much less effective sinks for gas atoms, since it decreases their growth (i.e. approaching a balance among absorbed and re-emitted atoms) and thus continuously increases a fraction of the source term flux (i.e. diffusion flux from grain to GF) eventually transported to grain edges. In particular, this allows prediction of a noticeable gas release from fuel when the grain-face coverage is far below the saturation value ≈50% determining onset of grain-face-bubble interlinking (see Eq. (41)). The surface-diffusion model is formulated in terms of the fluxes to the GF, GE and GC bubbles: f Ff = F␦ − Ff→ec − ωrsi Yf , e Fe = ηe Ff→ec − ωrsi Ye ,

c Fc = ηc Ff→ec − ωrsi Yc .

(44) where relative parts ηe and ηc of the total flux out of grain-face to edge and corner bubbles, respectively, for simplicity are assumed to be equal to the relative numbers of corresponding bubbles on an edge: ηe = Nbpe /(Nbpe + 1) and ηc = 1 − ηe . The term F␦ describes an overall flux of xenon atoms deposited on grain boundaries and takes the form: (Xe)

F␦ = Fgb

+ Frsi



where Fgb is the total xenon-atoms flux out of the fuel grain (including diffusion and biased-migration fluxes of gas atoms and bubbles, fluxes due to grain-boundary sweeping, etc.) and Frsi is the sum of re-solution fluxes from all extra-granular bubbles: (Xe)

f Yf Frsi = ωrsi

e (Xe) c (Xe) + ωrsi Ye + ωrsi Yc ,


ωrsi being the irradiation induced re-solution rates. The flux matches, Eq. (45), self-consistently determines the build-up concentration barrier C␦ of the re-solution


layer, which reduces the flux from the grain Fgb , on the one hand, and determines the net flux deposited on the grain boundary F␦ ≈ DC␦ /δ, on the other: C␦ =

dgr δ Frsi , 6Dg


where δ is the effective thickness of the re-solution layer estimated as δ ≈ 10 nm. The flux of Xe atoms from grain faces to grain edges and corners, Ff→ec , is calculated in the meanfield approximation: Ff→ec = F␦ (1 − Af )


ke2 , + ke2


where kf2 and ke2 are the total sink strengths of face bubbles and face periphery (edges), Af is determined in Eq. (40). Finally, transport of Xe atoms from the grain to the open porosity is described in accordance with Eq. (43) by (replacing the index i by Xe): d (Xe) Y = [Fe + Fc + Ff (1 − θf )](1 − θe ). dt out


The model was validated against the tests (Une and Kashibe, 1990) and (Kashibe and Une, 1991). In these tests, 3 and 4 BWR cycle samples with burnup ≈2.4 and 2.9% FIMA (≈23 and 28 GWd/t), respectively, were taken from the outer-pellet region (between rim and mid-point), and the fractional coverage of grain faces by bubbles was evaluated from SEM photographs as ≈6 and 10% (Kashibe and Une, 1991), respectively. Despite such low values of the grain-face coverage, significant fractional fission-gas release (up to 20–30%) during their base irradiation was measured by pin puncture tests from these samples. Therefore, a noticeable gas release from these fuel samples occurred at coverage far below the saturation value ≈50% and without visible bubble interlinking on the grain faces. Results of the model calculations are compared with the measured values in Table 1 and confirm the important conclusion on the onset of gas release at grain face coverage far below the saturation value. In particular, this conclusion allows application of the new model to analysis of MOX fuel irradiated in PWR rods, in which fission gas release is relatively high under normal-irradiation conditions despite there is no evidence of interconnected intergranular bubbles.

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200 Table 1 Modelling of experiments (Une and Kashibe, 1990; Kashibe and Une, 1991) Face bubble diameter (nm) MFPR calculations 216 Experiment 229

Fractional coverage (%)

Gas release (%)

5.9 6.6–10.1

10 10–20

2.5.3. Intergranular-swelling model Coalescence of face bubbles due to their random migration is considered as the main mechanism of grain-face bubble relaxation. The model of this phenomenon proposed in (Zacharie et al., 1998) for annealing conditions is refined and extended to the general case of fuel-operation conditions (e.g. steady-state irradiation and temperature transients) (Berdyshev and Veshchunov, 2002), on the basis of the general theory (Geguzin and Krivoglaz, 1971) of grain-face-bubble migration and coalescence. The general equation for the evolution of bubblesurface concentration is written in the form: dXf = ωnuc Xg2 − ωcls Xf2 , dt


where ωnuc is the frequency of bubble nucleation, ωcls the frequency of bubble coalescence, and Xg and Xf are the surface concentrations of atoms and bubbles on grain faces, respectively. After cessation of fuel irradiation, the surface concentration of atoms considerably drops due to their rapid sinking into bubbles, and the nucleation term in Eq. (50) tends to zero: dXf = −ωcls Xf2 . dt

8πDb ≈ 8παDb ln((Db /2R2f )τ0 )

Bubble diffusion coefficient Db is generally represented in the form:    1/3 q Q Ω Db = D0 exp − , (53) RT Rf where the exponent q depends on the dominating microscopic mechanism of bubble migration. By fitting procedure of the measured and calculated results, the authors of (Zacharie et al., 1998) evaluated the main parameters in Eq. (53) as Q = 310 kJ/mol and q = 3.4, which corresponds to the mixed mechanism (between surface, q = 4, and volume, q = 3) of bubble movement on the grain surface. After attainment of the grain-face-saturation coverage, the face-bubble concentration additionally obeys the limiting coverage condition, Eq. (41). Under isothermal-annealing conditions, solution of the coupled Eqs. (51) and (40), (41) has an asymptotic form:  4(q + 2)A∗ αD0 Ωq/3 (0) R f ≈ Rf + sin2 θf   1/q+2 Q × exp − . (54) t RT The pre-exponential factor in Eq. (53) is determined by fitting of MFPR calculation results for Rf (t) to measurements (Zacharie et al., 1998): αD0 ≈ 30 m2 /s. Validation of the new model against the fuel grain face microstructure and swelling measurements (Zacharie et al., 1998) allows reasonable agreement of the code predictions with measurements (Figs. 11 and 12). 2.6. Fuel oxidation in steam/hydrogen mixtures


The coalescence frequency of bubbles randomly moving on a surface can be represented by the formula derived in (Geguzin and Krivoglaz, 1971): ωcls =



where τ 0 is the characteristic time of the two-fold increase of the mean bubble radius Rf . Being a weak function of its argument (and thus, of Db ), logarithm can be approximated by a constant value, i.e. α ≈ const.

The oxygen/uranium-equilibrium phase diagram (SCDAP/RELAP5 MOD2 CODE, 1990) shows that at temperatures above 1200 K a single phase of nonstoichiometric UO2+x exists in a wide range of the composition. The extent of urania oxidation significantly affects both diffusivities of FP elements and distribution of FP elements between chemical states. For this reason, the MFPR code includes description of interactions between oxide fuel and oxidizing/reducing gas. The model implemented in the current version of MFPR describes fuel oxidation/vaporization in steam/hydrogen/inert mixtures by self-consistent con-


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

where R is the radius of the fuel sample assumed to be a cylinder, jO(sol) the boundary oxygen flux, and xs is the surface stoichiometric deviation. To simplify the problem, this equation is supplemented by an approximate algebraic relationship between x¯ , xs and jO(sol) , which is derived by joining of the two asymptotic solutions of the solid-state-diffusion problem (Dobrov et al., 1998). In a steady-state regime, the surface stoichiometry is determined by the oxygen partial pressure at the gas/solid interface,POs 2 , according to the equation: POs 2 = P¯ O2 (xeq ),

Fig. 11. Calculated number of face bubbles per millimetre of grain boundary as a function of annealing time at different temperatures. Dots are experimental data (Zacharie et al., 1998).

sideration of oxygen diffusion within UO2 matrix, effects of oxygen exchange at the gas/solid interface, fuel vaporization and mass transfer in the multicomponent gas phase. The process of fuel oxidation/reduction is described in terms of average stoichiometric deviation, x¯ . For the case of slow vaporization kinetics, integration of the solid-state-diffusion equation for the oxygen concentration in the fuel matrix yields: 1 d¯x R = jO(sol) (xs ), 2 dt



where P¯ O2 (x) is the equilibrium oxygen pressure over uranium oxide with the stoichiometry x provides by the Lindemer–Besmann correlation (Lindemer and Besmann, 1985). In the case when the actual surface stoichiometry xs differs from the equilibrium one, the surface oxidation kinetics provides the boundary condition for the solidstate diffusion represented in the form: jO(sol) = αF (xs , xeq POs 2 ),


where xeq (POs ) is the solution of Eq. (56), α is the specific constant characterizing the rate of oxygen exchange at the solid/gas interface, the function F(xs , xeq ) is determined by the surface-exchange mechanism, and F = 0 for xs = xeq . The MFPR code includes two models for the oxygen surface exchange effects. In the first, phenomenological model proposed initially by (Carter and Lay, 1970) and modified by (Abrefah et al., 1994), the driving force of equilibration has the linear form:  F (xs , xeq ) = PHs 2 O (xeq − xs ). (58) The second, semi-phenomenological model is based on the results by Gala and Grabke applied to the case of steam-hydrogen atmospheres by (Abrefah et al., 1994):     ¯ O2 (xs ) P F (xs , xeq ) = PHs 2 O 1 − . (59) POs 2

Fig. 12. Calculated inter-granular swelling as a function of treatment time at different temperatures. Dots are experimental data (Zacharie et al., 1998).

For these two models, the rate constant is represented in the Arrhenius form:   Kox E (mol/cm2 /s), α= , Kox = exp − ρUO2 T + E1 (60)

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

where (E = 25 800 K, E1 = 2.523) for the Carter–Lay model, and (E = 21 200 K, E1 = −3.923) for the Gala–Grabke model (Dobrov et al., 1998). In many cases, oxygen transport in the gas phase can be the rate-controlling process of fuel oxidation. Thus, the final step in the formulation of the problem is an approximate solution of the mass transfer problem in the multi-component gas mixture, giving an expression for the gas-phase oxygen flux, jO(gas) , in terms of the surface, POs 2 , and bulk, POb 2 , pressures of oxygen and other gas components. Matching of the oxygen fluxes at the gas/solid interface yields the equation: jO(gas) (POs 2 , POb 2 ) = α F (xs , xeq (POs 2 )),


which in combination with Eqs. (55) and (57), describes the fuel oxidation processes. As a result, the model self-consistently describes either fuel oxidation when POs 2 < POb 2 , or reduction when POs 2 > POb 2 . The model was validated against various tests with laminar regimes of the gas flow in the channel (Cox et al., 1991; Lewis et al., 1995, 1998; Abrefah et al., 1994), and applied for modelling of integral tests.

3. Application to integral tests VERCORS The VERCORS program was conducted for IRSN (Institut de Radioprotection et de Sˆuret´e Nucl´eaire) in the LAMA facility in CEA Grenoble (Ducros et al., 2001).


It aims to reproduce, in an annealing regime, kinetics of FP release in different atmospheres and for different burnups. Two of the six tests performed, VERCORS 4 and 5, are considered owing to their especially reliable results and extreme atmosphere conditions, reducing for VERCORS 4 and pure steam for VERCORS 5. The samples are constituted by three irradiated pellets (∼40 GWd/tU) in their original cladding and two depleted half pellets bounding the irradiated pellets. These samples were stored for 10 years, and then reirradiated at low power in the SILOE facility. The experimental sequence begins less than 40 h after reirradiation. The samples are placed in a zirconia crucible within a furnace, which contains two concentric channels sealed by dense zirconia sleeves. Steam and/or hydrogen is injected into the internal channel containing the sample, whereas helium is injected in the external one containing a tungsten susceptor. Two γ spectrometers give the kinetics of fission-product release. One is focused on the fuel sample and permits on-line FPrelease measurements as well as global release by comparison of pre- and post-test γ measurement. Another is focused on the gas container to follow on-line release of noble gases. In both tests after initial annealing at low temperatures (between 700 and 1300 K in He or mixed H2 /H2 O atmosphere), temperature is increased (1 K s−1 ) up to a nearly 1 h plateau at intermediate temperature (∼1573 K) in a mixed atmosphere (H2 O/H2 flow molar

Fig. 13. (a) Calculated Xe release in VERCORS 4 as compared with experimental data (135 Xe) (Ducros et al., 2001). (b) Calculated Xe relative fractions evolution in grain (atoms in grain, Xe in grain bubbles, Xe in bubbles on dislocations) in VERCORS 4 experiment.


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

Fig. 14. (a) Calculated Cs behaviour in VERCORS 4 as compared with experimental data (137 Cs) (Ducros et al., 2001). (b) Calculated Cs behaviour in VERCORS 5 as compared with experimental data (137 Cs) (Ducros et al., 2001).

ratio near 0.85–0.9) to obtain complete clad oxidation. Then, the temperature is increased again to a high-temperature plateau maintained for half an hour. The main difference between the two experiments is connected with the high-temperature plateau in He/H2 atmosphere in VERCORS 4 and pure steam atmosphere in VERCORS 5 (this gas phase composition is injected just after the oxidation plateau). Examples of MFPR results for gas release in VERCORS 4 are presented in Fig. 13a and b, and for caesium release in VERCORS 4 and 5-in Fig. 14a and b.

Release of the gas intergranular content accumulated during reactor operation takes place during the first step prior to the clad oxidation. This step is not seen in the experimental data, since only the short half-life element (generated during re-irradiation) is measured. Gas behaviour during the temperature ramp and high temperature plateau is determined by physical mechanisms discussed in Section 2.4.4. The calculated Xe release kinetics presented in Fig. 13a, is in a rather good agreement with measurements. This kinetic behaviour is accompanied with Xe redistribution within the grain

Fig. 15. (a) Calculated behaviour of Cs ternary compounds in VERCORS 4. (b) Calculated behaviour of Cs ternary compounds in VERCORS 5.

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

(among dissolved atoms, intragranular bubbles and bubbles on dislocations), as calculated by MFPR in Fig. 13b. The calculated caesium release kinetics in VERCORS 4 and 5 experiments presented in Fig. 14a and b, respectively, are also in a good agreement with experimental data. Additionally, it is shown in these figures that the behaviour of Cs during the experiments is mainly associated with stability of the ternary compound (grey) phase. In particular, as calculated by MFPR and shown in Fig. 15a for VERCORS 4, during the clad oxidation plateau (∼6000–10 000 s) some of the produced hydrogen reacts with the fuel, resulting in a partial Cs2 UO4 destruction and Cs release. Then after complete clad oxidation the remaining injected steam interacts with fuel, producing some Cs2 MoO4 . In the final step of experiment, in a pure hydrogen atmosphere Cs2 MoO4 and the remaining part of Cs2 UO4 become unstable, giving the final release. In VERCORS 5, a rather similar scheme shown in Fig. 15b is realised. However, in this case with pure steam atmosphere, the transfer of Cs from Cs2 UO4 to Cs2 MoO4 is nearly complete, and the first Cs release step is followed by vaporisation of Cs2 MoO4 near 2300 K, which gives the final release.


multi-phase and multi-component equilibrium at the grain boundary in association with accurate calculation of fuel oxidation. Examples of application to VERCORS 4 and 5 experiments show the possibilities of the code in the frame of severe accident interpretation like the PHEBUS experiments (Dubourg et al., 2005). Further development and validation of the MFPR code are foreseen, which will concern association of the fission-product behaviour with fuel-structure evolution, exhaustive description of chemically active elements and application to LOCA and severe-accident interpretation.

Acknowledgements The authors thank Dr. M.P. Kissane (IRSN) for supporting the code development for many years. Support and collaboration of Dr. P. Giordano (IRSN) are highly appreciated.

Appendix A. Thermo-chemical equilibrium in the multi-phase system The principal phases that appeared in the irradiated UO2 and considered in the model are the following:

4. Conclusions The advanced mechanistic code MFPR for simulation of fission-product release from irradiated UO2 fuel is currently under development in collaboration between IBRAE (Moscow, Russia) and IRSN (Cadarache, France). The present paper includes the main characteristics of the MFPR code, brief description of the new models and their validation against analytical tests and integral VERCORS experiments. The code is based on self-consistent consideration of evolution of various point defects (gas atoms, vacancies and interstitials) and extended defects (gas bubbles, dislocations, vacancy clusters and pores) and their mutual interactions under various irradiation and annealing regimes of UO2 fuel operation. Results of the validation have confirmed the enhanced predictive power of the self-consistent mechanistic approach. In the field of chemically active elements, the mechanistic approach of MFPR is based on complex association of a diffusion–vaporisation mechanism involving

• The fuel-FP oxide solid-solution that includes, apart from UO2 and dissolved oxygen: caesium Cs(c), alkaline earth metals Ba(c) and Sr(c) and their oxides BaO(c), SrO(c), zirconium and niobium in the forms Zr(c), ZrO2 (c), Nb(c), NbO(c), NbO2 (c), rare earth elements La(c), Ce(c), Eu(c), Nd(c), and their oxides LaO3/2 (c), CeO3/2 (c), CeO2 (c), EuO3/2 (c), EuO(c), NdO3/2 (c), metalloid Sb(s) and SbO3/2 (c), noble metals Mo(s), Ru(s) and their oxides MoO2 (c), RuO2 (c). • The metal phase composed of Mo(c) and Ru(c). • The phase of complex ternary-compounds (grey phase) including molybdates, zirconates and uranates of Ba, Sr and Cs in the form: BaUO4 (c), SrUO4 (c), Cs2 UO4 (c), BaMoO4 (c), SrMoO4 (c), Cs2 MoO4 (c), BaZrO3 (c), SrZrO3 (c) and Cs2 ZrO3 (c). • The separate solid-phase of CsI(c). The gas phase with the main components: Xe(g), Te(g), I(g), Cs(g), CsI(g), Cs2 MoO4 (g), MoO3 (g),

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200


(MoO3 )2 (g), (MoO3 )3 (g), RuO(g), RuO2 (g), RuO3 (g), Ba(g), Sr(g), ZrO(g), LaO(g), CeO(g), NdO(g), NbO(g), O2 (g). Thermo-chemical equilibrium in the system determines partitioning of elements between allowed compounds, solid phases and the gas phase. In the MFPR model, chemical equilibrium is treated in the framework of the equilibrium constant method based on the mass action law and discussed, for instance, in Imoto (1986), Ball et al. (1989), Cordfunke and Konings (1993), and Baibuz et al. (1985). A general approach to description of the ‘UO2 – FP–O’ system used in the MFPR model is formulated as follows (Baibuz et al., 1985). The thermo-chemical equilibrium in the multiphase system with np phases that consist of ns species composed of ne elements is considered. It is supposed that the system includes all simple (mono-atomic) substances and ns ≥ ne . Each species is to be a component of only one phase. Therefore, for atomic molybdenum, ruthenium and antimony, two condensed states are considered, that is Mo(s) and Mo(c), Ru(s) and Ru(c), Sb(s) and Sb(c), which are the constituents of solid solution (SS) and metal phase (MPh), respectively. The reaction of formation of each species Si from the ‘basic’ components can be written as Si =


bij Sj , i = 1, . . . ns ,



where bij are the stoichiometric coefficients. For equilibrium in the reactions of formation, Eq. (A.1), the mass action law yields: A i = Ki


(Aj )bij ,



where Ai is the activity of the i-th species. For the gaseous species Ai coincides with the partial pressure in atm: Ai = Pi .


It is assumed that the partial pressures are related to moles by the ideal gas law: Pi = Ni

RT , p0 Vgas


where Ni is the number of moles of the species ‘i’, T the temperature in K, R the universal gas constant, and Vgas is the volume of the gas phase, p0 = 1.01325 × 105 kg/(m s2 atm). In the case of condensed species, except dissolved oxygen, the activity is given by Ai = γi Xi .


where γ i is the activity coefficient, Xi is the mole fraction of the species ‘i’ defined as

Ni Nk , (A.6) Xi = (p) , M (p) = M k

where M(p) is the sum of moles in the condensed phase ‘p’ that includes this species. There are no reliable experimental data and well established models for the activities of the majority of species in the considered system. In the current MFPR model the activity coefficients γ i are set equal to 1 for almost all species. Exception is made for the following low soluble elements: condensed caesium, Cs(c), metals Mo, Ru and Sb, and alkaline earth elements Ba and Sr. The noble metals are known to be practically immiscible with UO2 . On the other hand, in the presence of the metallic phase the activity coefficients of these species in the phase SS can be estimated as an inverse maximal solid solubility. Therefore, in accordance with estimates presented in (Kleykamp, 1993) the activity coefficients for Mo(s) and MoO2 (c) are estimated as 105 . The same values of the activity coefficients are ascribed to the pairs Ru(s) and RuO2 (c), Sb(s) and SbO3/2 (c). In analogy to metals and basing on data for maximal solid solubilities (Kleykamp, 1993), the activity coefficient of Cs(c) is set equal to 103 , for barium (and BaO(c)) γ = 103 , and for strontium (and SrO(c)) γ = 102 . Additionally, the activity coefficient for the caesium molybdate is γ(Cs2 MoO4 ) = 102 . Although enhancing of the caesium molybdate activity cannot yet be well grounded by theoretical or experimental arguments, the reasons for such setting follow from the results of modelling of integral tests with a wide range of physical conditions including temperature regimes and oxidation/reduction environment (see Section 2.6). In the MFPR chemistry model, the stoichiometry deviation x in UO2+x is interpreted as a ‘mole frac-

M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

tion’ (negative for hypostoichiometry) of the oxygen dissolved in the matrix, while the number of ‘moles’ of the dissolved oxygen is related to the stoichiometric deviation by NO(c) = xNUO2 .


In the equilibrium state, the partial pressure of the gaseous-molecular oxygen is related to the stoichiometry deviation by the equation: PO2 = P¯ O2 (x, T ),


while the activity of the dissolved oxygen is defined as  (A.9) AO(c) = P¯ O2 (x, T ). The function P¯ O2 (x, T ) is described in literature in detail. In the current version of the MFPR module, the well-known correlation proposed by Lindemer and Besmann (1985) is used. The equilibrium constant Ki in the mass action Eq. (A.2) is related to the Gibbs free energy of formation of the species ‘i’ by   f G◦i Ki = exp − , (A.10) RT f G◦i = G◦i −


bij G◦j ,



where G◦j is the chemical potential or the Gibbs free energy per mole of pure species in the standard state. The equilibrium state of the considered system is determined by the solution of the mass action equations which satisfy the mass balance constraints given by Eitot



Nj bji ,



and obeyed the additional condition: Ni ≥ 0.


Here, Eitot is the total number of moles of the element ‘i’ in the system, bij are stoichiometry coefficients introduced in Eq. (A.2). As mentioned above, the condition, Eq. (A.13), is allowed to be violated only for oxygen dissolved in the matrix because the oxygen ‘moles’


defined by Eq. (A.7) can be negative in the case of hypostoichiometry. References Abrefah, J., de Aguiar Briad, A., Wang, W., Khalil, Y., Olander, D.R., 1994. J. Nucl. Matter 208, 98–110. Baibuz, A., Zitserman, V.U., Golubushkin, L.M., Chernov, U.G., 1985. Chemical Equilibrium in Non-ideal Systems. Institute of High Temperatures, Moscow (in Russian). Baker, C., 1977. J. Nucl. Mater. 66, 283. Baker, C., Kileen, C., 1987. Fission gas release during post irradiation annealing of UO2 . Materials for Nuclear Reactor Core Applications. BNES, London, Paper 24, pp. 183. Ball, R.G.J., Burns, W.G., Henshaw, J., Mignanelli, M.A., Potter, P.E., 1989. J. Nucl. Mater. 167, 191–204. Berdyshev, A.V., Veshchunov, M.S. Modelling of grain face diffusion transport and swelling in UO2 fuel,-Preprint IBRAE-2002-14, Moscow, 2002. Brailsford, A.D., Bullough, R., 1981. Philos. Trans. R. Soc. A 302, 87. Bullough, R., Eyree, B.L., Krishan, K., 1975. Proc. R. Soc. London A 346, 81. Carter, C.E., Lay, K.W., 1970. J. Nucl. Mater. 36, 77–86. Cordfunke, E.H.P., Konings, R.J.M., 1988. J. Nucl. Mater. 152, 301–309. Cordfunke, E.H.P., Konings, R.J.M., 1993. J. Nucl. Mater. 201, 57–69. Cox, D.S., Hunt, C.E.L., Liu, Z., Iglesias, F.C., Keller, N.A., Barrand, R.D., O’Connor, R.F. A model for the release of lowvolatility fission products in oxidising conditions, Proceedings of the 12th Annual Conference of the Canadian Nuclear Society, Saskatchewan, Canada, June 9–12, 1991, AECL Report, AECL10440. Dobrov, B.V., Likhanskii, V.V., Ozrin, V.D., Solodov, A.A., Kissane, M.P., Manenc, H., 1998. J. Nucl. Mater. 255, 59–66. Dollins, C.C., Nichols, F.A., 1978. J. Nucl. Mater. 78, 326. Dubourg, R., Faure-Geors, H., Nicaise, G., Barrachin, M., 2005. Fission product release in the first two PHEBUS tests FPT0 and FPT1. Nucl. Eng. Des. 235 (20), 2183–2208. Ducros, G., Malgouyres, P.P., Kissane, M., Boulaud, D., Durin, M., 2001. Nucl. Eng. Des. 208, 191. Evans, J.H., 1997. J. Nucl. Mater. 246, 121. Geguzin, E.Ya., Krivoglaz, M.A., 1971. The Movement of Microscopic Inclusions in Solid State. Metallurgiya, Moscow (in Russian). Grimmes, R.W., Catlow, C.R.A., 1991. Phil. Trans. R. Soc. London A 335, 609–634. Hayns, M.R., 1975. J. Nucl. Mater. 56, 267. Heames et al., T.J. VICTORIA: A Mechanistic Model of Radionuclide Behavior in the Reactor Coolant System Under Severe Accident Conditions, NUREG/CR-5545 SAND90-0756 Rev 1 R3, R4, 1992. Imoto, S., 1986. J. Nucl. Mater. 140, 1927. Kashibe, S., Une, K., 1991. J. Nucl. Sci. Technol. 28, 1090. Kashibe, S., Une, K., 1993. J. Nucl. Mater. 206, 22.


M.S. Veshchunov et al. / Nuclear Engineering and Design 236 (2006) 179–200

Kleykamp, H., 1985. J. Nucl. Mater. 131, 221–246. Kleykamp, H., 1993. J. Nucl. Mater. 206, 82–86. Lewis, B.J., Andre, B., Morel, B., Dehaudt, P., Maro, D., Purdy, P.L., Cox, D.S., Iglesias, F.C., Osborne, M.F., Lorenz, R.A., 1995. J. Nucl. Mater. 227, 83–109. Lewis, B.J., Corse, B.J., Thompson, W.T., Kaye, M.H., Iglesias, F.C., Elder, P., Dickson, R., Liu, Z., 1998. J. Nucl. Mater. 252, 235–256. Lindemer, T.B., Besmann, T.M., 1985. J. Nucl. Mater. 130, 473–488. Matzke, H., 1986. Adv. Ceram. 17, 1. Mikhlin, E.Ya., 1979. J. Nucl. Mater. 87, 405. Moriyama, K., Furuya, H., 1997. J. Nucl. Sci. Technol. 34, 900–908. Nelson, R.S., 1969. J. Nucl. Mater. 31, 153. Nogita, K., Une, K., 1994. Nucl. Instrum. Meth. Phys. Res. B91, 301. Osborne, M.F., Lorenz, R.A., 1992. Nucl. Safety 33, 344. Ozrin, V.D., Shestak, V.E., Tarasov, V.I., Veshchunov, M.S. Modelling of fission gas release during high-temperature annealing of irradiated UO2 fuel”, Preprint IBRAE-2002-19, Moscow, 2002. Ray, I.L.F., Thiele, H., Matzke, Hj., 1992. J. Nucl. Mater. 188, 90. Rest, J., Zawadzki, S.A. FASTGRASS, A Mechanistic Model for the Prediction of Xe, I, Cs, Te, Ba and Sr release from Nuclear Fuel under Normal and Severe-Accident Conditions, NUREG/CR5840 TI92 040783, 1994. SCDAP/RELAP5 MOD2 CODE, Manual, vol. 4: MATPRO—A Library of Material Properties for Light Water Reactor Accident Analysis. NUREG/GR-5273, EGG-255, vol. 4, 1990.

Shestak, V.E., Tarasov, V.I., Veshchunov, M.S. Modelling of Defect Structure Evolution in Irradiated UO2 Fuel in Framework of the MFPR Code, Preprint IBRAE-2004-02, Moscow, Russia, 2004. Stehle, H., Assmann, H., 1974. J. Nucl. Mater. 52, 303. Turnbull, J.A., Cornell, R.M., 1970. J. Nucl. Mater. 36, 161. Une, K., Kashibe, S., 1990. J. Nucl. Sci. Technol. 27 (11), 1002–10016. Veshchunov, M.S., 2000. J. Nucl. Mater. 227, 67–81. Veshchunov, M.S., Berdyshev, A.V., Tarasov, V.I. Development of Fission Gas Bubble Models for UO2 Fuel in Framework of MFPR Code, Preprint IBRAE-2000-08, Moscow, 2000. Veshchunov, M.S., Tarasov, V.I. Modelling of grain growth kinetics in irradiated and non-irradiated UO2 fuel pellets in framework of the MFPR code, Preprint IBRAE-2002-24, Moscow, Russia, 2002. Whapham, A.D., Sheldon, B.E., Electron Microscope Observation of the Fission-Gas Bubble Distribution in UO2 AERE-R-4970, 1965. Whapham, A.D., 1966. Electron microscope observation of the fission gas bubble distribution in UO2 . Nucl. Appl. 2, 123. White, R.J., Tucker, M.O., 1983. J. Nucl. Mater. 118, 1–38. Yoo, M.H., 1977. J. Nucl. Mater. 68, 193. Zacharie, I., Lansiart, S., Combette, P., Trotabas, M., Coster, M., Groos, M., 1998. J. Nucl. Mater. 255, 85–91.