© T. Ikonen et al., published by EDP Sciences, 2016

1 Introduction

In light water reactors, the thermal and mechanical behaviour of the fuel rods strongly influences the behaviour of the reactor in both steady state and transient conditions. For example, the power of the reactor is sharply affected by the fuel temperature due to the absorption of neutrons by Doppler-broadened cross sections. This coupling is important both in the steady state and, even more so, in transients. Similarly, transient heat transfer to the coolant and avoiding departure from nucleate boiling is dependent on the heat conductance of the pellet-cladding gap. The gap conductance is a notoriously complicated function of both thermal and mechanical properties of the fuel rod. Therefore dedicated fuel behaviour codes are often used in multiphysics simulations to solve the heat transfer in the rod self-consistently with, for example, the reactor power calculated by a neutronics code.

However, the fuel performance code or the expertise for its use is not always available for the multiphysics modeller. Therefore a somewhat simpler approach has been studied at VTT, where the fuel behaviour module FINIX has been recently developed. The FINIX module adopts a middle-ground between elaborate fuel performance codes and thermal elements in order to give a “simple but sufficient” description of the fuel rod's thermal and mechanical behaviour. The aim is to make the fuel behaviour modelling more easily approachable to multiphysics modellers without imposing on them the complete fuel behaviour phenomenology, but still provide significant improvements to stand-alone simulation methods.

As examples of multiphysics simulations, we use the FINIX fuel behaviour module integrated into Serpent 2 Monte Carlo reactor physics code and HEXTRAN reactor dynamics code. As a demonstration of the dynamical capabilities of the Serpent-FINIX system we simulate a prompt transient, where the pin power and fuel temperature are solved self-consistently, leading to termination of the transient due to Doppler-broadening of the cross sections. With HEXTRAN-FINIX, we simulate a control rod ejection accident in a full-core geometry, where we use FINIX to evaluate the effect of burnup on the outcome of the transient. We also compare the results to stand-alone HEXTRAN simulations.

2 Description of the FINIX module

2.1 Role in a multiphysics code package

FINIX is a fuel behaviour module designed to be integrated as a subprogram into a larger simulation code, where FINIX replaces the existing fuel model. The main design philosophy of FINIX is to provide a simple but sufficient model of the fuel rod that can be used in different types of simulation codes, including neutronics, reactor dynamics, thermal hydraulics, and system codes. While multiphysics capabilities can also be achieved with direct code-to-code coupling (see, e.g., Refs. [1,2]), such approach is often laborious when highly specialized software is involved. Thus, in the design of the FINIX code, flexibility of the interface between FINIX and the main simulation code has been prioritized. This facilitates integration into a wide range of different simulation codes.

FINIX is integrated into the main simulation code, the host code, at the source code level. This has the benefits of reduced data transfer between the codes, and allows the host code to have direct access to FINIX's functions and data structures. However, FINIX also has a high-level interface, through which the most common functionalities can be used without detailed knowledge of the FINIX data structures. In addition, to reduce the user's need for fuel-specific knowledge, FINIX has an internal database for different fuel types, from which the required fuel simulation parameters can be loaded in by just specifying the desired fuel rod type.

2.2 Overview of models

FINIX is primarily designed as a transient simulation code. Compared full-fledged fuel performance codes such as FRAPTRAN [3], TRANSURANUS [4] and BISON [5], FINIX emphasizes flexibility and the needs of the host code in favour of the fuel performance modelling. For example for neutronics simulations of fast transients, the most important quantity to solve is the time-dependent temperature distribution. On the other hand, many of the phenomena modelled in fuel performance codes such as cladding irradiation growth are less important in this case and are currently not considered in FINIX. However, similar to the work of reference [6] and for example the fuel model of RELAP [7], the coupling between the mechanical behaviour of the gap and the radial heat transfer is modelled. Here FINIX takes an approach similar to FRAPTRAN, although some simplifications are done as described below.

The FINIX model itself involves solving both the thermal and mechanical behaviour of the fuel rod, allowing not only thermal effects but also changes in rod geometry to be taken into account in the host code. The thermal and mechanical models are coupled by the gap pressure and conductance, which are functions of both the rod temperature and mechanical dimensions. Both the heat equation and the mechanical behaviour are solved radially in one dimensional, cylindrical and axisymmetric geometry, independently for several axial nodes. The solutions of the different axial nodes are coupled via the gap pressure, which is solved simultaneously for the whole rod. This scheme constitutes what is generally referred to as the 1.5-dimensional model. The main modules of FINIX and their interrelationships are shown in Figure 1. Properties such as thermal conductivity, thermal expansion, Young's moduli, coolant heat transfer, etc., are solved using publicly available correlations. For instance, the gap conductance is modelled similar to FRAPTRAN [3], taking into account the gas thermal properties and the gap thickness, and the heat transfer to coolant using Dittus-Boelter and Thom correlations (see, e.g., Ref. [11] for details). In coupled thermal hydraulic simulations, it is recommended to use the more sophisticated coolant heat transfer model of the host code.

Currently, FINIX is not equipped with modules to describe long-term evolution of phenomena such as cladding creep and oxidation, fuel swelling, and accumulation of fission products. FINIX thus lacks most of the models needed for simulating the effects of burnup accumulation over a long steady state irradiation. Therefore, for non-fresh fuel, the initial state of the fuel rod prior to the transient should be obtained by other means, for example from a FRAPCON [8] simulation.

For a comprehensive description of the FINIX module and its models, the interested reader is referred to references [9–11].

Fig. 1 The models of the FINIX module and its role in a multiphysics simulation. The iteration of the thermal and mechanical solutions is indicated by the flowchart. The convergence checks are assumed to automatically fail on the first iteration.

3 Serpent 2 simulation of VVER-1000 prompt transient

FINIX has been integrated into Serpent 2, a 3D continuous-energy Monte Carlo reactor physics burnup calculation code developed at VTT Technical Research Centre of Finland [12]. The development of Serpent 2 has a major focus on multi-physics applications. A universal multi-physics interface for code coupling is complemented with new methodology for the treatment of continuous temperature [13–15] and density [16] distributions. The recently implemented time dependent simulation mode [17] extends the applications of Serpent 2 even further. For further information on Serpent 2 see [18] and for the most recent multi-physics advances in Serpent 2 see [19].

3.1 Integration of FINIX into Serpent

The coupling of Serpent 2 and FINIX is done at the source code level. Serpent is responsible for solving the power distribution in the system while FINIX models the thermal and mechanical response of the fuel rod in transient or steady state conditions. The solution transfer between FINIX and the neutron transport part of Serpent is handled by a set of internal routines that form the fuel behaviour multi-physics interface in Serpent 2.

The fission power in fuel rod is tallied by Serpent and provided to FINIX nodes while conserving the total power generation as well as the local power generation in each node volume. The temperature solution calculated by FINIX can be used in Serpent as is with linear interpolation between the node points, without any mesh transformation. The changes in geometry obtained by FINIX can also be used in the neutron tracking without any transformation.

The multi-physics routines in Serpent 2 provide the neutron transport routines with the correct temperature and density distributions at different points in time and space so that the effect of the realistic temperature and density distributions can be accounted for in the interaction physics. The data transfer between Serpent and FINIX is done internally without disk operations.

For transient scenarios, before calculating the time-dependent solution, the steady state solution is obtained with the coupled system for the power level in question. This fuel behaviour solution is then used as the initial state of the fuel for the transient analysis. The steady state simulation is also used to create the initial neutron source for the transient simulation.

The time-dependent coupled solution is obtained by a sequential and iterative solving of the fission power distribution by Serpent and the temperature and strain distributions by FINIX. At the first iteration of a time step the temperature and strain distributions are taken from the end of the previous step and at the later iterations they are interpolated between the beginning of step and end of step distributions yielding a semi-implicit scheme. First the neutronics solution is obtained for the new time interval, after which the temperature distribution at the end of the interval is calculated by FINIX. A convergence criterion is applied after this and if the convergence of the coupled solution is not deemed sufficient, the neutronics solution for the next iteration can be obtained. The power distribution is relaxed over the different iterations.

3.2 Results and discussion

The simulation presented here is a prompt super-critical transient for a 2D VVER-1000 pin-cell. The geometry and material properties as well as the fuel rod specifications are taken from the UAM benchmark [20]. The goal of the simulation is to test and present the coupled calculation capabilities of the Serpent 2 – FINIX code system and should not be viewed as a realistic physical transient.

The pellet inner and outer radii were 0.070 and 0.378 cm, respectively. The cladding inner and outer radii were 0.386 and 0.455 cm. The lattice pitch for the hexagonal unit cell was 1.275 cm. The fuel was pure UO₂ with the enrichment of 3.3%. The cladding material was Zr + 1% Nb. The coolant temperature was set to 560 K.

The fuel rod was depleted until the average burnup of 10 MWd/kgU with 10 radial depletion zones to yield the radial burnup distribution shown in Figure 2a. The depletion calculation was done without fuel behaviour feedback. The realistic radial burnup distribution affects the radial power density distribution as well as the fuel thermal conductivity in the FINIX thermal model.

For the coupled calculation, the FINIX model of the rod used 101 equally spaced nodes in the pellet and 51 in the cladding. The temperature and strain distributions were brought into Serpent at these nodes. The power distribution was tallied in 20 radial zones with equal area. The cladding outer temperature was set to 570 K as a boundary condition for the thermal model of FINIX. In this simulation, this approximation may lead to slightly increased heat transfer from the fuel after pellet-cladding contact has been established. In general, setting a fixed coolant temperature boundary condition is preferable. This is also possible in FINIX [10], but as the current thermal hydraulic model lacks correlations for rapid transients, the simple approach of fixing the cladding surface temperature was used instead.

The power distribution tallied on time interval i was used to calculate the end of the interval (EOI) temperature distribution for interval i as well as an initial guess for the EOI temperature distribution for interval i + 1. A pointwise convergence criterion of 1 K was applied to the absolute difference between the temperature fields of subsequent iterations.

An initial steady state coupled calculation was conducted for the depleted system. The system was held critical at 233 W/cm linear power with soluble absorber in the coolant. The resulting radial power density distribution in steady state at 233 W/cm linear power is shown in Figure 2b. The critical steady state fuel behaviour solution produces the radial temperature and strain distributions shown in Figure 3. These were used as the initial conditions for the following time dependent simulation.

For the onset of the transient, the coolant boron concentration was decreased from 413 to 345 ppm. This corresponds to an instantaneous reactivity insertion of 2575 pcm (4.4 $) and makes the system prompt super-critical. After this initial modification, the system was allowed to evolve freely from the initial conditions for 30 ms. Delayed neutron emission was not included in the transient. The prompt super-criticality of the system and the short timescale of the transient mean that the effect of delayed neutrons would be negligible, however.

Figure 4a shows the evolution of the system linear power over time. The time behaviour is an initial exponential increase in system power until the increase in fuel temperature shuts down the transient. The development of the centreline and pellet surface temperatures can be seen in Figure 4b. The temperature at the pellet inner surface reaches its maximum value of 2782 K at 20.9 ms. This leaves a margin of approximately 325 K to the fuel melting temperature. The temperature of the pellet surface drops sharply after the thermal expansion of the pellet brings it into contact with the inner surface of the cladding greatly decreasing the thermal resistance of the gas gap. The total deposited energy can be integrated from the linear power giving 3.49 kJ/cm.

This calculation demonstrates the capabilities of the coupled code system to model a prompt power pulse starting from known initial conditions. The fission power and fuel behaviour are solved self-consistently and the code system captures the qualitative behaviour of the system very well. In a more realistic scenario, this simulation technique can solve both the size and shape of the power pulse and the values for interesting safety parameters such as pellet enthalpy and fuel and cladding maximum temperatures.

One of the main limitations of the Serpent – FINIX coupled code system at the moment is the omission of a delayed neutron emission model in the time dependent simulation mode of Serpent, limiting the applicability to fast transients. In longer transients, the coolant behaviour will also pay a greater role and has to be solved with a separate tool as has been done with e.g. the Serpent – SUBCHANFLOW coupled code system [21].

Fig. 2 The radial burnup distribution (a) and the resulting power distribution (b) for the VVER1000 pin-cell at 10 MWd/kgU average burnup calculated by Serpent.

Fig. 3 The steady state temperature (a) and strain (b) distributions in the VVER-1000 fuel rod.

Fig. 4 Evolution of system linear power (a) and rod centre line (green) and pellet surface (red) temperatures (b).

4 HEXTRAN simulation of VVER-440 rod ejection accident

4.1 Integration of FINIX into HEXTRAN

FINIX has also been integrated into the reactor dynamics code HEXTRAN [22] that is a coupled neutronics – thermal hydraulics code developed at VTT for transient and accident analysis of VVER reactors. HEXTRAN solves the two-group neutron diffusion equations with a nodal expansion method in a three-dimensional hexagonal fuel assembly geometry. Thermal-hydraulics of the reactor core is solved in separate one-dimensional hydraulic channels, which can be further divided into axial sub-regions. Usually each channel is coupled with one fuel assembly. Channel hydraulics is based on conservation equations for steam and water mass, total enthalpy and total momentum, and on a selection of optional correlations. During the hydraulics iterations, a one-dimensional heat transfer calculation is done for an average fuel rod of each assembly. HEXTRAN solves the heat transfer in a fuel rod on the basis of one-dimensional axially uncoupled equations using the theta method for time discretization. The heat capacity of the pellet and cladding and the conductivity of the cladding are given in input with temperature dependent correlations, whereas the heat conductivity of the pellet can also depend on burnup. Conductance of the gas gap is in practice always modelled using linear interpolation from a simple temperature dependent table. It is possible to define several fuel rod types in HEXTRAN and provide separate values for each fuel rod type. In practice, the lack of proper data diminishes the reliability and feasibility of this kind of approach.

In the new coupling with FINIX, HEXTRAN's own fuel heat transfer solution can be replaced with the FINIX module. The choice between FINIX and HEXTRAN's own models is done in the HEXTRAN input. HEXTRAN solves the power distribution in the reactor core and the distribution is transferred to the FINIX module. Also the bulk coolant temperature and the heat transfer coefficient between the cladding outer surface and the coolant are calculated in HEXTRAN and are used as boundary conditions for FINIX. The FINIX calculation is done for one average fuel rod of each fuel assembly, in total 150–500 fuel rods depending on the reactor type. HEXTRAN controls the data for each fuel rod and calls FINIX consecutively with the current state parameters of one fuel rod. FINIX returns the temperature distribution of the fuel rod for the neutronics calculation and other data that has to be stored for the next time step. Deformation of the fuel rod is taken into account in the heat transfer solution but at the moment deformation is ignored in the flow channel modelling. The same axial discretization is used in FINIX and in HEXTRAN. A typical number of axial levels is approximately 20. Radially typically 6 to 10 nodes are used for the pellet and one node for the cladding. The same time step is used in FINIX and HEXTRAN.

4.2 Results and discussion

Here the FINIX module's capability for reactor dynamical simulations is demonstrated with the recalculation of 3rd dynamic AER benchmark that was originally specified and calculated in 1994–1995 [23,24]. The benchmark concerns a control rod ejection in a VVER-440 reactor. Purpose of the benchmark was to test reactor dynamics codes capability to model asymmetric reactivity accident in hexagonal core geometry and to perform code-to-code comparisons. Typically in safety analyses of LWRs large amount of control rod ejection transients are analysed, starting from different initial states. In this benchmark the reactor is at the end of its first cycle. Here we study the differences in calculating such a benchmark case with and without explicit thermomechanical coupling in the fuel behaviour model.

In the scenario, the inlet temperature, inlet and outlet pressures of the core and coolant flow through the whole core are given and kept constant during the transient calculation. Coolant flow corresponds to conditions in which three of the six main coolant pumps are on. The initiating event of the benchmark is the ejection of the follower-type control rod from hot zero power state of the reactor when all seven rods of one rod bank are inserted 200 cm and all other rods are fully out of the core. Height of the core is 250 cm and the control rod is fully out from the core in 0.16 s. No reactor trip is modelled. One average fuel rod from each of the 349 fuel assemblies hasbeen modelled with FINIX. Fuel rods are divided axially to ten layers. Radially the fuel rod is modelled using five nodes for the pellet and one node for the cladding.

The reference case has been simulated using HEXTRAN's own heat transfer model using for the pellet and cladding temperature-dependent thermal conductivities and heat capacities that were defined in the benchmark specification. In these specifications, the thermal conductance was defined to depend on the average fuel temperature. Gap conductance has a constant value of 4 kW/K m² below 800 K, and increases linearly to 20 kW/K m² at 1500 K, above which the value is constant. Such an approach is often taken also in safety analyses of transients and accidents, which are typically simulated using very simple models for the gas gap. Radially uniform heat generation in the pellet is assumed, which is reasonable because the burnup is relatively low during the first cycle of the reactor.

The control rod ejection induces a strong power increase, which is cut off due to the Doppler effect. Fission power behaves similarly with FINIX and with the reference correlations, as shown in Figure 5. This is because at the time scale of the power transient, the relevant quantity is the fuel heat capacity, which is very similar in both models. However, the maximum void fraction at the core outlet is higher with the reference model, indicating a difference in heat transfer from the rods to the coolant. The gas gap conductance computed by FINIX is lower, and for that reason heat is transferred more slowly from the rods to the coolant. Consequently, the fuel rod temperature decreases more slowly as shown in Figure 6. Maximum fuel temperatures increase approx. 670 °C during the transient, but remain modest because the initial temperature is low.

The power transient is strongly asymmetric in the reactor core and for that reason also the fuel properties vary from one bundle to another. As an example, the maximum gap conductance of each fuel assembly 0.3 s after the control rod ejection is shown in Figure 7a. Dimensional changes of the gas gap in the fuel assembly in the vicinity of the ejected rod are also shown in Figure 7b. The dimensions change mainly due to thermal expansion. In the reference calculation the fuel rod dimensions do not change.

Fig. 6 Fuel maximum temperature and the maximum of radially averaged fuel pellet temperature (left) and the maximum and minimum gas gap conductances (right) during the control rod ejection transient with the HEXTRAN code.

Fig. 5 Fission power (left) and maximum void fraction at core outlet (right) during control rod ejection transient with the HEXTRAN code.

Fig. 7 Maximum gas gap conductance in each fuel assembly 0.3 s after beginning of the control rod ejection (left) and time evolution of gas gap dimensions in a fuel assembly located near to the ejected control rod (right). The location of the fuel rod is shown on the left panel.

5 Summary

A light-weight fuel behaviour module FINIX has been developed. FINIX is aimed especially for multiphysics simulations, where it takes the role of the simulation's fuel behaviour model. FINIX has been designed to be integrated into a wide array of simulation codes, and to provide an identical description of the fuel thermal behaviour across different disciplines such as reactor physics and thermal hydraulics. In cases where the role of fuel performance simulations has been taken by simple correlations, calculation of the thermal response can be improved by including mechanical feedback and power history dependence.

As a demonstration of the applicability of FINIX to multiphysics simulations, we have integrated FINIX into Serpent 2, a Monte Carlo reactor physics code, and to HEXTRAN, VTT's in-house reactor dynamics code. With Serpent, we model a test case of an initially supercritical VVER-1000 pin-cell, where the exponential increase of power is terminated by the negative reactivity feedback from the increasing fuel temperature. With HEXTRAN, a control rod ejection accident in VVER-440 reactor is simulated combining neutronics, thermal hydraulics and fuel performance simulations. The case involves running multiple (several hundred) instances of FINIX in one simulation, and showcases how FINIX can be used to model the fuel thermal behaviour of a whole reactor.

Development of FINIX is an on-going work. Future developments include enhancing capabilities in fuel modelling such as cladding mechanics, and fission gas release modelling to facilitate modelling of loss of coolant accidents and steady state behaviour in addition to the fast transients discussed in this work. In addition, the development of the FINIX interface to multiphysics simulation packages remains a topic of high interest [25].

Acknowledgments

This work was funded by the Finnish research programmes on nuclear power plant safety SAFIR2014 and SAFIR2018, and the NUMPS project of the Academy of Finland.

References

All Figures

	Fig. 1 The models of the FINIX module and its role in a multiphysics simulation. The iteration of the thermal and mechanical solutions is indicated by the flowchart. The convergence checks are assumed to automatically fail on the first iteration.
In the text

	Fig. 2 The radial burnup distribution (a) and the resulting power distribution (b) for the VVER1000 pin-cell at 10 MWd/kgU average burnup calculated by Serpent.
In the text

	Fig. 3 The steady state temperature (a) and strain (b) distributions in the VVER-1000 fuel rod.
In the text

	Fig. 4 Evolution of system linear power (a) and rod centre line (green) and pellet surface (red) temperatures (b).
In the text

	Fig. 6 Fuel maximum temperature and the maximum of radially averaged fuel pellet temperature (left) and the maximum and minimum gas gap conductances (right) during the control rod ejection transient with the HEXTRAN code.
In the text

	Fig. 5 Fission power (left) and maximum void fraction at core outlet (right) during control rod ejection transient with the HEXTRAN code.
In the text

	Fig. 7 Maximum gas gap conductance in each fuel assembly 0.3 s after beginning of the control rod ejection (left) and time evolution of gas gap dimensions in a fuel assembly located near to the ejected control rod (right). The location of the fuel rod is shown on the left panel.
In the text