Multiphysics simulation of VVER-1200 fuel performance during normal operating conditions

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Multiphysics simulation of VVER-1200 fuel performance during normal operating conditions

Khaled M. Yassin ，

Mohamed H. Hassan，

Mohammad M. Ghoneim，

Mostafa S. Elkolil，

Adel Alyan，

Said A. Agamy

Nuclear Science and Techniques

Vol.34, No.2

Article number 28

Published in print Feb 2023

Available online 25 Feb 2023

DOI：10.1007/s41365-023-01172-9

34600

Nuclear fuel performance modeling and simulation are critical tasks for nuclear fuel design optimization and safety analysis under normal and transient conditions. Fuel performance is a complicated phenomenon that involves thermal, mechanical, and irradiation mechanisms and requires special multiphysics modules. In this study, a fuel performance model was developed using the COMSOL Multiphysics platform. The modeling was performed for a 2D axis-symmetric geometry of a UO₂ fuel pellet in the E110 clad for VVER-1200 fuel. The modeling considers all relevant phenomena, including heat generation and conduction, gap heat transfer, elastic strain, mechanical contact, thermal expansion, grain growth, densification, fission gas generation and release, fission product swelling, gap/plenum pressure, and cladding thermal and irradiation creep. The model was validated using a code-to-code evaluation of the fuel pellet centerline and surface temperatures in the case of constant power, in addition to validation of fission gas release (FGR) predictions. This prediction proved that the model could perform according to previously published VVER nuclear fuel performance parameters. A sensitivity study was also conducted to assess the effects of uncertainty on some of the model parameters. The model was then used to predict the VVER-1200 fuel performance parameters as a function of burnup, including the temperature profiles, gap width, fission gas release, and plenum pressure. A compilation of related material and thermomechanical models was conducted and included in the modeling to allow the user to investigate different material/performance models. Although the model was developed for normal operating conditions, it can be modified to include off-normal operating conditions.

VVER-1200Fuel performanceCOMSOL codeZr-1%Nb claddingUO2 fuel rod

Introduction

Nuclear fuels are exposed to highly challenging harsh operational conditions in the reactor core, where corrosive media, mechanical stresses, and high temperatures are combined with intensive radiation effects on fuel elements [1]. Studies that have attempted to understand the irradiation behavior of fuels in nuclear reactors have found considerable alterations in the geometry, dimensions, composition, and microstructure of fuels during and after irradiation [2].

Understanding and predicting the evolution of nuclear fuel properties are necessary for fuel design, operational behavior, and long-term storage. Nuclear fuel performance codes are necessary because of the difficulty and high cost of nuclear experimental measurements. These multiphysics performances inextricably link the mechanical, thermal, and chemical phenomena, which are linked with many control parameters and associated uncertainties.

Early fuel modeling codes used a one-dimensional or asymmetrical-two-dimensional description for the fuel geometry, such as GAPCON [3], FEMAXI [4], FALCON [5], FRAPCON[6], and FRAPTRAN[7]. These codes usually apply several correlations with fixed coupling among the different phenomena. Fuel modeling codes are commonly built as isolated computer software. The limitation of this path is that codes are often restricted to the end user because they have hardcoded functionality. This requires considerable source code modification for various purposes [8].

For VVER Fuel performance modeling, some western PWRs fuel performance codes have been modified to predict fuel performance, including TRANSURANUS [9], FRAPCON [6], and PIN-Micro (based on GAPCON) [3]. Additionally, START-3, a Russian-developed code, was used and validated after feedback modifications from the FUMEX program [10].

Due to the complicated nature of fuel performance, multiphysics modeling has been developed to predict fuel performance. BISON is a code developed by the Idaho National Laboratory (INL) [11] and can be applied to various multidimensional abilities. This code is strongly upgradable and used for treating1D, 2D(axis-symmetric), and 3D configuration nuclear fuel behavior scenarios [12]. This is achieved by using the numerical solving abilities of the Multi-physics Object Oriented Simulation Environment (MOOSE), a commonly used objective, open-source numerical modeling platform that has applications outside the nuclear field. Notably, BISON is open-source code. However, it has a restricted distribution and should be provided by INL.

The COMSOL^® Multiphysics framework is adaptable commercial finite element software that can solve a broad spectrum of ordinary and partial differential correlations with a type of coupling for its dependent variables. It offers flexibility when managing physical problems. The greatest of them is the ability to define the system’s physics in the form of PDEs for any given problem by using a built-in mathematical module. This ability can be used to simulate irradiation effects, which are not included in other Multiphysics platforms such as ANSYS or ABAQUS, which are more suitable for studying Computational Fluid Dynamics cases, as well as solid mechanics.

COMSOL was used to model nuclear fuel performance for PWRs [13] and CANDU reactors [8, 14]. Prudil et al. [8, 14] developed the FAST code using the COMSOL Multiphysics platform to simulate the fuel performance of a CANDU reactor under normal and off-normal conditions. The FAST code postulates axis-symmetry, permitting the use of a 2D geometry illustration of the fuel element, considering pellet dishing and chamfering.

In this study, the development of an advanced VVER-1200 nuclear fuel performance modeling was undertaken using the COMSOL^® Multiphysics platform (version 5.3a). This approach helped bypass any difficulties or restricted distribution of codes developed by other research/commercial organizations. In addition, using multiphysics modules enables the inclusion of irradiation models describing the irradiation behavior of fuel performance. Finally, developing the code from scratch enables the future incorporation of different material properties and irradiation models to investigate different performance phenomena and compare the results with available experimental measurements. This capacity-building step is necessary for operating and regulatory organizations.

Model development

The COMSOL Multiphysics platform was used to develop the current model and predict fuel performance. In this context, the material properties and behavior were modeled, analyzed, and executed. Diverse physical phenomena are nonlinear with nonlinear material properties and highly coupled. The geometry simulation involved fuel pellets, clads, and gaps. The material properties are defined as temperature- or burnup-dependent parameters or other function-based variables. Most models were implemented in the simulation in a completely coupled manner. The modeling of some COMSOL coupling operators (e.g., mapping, integration) used in the FAST code was applied in this study.

2.1

Modeling geometry

The fuel rod modeled in this research has 10 individual uncracked UO₂ pellets with an axial centerline hole of radius that is approximately 0.6 mm and Zr-1%Nb cladding with an initial 60 μm radial gap and short gas plenum at the top. Modeling was performed on a 2D axisymmetric pellet [15].

Figure 1 presents the geometry of the pellet and meshing implemented in the modeling. Variations in the axial direction were not considered, owing to the power level and coolant temperature. As such, periodic boundary conditions were used in the modeling in the z-coordinate to establish an endless pellet stack [8, 14]. This methodology saves computational resources and time by reducing the number of degrees of freedom in the modeling. The methodology has been successfully applied [8, 12, 14, 16, 17]. The validity of this methodology was proved by contrasting the outcomes acquired using one pellet geometry and many pellet geometries, which revealed a minor deviation in the prediction of the temperature and amount of fission gas released.

Fig. 1

2D-axisymmetric single pellet (domain 1) with clad (domain 2) geometry and meshing applied in the modeling.

The dimensions used in this study represented a typical design of the VVER-1000 fuel element, used in the validation analysis, and the VVER-1200 fuel element, the focus of this study, and were as follows: cladding outer diameter, 9.1 mm; cladding thickness, 0.69 mm; pellet diameter, 7.6 mm; and axial hole diameter, 1.2 mm [18-21].

A symmetrical convective boundary condition on the clad outer surface was applied to simulate the heat flux from the clad to the coolant bulk. The time evolution of the linear power was as follows: over 4 h, it increased linearly and then remained constant at 420 W/cm for 4.0 years to reach a maximum burnup of approximately 70 MWd/kgU. Notably, the nominal burnup for the VVER-1200 fuel is 60 MWd/kgU in the IAEA Status report [22]- VVER-1200 (V-491) and 65 MWd/kgU in [19], whereas in [20], it is 70 MWd/kgU. As such, in this study, the highest burnup is considered to predict the worst damage to force the fuel to its design limits. These design parameters and operating conditions were used for a typical VVER-1200, as described in [19, 20].

2.2

Materials properties library

Table 1 Highlights the models of material properties employed in this work, which are primarily from MATPRO (Material Properties for Light Water Reactor) versions, updates (for VVER-1200 cladding), and open literature references [23-25].

Main material property models implemented in the current work

Property	UO₂ models	Zr-1%Nb (E110)
Thermal conductivity (W/mK)	Function of temperature, stoichiometry, porosity, fission product buildup, radiation damage [23]	Function of temperature [24]
Density (kg/m³)	Function of initial density, material strain [26]	Function of initial density, material strain [25]
Heat capacity(J/kg $\cdot$ K)	Function of temperature [24]	Function of temperature, heating rate [25]
Thermal expansion	Function of temperature [27]	Function of temperature [25]
Young’s modulus (N/m²)	Function of temperature [23]	Function of temperature [25]
Poisson ratio	0.316 [23]	Function of temperature [25]
Emissivity	Function of temperature [23]	Function of temperature [23]

2.3

Heat production and conductance

The heat transport model was conducted by utilizing the COMSOL built-in module “Heat transfer in Solids”. The module calculates the temperature T (Kelvin) as a dependent variable, using the heat conduction equation: $ρ C_{P} \frac{\partial T}{\partial t} = \nabla . (K \nabla T) + Q_{prod},$ (1) where ρ, K, and C_P are the density, thermal conductivity, and specific heat capacity, respectively, as material properties. The heating term, Q_prod, defines the volumetric heat generated inside the fuel. In the case of a single pellet, the heat-source term is assumed to be uniformly distributed, which is given by [28, 29] $Q = \frac{P_{lin}}{π . a_{pellet}^{2}},$ (2) where P_1in denotes linear power, and a_pel1et denotes pellet radius.

Heat is transferred through the pellet-to-clad gap and is modeled assuming 1-dimensional stationary heat transfer by utilizing the formula of radial heat flux as follows: $Q_{r} = (h_{gas} + h_{solid} + h_{rad}) (T_{fuel} - T_{Clad}),$ (3) where the overall heat conduction through the gap is the summation of the conductance caused by solid-to-solid contact (h_solid), gas conductance (h_gas), and conductance owing to radiative heat flux (h_rad). This formula is utilized between the fuel and cladding, following the modeling path implemented in the FAST and BISON codes [8, 12]. Gas conduction h_gas is represented by the formula suggested by Ross and Stoute [30]: $h_{gas} = \frac{k_{g} (T_{gas})}{C_{r} (R_{f} + R_{c}) + d_{gap} + g_{1} + g_{2}},$ (4) where k_gas is the gas conductivity in the gap and is represented by the formula used in [17, 31]; d_gap refers to the gap thickness calculated dynamically in the solid mechanics module; C_r refers to the surface roughness coefficient; r₁ and r₂ are the roughness of the fuel and clad surfaces, respectively; and g₁ and g₂ refer to the temperature jump distances for the two surfaces.

The heat conduction caused by the solid-solid contact, h_solid, is calculated using the empirical model proposed by Olander [32]: $h_{solid} = C_{s} (\frac{2 k_{1} k_{2}}{k_{f} + k_{s}}) (\frac{P_{c}}{δ^{0.5} H}),$ (5) where C_s is a constant (= 1.0), k₁ and k₂ are the thermal conductivities of the fuel and cladding when in contact; P_c is the contact pressure, which is simulated using a penalty approach, as explained in the mechanical deformation section; δ is the gas film average thickness; and H is the Meyer’s hardness of the cladding material.

The conductance resulting from radiant heat transfer, h_rad, was computed using a 1-D infinite-parallel plates approach (i.e., considering the view factor = 1). According to the Stefan-Boltzmann equation $q_{rad} = \frac{σ}{\frac{1}{ε_{1}} + \frac{1}{ε_{2}} - 1} (T_{1}^{4} - T_{2}^{4}) = h_{rad} (T_{1} - T_{2}),$ (6) where ε₁ and ε₂ are the emissivities of the radiating surfaces; σ is the Stephan-Boltzmann constant; and T₁ and T₂ are the fuel and cladding surface temperatures, respectively. Radiative conduction can be expressed as $h_{gab, rad} = \frac{σ}{\frac{1}{ε_{e, f}} + \frac{1}{ε_{e, c}} - 1} (T_{f}^{2} + T_{C}^{2}) {(T_{f} + T_{C})}^{} .$ (7) The heat transport from the cladding to the coolant was calculated using a thin-film heat transfer approach. The peripheral heat flux transfer from the cladding to the coolant Q_r can be computed as $Q_{r} = \frac{(T_{clad}^{} - T_{Coolant}^{})}{\frac{t_{oxide}^{}}{k_{oxide}^{}} + \frac{1}{h_{film}^{}}} - Q_{oxide},$ (8) where T_coolant and T_cladding are the coolant and cladding temperatures, respectively. k_oxide is the therma1 conductivity of the zirconium-oxide (ZrO₂) layer, and t_oxide refers to the layer thickness. h_film is the effective clad heat transfer coefficient, and Q_oxide defines the heat flux resulting from the oxidation reaction of cladding, as proposed by Leistikow and Schanz [33]: $Q_{oxide} = \frac{42.426 \times 10^{9}}{1.56} \frac{d t_{oxide}}{d t},$ (9) where the ZrO₂ layer depth, t_oxide, depends on the corrosion rate. The formula in MATPRO-A for clad oxidation was adopted for the outer surface of the cladding. This formula is recommended for E110 cladding with different activation energies and rate constants [25]. Cladding oxidation during normal LWR operation occurs in two stages, depending on the oxide thickness and, to some extent, on the temperature of the oxide. For thin oxides, the oxidation rate is controlled by the entire oxide layer. When the oxide layer thickens, a change in the outer portion occurs, and further oxidation is controlled by the intact linker layer [23]. In the pre-transition stage, the corrosion rate was suggested by the Arrhenius correlation: $\frac{d t_{oxide}}{d t} = C_{1} (\frac{- Q_{1}}{R T_{i}}), for t_{oxide} \leq t_{trans},$ (10) where Q₁ and C₁ are the activation energy in the pre-transition oxidation stage and rate constant, respectively. From [24], the values of activation energy Q₁=14680 and rate constant C₁= 5.19×10^-7 for the E110 alloy were obtained.

For the post-transition period, the oxidation growth rate is given by $\frac{d t_{oxide}}{d t} = C_{2} (\frac{- Q_{2}}{R T_{i}}), for t_{oxide} \geq t_{trans} .$ (11)

In addition, from [24], the values of activation energy Q₂=15355 and rate constant C₂= 17.72×10^-7 for the E110 alloy were obtained.

2.4

Deformation mechanics

2.4.1

Elastic deformation

In a nuclear reactor, the fuel element deformation is described using Cauchy's equation [8, 14, 16], which is included in the structural mechanics module in COMSOL in the following form: $- \nabla \cdot σ = F_{v},$ (12) where σ refers to Cauchy’s stress tensor, and F_v is the body force per unit volume, which is caused by applied forces, fission gas swelling, fuel densification, clad creeping, and thermal expansion. The stress was computed using the following linear elastic intrinsic model: $σ = [C] [ε],$ (13) where C and ε are the material property matrix and the elastic strain vector, respectively. The elastic strain was computed as the total displacement. In the case of 2D axis-symmetry, which is the case in this work, the strain components are represented by $ε_{r} = \frac{\partial u}{\partial r}, ε_{z} = \frac{\partial v}{\partial z}, ε_{r z} = \frac{1}{2} (\frac{\partial v}{\partial r} + \frac{\partial u}{\partial z}), ε_{r θ} = ε_{θ z} = 0.$ (14)

In this study, a modified penalty method used in the FAST code [8] was used to simulate the pellet-clad mechanical contact with an applied force perpendicular to the clad inner surface. The penalty function uses the interface pressure, P_i, which increases exponentially and reaches P_est, the estimated contact pressure. Once the fuel and clad surfaces contact each other, P_i linearly increases with the penetration distance (d_pen) and varies at a specific rate depending on the value of P_f, that is, the penalty factor, to stop the high stiffness that might occur due to the exponential term [8, 14]. $P_{i}^{} = {\begin{array}{l} P_{est} \exp (\frac{P_{f} d_{pen}}{P_{est}}) & d_{pen} < 0 \\ P_{est} + P_{f} d_{pen} & d_{pen} \geq 0 \end{array}$ (15)

The penalty factor and estimated contact pressure are parameters that can be adapted to the system under consideration. For the LWR fuel elements, P_est = 1 MPa and P_f =5×10¹³ Pa/m were applied for the pellet-clad mechanical contact [16]. Penalty contact was applied in the COMSOL built-in contact model, which is a part of the solid mechanics module. The contact tolerance was set to 2 µm, which is close to the fuel and cladding surface roughness. Using the penalty function and contact pressure allows for simple management of pellet cladding mechanical interaction (PCMI) because the phenomena are complicated but cannot be neglected [13].

2.4.2

Thermal strain/expansion

Fink’s thermal expansion model [27] was used. It assumed an isotropic UO₂ fuel pellet. The thermal strain ε_th is computed as $ε_{th} = {\begin{array}{l} 0.9973 + 9.082 \times 10^{- 6} T - 2.705 \times 10^{- 10} T^{2} + 4.391 \times 10^{- 13} T^{3}, 273 K < T \leq 923 K \\ 0.9672 + 1.179 \times 10^{- 5} T - 2.429 \times 10^{- 9} T^{2} + 1.219 \times 10^{- 12} T^{3}, 923 K < T \leq 3120 K \end{array}$ (16) using the assumption that the strain is zero at the reference temperature of 273 K.

2.4.3

Fuel densification and fission products’ swelling

Fuel densification can be described by employing the empirical correlation developed by Hastings and Evans[34]: $F_{poreR} = 0.6 - exp (- 0.506 - 8.67 \times 10^{- 10} T^{3} (1 - \exp (- 2.867 \times 10^{- 2} B u)))$ (17) $ε_{vol, dens} = \frac{Δ V_{dens}}{V_{0}} = \frac{1 - P_{o}}{1 - P_{o} (1 - F_{dens})}$ (18) where F_dens is the fuel porosity fraction removed by the densification process, and p₀ is the initial porosity.

The swelling caused by the effect of fission products (solids and gases) is computed by applying empirical correlations given in MATPRO [23].

Solid fission product swelling results in a volumetric strain, computed as a linear function of burnup: $\frac{δ ε_{sw - s}}{δ t} = 5.577 \times 10^{- 5} ρ \frac{δ B u}{δ t}$ (19) where ε_sw-s refers to the solid fission product swelling, ρ is the density (kg/m³), and Bu is the local burnup in fissions/atoms U.

The swelling caused by fission gases is computed using a semi-empirical correlation in MATPRO [23]: $\frac{δ ε_{sw - g}}{δ t} = 1.96 \times 10^{- 31} ρ \frac{δ Buf}{δ t} {(2800 - T)}^{- 11.73} e^{[- 0.0162 (2800 - T)]}^{} e^{- 0.0178 ρ B u}$ (20) where ε_sw-g is the volumetric gas swelling, and Bu and δBu are the burnup and burnup increments (fissions/atoms-U), respectively.

2.4.4

Clad thermal and irradiation creep.

The creep rate correlation is provided in MATPRO-11[35]. This relation computes the thermal and irradiation creeps in the radial direction. ε_ϕ is the thermal creep rate, which is calculated as $\frac{\partial ε_{\emptyset}}{\partial t} = 5 \times 10^{- 23} σ_{ϕ}^{2} (3.47 \times 10^{- 23} \frac{σ_{ϕ}^{3}}{| σ_{ϕ} |} \exp (\frac{- U}{T}) - ε_{ϕ}) \exp (\frac{- U}{T}) .$ (21) where σ_ϕ is the stress component in the circumferential direction; T is the temperature in Kelvin; and U is the activation energy divided by R, which is the ideal gas constant. $U = 2 l 2.7 - 0.5324 T + 1.17889 \times 10^{- 4} T^{2} + 3.3486 \times 10^{- 7} T^{3}$ (22)

The creep is driven by the effects of fast neutrons. We assumed that fast neutrons had the largest irradiation effect because they increased the concentration of point defects. The strain rate due to irradiation creep is described in [36]. $\frac{\partial ε_{\emptyset}}{\partial t} = \frac{2.2 \times 10^{- 7} σ_{\emptyset} exp {(\frac{- 5000}{T})}^{} \emptyset_{Fflux}^{- 0.65}_{}}{T^{7}},$ (23) where ϕ_Fflux is the flux of the fast neutrons (1/m^2· $\cdot$ s) with a kinetic energy of more than 1 MeV.

2.5

Fission gases production and release

The generation of fission gases in nuclear fuel has a strong impact on the thermomechanical performance of the fuel rods. In the early stage, fission gases tend to accumulate into bubbles, resulting in fuel swelling, which promotes pellet-cladding gap closure and leads to PCMI. At a later stage, fission gas release (FGR) to the fuel element free volume causes pressure rise-up and thermal conductivity degradation of the rod filling gas.

In the current work, the applied FGR model is presented in [37]. The release process is assumed to occur in two stages. In the first stage, atoms of fission gases are generated inside the fuel grains and then move toward the grain edges, where they pile up to form intergranular bubbles. In the second stage, the intergranular bubbles continue to grow until they connect (saturate) and discharge the gas to the fuel element free volume.

2.5.1

Release of fission gases to the grain boundaries

The release of fission gases to fuel grain boundaries can be simulated using Booth’s diffusion model [38, 39]. In this model, the fuel grains are considered to be idealized and homogenous spheres in which fission gases are generated symmetrically and demonstrate Fickian diffusion. The grain boundaries were considered to be ideal sinks, with atoms diffusing toward the grain surface to enter the intergranular bubbles. $\frac{\partial C}{\partial t} = - D \nabla^{2} C + P_{fg},$ (24) where C is the fission gas atom concentration in the fuel grains, P_fg is the volumetric rate of generation of fission gas atoms, and D is the diffusion coefficient of fission gases within the fuel grain.

The average local UO₂ grain diameter was computed using the grain growth correlation given in [40]. The grain growth rate (m/s) was calculated as follows: $\frac{d g_{d}^{}}{\partial t} = k_{g} (\frac{1}{g_{d}} - \frac{1}{g_{max}} - \frac{1}{g_{ir}}),$ (25) where g_d is the average grain diameter (2g_r), and k_g (m²/s) is the rate constant, which is computed as follows: $k_{g} = 1.46 \times 10^{- 8} exp (\frac{- 321001}{T}),$ (26) where T is the temperature (K), and g_max is the restricting grain size (m), which is a function of temperature given by $g_{max} = 2.23 \times 10^{- 3} exp (\frac{- 7620}{T}),$ (27) and g_ir accounts for the irradiation effects on grains, given by $g_{ir} = \frac{6.71 \times 10^{18} exp (\frac{- 5620}{T})}{F_{rate} T},$ (28) F_rate refers to the fission rate density (m⁻³).

This1model does not consider grain size distribution within a region; it only describes average grain size. This approximation provides accurate findings forFGR analysis, regardless of the effect of wide variations in grain size distribution [40].

2.5.2

Fission gas release to the free volume of the fuel element

After fission gases have been released to the grain borders, they are confined as intergranular bubbles between fuel grains. When the gas atoms heaped up on the grain boundaries and achieved G_bsat, the saturated gas concentration, fission gases were liberated to the fuel element free volume. G_bsat(atoms/m³) was computed as described by White and Tucker [39]: $G_{bsat} = \frac{4 r_{f} f (θ_{fg}) f_{B}}{3 k_{B} T {sin}^{2} (θ_{fg})} (P_{ext} + \frac{2 γ_{se}}{r_{f}}) (\frac{3}{g_{r}}),$ (29) $f (θ_{}) = 1 - 1.5 \cos (θ) + 0.5 {cos}^{3} (θ),$ (30) where r_f = 10^-7m is the radius of curvature of the fission gas bubbles, f(θ_fg)1is a function that describes bubble shape, f_B=0.5 is the fraction of the grain surface covered in bubbles when interlinkage occurs, θ_fg is the half-angle between bubble surfaces, T is the temperature in K, k_B = 1.3806 $\times$ l0^-23 J K^-1 is the Boltzmann constant, γ_Se = 0.6261 J · m^-2 is the surface energy of the bubbles, P_ext is the externally applied hydrostatic pressure in Pa, and g_r is the radius of the grains in m. By using these values, the grain boundary saturation becomes $G_{bsat} = \frac{2.0811 \times 10^{16}}{T g_{r}} (P_{ext} + 2.504 \times 10^{6}) .$ (31)

The hydrostatic pressure in the vicinity of the bubble effectively increased the pressure of the gas inside the bubbles. This increased the quantity of gas required for the bubbles to be sufficiently large to be connected. Therefore, a higher hydrostatic pressure decreased the fission gas released from the fuel. For simplicity, the external pressure was assumed to be zero, similar to the approximations made in [13].

Model validation

The validation step is necessary to test the underlying models/assumptions, identify limitations, and help model debugging. In addition, this step is intended to provide a proof of concept for the developed model. The open literature on fuel performance experimental data for VVER under relevant conditions is limited; thus, validation was conducted by comparing the results of this study with those in the literature that use other similar codes. Other studies have discussed either thermal and mechanical behavior or irradiation consequences, such as fission gas release and swelling. For this reason, the validation has been divided into two steps: thermal and mechanical performance and Fission Gas release. These two steps are discussed in detail in the following subsections.

3.1

Thermal and Mechanical Performance

FINIX is a fuel behavior module designed to be integrated as a subroutine into a larger simulation code, where FINIX replaces the corresponding fuel model [41]. FINIX was integrated into Serpent 2 in [42] to study a VVER-1000 fuel rod. The geometry, material properties, and fuel rod specifications were obtained from the UAM benchmark [43]. The inner and outer radii of the pellets were 0.070 cm and 0.378 cm, respectively. The inner and outer radii of the cladding were 0.386 and 0.455 cm, respectively. The fuel was pure UO₂ with an enrichment of 3.3%. The cladding material was Zr+1% Nb. The coolant temperature was 560 K. The fuel rod was depleted until an average burnup of 10 MWd/kgU was achieved. The system was maintained at a linear power of 233 W/cm. Figures 2 and 3 compare the results obtained using FINIX and our COMSOL model, illustrating satisfactory agreement between them. In Figure.2, the COMSOL-based modeling results show a higher temperature in the pellet domain because it includes irradiation effects (i.e, densification, swelling, cladding creep, and thermal expansion), which are not included in this version of the FINIX code, as explained in [42].

Fig. 2

Temperature distribution through fuel element in radial coordinates in the COMSOL model and FINIX module

Fig. 3

Total displacement comparison between COMSOL model and FINIX module

Figure 3 shows the total displacements computed by the developed COMSOL model due to densification, cladding creep, fuel swelling, and thermal expansion, which are temperature dependent. Consequently, the results of COMSOL and FINIX differed slightly.

3.2

Fission Gas Release

Operational and (post-irradiation examination) PIE data for the Zaporozhye NPP, FA-E0325, and VVER-1000 fuel rods were available in the OECD NEA IFPE Database and used to perform comparative calculations among several fuel performance codes [44].

The calculations were performed using three computer codes for the VVER-1000 fuel rod performance analysis: PINw99, TOPRA-2, and TRANSURANUS (V1M1J03). All these results were compared with those obtained using the COMSOL-developed model. As illustrated in Fig. 4, the results of COMSOL show good agreement with the other codes and are closer to the final measurement values.

Fig. 4

Modeling results comparison of fission gas released from VVER-1000 fuel.

Sensitivity analysis

Sensitivity analysis was performed to assess the effects of uncertainty on some of the model parameters. This step was performed by systematically varying the model parameters one at a time and comparing the results to a baseline case. This analysis helps identify parameters that contribute significantly to the uncertainty of the results and therefore require the most attention when modeling. Notably, this is not a comprehensive error analysis because it does not consider the combined effect of uncertainty from multiple parameters simultaneously.

In this work, the sensitivity analysis focused on quantify the changes in the mid-pellet clad strain, circumferential ridge strain, and fission gas release to perturbations in various modeling parameters. This analysis was performed to help manage the large volume of generated data (i.e., several values per case instead of several time/burnup-dependent functions). An uncracked pellet is used in the model. Apparently, it produced results comparable to other codes, as shown in the validation cases. A sensitivity analysis was performed on the six model parameters. These were three code input parameters (i.e., operating conditions of the nuclear fuel)—linear power, coolant temperature, and clad-to-coolant heat transfer coefficient—and three material properties: UO₂ thermal conductivity, UO₂ thermal expansion, and clad thermal expansion. Additionally, these cases were examined for three linear power/temperature levels: 30, 45, and 60 kW/m. This covers a wide range of operating power levels, reaching the maximum allowable linear power in VVER-1200, as well as a simple power history that helps manage calculation times and isolate power transient effects. Because of the length of the sensitivity data, only the representative results obtained for 45 and 60 kW/m are given. This is because they have a higher significant effect than the 30 KW/m case.

The impact of the six model parameters on the FGR is presented in Sects. 4.1 and 4.2. Because the mid-pellet clad strain, and the circumferential ridge clad strain, were not sensitive to the input parameters, their sensitivity relative to the material parameters is outlined in Sects. 4.3 and 4.4, and meshing sensitivity is outlined in 4.5. These results are expected to be valid for the range of variations attempted in this study.

4.1

Sensitivity of FGR to Model Input Parameters

The most sensitive FGR model input parameters were linear power and coolant temperature. High sensitivity to linear power and coolant temperature occurs because they significantly influence the temperature of UO₂. When the temperature increased, the speed of fission gas diffusion increased, and the concentrations at which grain boundary saturation decreased. An increase in linear power leads to larger thermal gradients and higher temperatures. An increase in the coolant temperature effectively increases the temperature of the heat sink, leading to a corresponding increase in all predicted temperatures, assuming a constant thermal conductivity. Both cases also have additional positive feedback because the thermal conductivity of UO₂ decreases as the temperature increases.

The 45 kW/m case resulted in a change from –18.3% to 39% in fission gas release from a ±5% change in the linear power. At 60 kW/m, the gas release changed from –6.5% to 16% under the same conditions. The percentage is high for the 45 kW/m case, partly due to the small initial baseline value. The coolant temperature perturbation resulted in a change of –17.6% to 12.4% at 45 kW/m and a change of –9.5% to 13.2% at 60 kW/m from a ±10% change. The effect is more pronounced in the 45 KW/m case because ±10% is a relatively high change in coolant temperature.

4.2

Sensitivity of FGR to Material Parameters

The FGR model was sensitive to the thermal conductivity of UO₂ because the aforementioned gas release parameters (section 2.5) are temperature dependent. A change of ±20% W/m.K in the thermal conductivity of the UO2 was found to change the fission gas volume by +2.3 mL and -0.6 mL STP (standard temperature and pressure) for the 45 kW/m case. For the 60 kW/m case, the same thermal conductivity change resulted in a change of +4 mL and –3.2 mL at STP.

For the fuel thermal expansion, a change of ±2% was found to change the fission gas volume by +0.08 mL and –0.006 mL at STP for the 45kW/m case. At 60 kW/m, this perturbation led to +0.47 ml and –0.273 ml. Thus, fission gas quantity is less sensitive to changes in UO₂ thermal expansion.

4.3

Sensitivity of Mid-Pellet Cladding Strain to Material Parameters

As we have outlined, the mid-pellet clad strain is not sensitive to the model input parameters. The material parameters are sensitive to only one input parameter: cladding thermal expansion. As expected, an increase in the thermal expansion strain resulted in an increased cladding strain in all cases (and vice versa). Similarly, increasing the thermal conductivity decreased the sheath strain in all cases. At 30 and 45 kW/m, the perturbations in the thermal expansion strain led to a minor change of ±0.01%; however, at 60 kW/m, the thermal expansion strain perturbation produced an increase in the mid-pellet clad strain of +0.04% for a 10% reduction in cladding thermal expansion.

4.4

Sensitivity of Cladding Strain at Circumferential Ridge to Material Parameters

The clad strain at the circumferential ridge was dependent on the same factors, which influence the mid-pellet sheath strain because of the reasons we have discussed. As a result, the sensitivity of the circumferential ridge strainexhibited the same trend as the mid-pellet clad strain.

4.5

Sensitivity of Results to meshing density

The finite element mesh used in COMSOL in this study was generated according to a defined meshing sequence, such as the geometry definition of the model (Fig.1). In this case, a mapped quadrilateral mesh was produced for each subdomain according to the elements defined on the edges.

To study how the code results may be affected by the meshing, the baseline case mesh density used was multiplied by 4 and 8. The baseline results appear to converge well with the finite element mesh because the results show very small changes when the mesh density is increased.

Results and discussions

5.1

Temperature Profiles

Figure 5 presents the estimated temperature profiles at the fuel pellet centerline, fuel outer surface, and clad inner surface throughout the power up and steady state operation by applying the developed model.

Fig. 5

Temperature profile predictions of the developed model

The fuel centerline temperature, pellet outer surface temperature, and clad inner surface are predicted to be approximately 1720 K, 800 K, and 680 K, respectively, during the reactor operation. The centerline temperature continues to increase because of the degradation of the thermal conductivity of UO₂ with the rise in temperature of the fuel as fuel burnup continues [45] (Fig. 6). In addition, the dramatic decrease in temperature between the fuel centerline and the outer surface is due to structural changes, such as fuel densification. The clad creeping down and fuel swelling led to progressive gap closure, followed by a decrease in fuel temperature. According to the obtained results, the gap was closed at a burnup of approximately 48 MWd/kgU (Fig. 7).

Fig. 6

Temperature profile across the pellet radius.

Fig. 7

Calculated gap size evolution

5.2

Gap width

The gap width is a result of the combined effects of swelling, fuel densification, thermal expansion, and clad creeping. As shown in Fig. 7, the primary gap width was approximately 60 μm at the beginning of the operation. The gap thickness was reducedsubstantially because of the thermal expansion of the fuel pellet, as it achieved the operational temperature. Subsequently, the gap width increases for a short time because of fuel densification. However, this phenomenon finally saturates, and fuel swelling causes the gap to shrink until the fuel and clad contact each other. The gap size from the developed model was found to decrease dramatically and reach 3.04 μm (i.e., approaching the surface roughness of the pellet and cladding inner surfaces) at 48 MWd/kgU, showing an earlier gap closure. Once the gap collapsed, the fuel surface temperature remained approximately constant.

5.3

Fission gas release

The fission gas release against the fuel burnup is depicted in Fig. 8. Maximum FGR at the end of calculations was approximately 3.52%, which provides a good comparison with that in the validation case (2.5% in VVER-1000 fuel rod) [44]. This result is predictable because the fission gas release process is influenced by fuel temperature. In the case of an uncracked pellet, the model predicts the release of fission gases from the beginning of operation because tension causes extremely low grain boundary saturation near the fuel surface.

Fig. 8

Calculated fission gas release evolution

5.4

Gap/Plenum pressure

Fig. 9 illustrates the gap/plenum pressure evolution versus fuel burnup. It shows a high gap/plenum pressure because of the poor thermal conductivity of the gap/plenum, which increases the fuel temperature and results in additional fission gas release. The maximum gas pressure in the current state was approximately 13.8 MPa and less than the coolant pressure of 16.2 MPa. According to A. Medvedev et al., “the minimum margin for the gas pressure limit; to ensure that it is not exceeding coolant pressure; is 1.5”. Thus, the failure of pellet-cladding heat transfer as a result of “lift-off” effect is eliminated [46].

Fig. 9

Calculated gap/plenum pressure evolution

Conclusion

A fuel performance model was constructed using the COMSOL Multiphysics platform. The modeling was performed for a 2D axis-symmetric geometry of the UO2 fuel pellet with E110 (Zr-1%Nb) cladding for VVER-1200 fuel. The model considers all relevant phenomena, including heat generation and conduction, gap heat transfer, elastic strain, mechanical contact, thermal expansion, grain growth, densification, fission gas generation and release, fission product swelling, gap/plenum pressure, and cladding thermal and irradiation creep. The model was validated using a code-to-code comparison of fuel surface and centerline temperatures for stable power and mechanical strain in addition to validation of FGR predictions. A sensitivity study was also conducted to assess the effects of uncertainty on some of the model parameters. The model was then used to predict other fuel performance parameters as a function of burnup, such as temperature profiles, gap width, fission gas release, and plenum pressure.

The use of multiphysics modules enables the proper inclusion of irradiation models that describe the irradiation behavior of the fuel performance. In addition, a compilation of related material and thermomechanical models was conducted and included in the model to allow the user to investigate different materials/performance models to investigate different phenomena and compare the results with available experimental measurements. This capacity-building step is necessary for operating and regulatory organizations. Although the model was developed for normal operating conditions, it can be modified to include off-normal operating conditions. Future improvements of the model would include additional detailed treatment of the pellet-clad heat transfer and mechanical interaction and enhanced models, such as fission gas release, pellet crack propagation, and fuel irradiation creep.

References

[1]

IAEA, Quality and Reliability Aspects in Nuclear Power Reactor Fuel Engineering. Vienna, 2015.

[2]

P.R. Roy, D.N. Sah,

Irradiation behaviour of nuclear fuels

. Pramana-J. Phys 24, 397–421 (1985), doi: 10.1007/BF02894841