COUPLING OF nTRACER TO COBRA-TF FOR HIGH-FIDELITY ANALYSIS OF VVERs

Paul Scherrer Institut is developing a high-resolution multi-physics core solver for VVER analysis. This work presents the preliminary stages of the development, specifically the coupling of the 3D pin-by-pin neutronic solver nTRACER to the sub-channel thermalhydraulic code COBRA-TF for single assembly multi-physics steady state calculations. The coupling scheme and the modifications performed in the codes are described in details. The results of the coupled nTRACER/COBRA-TF calculations are compared to the ones of a standalone nTRACER calculation where the feedbacks are provided by a simplified 1D thermal-hydraulic solver. The agreement is very good with fuel temperature differences around 10 K which can be attributed to the different correlations used in the various solvers. The crosscomparison of the two multi-physics computational routes serves as a preliminary verification of the coupling scheme developed between nTRACER and COBRA-TF.


INTRODUCTION
The calculations performed for the design and operation of a NPP and its components are a key factor for the safety analyses of the plant. The conventional approach for the analysis of NPPs involves coarse mesh diffusion calculations with few energy groups coupled to thermal-hydraulic codes capable of assemblywise analysis of the core and coarse mesh modelling of the Reactor Pressure Vessel. In recent years, the evolution of computing clusters and computational techniques lead to the development of codes capable of multi-physics simulations of high-resolution models. These high-fidelity high-resolution codes allow the prediction of safety parameters, like the fuel temperature, on a local scale, simulating heterogeneities and local effects which are approximated or excluded in coarse mesh calculations. Due to their high computational cost, high-resolution methods can only be used as an audit tool of the conventional approach. Currently, VVER technology is expanding with new installations worldwide. Nonetheless, there is no highresolution tool available for VVER steady state and transient analysis, at the time of writing. Such a refined multi-physics code system is necessary for hexagonal geometries. Paul Scherrer Institut has already taken some steps on that direction with the ongoing validation of the 3D pin-by-pin neutron transport code nTRACER for VVER geometries [1]. The goal of this study is the development of a high-resolution core solver based on nTRACER, coupled to the sub-channel code COBRA-TF for core steady state and transient VVER analysis. nTRACER is developed by Seoul National University (SNU) and was recently expanded to the hexagonal geometry of VVERs. It is referred to as nTF in the text. On the other hand, COBRA-TF (CTF) has been verified for 3D VVER full core calculations on a coarse mesh [2] and high-resolution hexagonal assembly analysis [3]. The two codes have been coupled before for steady state calculations in Cartesian geometry through the SALOME platform on the premise of the NURESIM project [4]. The coupling of CTF to nTF, as well as the use of CTF for full core VVER sub-channel safety analysis are, to the extent of the writer's knowledge, a novel concept. This paper focuses on the initial steps of the development of the coupled code system. The latest version of nTF is coupled to COBRA-TF v3.6. The coupling scheme, which can be applied only for single assembly calculations presently, is described in detail in the following sections. The single 3D assembly was preferred as a test-case for the initial stage of this work in order to verify the coupled code system with simple VVER configurations before extending it to full core calculations. In addition, the simple 1D thermal-hydraulic (1D-TH) solver, embedded in nTRACER for Cartesian geometries, is also expanded to hexagonal geometries. The purpose is to cross-examine the results of the two multi-physics solvers, nTF/CTF and nTF with 1D-TH. As one might expect, the crude 1D-TH solver cannot be used to assess the accuracy of a subchannel code like CTF. However, the cross-comparison between nTF/CTF and nTF with 1D-TH allows uncovering coding errors in the coupling scheme. The comparison is performed in a step by step manner, increasing gradually the complexity of the problem and its computational cost. Namely, the examined configurations range from a VVER assembly of 37 rods 40cm in height, to a 3D VVER-440 assembly (127 rods) of 2.5m with axial reflector.

CODE DESCRIPTION
This section is focused on the description of the codes of the high resolution core solver, nTF and CTF, and of the Monte Carlo code Serpent.

nTRACER
nTRACER [5] is a direct whole core transport calculation code. In contrast to the conventional approach, the direct whole core calculation requires the explicit modeling of the nuclide and temperature distributions inside each fuel pellet. nTRACER solves the neutron transport equation at sub-pin level with the use of the planar (2D) Method of Characteristics (MOC) based on 3D CMFD formulation. The axial coupling is resolved by a 1D two-node SENM (Source Expansion Nodal Method) SP3 solver, embedded in the CMFD scheme. From a modeling point of view, Cartesian and hexagonal geometries only differ in the radial direction, thus the introduction of a hexagon-based modular ray tracing scheme for the MOC and a hexagonal CMFD formulation allowed the use of nTF for VVER analysis [6]. The 1D TH solver is calculating the energy balance with the heat convention equation, ܳ = ݉̇ܿ (ܶ ௨௧ − ܶ ), on the axial direction for each pincell to define the coolant temperature for each axial layer. This temperature is used as a boundary condition to calculate the fuel temperature profile by a 1D radial fuel conduction model.

COBRA-TF
CTF [7] is a ΤΗ simulation code designed for Light Water Reactor (LWR) vessel and core analysis. As a sub-channel code, CTF models the detailed intra-assembly flow and sub-pin heat conduction, resulting in higher resolution TH calculations. The COBRA TH analysis is carried out in an array of parallel channels delimited by fuel rods and open gaps, and divided into axial intervals. The volumes bounded by axial planes and channel lateral borders make up the three-dimensional grid of computational cells (control volumes) where the flow differential equations for mass, energy, axial and lateral momentum are solved. CTF uses a EPJ Web of Conferences 247, 02008 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124702008 two-fluid three-field modeling approach. To obtain TH parameters such as pressure gradients and fuel rod temperature in each control volume, one needs to couple continuity, energy and pressure equations.

MULTI-PHYSICS SOLVERS
The coupling scheme is presented in detail, as well as the modifications performed in nTF for the extension of the embedded TH solver to hexagonal geometries.

nTF/CTF
As it was mentioned in the previous section nTRACER was coupled in the past to CTF for Cartesian geometries [4]. A coupling interface was introduced in nTRACER, specifically for data exchange with CTF. In the previous coupling scheme, nTRACER functions as a master code which calls pre-existing CTF subroutines through the coupling interface, to perform the TH calculation when necessary. nTRACER provides the basic geometric description and boundary conditions of the TH model to CTF, as well as the pin power distribution for each iteration of feedback. CTF returns the moderator temperature and density for each pincell. All variables are exchanged directly though the coupling interface allowing an internal coupling. The fuel rod temperature is calculated by the fuel conduction model embedded in nTRACER. CTF is called iteratively by nTRACER until the water temperature calculation is converged. In the new coupling scheme, nTF remains the master code, however the coupling interface is eliminated and CTF is embedded as a library in the neutronic code. This minimizes the time required for the exchange of data between the codes. The internal coupling algorithm is presented in fig.1. nTF creates directly the CTF input with the basic parameters for the TH calculation and the geometry model. This input initializes the CTF calculation. CTF is provided with the detailed pin power distribution by nTF and returns the moderator temperature and density for each pincell, as well as the temperature profile for the fuel rod (fuel, gap and cladding). nTF uses the TH parameters to calculate new cross-sections with the subgroup method and then proceeds to the CMFD and planar MOC calculations. TH feedback is not required for every outer iteration of nTF. As the neutronic solution is converging the pin power distribution is changing. However, the neutronic solution requires a higher accuracy (10 -6 ) than the TH calculation (10 -4 on moderator temperature). Thus, after a number of iterations the changes in the power distribution do not affect significantly the TH feedback which is considered converged; the calculation proceeds in standalone mode.

Domain Decomposition
nTF and CTF use different meshes for the 3D domain decomposition in hexagonal geometries. The basic radial decomposition for the two codes is presented in fig. 2, together with the superimposition of the relative meshes. The smallest structural element for the nTF models is the pincell, which consists of a fuel rod segment surrounded by coolant. The pincell is divided radially in several regions of equal volume with spatially flat cross-sections, which are defined by a unique set of parameters (e.g. temperature) for each region (FXR). As can be seen in fig. 2, the shape of the pincell depends on its position in the fuel assembly. This irregularity is particular for hexagonal assemblies due to their geometry. The inter-assembly water gap is modeled by individual water cells. Axially the solution domain is divided into several layers on a coarse mesh, since most reactor cores are radially more heterogeneous than axially. As it was mentioned also in the previous section, CTF decomposes the domain into an array of channels which are delimited by fuel rods and the open gaps, defined on the vertical planes crossing the centers of neighboring rods. For the treatment of hexagonal geometry CTF combines the use of triangular channels in the assembly and square channels at the assembly boundary, which include the inter-assembly gap. Axially, the decomposition is the same as nTF. The methodology used by CTF allows the modeling of any possible configuration of fuel rods and assemblies as long as the neighboring channels are defined in the input data. For the coupled code system the temperature profile of each pincell (fuel, gap, cladding, moderator) is transferred to nTF for the cross-section generation. The water density and temperature in the pincell is calculated by averaging the corresponding parameters of the neighboring channels, according to the fraction of the wetted perimeter of the pin that correspond to each channel ( fig.2). This means that all the FXR regions of the moderator in a pincell have the same water density and temperature.

Modifications in the codes
The latest version of nTF is extensively modified in order to perform coupled TH calculations, with CTF, for single VVER fuel assemblies. In the previous coupling scheme, a pre-processor was used, provided by the COBRA-TF coding team, to create the CTF input for Cartesian geometries. Thus an extra independent executable was included in the coupling. nTRACER fed all relevant TH parameters to the pre-processor, which then produced the CTF input. In the new coupling scheme the use of a pre-processor is eliminated.
A new subroutine is included in nTF based on the work done by KIT on the modeling of VVER assemblies nTF CTF EPJ Web of Conferences 247, 02008 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124702008 with CTF [3]. This subroutine creates the CTF input directly, without the use of an extra executable, simplifying thus the execution of the coupled code system.
Besides the structural changes of the coupling scheme, the existing relevant subroutines are extended to treat the hexagonal geometry. Several geometry dependent functions are modified to match the architecture of the VVER fuel assembly. This is necessary due to the differences in the domain decomposition used by the codes for Cartesian and hexagonal geometry. One of the main modifications is the introduction of a new subroutine for the correlation of the different fuel rod numbering schemes used in the neutronic and the TH code, as depicted in fig. 2. Moreover, due to the difference in the modeling of the inter-assembly water gap in Cartesian and hexagonal geometries a new subroutine for water-gap treatment is necessary. For Cartesian geometries the water gap was divided into several segments, aligned to the pincells on the assembly side.
On the other hand, for the hexagonal assembly the water gap is divided in twelve segments, two for each side of the hexagon, independently of the assembly pitch. For this work the moderator temperature and density in a gap segment is calculated as a weighted average of the neighboring cells, according the length of the boundary of the cells that is connected to the gap. Finally, the treatment of guide tubes (GTs) and axial reflectors is also included. A GT is modeled explicitly with CTF, whereas the axial reflector is defined as an additional axial layer of the fuel rods with zero power. This approximation is considered sufficient at this preliminary stage of the coupling scheme. nTF CTF Figure 3. Radial discretization of the fuel pin in nTF and CTF and for the fuel conduction solver.
The previous coupling scheme created for Cartesian geometries required only the moderator temperature and density for each pincell as TH feedback from CTF. The temperature profile of the fuel pin, the gap and the cladding was calculated by the fuel conduction solver of nTRACER for each pincell. However, in the high-resolution coupled code system the temperature of the fuel rod materials should be provided by the TH code, since the models and correlations used for the fuel conduction by CTF allow higher accuracy than the simplified solver which is included in nTF. The discretization of the fuel rod in nTF and CTF for the solution of the conduction equation is presented in fig.3. On the fig. the calculation points of the temperature are presented as well for the two codes. nTF divides the fuel pin in equispaced nodes. The fuel temperature is initially calculated on the surface of the nodes and subsequently on the node centerline with interpolation. The TH nodes in the fuel pin differ from the FXRs used in the neutronic calculation. The corresponding FXR temperature is calculated as the volume average value of the TH nodes or segments of the nodes contained in the FXR. The cladding is divided in two nodes. The cladding temperature used for the crosssection generation is the volume averaged value at the centerline of the cladding. The gap temperature is defined as the mean value of the temperature at the fuel outer surface and the inner surface of the cladding. CTF on the other hand, further discretizes the TH radial nodes in azimuthal segments and defines the fuel EPJ Web of Conferences 247, 02008 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124702008 temperature at the center of the segment for all inner subdivisions and at the surface of the fuel pin for the last node which is half the width of the inner nodes. The cladding is divided in two nodes and the temperature is defined at the inner and outer surface. The gap temperature is not calculated by CTF, which assumes constant gap conductance. In order to use the CTF fuel temperature in nTF, the CTF fuel pin discretization was introduced in nTF. A subroutine was added in CTF to calculate the azimuthal average of the temperature for each radial node and transfer the data to nTF. In the neutronic calculation the temperature for the FXR regions is calculated with the same method, nonetheless the pin surface temperature is used for the outer TH node. In general, the fuel conduction calculation requires a large number of TH nodes for the calculation of the fuel temperature (~30), significantly larger than the number of FXRs used (~6). Thus, the use of the surface temperature rather than the node average has a small impact in the cross-section generation. The cladding temperature is defined as the volume average of the surface temperatures. Finally, the gap temperature is calculated as the mean of the temperatures at the fuel outer surface and the cladding inner surface.

nTF with 1D-TH
The 1D TH solver of nTRACER is also extended to model VVER assemblies, including GTs and the axial reflector. The water gap is treated the same way as for nTF/CTF. Several approximations are involved in the simplified solver. As it was already mentioned the 1D heat convection equation is solved independently for each pincell. The coolant mixing in the assembly is not taken into account. Despite the different shapes of the pincells ( fig.2) the hydraulic diameter used is calculated with the area of the inner pincells and the average diameter of fuel rods and GTs. The moderator temperature in GT channels is calculated as the average moderator temperature of the assembly. The temperature of the channels in the top reflector is the outlet temperature of the last active layer for each channel. In order to define the fuel temperature profile, each pin is discretized ( fig.3) and the fuel conduction equation is solved with the volumetric heat produced in the pin as a scaling factor and the pellet surface temperature as a boundary condition.

RESULTS
The cross-comparison of the two multi-physics solvers ranges from a simplified VVER assembly of 37 rods 40cm in height, to a 3D VVER-440 assembly (127 rods) of 2.5m with axial reflector. In order to test the coupling scheme, the outcome of the calculations is presented for several configurations of different complexity. Fig. 4 presents the radial distribution of the difference between the pin power, coolant temperature (T) and fuel T as calculated by nTF/CTF and nTF with 1D-TH for an assembly of 37 pins and a VVER-440 assembly, both 2 axial layers thick with reflective boundary conditions 1 . Fig. 5 depicts the profile of the difference between the average pin power, coolant T and fuel T of each axial layer of an assembly of 37 pins with 2m height (reflective b. c.) and a 2.5m VVER-440 assembly with an axial reflector.
The calculations for the small assembly of the 2 axial layers result in minor differences between nTF/CTF and nTF/1D-TH for the pin power and the coolant T. The fuel T difference reaches 10K so further investigation was necessary. The same configuration was simulated using the nTF/CTF coupled system, but this time with the fuel T feedback provided by the fuel conduction model in nTF. Similar discrepancies were observed when cross-comparing the result of the CTF and nTF fuel conduction model for the same moderator T. The cause of the discrepancies can then be attributed to the different correlations used by CTF and nTF for the calculations of key parameters like the fuel conductivity. It is obvious by figs 4 & 5 that the magnitude of the difference between the results of the two multi-physics solvers is increasing as the size of the assembly increases axially and radially. Specifically, for the 3D VVER-440 assembly the pin power difference reaches 2% and the difference in the fuel T 14K.   The approximations involved in the 1D TH solver are the cause of the discrepancies. As the size of the assembly increases radially, the power, and subsequently fuel T, distribution is shifting to a sharper cosinus profile. The temperature heterogeneity results in increased coolant mixing. At the same time as the height of the assembly increases the temperature heterogeneity increases as well, thus the phenomena are accentuated. Fig 5 shows that the larger discrepancies for the coolant and fuel T occur on the axial levels where the radial power profile is more heterogeneous. Nonetheless, the differences remain in acceptable levels which gives confidence that no large coding error exists in the coupling scheme nTF/CTF for hexagonal geometries. Finally, it must be pointed out that the total calculation time required for TH analysis of the VVER-440 3D assembly is 265sec for the coupled code system whereas the simplified 1D solver requires ~3sec. This is expected, considering that CTF is a separate high fidelity code. Nonetheless, in both cases the time required for the TH calculation is <1% of the time needed by nTF.

CONCLUSIONS
The preliminary steps for the development of a high resolution core solver with nTF and CTF for VVER safety analysis are presented in this work. The existing coupling scheme between the two codes is improved and extended to hexagonal geometries. nTF functions as a master code, internally coupled to CTF, which is embedded in the source of the neutronic code as libraries. Several modifications are necessary due to the different domain decomposition and geometry models used by nTF for hexagonal assemblies in comparison to Cartesian geometry (e.g. water gap). In addition, both codes are modified in order to allow the exchange of the fuel temperature profile which in the previous scheme was calculated by the fuel conduction model in nTF. The newly developed coupled code system is used for the simulation of single 3D VVER assemblies ranging in size axially and radially. The outcome of the simulations is cross-compared with standalone nTF, with TH feedback provided by the embedded 1D TH solver, which was also extended to hexagonal geometries on the course of this work. The results verify the coupling scheme which will now be expanded for full core calculations and validated by experimental data (Rostov-2 OECD benchmark).