SERPENT/SUBCHANFLOW COUPLED CALCULATIONS FOR A VVER CORE AT HOT FULL POWER

An increasing interest on the development of highly accurate methodologies in reactor physics is nowadays observed, mainly stimulated by the availability of vast computational resources. As a result, an on-going development of a wide range of coupled calculation tools is observed within diverse projects worldwide. Under this framework, the McSAFE European Union project is a coordinated effort aimed to develop multiphysics tools based on Monte Carlo neutron transport and subchannel thermal-hydraulics codes. These tools are aimed to be suitable for high-fidelity calculations both for PWR and VVER reactors, with the final goal of performing pin-by-pin coupled calculations at full core scope including burnup. Several intermediate steps are to be analyzed in-depth before jumping into this final goal in order to provide insights and to identify resources requirements. As part of this process, this work presents the results for a pin-by-pin coupling calculation using the Serpent 2 code (developed by VTT, Finland) and the subchannel code SUBCHANFLOW (SCF, developed by KIT, Germany) for a full-core VVER model. For such purpose, a recently refurbished master-slave coupling scheme is used within a High Performance Computing architecture. A full-core benchmark for a VVER-1000 that provides experimental data is considered, where the first burnup step (i.e. fresh core at hot-full rated power state) is calculated. For such purpose a detailed (i.e. pin-by-pin) coupled Serpent-SCF model is developed, including a simplified equilibrium xenon distribution (i.e. by fuel assembly). Comparisons with main global reported results are presented and briefly discussed, together with a raw estimation of resources requirements and a brief demonstration of the inherent capabilities of the proposed approach. The results presented here provide valuable insights and pave the way to tackle the final goals of the on-going high-fidelity project.


INTRODUCTION
During the last decade, a global trend for the development of highly accurate methodologies in reactor physics has been observed in several projects worldwide [1,2], where the availability of vast computational resources allows the development of a wide range of coupled state-of-the-art neutronic and thermal-hydraulic (TH) calculation tools. These analyses, usually identified as high-fidelity, tend to diminish the level of approximations in the involved coupled physics through the full-scope coupling of codes with high detail modelling capabilities. Under this framework, the McSAFE European Union project [2] is a coordinated effort aimed to develop multiphysics tools based on Monte Carlo (MC) neutron transport and subchannel thermal-hydraulics codes, suitable for high-fidelity calculations for PWR and VVER reactors for steady-state, burnup and transient calculations. Regarding burnup calculations, one of the project stated goals is to develop full-core depletion with TH coupling at pin level, which represents a challenging target.
For this purpose, a series of tools have been developed, based on the coupling of MC codes and subchannel codes through diverse methods [6,7]. The proper verification and validation of these schemes in a wide range of geometries is compulsory, together with the performance analysis and the identification of potential bottlenecks. For practical reasons, this verification analysis is to be done step-by-step, where each stage provides valuable insights for the whole work path.
In this work, the coupling capability is tested for a full-core VVER case. For such purpose, a coupled master-slave approach between the MC Serpent 2 and the subchannel TH SUBCHANFLOW (SCF) codes is considered. A full scope pin-by-pin modelling is thus developed for a VVER benchmark that provides experimental results for diverse reactor states. As a result, two different reported states from first operating cycle [8] are calculated with this high-fidelity tool, for two cases with and without equilibrium xenon.
Despite no detailed (i.e. high-fidelity level) results are provided within this benchmark for these cases, the comparison of obtained results with these experimental values provide a global assessment of the consistency of the proposed scheme. Besides, in view of the stated high-fidelity goals, a series of detailed analysis are also presented to show the potential capabilities of this scheme, together with the associated computational requirements.

COUPLED SCHEME
A master-slave coupling scheme between Serpent MC [4] and SCF subchannel TH [3] codes is considered, recently rebuild from scratch to tackle burnup and transient calculations for a wide range of reactor geometries [7].

Serpent 2
Serpent 2 is a multi-purpose three-dimensional continuous energy Monte Carlo transport code, developed since 2004 at VTT Technical Research Centre of Finland Ltd. It represents a state-of-the-art code, aimed to perform static, burnup and dynamic 3-D calculations using standard ACE format Nuclear Data Libraries (NDL). It counts with several geometry definition alternatives, which allows to deal with almost all reactor geometries usually encountered. One of the main advantages of this code is its particle tracking, which relies on the combination of conventional based surface tracking and the delta-tracking method. In addition, the code counts with inherent multi-physics features, namely due to the combined capability to handle variable material densities fields (through methods based on rejection sampling techniques) and also manage temperature fields through a rejection sampling approach combined with Target Motion Sampling. In addition, these capabilities can be easily used with the definition of a mesh that superimpose such TH fields in the geometry (namely multiphysics Interface Files -IFC [11]).

SCF
SUBCHANFLOW is a subchannel-level thermal-hydraulic code for steady-state and transient calculations developed at KIT (Germany) [5]. The flow solver consists of four conservation equations, namely for mass, energy, and axial and lateral momentum. The main geometry is defined as a set of channels and rods with given hydraulic parameters and connectivities in a flexible manner, i.e. without assuming a particular geometry type such as square or hexagonal. In addition, a geometry preprocessor allows user to develop suitable inputs for PWR and VVER geometries through a high-level definition of lattices and geometry aspects [6].

Coupling Serpent and SCF considering equilibrium xenon
The coupling approach maintains all inherent Serpent capabilities [4,7], which allows to use the burnup capabilities, such as the Predictor-Corrector scheme (P-C). Moreover, from SCF point of view, each burnup step is done through a steady-state calculation, where the N-TH iteration is done according to Serpent selected scheme. A relaxation on TH fields between iteration steps is done in terms of Equation (1): where T i,j,k t represents the TH field in zone i,j,k for iteration t and w a relaxation factor. In addition a volume average of the fuel temperature within the pin is considered for the feedback, while a L2 norm is considered for the convergence of the TH parameters for each by i,j,k zone.

VVER-1000 BENCHMARK TEST CASE DESCRIPTION
The selected experimental benchmark [8] provides data for a full-core VVER-1000 reactor, whose main aspects are summarized in Table I. This benchmark is based on data from Khmelnitsky NPP 2 nd unit, where successive updates and corrections of reported data were held during last years [9,10]. The core is constituted from several TVSA-type fuel assemblies (FA) that include diverse enrichments and burnable poisons [8]. In addition, the benchmark provides Control Rod (CR) positions and critical boron concentrations for diverse reactor states during the first four operation cycles.
For the aim of comparison and assessment of the coupled scheme here proposed, two specific reactor states from the first cycle were selected, which are listed in Table II. For these cases, critical acid boric concentrations calculation results using standard tools are also reported within the benchmark [8,9] with a span between ~1.0 and 0.5 g/kg (i.e. ~175 to 87.5 ppm boron). Within this work, the values from Table II will be considered and the deviation from critical condition will be used as agreement parameter.

Serpent and SCF Models
A full-core model was developed for Serpent, as shown in in Figure 1, where the inner-core reflector and the upper and lower zones were modeled as homogenized compositions, as proposed in the benchmark specifications. The calculations were performed both using ENDF/B VII.0 and JEFF3.1.1 to estimate the NDL impact in the modelling. Besides, regular hexagonal fission power detectors [11] both at FA and pin levels were included to obtain coarse and detailed results.
a) x-y cut of the model-center of active length c) x-z cut in core center (not at scale) b) x-y cut detail in center of the core

Figure 1 Serpent model details
Regarding TH modelling in SCF, a fuel-centered channel level model [6] was developed, as show in Figure  2. Axially, 30 equispaced zones were considered. This model was coupled with Serpent one using ad-hoc mapping files that map one-to-one TH fields to proper nested hexagonal IFC [11]. Besides, the convergence of SCF was set to be at least two orders of magnitude below the N-TH coupling convergence of the fields. Table II, thus the inlet temperature, boron concentrations and total power were properly set for each case. In addition, to consider xenon concentration for case 2, two burnup steps up to 4.51 FPD were held using a P-C scheme with equilibrium xenon. For this case, as far as the detailed (i.e. pin-by-pin) division of the full-core depletion zones is not possible (due to excessive RAM memory requirements [12]) and the correct modeling of xenon distribution does not need such detail, a simplified material division was considered. Thus, for case 2, a division by type of pin on each fuel assembly with 30 axial zones (totalizing 8160 burnup zones) was held. Besides, the lower optimization mode from Serpent [11,13] was selected for this case to diminish RAM requirements.

Figure 2 SCF model details
Serpent runs were held with 2000 active cycles of 200000 particles each (totalizing 4e8 active histories per iteration step), obtaining a statistical convergence of ~6 pcm. Regarding N-TH iteration, the convergence of the TH fields is set to 30K and 2.5 K for fuel and coolant temperatures, while a 0.01 [g/cm3] value is set for the coolant density variation. Both cases were run in a mid-size High Performance Architecture cluster using a hybrid MPI-OpenMP compilation of the tool, usign 20 OMP threads and 40 MPI nodes (totalizing 800 CPU @ 2.60GHz Intel(R) Xeon(R) CPU E5-2660 v3).

Global results comparison
The converged coupled results for cases in Table II are presented in Table III, where it should be regarded that the experimental configurations are for critical states (i.e. keff = 1.0 is expected). A good agreement is observed in terms of global parameters, with differences in reactivity below 1200 pcm for both cases (with and without equilibrium xenon), where the impact of the considered NDL is found to be ~100 pcm. Besides, the obtained axial offset for both cases present a reasonable agreement regarding the axial details within the benchmark. Nevertheless, in-depth analysis of this issue is foreseen in further works.

Highly detailed results for case 2 (50% of total power and equilibrium xenon)
In addition to the global parameters presented above, the proposed high-fidelity approach offers highly detailed results for relevant fields. To show these capabilities, detailed results are presented for the case 2 (using JEFF 3.1.1 NDL). The temperatures profiles are presented in Figure 3 for a radial and axial cut. Besides, power density results both at FA and pin levels are presented in Figure 4 and Figure 5.  These results show the consistency of the approach, where the behavior of the interchanged fields is as expected. Besides, the Figure 5 demonstrates the capability to obtain the intended pin-by-pin results for this full-scope core VVER case (with an average statistical convergence of 6%).
Finally, in view of the results presented above, potential applications of these full scope high-fidelity analyses could be easily identified, such as the direct calculation of Safety parameters, f.e. the Departure from Nucleate Boiling Ratio (DNBR), which is a standard SCF output.

Potential drawbacks and limitations: main resources requirement analysis
As previously mentioned, the main drawback of these high-fidelity approaches is related to the amount of computational resources involved. A brief summary of average requirements is presented in Table IV for the case 2 (with JEFF 3.1.1 NDL), where it should be regarded that the RAM memory is moderated due both to the optimization mode chosen in Serpent and to the number of burnable materials considered.  [13] (lower memory usage). 2 Maximum stdev in fission for outer top zone of core. Values in center of core as for Figure 5. Reactivity statistical convergence is ~3 pcm in cases without xenon due to higher Serpent optimization mode. 3 For stated convergence criteria 4 For case with equilibrium xenon with lower optimization mode in Serpent (~8200 depletion zones).
The analysis of Table IV shows that the amount of computational resources is very high, directly related to the MC approach of the neutronics. In addition, the RAM requirement imposes a limitation for the full-core pin-by-pin depletion, to be tackled in further works through Domain Decomposition technique [12].

CONCLUSIONS
The McSAFE high-fidelity project developed several tools aimed to provide high-detailed solutions in reactor core physics, obtained with a low number of approximations for steady state, transients and burnup calculations. A proper verification and validation of these schemes in a wide range of geometries is compulsory, providing also valuable insights for the global performance of the approach and allowing the identification of potential bottlenecks for bigger cases, such as full-core scope.
The capabilities to develop full-core calculations with coupled N-TH for a VVER-1000 real geometry were assessed, using a refurbished version of the Serpent-SCF master-slave coupling. Results for coupled cases without and with a simplified equilibrium xenon distribution showed a good agreement with experimental data available, encouraging further calculations. Besides the global agreement, pin-level detailed (i.e. highfidelity) results for the interchanged fields were also analyzed, showing both the capabilities available within this approach and the consistency of the behavior of these parameters.
The results presented here provide valuable insights and pave the way to tackle the final goals of the McSAFE project. As a result, further works will deal with full-core VVER pin-by-pin depletion calculation for the whole burnup cycle, within a fully coupled MC-TH approach.