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

Multiphysics simulation of VVER-1200 fuel performance during normal operating conditions 1 corresp first-author Khaled M. Yassin eng.khaled_m@yahoo.com Khaled M. Yassin eng.khaled_m@yahoo.com eng.khaled_m@yahoo.com https://orcid.org/0000-0002-8087-5186 2 corresp Mohamed H. Hassan mohamed.hassan@alexu.edu.eg Mohamed H. Hassan mohamed.hassan@alexu.edu.eg mohamed.hassan@alexu.edu.eg 1 corresp Mohammad M. Ghoneim mohghon1@yahoo.com Mohammad M. Ghoneim mohghon1@yahoo.com mohghon1@yahoo.com 1 3 corresp Mostafa S. Elkolil moustafa_elkoliel@yahoo.com Mostafa S. Elkolil moustafa_elkoliel@yahoo.com moustafa_elkoliel@yahoo.com 1 corresp Adel Alyan adelalyan@yahoo.com Adel Alyan adelalyan@yahoo.com adelalyan@yahoo.com 2 corresp Said A. Agamy sa_agamy@yahoo.com Said A. Agamy sa_agamy@yahoo.com sa_agamy@yahoo.com Nuclear Research Center, Egyptian Atomic Energy Authority, Abo Zaabal, Cairo 13759, Egypt Nuclear Research Center, Egyptian Atomic Energy Authority, Abo Zaabal, Cairo 13759, Egypt Department of Nuclear and Radiological Engineering, Faculty of Engineering, Alexandria University, Alexandria 21544, Egypt Department of Nuclear and Radiological Engineering, Faculty of Engineering, Alexandria University, Alexandria 21544, Egypt Faculty of Engineering, Egyptian Russian University, Badr 11829, Egypt Faculty of Engineering, Egyptian Russian University, Badr 11829, Egypt 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 2 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-1200 Fuel performance COMSOL code Zr-1%Nb cladding UO 2 fuel rod Introduction 1 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 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 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 3 ], FEMAXI [ 4 4 ], FALCON [ 5 5 ], FRAPCON[ 6 6 ], and FRAPTRAN[ 7 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 8 ]. For VVER Fuel performance modeling, some western PWRs fuel performance codes have been modified to predict fuel performance, including TRANSURANUS [ 9 9 ], FRAPCON [ 6 6 ], and PIN-Micro (based on GAPCON) [ 3 3 ]. Additionally, START-3, a Russian-developed code, was used and validated after feedback modifications from the FUMEX program [ 10 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 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 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 13 ] and CANDU reactors [ 8 8 , 14 14 ]. Prudil et al. [ 8 8 , 14 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 1 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. Modeling geometry 2 The fuel rod modeled in this research has 10 individual uncracked UO 2 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 15 ]. Figure 1 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 8 , 14 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 8 , 12 12 , 14 14 , 16 16 , 17 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. 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 18 - 21 21 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 22 ]- VVER-1200 (V-491) and 65 MWd/kgU in [ 19 19 ], whereas in [ 20 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 19 , 20 20 ]. Materials properties library 2 Table 1 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 23 - 25 25 23-25 ]. 1 1  Property 1 1 UO 2 models 1 1 Zr-1%Nb (E110) 1 1 Thermal conductivity (W/mK) 1 1 Function of temperature, stoichiometry, porosity, fission product buildup, radiation damage [ 23 23 ] 1 1 Function of temperature [ 24 24 ] 1 1 Density (kg/m 3 ) 1 1 Function of initial density, material strain [ 26 26 ] 1 1 Function of initial density, material strain [ 25 25 ] 1 1 Heat capacity(J/kg <math id="M1"><mo>⋅</mo></math> K) 1 1 Function of temperature [ 24 24 ] 1 1 Function of temperature, heating rate [ 25 25 ] 1 1 Thermal expansion 1 1 Function of temperature [ 27 27 ] 1 1 Function of temperature [ 25 25 ] 1 1 Young’s modulus (N/m 2 ) 1 1 Function of temperature [ 23 23 ] 1 1 Function of temperature [ 25 25 ] 1 1 Poisson ratio 1 1 0.316 [ 23 23 ] 1 1 Function of temperature [ 25 25 ] 1 1 Emissivity 1 1 Function of temperature [ 23 23 ] 1 1 Function of temperature [ 23 23 ] Heat production and conductance 2 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: <math display="block" id="M2"><mrow><mi>ρ</mi><msub><mi>C</mi><mtext>P</mtext></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>∇</mo><mo>.</mo><mrow><mo>(</mo><mrow><mi>K</mi><mo>∇</mo><mi>T</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>Q</mi><mrow><mtext>prod</mtext></mrow></msub><mn>,</mn></mrow></math> 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 28 , 29 29 ] <math display="block" id="M3"><mrow><mi>Q</mi><mo>=</mo><mfrac><mrow><msub><mi>P</mi><mrow><mtext>lin</mtext></mrow></msub></mrow><mrow><mi>π</mi><mo>.</mo><msubsup><mi>a</mi><mrow><mtext>pellet</mtext></mrow><mn>2</mn></msubsup></mrow></mfrac><mn>,</mn></mrow></math> 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: <math display="block" id="M4"><mrow><msub><mi>Q</mi><mrow><mtext>r</mtext><mo> </mo></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>h</mi><mrow><mtext>gas</mtext></mrow></msub><mo>+</mo><msub><mi>h</mi><mrow><mtext>solid</mtext></mrow></msub><mo>+</mo><msub><mi>h</mi><mrow><mtext>rad</mtext></mrow></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>fuel</mtext></mrow></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>Clad</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo> </mo><mn>,</mn></mrow></math> 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 8 , 12 12 ]. Gas conduction h gas is represented by the formula suggested by Ross and Stoute [ 30 30 ]: <math display="block" id="M5"><mrow><msub><mi>h</mi><mrow><mtext>gas</mtext></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>k</mi><mtext>g</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>gas</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>C</mi><mtext>r</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>R</mi><mtext>f</mtext></msub><mo>+</mo><msub><mi>R</mi><mtext>c</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>d</mi><mrow><mtext>gap</mtext></mrow></msub><mo>+</mo><msub><mi>g</mi><mn>1</mn></msub><mo>+</mo><msub><mi>g</mi><mn>2</mn></msub></mrow></mfrac><mo>,</mo></mrow></math> where k gas is the gas conductivity in the gap and is represented by the formula used in [ 17 17 , 31 31 ]; d gap refers to the gap thickness calculated dynamically in the solid mechanics module; C r refers to the surface roughness coefficient; r 1 and r 2 are the roughness of the fuel and clad surfaces, respectively; and g 1 and g 2 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 32 ]: <math display="block" id="M6"><mrow><msub><mi>h</mi><mrow><mtext>solid</mtext></mrow></msub><mo>=</mo><msub><mi>C</mi><mtext>s</mtext></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mn>2</mn><msub><mi>k</mi><mn>1</mn></msub><msub><mi>k</mi><mn>2</mn></msub></mrow><mrow><msub><mi>k</mi><mtext>f</mtext></msub><mo>+</mo><mo> </mo><msub><mi>k</mi><mtext>s</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>P</mi><mtext>c</mtext></msub></mrow><mrow><msup><mi>δ</mi><mrow><mn>0.5</mn></mrow></msup><mi>H</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where C s is a constant (= 1.0), k 1 and k 2 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 <math display="block" id="M7"><mrow><msub><mi>q</mi><mrow><mtext>rad</mtext></mrow></msub><mo>=</mo><mfrac><mtext>σ</mtext><mrow><mfrac><mn>1</mn><mrow><msub><mtext>ε</mtext><mn>1</mn></msub></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mtext>ε</mtext><mn>2</mn></msub></mrow></mfrac><mo>−</mo><mn>1</mn></mrow></mfrac><mrow><mo>(</mo><mrow><msubsup><mi>T</mi><mn>1</mn><mn>4</mn></msubsup><mo>−</mo><msubsup><mi>T</mi><mn>2</mn><mn>4</mn></msubsup></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>h</mi><mrow><mtext>rad</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mn>1</mn></msub><mo>−</mo><msub><mi>T</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where ε 1 and ε 2 are the emissivities of the radiating surfaces; σ is the Stephan-Boltzmann constant; and T 1 and T 2 are the fuel and cladding surface temperatures, respectively. Radiative conduction can be expressed as <math display="block" id="M8"><mrow><msub><mi>h</mi><mrow><mtext>gab</mtext><mo>,</mo><mtext>rad</mtext></mrow></msub><mo>=</mo><mfrac><mi>σ</mi><mrow><mfrac><mn>1</mn><mrow><msub><mi>ε</mi><mrow><mtext>e</mtext><mo>,</mo><mtext>f</mtext></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>ε</mi><mrow><mtext>e</mtext><mo>,</mo><mtext>c</mtext></mrow></msub></mrow></mfrac><mo>−</mo><mn>1</mn></mrow></mfrac><mrow><mo>(</mo><mrow><msubsup><mi>T</mi><mtext>f</mtext><mn>2</mn></msubsup><mo>+</mo><msubsup><mi>T</mi><mtext>C</mtext><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>f</mtext></msub><mo>+</mo><msub><mi>T</mi><mtext>C</mtext></msub></mrow><mo>)</mo></mrow></mrow><mo> </mo></msup><mo>.</mo></mrow></math> 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 <math display="block" id="M9"><mrow><msub><mi>Q</mi><mtext>r</mtext></msub><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><msubsup><mi>T</mi><mrow><mtext>clad</mtext></mrow><mtext> </mtext></msubsup><mo>−</mo><msubsup><mi>T</mi><mrow><mtext>Coolant</mtext></mrow><mtext> </mtext></msubsup></mrow><mo>)</mo></mrow></mrow><mrow><mfrac><mrow><msubsup><mi>t</mi><mrow><mtext>oxide</mtext></mrow><mtext> </mtext></msubsup></mrow><mrow><msubsup><mi>k</mi><mrow><mtext>oxide</mtext></mrow><mtext> </mtext></msubsup></mrow></mfrac><mtext> </mtext><mo>+</mo><mtext> </mtext><mfrac><mn>1</mn><mrow><msubsup><mi>h</mi><mrow><mtext>film</mtext></mrow><mtext> </mtext></msubsup></mrow></mfrac></mrow></mfrac><mo>−</mo><msub><mi>Q</mi><mrow><mtext>oxide</mtext></mrow></msub><mo>,</mo></mrow></math> where T coolant and T cladding are the coolant and cladding temperatures, respectively. k oxide is the therma1 conductivity of the zirconium-oxide (ZrO 2 ) 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 33 ]: <math display="block" id="M10"><mrow><msub><mi>Q</mi><mrow><mtext>oxide</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mn>42.426</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mn>9</mn></msup></mrow><mrow><mn>1.56</mn></mrow></mfrac><mfrac><mrow><mtext>d</mtext><msub><mi>t</mi><mrow><mtext>oxide</mtext></mrow></msub></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>,</mo></mrow></math> where the ZrO 2 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 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 23 ] . In the pre-transition stage, the corrosion rate was suggested by the Arrhenius correlation: <math display="block" id="M11"><mrow><mfrac><mrow><mtext>d</mtext><msub><mi>t</mi><mrow><mtext>oxide</mtext></mrow></msub></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>C</mi><mn>1</mn></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><msub><mi>Q</mi><mn>1</mn></msub></mrow><mrow><mi>R</mi><msub><mi>T</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo><mtext> for </mtext><msub><mi>t</mi><mrow><mtext>oxide</mtext></mrow></msub><mo>≤</mo><msub><mi>t</mi><mrow><mtext>trans</mtext></mrow></msub><mo>,</mo></mrow></math> where Q 1 and C 1 are the activation energy in the pre-transition oxidation stage and rate constant, respectively. From [ 24 24 ], the values of activation energy Q 1 =14680 and rate constant C 1 = 5.19×10 -7 for the E110 alloy were obtained. For the post-transition period, the oxidation growth rate is given by <math display="block" id="M12"><mrow><mfrac><mrow><mtext>d</mtext><msub><mi>t</mi><mrow><mtext>oxide</mtext></mrow></msub></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>C</mi><mn>2</mn></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><msub><mi>Q</mi><mn>2</mn></msub></mrow><mrow><mi>R</mi><msub><mi>T</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo><mtext> for </mtext><msub><mi>t</mi><mrow><mtext>oxide</mtext></mrow></msub><mo>≥</mo><msub><mi>t</mi><mrow><mtext>trans</mtext></mrow></msub><mo>.</mo></mrow></math> In addition, from [ 24 24 ], the values of activation energy Q 2 =15355 and rate constant C 2 = 17.72×10 -7 for the E110 alloy were obtained. Deformation mechanics 2 Elastic deformation 3 In a nuclear reactor, the fuel element deformation is described using Cauchy's equation [ 8 8 , 14 14 , 16 16 ], which is included in the structural mechanics module in COMSOL in the following form: <math display="block" id="M13"><mrow><mo>−</mo><mo>∇</mo><mo>⋅</mo><mi>σ</mi><mo>=</mo><msub><mi>F</mi><mtext>v</mtext></msub><mo>,</mo></mrow></math> 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: <math display="block" id="M14"><mrow><mi>σ</mi><mo>=</mo><mrow><mo>[</mo><mi>C</mi><mo>]</mo></mrow><mrow><mo>[</mo><mi>ε</mi><mo>]</mo></mrow><mo>,</mo></mrow></math> 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 <math display="block" id="M15"><mrow><msub><mi>ε</mi><mi>r</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac><mo>,</mo><msub><mi>ε</mi><mi>z</mi></msub><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo> </mo><mo>,</mo><msub><mi>ε</mi><mrow><mi>r</mi><mi>z</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mtext> </mtext><mo>,</mo><msub><mi>ε</mi><mrow><mi>r</mi><mi>θ</mi></mrow></msub><mo>=</mo><msub><mi>ε</mi><mrow><mi>θ</mi><mi>z</mi></mrow></msub><mo>=</mo><mn>0.</mn></mrow></math> In this study, a modified penalty method used in the FAST code [ 8 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 8 , 14 14 ] . <math display="block" id="M16"><mrow><msubsup><mi>P</mi><mi>i</mi><mo> </mo></msubsup><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>P</mi><mrow><mtext>est</mtext></mrow></msub><mi>exp</mi><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>P</mi><mtext>f</mtext></msub><msub><mi>d</mi><mrow><mtext>pen</mtext></mrow></msub></mrow><mrow><msub><mi>P</mi><mrow><mtext>est</mtext></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>d</mi><mrow><mtext>pen</mtext></mrow></msub><mo><</mo><mn>0</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>P</mi><mrow><mtext>est</mtext></mrow></msub><mo>+</mo><msub><mi>P</mi><mtext>f</mtext></msub><msub><mi>d</mi><mrow><mtext>pen</mtext></mrow></msub></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>d</mi><mrow><mtext>pen</mtext></mrow></msub><mo>≥</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> 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 13 Pa/m were applied for the pellet-clad mechanical contact [ 16 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 13 ]. Thermal strain/expansion 3 Fink’s thermal expansion model [ 27 27 ] was used. It assumed an isotropic UO 2 fuel pellet. The thermal strain ε th is computed as <math display="block" id="M17"><mrow><msub><mtext>ε</mtext><mrow><mtext>th</mtext></mrow></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>0.9973</mn><mo>+</mo><mn>9.082</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>6</mn></mrow></msup><mi>T</mi><mo>−</mo><mn>2.705</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>10</mn></mrow></msup><msup><mi>T</mi><mn>2</mn></msup><mo>+</mo><mn>4.391</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>13</mn></mrow></msup><msup><mi>T</mi><mn>3</mn></msup><mo>,</mo><mo> </mo><mtext> </mtext><mn>273</mn><mo> </mo><mi>K</mi><mo><</mo><mi>T</mi><mo>≤</mo><mn>923</mn><mo> </mo><mi>K</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>0.9672</mn><mo>+</mo><mn>1.179</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><mi>T</mi><mo>−</mo><mn>2.429</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>9</mn></mrow></msup><msup><mi>T</mi><mn>2</mn></msup><mo>+</mo><mn>1.219</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>12</mn></mrow></msup><msup><mi>T</mi><mn>3</mn></msup><mo>,</mo><mtext> </mtext><mn>923</mn><mo> </mo><mi>K</mi><mo><</mo><mi>T</mi><mo>≤</mo><mn>3120</mn><mo> </mo><mi>K</mi></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> using the assumption that the strain is zero at the reference temperature of 273 K. Fuel densification and fission products’ swelling 3 Fuel densification can be described by employing the empirical correlation developed by Hastings and Evans[ 34 34 ]: <math display="block" id="M18"><mrow><msub><mi>F</mi><mrow><mtext>poreR</mtext></mrow></msub><mo>=</mo><mn>0.6</mn><mo>−</mo><mtext>exp</mtext><mrow><mo>(</mo><mrow><mo>−</mo><mn>0.506</mn><mo>−</mo><mn>8.67</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>10</mn></mrow></msup><msup><mi>T</mi><mn>3</mn></msup><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>exp</mi><mrow><mo>(</mo><mrow><mo>−</mo><mn>2.867</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>2</mn></mrow></msup><mi>B</mi><mi>u</mi></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> <math display="block" id="M19"><mrow><msub><mtext>ε</mtext><mrow><mtext>vol</mtext><mo>,</mo><mtext>dens</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mtext>Δ</mtext><msub><mi>V</mi><mrow><mtext>dens</mtext></mrow></msub></mrow><mrow><msub><mi>V</mi><mn>0</mn></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mn>1</mn><mo>−</mo><msub><mi>P</mi><mtext>o</mtext></msub></mrow><mrow><mn>1</mn><mo>−</mo><msub><mi>P</mi><mtext>o</mtext></msub><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msub><mi>F</mi><mrow><mtext>dens</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></math> where F dens is the fuel porosity fraction removed by the densification process, and p 0 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 23 ]. Solid fission product swelling results in a volumetric strain, computed as a linear function of burnup: <math display="block" id="M20"><mrow><mfrac><mrow><mi>δ</mi><msub><mi>ε</mi><mrow><mtext>sw</mtext><mo>−</mo><mtext>s</mtext></mrow></msub></mrow><mrow><mi>δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mn>5.577</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><mtext>ρ</mtext><mfrac><mrow><mtext>δ</mtext><mi>B</mi><mi>u</mi></mrow><mrow><mtext>δ</mtext><mi>t</mi></mrow></mfrac></mrow></math> where ε sw-s refers to the solid fission product swelling, ρ is the density (kg/m 3 ), 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 23 ]: <math display="block" id="M21"><mrow><mfrac><mrow><mi>δ</mi><msub><mi>ε</mi><mrow><mtext>sw</mtext><mo>−</mo><mtext>g</mtext></mrow></msub></mrow><mrow><mi>δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mn>1.96</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>31</mn></mrow></msup><mi>ρ</mi><mfrac><mrow><mi>δ</mi><mtext>Buf</mtext></mrow><mrow><mi>δ</mi><mi>t</mi></mrow></mfrac><msup><mrow><mrow><mo>(</mo><mrow><mn>2800</mn><mo>−</mo><mi>T</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>−</mo><mn>11.73</mn></mrow></msup><msup><mi>e</mi><mrow><mrow><mo>[</mo><mrow><mo>−</mo><mn>0.0162</mn><mrow><mo>(</mo><mrow><mn>2800</mn><mo>−</mo><mi>T</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow></mrow></msup><msup><mrow/><mo> </mo></msup><msup><mi>e</mi><mrow><mo>−</mo><mn>0.0178</mn><mi>ρ</mi><mi>B</mi><mi>u</mi></mrow></msup></mrow></math> where ε sw-g is the volumetric gas swelling, and B u and δ B u are the burnup and burnup increments (fissions/atoms-U), respectively. Clad thermal and irradiation creep. 3 The creep rate correlation is provided in MATPRO-11[ 35 35 ]. This relation computes the thermal and irradiation creeps in the radial direction. ε ϕ is the thermal creep rate, which is calculated as <math display="block" id="M22"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ε</mi><mo>∅</mo></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mn>5</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>23</mn></mrow></msup><msubsup><mi>σ</mi><mi>ϕ</mi><mn>2</mn></msubsup><mo> </mo><mrow><mo>(</mo><mrow><mn>3.47</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>23</mn></mrow></msup><mfrac><mrow><msubsup><mi>σ</mi><mi>ϕ</mi><mn>3</mn></msubsup></mrow><mrow><mrow><mo>|</mo><mrow><msub><mi>σ</mi><mi>ϕ</mi></msub></mrow><mo>|</mo></mrow></mrow></mfrac><mi>exp</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mi>U</mi></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><msub><mi>ε</mi><mi>ϕ</mi></msub></mrow><mo>)</mo></mrow><mi>exp</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mi>U</mi></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> 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. <math display="block" id="M23"><mrow><mtext>U</mtext><mo>=</mo><mtext> </mtext><mn>2</mn><mtext>l</mtext><mn>2.7</mn><mo>−</mo><mn>0.5324</mn><mtext> </mtext><mi>T</mi><mo>+</mo><mn>1.17889</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>4</mn></mrow></msup><msup><mi>T</mi><mn>2</mn></msup><mo>+</mo><mn>3.3486</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>7</mn></mrow></msup><msup><mi>T</mi><mn>3</mn></msup></mrow></math> 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 36 ]. <math display="block" id="M24"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ε</mi><mo>∅</mo></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mn>2.2</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>7</mn></mrow></msup><msub><mi>σ</mi><mo>∅</mo></msub><mo> </mo><mtext>exp</mtext><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mn>5000</mn></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mtext> </mtext></msup><msubsup><mo>∅</mo><mrow><mtext>Fflux</mtext></mrow><mrow><mo>−</mo><mn>0.65</mn></mrow></msubsup><msub><mrow/><mo> </mo></msub></mrow><mrow><msup><mi>T</mi><mn>7</mn></msup></mrow></mfrac><mo>,</mo></mrow></math> where ϕ Fflux is the flux of the fast neutrons (1/m 2· <math id="M25"><mo>⋅</mo></math> s) with a kinetic energy of more than 1 MeV. Fission gases production and release 2 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 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. Release of fission gases to the grain boundaries 3 The release of fission gases to fuel grain boundaries can be simulated using Booth’s diffusion model [ 38 38 , 39 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. <math display="block" id="M26"><mrow><mfrac><mrow><mo>∂</mo><mi>C</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mi>D</mi><msup><mo>∇</mo><mn>2</mn></msup><mi>C</mi><mo>+</mo><msub><mi>P</mi><mrow><mtext>fg</mtext></mrow></msub><mo>,</mo></mrow></math> 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 2 grain diameter was computed using the grain growth correlation given in [ 40 40 ]. The grain growth rate (m/s) was calculated as follows: <math display="block" id="M27"><mrow><mfrac><mrow><mtext>d</mtext><msub><mi>g</mi><mtext>d</mtext></msub><msup><mrow/><mo> </mo></msup></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>k</mi><mtext>g</mtext></msub><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mrow><msub><mi>g</mi><mtext>d</mtext></msub></mrow></mfrac><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>g</mi><mrow><mtext>max</mtext></mrow></msub></mrow></mfrac><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>g</mi><mrow><mtext>ir</mtext></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where g d is the average grain diameter (2 g r ), and k g (m 2 /s) is the rate constant, which is computed as follows: <math display="block" id="M28"><mrow><msub><mi>k</mi><mtext>g</mtext></msub><mo>=</mo><mn>1.46</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>8</mn></mrow></msup><mtext>exp</mtext><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mn>321001</mn></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where T is the temperature (K), and g max is the restricting grain size (m), which is a function of temperature given by <math display="block" id="M29"><mrow><msub><mi>g</mi><mrow><mtext>max</mtext></mrow></msub><mo>=</mo><mn>2.23</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>3</mn></mrow></msup><mtext>exp</mtext><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mn>7620</mn></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> and g ir accounts for the irradiation effects on grains, given by <math display="block" id="M30"><mrow><msub><mi>g</mi><mrow><mtext>ir</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mn>6.71</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mn>18</mn></mrow></msup><mtext>exp</mtext><mrow><mo>(</mo><mrow><mfrac><mrow><mo>−</mo><mn>5620</mn></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><msub><mi>F</mi><mrow><mtext>rate</mtext></mrow></msub><mi>T</mi></mrow></mfrac><mo>,</mo></mrow></math> F rate refers to the fission rate density (m −3 ). 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 40 ]. Fission gas release to the free volume of the fuel element 3 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 3 ) was computed as described by White and Tucker [ 39 39 ]: <math display="block" id="M31"><mrow><msub><mi>G</mi><mrow><mtext>bsat</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mn>4</mn><msub><mi>r</mi><mtext>f</mtext></msub><mi>f</mi><mrow><mo>(</mo><mrow><msub><mi>θ</mi><mrow><mtext>fg</mtext></mrow></msub></mrow><mo>)</mo></mrow><msub><mi>f</mi><mtext>B</mtext></msub></mrow><mrow><mn>3</mn><msub><mi>k</mi><mtext>B</mtext></msub><mi>T</mi><msup><mrow><mtext>sin</mtext></mrow><mn>2</mn></msup><mrow><mo>(</mo><mrow><msub><mi>θ</mi><mrow><mtext>fg</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mtext>P</mtext><mrow><mtext>ext</mtext></mrow></msub><mo>+</mo><mfrac><mrow><mn>2</mn><msub><mtext>γ</mtext><mrow><mtext>se</mtext></mrow></msub></mrow><mrow><msub><mi>r</mi><mtext>f</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mfrac><mn>3</mn><mrow><msub><mi>g</mi><mtext>r</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M32"><mrow><mi>f</mi><mrow><mo>(</mo><mrow><msub><mi>θ</mi><mo> </mo></msub></mrow><mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>−</mo><mn>1.5</mn><mi>cos</mi><mrow><mo>(</mo><mtext>θ</mtext><mo>)</mo></mrow><mo>+</mo><mn>0.5</mn><msup><mrow><mtext>cos</mtext></mrow><mn>3</mn></msup><mrow><mo>(</mo><mtext>θ</mtext><mo>)</mo></mrow><mo>,</mo></mrow></math> where r f = 10 -7 m is the radius of curvature of the fission gas bubbles, f ( θ fg ) 1 is 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 <math id="M33"><mo>×</mo></math> 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 <math display="block" id="M34"><mrow><msub><mi>G</mi><mrow><mtext>bsat</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mn>2.0811</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mn>16</mn></mrow></msup></mrow><mrow><mi>T</mi><msub><mi>g</mi><mtext>r</mtext></msub></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>P</mi><mrow><mtext>ext</mtext></mrow></msub><mo>+</mo><mn>2.504</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mn>6</mn></msup></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> 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 13 ]. Model validation 1 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. Thermal and Mechanical Performance 2 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 41 ]. FINIX was integrated into Serpent 2 in [ 42 42 ] to study a VVER-1000 fuel rod. The geometry, material properties, and fuel rod specifications were obtained from the UAM benchmark [ 43 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 2 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 Figures 2 and 3 3 compare the results obtained using FINIX and our COMSOL model, illustrating satisfactory agreement between them. In Figure.2 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 42 ]. Figure 3 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. Fission Gas Release 2 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 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 Fig. 4 , the results of COMSOL show good agreement with the other codes and are closer to the final measurement values. Sensitivity analysis 1 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 2 thermal conductivity, UO 2 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. Sensitivity of FGR to Model Input Parameters 2 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 2 . 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 2 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. Sensitivity of FGR to Material Parameters 2 The FGR model was sensitive to the thermal conductivity of UO 2 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 2 thermal expansion. Sensitivity of Mid-Pellet Cladding Strain to Material Parameters 2 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. Sensitivity of Cladding Strain at Circumferential Ridge to Material Parameters 2 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. Sensitivity of Results to meshing density 2 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 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 1 Temperature Profiles 2 Figure 5 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. 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 2 with the rise in temperature of the fuel as fuel burnup continues [ 45 45 ] ( Fig. 6 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. 7 ). Gap width 2 The gap width is a result of the combined effects of swelling, fuel densification, thermal expansion, and clad creeping. As shown in Fig. 7 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. Fission gas release 2 The fission gas release against the fuel burnup is depicted in Fig. 8 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 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. Gap / Plenum pressure 2 Fig. 9 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 46 ]. Conclusion 1 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 IAEA , Quality and Reliability Aspects in Nuclear Power Reactor Fuel Engineering . Vienna , 2015 . Quality and Reliability Aspects in Nuclear Power Reactor Fuel Engineering P.R. Roy , D.N. Sah , Irradiation behaviour of nuclear fuels . Pramana-J. Phys  24 , 397 – 421 ( 1985 ), doi: 10.1007/BF02894841 http://doi.org/10.1007/BF02894841 Irradiation behaviour of nuclear fuels C.R. Hann , C.E. Beyer , L. J. Parchen et al. , GAPCON-THERMAL-1: a computer program for calculating the gap conductance in oxide fuel pins . 1973 , Battelle Pacific Northwest Labs . https://digital.library.unt.edu/ark:/67531/metadc1019317/ https://digital.library.unt.edu/ark:/67531/metadc1019317/ GAPCON-THERMAL-1: a computer program for calculating the gap conductance in oxide fuel pins T. Nakajima , H. Saito , T. Osaka. , FEMAXI-IV: a computer code for the analysis of thermal and mechanical behavior of light water reactor fuel rods . Nucl. Eng. Des. 148 , 41 - 52 , ( 1994 ). doi: 10.1016/0029-5493(94)90240-2 http://doi.org/10.1016/0029-5493(94)90240-2 . FEMAXI-IV: a computer code for the analysis of thermal and mechanical behavior of light water reactor fuel rods W.F. Lyon , M.N. Jahingir , R.O. Montgomery et al. , Fuel analysis and licensing code: FALCON MOD01: Volume 3: Verification and validation . EPRI , Palo Alto, CA 1011309 ( 2004 ). Fuel analysis and licensing code: FALCON MOD01: Volume 3: Verification and validation G.A. Berna , C.E. Beyer , K.L. Davis et al. , FRAPCON-3: A computer code for the calculation of steady-state, thermal-mechanical behavior of oxide fuel rods for high burnup . US NRC, Office of Nuclear Regulatory Research NUREGKR-6534 , Volume 2 Washington , 1997 . FRAPCON-3: A computer code for the calculation of steady-state, thermal-mechanical behavior of oxide fuel rods for high burnup  J.M. Cuta , C.E. Beyer , K.J. Geelhood et al. , FRAPTRAN 1.4: a computer code for the transient analysis of oxide fuel rods . US Nuclear Regulatory Commission, Office of Nuclear Regulatory Research, NUREG/CR-7023 , 2011 . FRAPTRAN 1.4: a computer code for the transient analysis of oxide fuel rods A. Prudil , B.J. Lewis , P.K. Chan et al. , Development and testing of the FAST fuel performance code: Normal operating conditions (Part 1) , Nucl. Eng. Des. 282 , 158 - 168 , ( 2015 ). doi: 10.1016/j.nucengdes.2014.09.036 http://doi.org/10.1016/j.nucengdes.2014.09.036 . Development and testing of the FAST fuel performance code: Normal operating conditions (Part 1) P.V. Uffelen , P. Bliar , S. Boneva et al. , The Application of the TRANSURANUS Fuel Performance Code to WWER Fuel: An Overview . 13 International Conference on WWER Fuel Performance, Modeling and Experimental Support ( Nesebar Bulgaria ; 15-21 Sep 2019 ). The Application of the TRANSURANUS Fuel Performance Code to WWER Fuel: An Overview S. Stefanova , I.G. Kolev , P. Chantoin et al. , VVER Reactor Fuel Performance , International conference on WWER fuel performance, modelling and experimental support Proceedings ., ( St Constantine, Varna, Bulgaria , 7-11 November 1994 ). VVER Reactor Fuel Performance R.L. Williamson , J.D. Hales , S.R. Novascone et al. , Multidimensional multiphysics simulation of nuclear fuel behavior . J. Nucl. Mater. 423 , 149 - 163 ( 2012 ). doi: 10.1016/j.jnucmat.2012.01.012 http://doi.org/10.1016/j.jnucmat.2012.01.012 . Multidimensional multiphysics simulation of nuclear fuel behavior R.L. Williamson , J.D. Hales , S.R. Novascone et al. , BISON theory manual the equations behind nuclear fuel analysis. , Idaho National Lab.(INL) , Idaho Falls , 2016 , doi: 10.2172/1374503 http://doi.org/10.2172/1374503 . BISON theory manual the equations behind nuclear fuel analysis. R. Liu , W. Zhou , P. Shen et al. , Fully coupled multiphysics modeling of enhanced thermal conductivity UO 2 –BeO fuel performance in a light water reactor . Nucl. Eng. Des. 295 , 511 - 523 , ( 2015 ). doi: 10.1016/j.nucengdes.2015.10.019 http://doi.org/10.1016/j.nucengdes.2015.10.019 . Fully coupled multiphysics modeling of enhanced thermal conductivity UO2–BeO fuel performance in a light water reactor A. Prudil , B.J. Lewis , P.K. Chan et al. , Development and testing of the FAST fuel performance code: Transient conditions (Part 2) . Nucl. Eng. Des. 282 , 169 - 177 ( 2015 ). doi: 10.1016/j.nucengdes.2014.11.036 http://doi.org/10.1016/j.nucengdes.2014.11.036 . Development and testing of the FAST fuel performance code: Transient conditions (Part 2) L.T. Yanko , Nuclear fuel for VVER-1200 . Rosatom Seminar on Russian Nuclear Energy Technologies and Solutions , April 2-3, 2012 . Nuclear fuel for VVER-1200 R. Liu , A. Prudil , W.Z. Zhou et al. , Multiphysics coupled modeling of light water reactor fuel performance . Prog. Nucl. Energy 91 , 38 - 48 ( 2016 ). doi: 10.1016/j.pnucene.2016.03.030 http://doi.org/10.1016/j.pnucene.2016.03.030 . Multiphysics coupled modeling of light water reactor fuel performance D. Morgan , Dissertation, Royal Military College of Canada , Kingston, Ontario, Canada , 2007 . A.A. Galahom , Improvement of the VVER-1200 fuel cycle by introducing thorium with different fissile material in blanket-seed assembly . Nucl. Sci. Eng. 193 , 638 - 651 ( 2019 ). doi: 10.1080/00295639.2018.1560757 http://doi.org/10.1080/00295639.2018.1560757 . Improvement of the VVER-1200 fuel cycle by introducing thorium with different fissile material in blanket-seed assembly V. Molchanov , Nuclear fuel for VVER reactors. Actual state and trends . 8th International Conference on VVER Fuel Performance, Modeling and Experimental Support 27.09–02.10.2009 , Helena Resort, Bulgaria , 2009 . Nuclear fuel for VVER reactors. Actual state and trends Y. Semchenkov , Y. Styrin , Advancing of VVER Reactor Core . BULATOM International Nuclear Forum on Nuclear Energy - challenges and prospects , ( Varna, Bulgaria , 9-11 June, 2010 ). Advancing of VVER Reactor Core L.T. Yanko , Nuclear fuel for VVER-1200 . Rosatom Seminar on Russian Nuclear Energy Technologies and Solutions , ( Johannesburg, South Africa , 2-3 Apr, 2012 ). Nuclear fuel for VVER-1200 IAEA . Status report 108 - VVER-1200 (V-491) (VVER-1200 (V-491)) . 2011 ; Available from: https://aris.iaea.org/PDF/VVER-1200(V-491).pdf https://aris.iaea.org/PDF/VVER-1200(V-491).pdf . Status report 108 - VVER-1200 (V-491) (VVER-1200 (V-491)) C.M. Allison , D.T. Hagrman , G.A. Berna et al. , SCDAP/RELAP5/MOD3. 1 Code Manual Volume IV: MATPRO . U.S. Nuclear Regulatory Commission, Office of Nuclear Regulatory Research, NUREG/CR-6150 , Volume 4 . Washington , 1995 . doi: 10.2172/100327 http://doi.org/10.2172/100327 . SCDAP/RELAP5/MOD3. 1 Code Manual Volume IV: MATPRO A. Shestopalov , K. Lioutov , L. Yegorova , Modification of USNRC's FRAP-T6 Fuel Rod Transient Code for High Burnup VVER Fuel . U.S. Nuclear Regulatory Commission, NUREG/IA-0164 Washington , 1999 . Modification of USNRC's FRAP-T6 Fuel Rod Transient Code for High Burnup VVER Fuel A. Shestopalov , K. Lioutov , L. Yegorova , Adaptation of USNRC's FRAPTRAN and IRSN's SCANAIR Transient Codes and Updating of MATPRO Package for Modeling of LOCA and RIA Validation Cases with Zr-1% Nb (VVER Type) Cladding . U.S. NRC, NUREG/IA-0209 , Washington , 2003 . Adaptation of USNRC's FRAPTRAN and IRSN's SCANAIR Transient Codes and Updating of MATPRO Package for Modeling of LOCA and RIA Validation Cases with Zr-1% Nb (VVER Type) Cladding IAEA , Design and Performance of WWER Fuel . Technical Reports Series No. 379 , Vienna ( 1996 ). Design and Performance of WWER Fuel J.K. Fink , Thermophysical properties of uranium dioxide . J. Nucl. Mater. 279 , 1 - 18 ( 2000 ). doi: 10.1016/S0022-3115(99)00273-1 http://doi.org/10.1016/S0022-3115(99)00273-1 . Thermophysical properties of uranium dioxide J.R. Lamarsh , A.J. Baratta , Introduction to nuclear engineering . 3 rd edn, Prentice hall Upper Saddle River , NJ , 2001 , pp. 409 - 413 . Introduction to nuclear engineering N.E. Todreas , M.S. Kazimi , Nuclear Systems 1 Thermal Hydraulic Fundamentals, MIT, Hemisphere Publishing Corporation , 1990 , pp. 53 - 57 . Nuclear Systems 1 Thermal Hydraulic Fundamentals, MIT, Hemisphere Publishing Corporation A.M. Ross , R.L. Stoute , Heat transfer coefficient between UO 2 and Zircaloy-2 , Atomic Energy of Canada Limited , Canada , 1962 . Heat transfer coefficient between UO2 and Zircaloy-2 K. Shaheen , Dissertation ( Royal Military College of Canada , Kingston, Ontario, Canada , 2011 ). D.R. Olander , Fundamental aspects of nuclear reactor fuel elements. TID-26711-Pl ( Atomic Energy Commission , Mumbai ), 1985 . doi: 10.2172/7343826 http://doi.org/10.2172/7343826 . Fundamental aspects of nuclear reactor fuel elements. TID-26711-Pl S. Leistikow , G. Schanz , Oxidation kinetics and related phenomena of Zircaloy-4 fuel cladding exposed to high temperature steam and hydrogen-steam mixtures under PWR accident conditions . Nucl. Eng. Des. 103 , 65 - 84 ( 1987 ). doi: 10.1016/0029-5493(87)90286-X http://doi.org/10.1016/0029-5493(87)90286-X . Oxidation kinetics and related phenomena of Zircaloy-4 fuel cladding exposed to high temperature steam and hydrogen-steam mixtures under PWR accident conditions I.J. Hastings , P.J. Fehrenbach , R.R. Hosbons , Densification in irradiated UO 2 fuel . J. Am. Ceram. Soc. 67 ( 2 ), C-24 - C-25 , ( 1984 ), doi: 10.1111/j.1151-2916.1984.tb09613.x http://doi.org/10.1111/j.1151-2916.1984.tb09613.x . Densification in irradiated UO2 fuel D.L. Hagrman , G.A. Reymann , MATPRO-Version 11: A Handbook of Materials Properties for Use in the Analysis of Light Water Reactor Fuel and Behavior . Idaho National Lab (INL) , United States, NUREG/CR-0497 , 1979 doi: 10.2172/6442256 http://doi.org/10.2172/6442256 . MATPRO-Version 11: A Handbook of Materials Properties for Use in the Analysis of Light Water Reactor Fuel and Behavior N.E. Hoppe , Engineering model for zircaloy creep and growth . in International Topical Meeting on LWR Fuel Performance , Avignon (France) , 1991 . Engineering model for zircaloy creep and growth K. Forsberg , A.R. Massih , Diffusion theory of fission gas migration in irradiated nuclear fuel UO 2 . J. Nuc. Maters. 135 , 140 - 148 ( 1985 ). doi: 10.1016/0022-3115(85)90071-6 http://doi.org/10.1016/0022-3115(85)90071-6 Diffusion theory of fission gas migration in irradiated nuclear fuel UO2 S.D. Beck , The diffusion of radioactive fission products from porous fuel elements . Battelle Memorial Inst. , Columbus, Ohio , 1960 . doi: 10.2172/4158094 http://doi.org/10.2172/4158094 . The diffusion of radioactive fission products from porous fuel elements R.J. White , M.O. Tucker , A new fission-gas release model . J. Nucl. Maters. 118 ( 1 ), 1 - 38 ( 1983 ). doi: 10.1016/0022-3115(83)90176-9 http://doi.org/10.1016/0022-3115(83)90176-9 A new fission-gas release model O.V. Khoruzhii , S.Yu. Kourtchatov , V.V. Likhanskii , New model of equiaxed grain growth in irradiated UO 2 . J. Nucl. Maters. 265 ( 1-2 ), 112 - 116 ( 1999 ). doi: 10.1016/S0022-3115(98)00632-1 http://doi.org/10.1016/S0022-3115(98)00632-1 New model of equiaxed grain growth in irradiated UO2 T. Ikonen , H. Loukusa , E. Syrjälahti et al. , Module for thermomechanical modeling of LWR fuel in multiphysics simulations . Ann. Nucl. Energy 84 , 111 - 121 ( 2015 ) doi: 10.1016/j.anucene.2014.11.004 http://doi.org/10.1016/j.anucene.2014.11.004 Module for thermomechanical modeling of LWR fuel in multiphysics simulations E. Syrjälahti , V. Valtavirta , J. Kättö et al. , Multiphysics simulations of fast transients in VVER-1000 and VVER-440 reactors . 11th International Conference on WWER Fuel Performance, Modeling and Experimental Support , ( Varna, Bulgaria . 26 Sep - 3 Oct 2015 ). Multiphysics simulations of fast transients in VVER-1000 and VVER-440 reactors K. Ivanov , M. Avramova , T. Blyth et al. , Benchmark for uncertainty analysis in modeling (UAM) for design, operation and safety analysis of LWRs . Specification and Support Data for the Core Cases (Phase II) Version 1 , ( OECD NEA/NSC/DOC 2012 ) https://inis.iaea.org/collection/NCLCollectionStore/_Public/45/026/45026304.pdf https://inis.iaea.org/collection/NCLCollectionStore/_Public/45/026/45026304.pdf Benchmark for uncertainty analysis in modeling (UAM) for design, operation and safety analysis of LWRs G. Passage , A.S. Scheglov , V.N. Proselkov et al. , Comparative calculations and operation-to-PIE data juxtaposition of the Zaporozhye NPP, WWER-1000 FA-E0325 fuel rods after 4 years of operation up to 49 MWd/kgU burnup , 6 International conference on WWER fuel performance, modelling and experimental support ( Albena, Bulgaria , 19-23 Sep., 2005 ). https://inis.iaea.org/collection/NCLCollectionStore/_Public/37/098/37098340.pdf https://inis.iaea.org/collection/NCLCollectionStore/_Public/37/098/37098340.pdf Comparative calculations and operation-to-PIE data juxtaposition of the Zaporozhye NPP, WWER-1000 FA-E0325 fuel rods after 4 years of operation up to 49 MWd/kgU burnup IAEA , Thermal Conductivity of Uranium Dioxide . TECHNICAL REPORTS SERIES , Vienna , 1966 . Thermal Conductivity of Uranium Dioxide A. Medvedev , S. Bogatyr , V. Kouznetsov et al. , Fuel rod behaviour at high burnup WWER fuel cycles . In: Proceedings of the Fourth International Conference WWER Fuel Modelling and Experimental Support , ( Varna, Bulgaria , 29 Sep. - 3 Oct. 2003 ). Fuel rod behaviour at high burnup WWER fuel cycles All authors contributed to the study conception and design. Material preparation and data analysis were performed by Khaled M. Yassin, Mohamed H. Hassan and Mohammad M. Ghoneim. The first draft of the manuscript was written by Khaled M. Yassin, Answers to reviewer comments are prepared by Khaled M. Yassin, revised by Mohamed H. Hassan, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. The authors would like to acknowledge the support of Dr. Mahmoud Badawy, Mechanical Engineering Department, Alexandria University, for his support in using the COMSOL Multiphysics platform. 10.1007/s41365-023-01172-9 Funding Open access funding provided by The Science, Technology & Innovation Funding Authority (STDF) in cooperation with The Egyptian Knowledge Bank (EKB). 2022-08-31 2022-11-04 2022-11-22 2023-02-25 2023-03-20 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-023-01172-9 http://dx.doi.org/10.1007/s41365-023-01172-9 . This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4.0/. 2023 Nuclear Science and Techniques 2 34