Verification and validation of the high-resolution code nTF with VVER-1000 hot zero power monte carlo calculations and experimental data

Novel codes allow the prediction of parameters with increased resolution, in comparison to the conventional approach, sometimes using meshing smaller than the size of the fuel pin. Optimally, complex geometric structures, like the ones included in the VVER core e.g. corner stiffeners, would be modeled explicitly by novel neutronic codes. In this paper, the performance of the deterministic neutron transport code nTF is studied for VVER geometries. The code is verified standalone with Monte Carlo reference solutions for a range of Hot Zero Power (HZP) models, from 2D pincells to the 3D full core of a VVER-1000, described in the X2 benchmark. Several modeling options and approximations are studied to optimize the nTF X2 3D core model. The updated version of nTF, including all necessary modifications to improve its VVER modeling capabilities, presents good agreement with the reference solutions for the simple VVER configurations and the critical state of the X2 core. Specifically for the full core critical state, the RMS of the pin power difference is 1.5% and the eigenvalue difference 74 pcm with respect to a Monte Carlo reference solution. Any resulting discrepancies are analyzed, the more significant ones with the use of the lattice code CASMO5-VVER that employs similar methods as nTF. The optimized nTF X2 full core model is validated with the experimental data included in the X2 benchmark, produced by HZP start-up tests. The differences of the nTF results with the measured and simulated data of the X2 benchmark remain within the accuracy limits set by the utilities and the literature of deterministic code verification for most cases. The related discrepancies are discussed in detail.


Introduction
A key aspect for the safe operation of a Nuclear Power Plant is the accurate prediction of its core behavior through numerical tools, which are capable of simulating successfully all relevant physical phenomena, including neutronics.The conventional neutronic analysis aims to simplify the core geometry and material distribution by replacing it with a configuration of homogeneous nodes of one material with the same neutronic behavior as the section of the core they model.Then a simplified balance equation for the neutronic population of the full core is solved on this coarse mesh.The effect of local perturbations is handled through the coarse mesh neutronic parameters (cross-sections and discontinuity factors) determined for each homogeneous node.The local values of the neutron flux or power can be reconstructed by response functions, which also take into account the local geometry of the core (Bahadir and Georgieva, 2018).In order to achieve the desired level of accuracy, the neutronic code should be able to model explicitly all the geometric features of the reactor which have an impact on key parameters, or reproduce their effect with modeling techniques.Unlike the conventional approach, high-resolution solvers examine the core in a much more detailed level, sometimes using meshing smaller than the size of the fuel pin.All neutronic data are used directly by the code with no prior homogenization.Thus, a high-resolution neutronic solver needs to perform neutron balance calculations with few approximations.Ideally, all high-resolution solvers would have the capacity to model explicitly the core geometry in full detail, like Monte Carlo codes.However, deterministic codes are often limited in terms of geometry description, due to the discretization schemes involved in the solution of the transport equation.Therefore, case specific modeling approximations are necessary, together with corresponding numerical techniques (e.g.methods to increase stability) to allow the effective modeling of real reactor cores.High-resolution codes can be inefficient when strong geometry approximations are applied on elements of the core with a significant impact on the neutron flux (Papadionysiou et al., 2021).Thus, despite the complexity of the geometry, some structural elements need to be modeled explicitly, in order to successfully predict local parameters.However, complex geometries often result in computationally expensive calculations.Besides the actual reproduction of complex geometric structures, which may require the introduction of additional cells, the solver needs to build the appropriate computational grid to achieve the convergence of the solution.Obviously, the finer the computational mesh becomes, the more expensive is the calculation.The successful use of a high-resolution deterministic solver lies on the optimal balance between spatial resolution, modeling accuracy and computational requirements to ensure the accurate prediction of local and global key parameters, as well as cost-effective calculations.
Currently, the water-water energetic reactor (VVER) technology is under expansion, including new projects in emerging nuclear countries but also in established ones such as Finland.The VVER is a pressurized light water reactor.The main distinguishing feature of the VVER core is the hexagonal lattice of the fuel assemblies, which requires the development of specialized numerical methods to solve the transport equation.In addition, the core includes complex geometric structures, described thereafter, that need to be modeled effectively by deterministic codes, due to their impact on the neutron flux.A VVER-assembly can be built with corner stiffeners and spacer grids to prevent deformation.Another important feature of the VVER, is the heavy reflector.Concerning multi-physics calculations, the impact of the heavy reflector model on the accuracy of global and local parameters is very strong for the neutronic analysis (Sargeni et al., 2016).It can be attributed to the structure of the reflector, which contains localized moderator-rich cavities.This complex geometry poses an issue for deterministic codes (Petkov and Mittag, 2006).High-resolution solvers especially, must reproduce the effect of the reflector as realistically as possible, in order to grasp the local flux and power distribution accurately everywhere in the active core (Papadionysiou et al., 2021).
Considering the status of VVERs worldwide the Laboratory for Reactor Physics and Thermal-Hydraulics (LRT) in Paul Scherrer Institut (PSI) is currently developing a high fidelity multi-physics code system for VVER high-resolution steady state and cycle calculations.The neutronic analysis is performed with the deterministic direct whole core neutron transport code nTRACER-Fast (nTF), of Seoul National University (SNU) (Jung et al., 2013;Kim and Joo, 2020).The recently developed hexagonal module of nTRACER is capable of 3D sub-pin transport calculations.Older versions of the code were verified and validated (Papadionysiou et al., 2020;Papadionysiou et al., 2021) for 3D full core geometries in the course of the LRT VVER project.Besides verification and validation, the purpose of the standalone analysis is to discover the weak points of the neutronic code, before introducing the effect of the T/H code and the depletion library for multi-physics cycle analysis.It allows to distinguish the different sources of discrepancy in the calculations and optimize the use of the solver (e.g.computational mesh).Thus, it must also be determined, at this stage, if any modifications on the neutronic code are necessary to allow for better application of the computational methods or improved modeling of the core.This work is a continuation of (Papadionysiou et al., 2020;Papadionysiou et al., 2021).All observations and conclusions derived from (Papadionysiou et al., 2020;Papadionysiou et al., 2021) are taken into account to improve the modeling capabilities of nTF.The solver is compared with a Monte Carlo code, Serpent2, and a deterministic code of similar computational methods, CASMO5-VVER, in order to study its capacity to simulate the behavior of a VVER-1000 core and any accompanying discrepancies.The optimal way to perform this kind of study is the use of fresh Hot Zero Power (HZP) benchmarks; the material composition and temperature are fixed throughout the model, there is no need for T/H feedbacks.The factors which determine the outcome of the calculation are the accuracy of the geometry description, the cross-section libraries and finally the computational methods used.Similarly with (Papadionysiou et al., 2021), the reference Monte Carlo (Serpent2) solution for the HZP critical state and the HZP start-up tests of the X2 benchmark (Bilodid et al., 2020) are used to verify and validate the latest version of nTF.

Code description
This section is focused on the description of nTF and the codes used in this study for the verification of the neutronic solver.

nTF
nTF is a direct deterministic code capable of 3D sub-pin resolution (Jung et al., 2013) in hexagonal geometries.The code solves the neutron transport equation at the sub-pin level while considering the full extent of the reactor geometry with a 2D/1D scheme and 3D acceleration.Cross-sections are calculated on the fly by the subgroup method.nTF is using the nTRACER library, which is built by SNU and it is based on ENDF/B-VII.0.The explicit modeling of the nuclide and temperature distributions inside each fuel pellet is required.nTF is based on the planar (2D) Method of Characteristics (MOC) based on 3D Coarse Mesh Finite Differences (CMFD) formulation.It divides the 3D core into separate axial layers or planes, which are modeled explicitly.The planar transport equation is solved, and then the different planes are coupled through the neutron axial leakage source.The axial leakage is calculated by the CMFD formulation, which includes a 1D MOC axial solver in the latest version of nTF (Choi and Joo, 2020).The flux and cross-sections of the MOC planes are homogenized to create the 3D CMFD problem.The purpose of the CMFD formulation is to combine the 2D/1D scheme and accelerate the MOC calculations.For the hexagonal core, the planar MOC calculations use a hexagon-based modular ray-tracing scheme and the CMFD formulation is applied the unstructured geometry of the polygon cells creating the hexagon-based mesh (Yoo et al., 2015).nTF employs hybrid parallelization methods (MPI & OpenMP) with an axial domain decomposition.The default scattering treatment is Transport-Corrected P0.In the latest version of nTF there is the option for higher orders, which is provided and applied in some part of this study (see Section 4.1.1& 4.2.1).

Serpent2
Serpent2 (Leppanen, 2011) is a multi-purpose three-dimensional continuous-energy Monte Carlo particle transport code, developed at VTT Technical Research Centre of Finland, ltd.As a Monte Carlo code, Serpent2 can explicitly represent the phenomena occurring in the reactor core without the biases induced by the discretization of the mathematical models.Thus, it is commonly employed for the verification of deterministic transport codes.In this study, Serpent2 is used with ENDF/B.VII.0.

CASMO5-VVER
CASMO5 is a multi-group 2D transport lattice code for burnup calculations on BWR and PWR configurations.CASMO5-VVER (Ferrer et al., 20192019) is the extension of the code to hexagonal geometries for VVERs.The code can simulate cylindrical fuel rod bundles of varying composition in a hexagonal pitch array with allowance for water gaps and corner stiffeners in the regions separating fuel assemblies.The heavy reflector can be simulated with homogeneous pins of water and steel, mapping the corresponding complex geometric structures in the 2D full core model (Jang, 2021).Nuclear data for CASMO5 is collected in a library containing microscopic cross-sections in 586 energy groups, based on the ENDF/B-VII.0 nuclear data evaluation, which are collapsed to 19 or 35 groups.The energy-condensed neutronic data constitute the input to the 2D MOC transport calculation.CASMO5 uses predetermined ray spacing and azimuthal angles (0.05 cm & 72 angles) based on extensive testing to achieve high accuracy.The transport equation is solved with a Linear Source approximation to the MOC, which is available for both Transport-Corrected P0 isotropic and full anisotropic scattering treatments.For the hexagonal geometry, the code employs appropriate quadrature sets for the MOC and the Coarse Mesh Non-linear Diffusion acceleration.

Benchmark description
The aim of the X2 benchmark (Lötsch, et al., 2016) is to develop a VVER-1000 data platform for the verification and validation of reactor simulation tools.It is based on the operational data of the second VVER-1000 unit of the Ukrainian NPP Khmelnitsky.The benchmark consists of three stages including HZP experiments, cycle depletion and several transients.An effort to revise the X2 benchmark was initiated with an updated and refined description of the HZP experiments (Bilodid et al., 2020).The new benchmark specifications include detailed geometry models of the active core and the heavy reflector, with updated material compositions and operational data for the critical state (HZP) and all reactor states corresponding to the HZP experiments.A description of the HZP start-up tests together with the relevant measurements, like temperature reactivity coefficients (TRC), SCRAM worth and worth of specific rod banks, are also part of the revised specifications.In addition to experimental plant data, Monte Carlo reference numerical solutions are used to provide quantities not directly available from the measurements.The power distribution at the assembly and the pin level for the critical state is produced by the Monte Carlo code Serpent2 and provided for code verification.Fig. 1 presents the layout for the fresh core and the structure of the heavy reflector as depicted in the benchmark specifications.
The measurements of temperature reactivity coefficients (TRC) were conducted at zero-burnup state (no fission product poisoning, no Xe-135, no Sm-149) with varying inlet coolant temperature (Tin).The experiment was repeated for two different values (TRC1 & TRC2) of boron concentration (CB) and with control rods in banks 1 to 6 fully withdrawn (H = 100 %).Control rod movements in banks 7 to 10 compensated the reactivity insertion, caused by the coolant temperature change.The SCRAM worth start-up tests included two steps.First, from the nearly critical state all control rods, except the chosen ''stuck" cluster (see Fig. 2), were fully inserted.The SCRAM with stuck rod was measured.Secondly, the ''stuck" rod was dropped down and the full SCRAM worth was measured.The coolant temperature and boric acid concentration remained practically unchanged during the SCRAM test.Finally, the integral and differential worth of control rod bank #10 was measured with steps' length of 4 %-7%, starting from fully inserted (0 %).The core inlet coolant temperature and the pressure above the core were kept stable during the measurements.Table 1 lists the operational conditions, including control rod positions, for all reactor states studied in the benchmark.Fig. 2 presents the position of the control rod banks in the core, which are distributed with 120 • rotational symmetry.The "stuck" rod is also depicted in Fig. 2.Table 2..

nTF X2 model
The nTF X2 model used in this paper is built according to the revised benchmark specifications (Bilodid et al., 2020).Several approximations are necessary.Fig. 3 depicts all geometry approximations applied to the nTF model (Papadionysiou et al., 2021).Side (i) presents the X2 reference model built by Serpent2 and side (ii) & (iii) the equivalent nTF approximation.
A few but important modeling options for the geometry and the computational methods used in nTF are listed below.
• In the benchmark specifications, spacer grids are approximated as rings around the fuel cladding, with the same mass and thickness as the actual spacer grid (see Fig. 3.Ai).In nTF, the axial decomposition of the 3D core with nodes smaller than 10 cm in height can cause instabilities (Kelley, 2015).Thus, the spacer grids are modeled as an added layer of homogenized steel and water around the fuel cladding, extending to the same radius as the steel rings approximating the spacer grids in the Serpent2 model (see Fig. 3.Aii).• In the Serpent2 model, the axial reflectors are approximated with 4 layers of homogeneous materials, defined in the benchmark specifications (see Fig. 3.Ai).Due to the small thickness of these layers, the axial reflectors need to be further homogenized in the nTF model (Kelley, 2015).The segments T1 & T2 forming the top reflector are homogenized into a single region (see Fig.  reflector is divided in regions B1 & B2 below the lowerplug in the Serpent2 model.nTF simulates two homogeneous regions as well. The lower segment is built with the material B2.The top segment consists of a homogeneous mixture of the lowerplug, the 1st homogeneous region of the reference model (B1) and part of the 2nd homogeneous region of the reference model (B2) (see Fig. 3.Aii).• Fig. 3.Bi presents the explicit model of the corner stiffeners, surrounding the assemblies of the active core.The radial decomposition of a hexagonal assembly in nTF consists of a honeycomb mesh of fuel pins, surrounded by the water gap of the assembly subdivided in rectangular cells.The previous version of nTF, which was verified and validated in (Papadionysiou et al., 2021), lacked the capacity to define assembly gap cells of different material; thus, the corner stiffeners of the active core were not modeled explicitly in (Papadionysiou et al., 2021).Instead, the mass of the corner stiffeners, of one assembly, was homogenized with the mass of the moderator in the water gap (see Fig. 3.Bii).In (Papadionysiou et al., 2021), it became apparent that the pin power at the corners of each assembly was overestimated by nTF due to the homogenization of the corner stiffeners.Extra moderation at the corners of all the assemblies due to the coolant-steel mixture was creating anomalies to the power distribution.Hence, the code was further extended by SNU to model the corner stiffeners.The model remains semi-explicit, since the latest version of the code only allows to define an integer number of gap cells (see Fig. 3.Biii), starting from the assembly corners, with the corner stiffener material.This means that the length of the corner stiffener can only be approximated by a multiplier of the pincell size.Nonetheless, there is small difference between the actual length of the corner stiffener and the model.• According to the benchmark specifications, the central guide tube has a radius of 0.65 cm.In Serpent2, it can be modeled explicitly; however, the radius is reduced in the nTF model since it exceeds the size of the pincell (flat-to-flat is 1.275 cm).• The reference Serpent2 calculation is performed with the ENDF/B-VII.0 library.The nTF cross-section library is also based on ENDF/ B-VII.0.The X2 solutions are produced with transport corrected P0 scattering.The HZP critical case is also simulated with P2 scattering to study the effect on a 3D model.The latest version of nTF was developed by SNU to allow the use of the P2 scattering approximation for hexagonal geometry.• The heavy reflector is modeled explicitly in Serpent2, except the groove region (Fig. 1, right), which is reproduced as a homogeneous mixture of coolant and steel, equivalent of the real geometry at the edge of the core basket (cylindrical with horizontal grooves).nTF allows the modeling of the water holes, the water liner between the core basket and barrel, the downcomer and the groove with homogeneous pincells of water and steel that map approximately their shape in the reflector assemblies (Fig. 3.Cii).In (Papadionysiou et al., 2021), the water gap surrounding the active care is modeled also with homogeneous pincells, as a volume averaged mixture of water and steel (Fig. 3.Cii).After modifications by SNU (Kim and Joo, 2020), nTF is also capable of reproducing the water gap explicitly, except the corner cell (Fig. 3.Ciii) (Papadionysiou et al., 2021).• Control rods (CRs) are modeled explicitly in nTF by adjusting the size of the axial layers to match the insertion level of the CRs in the core and the change of CR material (30 cm of B4C and the rest Dy2O3•TiO2).As a result, the nTF axial mesh is differing from the Serpent2 model (20 layers × 17.75 cm).Nonetheless, this does not pose any obstacle for the comparison of the results and removes a possible source of discrepancy, which would be the homogenization of the CR material in partially rodded nodes.

Results
As a first step, in order to examine thoroughly the performance of the novel neutronic code, a sensitivity study is realized which includes the comparison of nTF deterministic calculations with Serpent2, for simple VVER configurations.Subsequently, the study moves to the 3D full core X2 benchmark, where the performance of the code is verified against Serpent2 and validated against experimental data in HZP state.Different options for the modeling of the heavy reflector are examined, including what was presented in (Papadionysiou et al., 2021).CASMO5-VVER is also compared with Serpent2 for 2D VVER-1000 full core calculations, in order to evaluate the effect of the heavy reflector model and the scattering treatment with another MOC code besides nTF.

nTF vs Serpent2 for simple VVER configurations
The goal of this study is to analyze separately the various models required to simulate successfully elements of the VVER core with a deterministic code.The cross-examination is performed for a range of models, from pincells to 3D fuel assemblies with Gd pins, increasing gradually the complexity of the problem.Each step aims at the evaluation of different nTF features, which could be sources of discrepancies.Several 2D pincell models of different enrichment (with and without Gd) are cross-compared in order to assess the nTRACER library, e.g.resonance self-shielding methods.Pincell models are appropriate for this purpose since they are less affected by heterogeneous transport.This is studied with the simulation of 2D VVER-1000 assemblies with reflective boundary conditions.One assembly model contains pins of the same low enrichment whereas the other assembly presents strong radial heterogeneity through the use of Gd pins and fuel pins of different enrichment.The same assemblies are also modeled in 3D, with top & bottom reflectors, to determine the effect of the 1D MOC axial solver and void boundary conditions.

Pincell models
The 2D pincell models of different enrichment correspond to the fuel compositions described in the X2 benchmark (Bilodid et al., 2020).pins that are modeled with 10.
Overall, the eigenvalue comparison suggests a good agreement of the nTF and Serpent2 libraries since it remains well below the target accuracy (<200 pcm) set by literature for the verification of deterministic codes (Kim, et al., 2017) for all models.With P0, nTF overestimates the eigenvalue in all cases, except the Gd pins and the 1.3 % enriched fuel.The difference is increasing with enrichment, except for the 4.0 % pin.On the other hand, the effect of resonance self-shielding is obvious on the Gd pincells, even when using 10 radial subdivisions.nTF underestimates the eigenvalue in comparison to Serpent2, in contrast with the pincells without Gd.When P2 is used nTF overestimates the eigenvalue for all pincells.The difference is increased; however, it follows the same trend as observed for P0, increasing with enrichment.Again, the Gd pins present a lower discrepancy in comparison to the pins without Gd, nonetheless the difference between the discrepancies (with and without Gd) is decreased.

Assembly models
The assembly models studied in this section correspond to two X2 fresh fuel assemblies described in the benchmark specifications.As illustrated in Fig. 5, assembly 13AU is built with a single type of low enrichment fuel (1.3 %) and assembly 390GO (4.0 % inside) contains a ring of different enrichment pins at the outer boundary (3.6 %), together with several Gd pins (3.3 % with 5.0 % Gd) close to the center.These specific models were selected to illustrate the effect of MOC transport for various "levels" of heterogeneity in the lattice.Both assemblies are simulated at 600 K, with all material compositions corresponding to the X2 critical state.The P0 transport corrected scattering treatment is used.As it was discussed in Section 3.1, the outer diameter of the central guide tube is larger than the pincell size in nTF and it cannot be modeled according to the specifications.Thus, the diameter is decreased both in the nTF models of the two assemblies and in Serpent2.
Fig. 6 presents the relative difference of the normalized power, [(Nor.Pow.Serpent2 -Nor.Pow.nTF ) / Nor.Pow.Serpent2 ], on the pin level for the 2D assembly.The uncertainty of the Serpent2 calculations remains < 0.03 % for the pin power and 2 pcm for the eigenvalue (6,000 active batches and 140 inactive batches of 500,000 histories each).The difference of the eigenvalue for the 2D assembly calculations is respectively 59 pcm and 25 pcm, remaining in acceptable levels (<200 pcm).The effect of the 2D calculation on the eigenvalue is not consistent with the pincell results.nTF underestimates the eigenvalue for the 2D assemblies, where the maximum absolute eigenvalue difference corresponds to the low enrichment assembly.For the 13AU assembly (1.3 % enriched fuel), the relative power difference remains < 0.2 % with an RMS of 0.06 %.The effect of heterogeneous transport from the MOC calculation is minimal.The heterogeneous 390GO assembly presents higher discrepancies with a maximum at 0.33 % and RMS at 0.08 %.The RMS is very close to that of the homogeneous assembly since the maximum discrepancies appear only for the Gd pins.For both cases nTF presents a tendency to flatten the power profile, overestimating the power at the edge and underestimating it in the center of the assembly.This is a common occurrence with MOC codes (Jang, 2021) and is studied further in this work for the full core.It is clear that the flux dip of the Gd pins has an effect on the nTF pin power prediction, however there is still very good agreement with Serpent2.The heterogeneous 390GO assembly is simulated also with 4 subdivisions in fuel pins and 7 subdivisions in the Gd pins to study the effect of a coarser mesh.The pin power discrepancies (maximum and RMS) of the Gd pins remain the same.The eigenvalue difference increases by 10 pcm.
Finally, in order to study the effect of the 1D MOC axial solver the two assemblies are modeled in 3D.The same pin layout, as depicted in Fig. 5, is extended to 3.55 m, the height of the active core for the X2 benchmark.A coolant bottom and top reflector of 17.75 cm is added to reduce the axial leakage due to the void boundary conditions (22 axial nodes in total).The uncertainty of the Serpent2 calculations remains < 0.08 % for the axial power profile and 4 pcm for the eigenvalue.For the 3D 13AU assembly, the maximum relative power difference in the axial profile is 0.37 % and the RMS 0.23 %.The difference in the eigenvalue is 40 pcm.Similarly, the 390GO 3D assembly has a maximum relative power difference of 0.36 % and RMS of 0.15 % and a difference in the eigenvalue of − 2 pcm.The axial solver seems to have some effect on the power distribution; however, it is independent of the assembly heterogeneity and it remains at an acceptable level.The maximum absolute eigenvalue difference corresponds to the low enrichment 3D assembly, following the trend observed for the 2D assembly models.However for both 3D assemblies the difference is decreased in comparison to the 2D cases, which can be attributed to the axial solver.
Overall, the comparison of nTF vs Serpent2 for small VVER configurations proved the capacity of the code to deal with computational features, which can be sources of inaccuracy in deterministic transport codes.The effect of some said features (e.g.flux dip of Gd pins) is observed in the nTF results and can be expected to emerge in the calculations that follow.However, the corresponding discrepancies are small enough to verify that nTF is capable to handle basic functions of simulating hexagonal fuel elements.

Verification of nTF for the X2 HZP critical state
The Monte Carlo reference solution, included in the benchmark specifications, corresponds to the critical state of the reactor.For the critical conditions, the eigenvalue calculated by Serpent2 is equal to 1.00012 ± 0.00001.The results of five independent Serpent2 runs, each simulating 48 × 10 9 neutron histories, are averaged.The standard deviation for the assembly and pin power values is < 0.1 %.The full core is simulated with both codes, however all results are provided for the 1/ 6th core.nTF employs a ray spacing of 0.01 cm, with 12 azimuthal angles and 4 polar angles, similarly to (Papadionysiou et al., 2021).The fuel pins are modeled with 4 radial subdivisions, except Gd pins that are modeled with 7, following the observations of Section 4.1.2.3.Cii) and the smearing of the corner stiffener material in the water gap of each fuel assembly (see Fig. 3.Bii).In model II (Papadionysiou et al., 2021) the corner stiffeners are again homogenized in the assembly gap, but the water liner is modeled semi-explicitly (see Fig. 3.Ciii).Model III corresponds to the latest version of nTF which can model the water liner and the corner stiffeners semi-explicitly (see Fig. 8.Biii).In this case, a single full core calculation requires 300 cpus (12 cpus × 25 axial layers) for 15,256 secs (~4 hrs and 14 mins) and 7.44 TB (88 % of the total memory corresponding to 25 nodes).All solutions are produced with transport corrected P0 scattering.As it was discussed in (Papadionysiou et al., 2021), the map of the pin power difference for model I produces its maximum at 7.6 % on the pins neighboring the reflector.This is a significant discrepancy (>3%) according to literature for the verification of deterministic codes (Kim, et al., 2017).However, the corresponding RMS equals 1.1 % remaining within target accuracy (<1.5 %).The eigenvalue difference between nTF and the reference solution, [keff Ser- pent2 -keff nTF ], is − 36 pcm, well below the target accuracy (<200 pcm) (Kim, et al., 2017).For the calculation involving the approximated model only of the corner stiffeners (model II) the eigenvalue difference is − 37 pcm.The map of the pin power difference produces a maximum of 5 %, which is still above target accuracy (>3%) (Kim, et al., 2017).The corresponding RMS equals 0.9 %.In the nTF calculation with the semiexplicitly modeled corner stiffeners and water liner (model III), the maximum pin power difference drops to 3.4 % and the RMS equals 1.5 %.The eigenvalue difference is increasing to -74 pcm.Nonetheless, the eigenvalue and the RMS are still in acceptable levels (Kim, et al., 2017).The maximum pin power difference exceeds the target accuracy also in this case, but not significantly.
The strong heterogeneities presented in the nTF solutions with the approximations, at the corners of the assemblies and at the layer of pins neighboring the heavy reflector, have disappeared in the solution of the new model.In Fig. 7 (model III), there is a decrease of the absolute pin power difference at the edge of the assemblies next to the reflector, caused by the model of the corners of the water gap surrounding the active core (see Fig. 8.Ciii).In addition, the pincells around the central guide tube exhibit stronger discrepancies than the surroundings.This is caused by the reduced central guide tube diameter, described in Section 3.1.Finally, it must be pointed out that the power of all Gd pins is overestimated by nTF, similarly to what was observed in Section 4.1.2.In the outer region where nTF generally overestimates the power, the Gd pins appear to have a larger difference.Similarly, in the center where nTF underestimates the power, the Gd pins are not matching the profile of their surroundings and their difference is smaller.However, all these discrepancies are not comparable in magnitude to the effect of the approximated corner stiffeners and water liner, presented in Fig. 7 for model I and model II.
All profiles of Fig. 7 present a tilt, which is not significant in the  solution of model I but is substantially pronounced in the solution of model III.nTF underestimates the power in the center of the core in comparison to Serpent2 and overestimates it in the outer assemblies, flattening thus the power profile, similarly with the single assembly models of Section 4.1.2.The tilt is partially caused by the approximations involved in the nTF reflector modeling approach.CASMO5-VVER employs a similar approach for the modeling of heavy reflectors in full core models.In the framework of the LRT VVER project, CASMO5-VVER was used to simulate a VVER-1000 2D full core (Rostov-2 (OECD/NEA), which was compared with a Serpent2 reference solution (Jang, 2021).
The geometry of the Rostov-2 2D core is identical to the X2 reference model, except the groove region in the heavy reflector (see Fig. 1).Fig. 8 presents the CASMO5-VVER model of Rostov-2 and the pin power comparison of CASMO5-VVER vs Serpent2 which results also in power tilts (Jang, 2021).nTF uses one spectrum for the condensation of cross-sections that corresponds to the fuel region and not the reflector material.This means that the self-shielding in steel is not taken into account accurately, which could be an additional source of discrepancy.The tilt could also be partially attributed to the treatment of anisotropic scattering by the transport correction.This hypothesis is tested in Section 4.2.1.Besides the effect of the reflector model and P0, the tilt in the power profile is a common outcome when comparing deterministic solvers with Monte Carlo codes (Gunow et al., 2019).With the older nTF version, the tilt is partially (model II) or almost fully (model I) compensated due to the heterogeneities caused by the approximation of the water liner and the corner stiffeners, resulting in a lower pin power difference RMS in comparison to the solution of the new model (model III) (see also (Papadionysiou et al., 2021).Despite the observed discrepancies, the new X2 model of nTF presents very good agreement with the reference solution on all accounts.
In order to examine the axial evolution of the core power map, Fig. 9 illustrates the relative difference of the normalized assembly power profile of the 1/6th core as a function of height, together with the relative difference of the axial profile at the same height.The last few axial nodes of the nTF X2 model do not have the same thickness as the reference Serpent2 model thus the assembly power is compared up to ~ 2.6 m.The assembly power profiles of nTF and the reference solution are normalized to 1 for each axial layer.Thus, Fig. 9 presents changes in the power tilt as a function of height, which could signify if the nTF pin power difference from the reference solution, presented in Fig. 7 (model III), is compensated due to axial averaging.The relative difference of the normalized assembly power ranges from − 2.4 % to 3.2 %.It is apparent in Fig. 9 that for each assembly the difference remains almost steady throughout the height of the core, which means that the pin power difference of the new model (see Fig. 7 (model III)) does not result from an error compensation due to axial averaging.The axial profile of the new nTF solution is almost the same with what was presented in (Papadionysiou et al., 2021).The peak power is underestimated in nTF (<1.7 %).The highest relative difference is found in the layer next to the bottom reflector (~5%).This can be attributed to different axial leakage, due to the approximations applied to the nTF model, especially on the axial reflectors and the spacer grids.

nTF X2 HZP solution with P2 anisotropic scattering
In order to test the effect of the P0 transport-corrected scattering treatment, the new nTF model (model III) is simulated with P2 anisotropic scattering.Fig. 10 presents the relative pin power difference for the new nTF model with P2 scattering vs Serpent2.All other parameters of the nTF calculation are identical.In the same figure, the pin power difference of nTF-P0 vs nTF-P2 is also presented.The color scale is the same for both plots.Table 3 summarizes the results of the all nTF X2 critical state calculations presented in this paper.
The eigenvalue difference increases to − 111 pcm.Also, the maximum pin power discrepancy increases to 3.6 %.However, the RMS slightly drops to 1.4 %.The subplot illustrates the effect of P2 on the nTF pin power map.It is obvious in Fig. 10 that the differences close to the reflector are generally smaller in comparison to Fig. 7 (model III).Nonetheless, it must be pointed out that the discrepancies of the Gd pins are more pronounced in comparison to the surrounding profile.The RMS of the relative difference between the nTF solution with P0 and P2 is 0.2 %.The use of the P2 scattering approximation leads to reduction of the pin power discrepancies with respect to Serpent2; however not enough to reduce the tilt significantly.
In order to study the effect of the scattering treatment with another MOC code, CASMO5-VVER is used to simulate Rostov-2 in 2D, based on the model developed in LRT in the framework of (Jang, 2021), with P0, P3 and P5 approximations for the anisotropic scattering.As it was mentioned already in this Section, the geometry of the Rostov-2 2D core is identical to the X2 reference model, except the groove region in the heavy reflector (see Fig. 1).The results are compared with a Serpent2 calculation (100,000 batches × 1,000,000 histories, 1,000 inactive cycles).The reflector model in Serpent2 is built, following the approximations of CASMO5-VVER (see Fig. 8 (Left)), with homogenized pins of water and steel, preserving the real volume of the corresponding material in the actual reflector.The goal of using an approximated reflector model in Serpent2 is to remove the bias introduced in CASMO5-VVER and quantify the effect of the scattering treatment.Fig. 12 illustrates the relative difference of the pin power distribution as calculated by CASMO5-VVER vs Serpent2, [(Nor.Pow.Serpent2 -Nor.Pow.CASMO5-VVER ) / Nor.Pow.Serpent2 ].The pin power uncertainty reaches 0.3 % in Ser-pent2 and the eigenvalue uncertainty 1 pcm.Table 4 presents the Fig. 9. Relative assembly power difference between nTF and Serpent2 vs core height and relative difference of the axial power profile at the same level.
The discrepancies presented between CASMO5-VVER and Serpent2 remain in acceptable levels for all 2D full core cases; however, the maximum difference is higher than the target accuracy with P3 and P5.
When the order of the scattering treatment is increased, the eigenvalue difference is decreased significantly.Nonetheless, the opposite occurs with pin power.This is not the case for the verification of the nTF X2 3D full core model in terms of eigenvalue and RMS (see Table 3), even though similar behavior is observed for the maximum pin power difference.Nonetheless, both MOC solvers present a tilt in the pin power distribution, where the relative power is underestimated in the center of the core and overestimated in the edge, making the overall radial profile more flat.Fig. 11 (P0) proves that the simulated geometry of the heavy reflector is not the main source of the tilt, since the model of the heavy reflector is identical in CASMO5-VVER and Serpent2.However, the generation of cross-sections for the heavy reflector can still be a source of discrepancy.In addition, as it was mentioned before the tilt is an inherent feature of deterministic codes.For both nTF and CASMO5-   VVER, the tilt is not eliminated when the order of the scattering treatment is increased.This finding is consistent with previous studies on the topic (Ryu et al., 2015).Concerning nTF, even though P2 has a positive impact to the solution, the improvement is not enough to justify the resources required for such a calculation.P2 is so computationally expensive that an alternative less memory consuming ray-tracing scheme (Group Major instead of Node Major) was used by nTF, in order to fit the calculation in the cluster.This resulted in severe penalties on the running time, increasing to ~ 22 hrs (300 CPUs for both calculations).

Validation of nTF for the X2 HZP start-up tests
At this stage, the nTF HZP model is optimized and verified, allowing the study to proceed with the validation of said model and nTF against experimental data.As it was mentioned in Section 3, the X2 benchmark specifications contain measurements done during start-up tests: two temperature reactivity coefficients (TRC), the SCRAM worth with a "stuck" rod and the measurement of the differential worth (s-curve) for CR bank #10 (see Fig. 2).Starting with the temperature reactivity coefficients, the different core states of the HZP experiments, described in Table 1, are modeled in nTF using the new X2 model, verified in the previous section, with the same calculation options as for the critical state.The full core is simulated for each case; hence, all calculation costs are similar to the ones described in the previous section.Table 4 presents the temperature reactivity coefficients as simulated by nTF with the new model and the model with approximated corner stiffeners (Papadionysiou et al., 2021), Serpent2 and measured.
According to the benchmark specifications, the reactivity coefficients are calculated as a reactivity difference between two corresponding calculations: where k eff1 and k eff2 are the eigenvalues of the two reactor states before and after the perturbation respectively.Concerning the first temperature reactivity coefficient, the nTF TRC is within 1σ of the measured value and 2σ of the Serpent2 result.For the second TRC, nTF underestimates the value>2σ in comparison to both measured and simulated data.The same discrepancies are observed in (Papadionysiou et al., 2021) with the previous version of nTF, which could not model the corner stiffeners semi-explicitly.Nonetheless, the value calculated with the new model, using semiexplicitly modeled corner stiffeners, for the first TRC is closer to the experimental and Serpent2 outcome.In both cases, nTF under-predicts the reactivity change due to the isothermal change in reactor temperature, as was discussed in (Papadionysiou et al., 2021).The eigenvalues of the Serpent2 calculations for the four TRC points, obtained by the benchmark team, are compared with the corresponding nTF solutions of the new model, repeating the study presented in (Papadionysiou et al., 2021).As depicted in Table 5, the eigenvalues predicted by nTF are in good agreement with the Monte Carlo results (<200 pcm) (Kim, et al., 2017).Such discrepancies prove the capacity of the latest nTF version to model the VVER core in different operational states.Similarly to what was observed with the old model, the inconsistency on the TRC values is caused by the small temperature difference between operational points.The eigenvalue difference between the two points of each TRC is only 2 or 3 pcm.However, when divided by the temperature difference between TRC points (~4.5 • C), the minor 3 pcm discrepancy can shift the TRC value>0.67 pcm/ • C, which is larger than the uncertainty of the experiment and the Serpent2 calculation.
The worth of CR #10 is measured starting with the rod fully inserted (0 %) while all other rods are out (100 %).The stepwise extraction of CR #10 is compensated by adjusting the boron concentration.Temperature and pressure remain the same.There is no experimental uncertainty associated to these measurements in the benchmark specifications.The uncertainty of Serpent2 reaches a maximum of 4 pcm for integral rod worth and ranges from 0.2 to 0.4 pcm/% for differential rod worth (pcm/% means pcm per 1 % of insertion, the same unit as the differential worth).Fig. 12 illustrates the experimental s-curve (sΔρ) and the predictions of Serpent2 and nTF, together with the differential worth (dρ/ dh) of CR #10 as a function of end-of-step position.The comparison of nTF with the experimental and calculated differential worth demonstrates a good overall agreement except for few points in the central part, where the Serpent2 results are significantly lower than the measurement; nTF is even lower.At the same time, the measured data points do not fit into the expected smooth curve, which suggests a large experimental uncertainty, as is stated in the benchmark specifications.The specifications also mention that apart from measurement uncertainties, the possible reason for deviation can be the specified density of the CR materials, which can be bigger in reality.The maximum discrepancy of differential rod worth of nTF from the measurement is 2.0 pcm/% and from Serpent2 1.2 pcm/%, which surpasses the 2σ uncertainty in Serpent2.The total CR worth, as it is calculated by the s-curve, is underestimated by 50 pcm by nTF in comparison to the measurement and 33 pcm in comparison to Serpent2.This is respectively ~ 11 % and 8 % of the total CR worth, which according to the utility perspective is above the limit when nTF is compared to the measurement (>10 %) (Tractabel, , 2021).In order to get a better understanding of the discrepancy, a 2D model of the assembly with CR bank #10 (assembly AU22 in the specifications (Bilodid et al., 2020), with 2.2 % enrichment in all pins) is simulated in nTF and Serpent2, rodded and unrodded, with P0 & P2.The results are summarized in Table 5, [keff Serpent2 -keff nTF ]. nTF underestimates significantly the rod worth with P0, for both materials of the CR, B4C and Dy2O3•TiO2.When P2 is used the discrepancy is reduced to acceptable levels.The underestimation of the rod worth agrees with what is observed in Fig. 12.The total CR worth of bank #10 is also calculated with the P2 approximation for the scattering source.However, there is no significant improvement to the solution (underestimation of 48 pcm from the measurement and 31 pcm from Serpent2).
According to the benchmark specifications, the SCRAM worth with the "stuck" rod and the full SCRAM worth was measured as follows: "First, from the nearly critical state all control rods, except the chosen ''stuck" cluster, were dropped to the lowest position.The worth of SCRAM with one cluster ''stuck" was measured.Second, after about 1 min the ''stuck" rod was dropped down and the full SCRAM worth was measured."(Bilodid et al., 2020).The measured values and the calculated ones by Serpent2 and nTF are presented in Table 6.The nTF SCRAM worth is closer to the experiment than the Serpent2 values.In the first case, nTF overestimates the reactivity by 6.7 %, exceeding the 1σ uncertainty of the experiment, and in the second case by 5.4 %, remaining within 1σ.From the utility perspective this is acceptable since both differences are<10 % (Tractabel, 2021).The results of nTF are closer to the Serpent2 prediction than the experiment, nonetheless far from the 1σ uncertainty of the Monte Carlo code (~180 pcm).According to the benchmark specifications, the systematic deviation between predicted and experimental values of SCRAM worth is a known issue for VVER-1000 units (Bikeev et al., 2018) due to the applied measurement technic (Afanasiev et al., 2014).In addition, the experiment was performed in a dynamic way, whereas the calculations are static, which can also cause discrepancies (for more information refer to the benchmark specifications (Bilodid et al., 2020).The SCRAM cases are also attempted with P2 treatment for the anisotropic scattering.However, the calculations cannot converge.

Conclusions
The first stage of the development of the LRT high fidelity highresolution core solver for VVER geometries includes the verification and validation of the neutronic code and is presented in this paper.The ability of the deterministic code nTF to model the behavior and complex geometric structure of a VVER-1000 is studied standalone, before introducing the effect of the T/H code and the depletion library.To that end, the benchmark specifications for the HZP conditions of a VVER-1000 (X2) are employed.X2 includes experimental and simulated data produced by the benchmark team with Serpent2.In order to analyze separately the various nTF functions required to simulate successfully elements of the VVER core, the verification of the neutronic code is performed for simple VVER configurations, before the 3D full core.nTF is compared with Serpent2 for a range of models, including 2D pincells with P0 and P2 scattering treatment and 2D & 3D assemblies.The code presents good agreement with Serpent2 in all cases; however, specific trends are observed, like the discrepancies for Gd pins and a tilt in the pin power difference of the assembly models, with nTF flattening the power profile.The optimization of the nTF model for the 3D X2 core, presented in this work, is the continuation of (Papadionysiou et al., 2021).The code was further developed by SNU to achieve better accuracy according the observations of (Papadionysiou et al., 2021).The whole study is included in this paper.Different modeling options are tested and the impact of approximations is quantified, specifically for the heavy reflector, which is always a problematic region for VVERs, and the modeling of corner stiffeners.It becomes apparent that, highresolution solvers must allow very limited approximations to achieve target accuracy in local predictions.The optimized nTF 3D full core model is verified for steady state neutronic analysis with the X2 critical state Serpent2 solution.The resulting discrepancies are analyzed and attributed either to the nTF X2 model or to the methods employed in the neutronic code, according also to the trends observed for the simple VVER configurations.Such a discrepancy is the tilt of the pin power difference, which is partly caused by the modeling of the heavy reflector but is also an inherent feature of deterministic codes, as is proven by the comparison of CASMO5-VVER with a Serpent2 2D VVER-1000 full core solution.The assumption of transport corrected isotropic scattering (P0) is shown to have a negligible effect on the performance of MOC solvers compared to P2, P3 or P5 approximations of the scattering source.Finally, the optimized nTF 3D full core model is validated with the HZP start-up tests included in the X2 benchmark, like TRCs, SCRAM tests and CR insertion.The TRCs calculated by nTF are lower than the experimental value.However, this is attributed to the sensitivity of the TRC parameters rather than the accuracy of nTF, as proven by eigenvalue comparisons with the relevant simulated core states.The total CR worth, as it is calculated by nTF, is underestimated by ~ 11 % from the measurement and 8 % from Serpent2, which according to the utility perspective is above the limit when nTF is compared to the measurement (>10 %) The SCRAM worth as calculated by nTF is overestimated by 6.7 % in the case of the "stuck" rod, exceeding the 1σ uncertainty of the experiment, and for the full SCRAM by 5.4 %, remaining within 1σ.From the utility perspective this is acceptable since both differences are<10 %.The same applies for the comparison of the nTF SCRAM worth with the Serpent2 results.The verification and validation of the nTF for an actual VVER-1000 core proves that deterministic codes can achieve high-resolution predictions of equivalent accuracy with Monte Carlo codes.This work continues with the verification of the multiphysics LRT core solver with a coupled code system of similar capabilities.

Fig. 1 .
Fig. 1.Fuel designs, core layout and Serpent2 model of the heavy reflector for the X2 benchmark.

Fig. 2 .
Fig. 2. Positions of the control rod banks in 1/3rd of the reactor core, together with the "stuck" rod.

Fig. 4
Fig.4presents the eigenvalue difference of nTF vs Serpent2, [keff Serpent2 -keff nTF ], for different enriched states with Gd and without, with P0 transport corrected and P2 scattering treatment.The coolant density is set at 0.7629 g/cc, including 1207 ppm of natural Boron (X2 critical state).The temperature of all elements is set at 600 K, and not 554.15K,

Fig. 7
Fig. 7 illustrates the relative difference of the axially averaged normalized pin power map of nTF in comparison to Serpent2, [(Nor.Pow.Serpent2 -Nor.Pow.nTF ) / Nor.Pow.Serpent2 ], as presented in (Papadionysiou et al., 2021) and with the latest version of nTF.Model I(Papadionysiou et al., 2021) employs the approximation of the water liner around the active core with homogeneous pincells (see Fig.3.Cii) and the smearing of the corner stiffener material in the water gap of each fuel assembly (see Fig.3.Bii).In model II(Papadionysiou et al., 2021) the corner stiffeners are again homogenized in the assembly gap, but the water liner is modeled semi-explicitly (see Fig.3.Ciii).Model III corresponds to the latest version of nTF which can model the water liner and the corner stiffeners semi-explicitly (see Fig.8.Biii).In this case, a single full core calculation requires 300 cpus (12 cpus × 25 axial layers) for 15,256 secs (~4 hrs and 14 mins) and 7.44 TB (88 % of the total memory corresponding to 25 nodes).All solutions are produced with transport corrected P0 scattering.As it was discussed in(Papadionysiou et al., 2021), the map of the pin power difference for model I produces its maximum at 7.6 % on the pins neighboring the reflector.This is a significant discrepancy (>3%) according to literature for the verification of deterministic codes(Kim, et al., 2017).However, the corresponding RMS equals 1.1 % remaining within target accuracy (<1.5 %).The eigenvalue difference between nTF and the reference solution, [keff Ser-

Table 1
Operational conditions for all HZP simulated states of the X2 core.

Table 2
Pin power and eigenvalue difference of nTF vs Serpent2 for the X2 3D core.

Table 4
Temperature reactivity coefficients simulated and experimentally measured.

Table 5
Eigenvalue difference of nTF from Serpent2 for the 2D assembly model of CR bank #10.

Table 6
Experimentally measured and calculated SCRAM worth.