VERIFICATION OF THE OpenMC HOMOGENIZED MYRRHA-1.6 CORE MODEL

The OpenMC code is being employed both as a multi-group nodal macroscopic cross-section generator and a reference multi-group Monte Carlo (MGMC) solution. The aim is to do a neutronic benchmark verification study versus a deterministic model (based on the MYRRHA-1.6 core) performed by the PHISICS simulator. MYRRHA, a novel research accelerator driven system concept that is also foreseen to work as a critical configuration, offers a rich opportunity of testing state-of-the art methods for reactor physics analysis due to its strong heterogeneous configuration utilized for both thermal and fast spectra irradiation purposes. The original core configuration representing MYRRHA-1.6 and formed by 169 assemblies, was launched in OpenMC for producing a homogenous nodal model that, when executed in its multi-group Monte Carlo mode, it produced a that differs in almost 500 pcm from the original case. This means that in the future, such approximation should correct the nodal cross-sections to preserve the reaction rates in order to match those ones from the heterogeneous model. Nevertheless, such MGMC mode of operation offered by OpenMC could be exploited in order to verify deterministic core simulators. By inputting the same nodal multi-group cross-section model into the transport solver of the PHISICS toolkit, the neutronic benchmark showed a difference of 171 pcm in eigenvalue while comparing it to its OpenMC MGMC counterpart. Also, local multi-group and energy-integrated nodal profiles of the neutron flux showed a maximum relative difference between methodologies of 15% and 1%, respectively. This means that the MGMC capabilities offered by OpenMC can be employed to verify other deterministic methodologies.


INTRODUCTION
MYRRHA (Multi-purpose hYbrid Research Reactor) is a multi-functional irradiation facility for innovative applications, and is intended to be ready by 2035 [1]. Although MYRRHA is a one-of-a kind subcritical reactor prototype, cooled by a lead-bismuth eutectic (LBE) mixture and driven by a particle accelerator intended mainly for the transmutation proof-of-concept, its 1.6 core version has also been foreseen to work as a critical facility for material irradiation purposes, and even for the production of radioisotopes inside In-Pile Test sections that are connected to other water coolant loops [2]. Thus, this core design brings interesting neutronic characteristics (such as strong flux-spectra gradients) among its different components such as fuel assembly batches, irradiation test sections and control rod banks. At first, a neutron transport analysis of the MYRRHA core 1.6 was based with well-established Monte Carlo (MC) codes such as MCNPX 2.7.0 [3], where a criticality computation of the domain was carried out along with neutronic observables such as, for example, irradiation dosage (i.e. DPA), assembly and pin peaking factors, axial and radial core flux distributions, among others. These characteristics were all computed as a function of irradiation by means of the ALEPH2 code [4], and as stated in references [2,5].
Nowadays, the use of other Monte Carlo codes that are very much reactor-physics oriented and more dedicated towards neutron transport have emerged in the academic community, as it is the case of OpenMC [6]. The 0.10.0 version of the code (which corresponds to the last released stable version) offers, among other modules, the possibility of computing both energy-collapsed (in terms of multi-group bins) and homogenized macroscopic cross-sections over pre-defined volumetric domains, including a full computation of multi-group scattering matrices along the different Legendre order-polynomials that can be defined by the user (OpenMC is also capable of defining scattering matrices by their angular dependence in terms of polar and azimuthal bins). Studies have already been conducted about employing OpenMC for the creation of so-called nodal databases of homogenized and few-group cross-sections (for the most recent studies see references [7][8]), with the aim of being utilized as input parameters in deterministic transport codes that can be employed at low scales (i.e. from the pin cell up to a fuel assembly), or even larger scales (i.e. a group of assemblies configuration or a part of a core).
In this work, OpenMC has been employed as a multi-group cross-section generator at the scale of the assembly level for a configuration that matches the MYRRHA 1.6 core. This means that a global homogenization was carried out for each assembly utilizing the fluxes from the original full core heterogeneous configuration (i.e. no reflective boundary condition assumption for each separate assembly was required). In the end, the aim was to create an input database for the deterministic transport code known as PHISICS [9], so a coarse-mesh model in space could be fastly solved. Regarding the energy mesh, it was found out that due to the unique characteristics of spectral changes in the core, 20 groups were optimal for reducing the statistical uncertainty while tallying the different nodal reaction rates.
Another interesting capability that OpenMC has to offer, is the capacity to perform multi-group Monte Carlo (MGMC) calculations by reading homogenized and multi-group macroscopic cross-sections. Therefore, the objective of this study is two-fold: on the one hand, to demonstrate (once more) the capabilities that OpenMC has in computing MC-based homogenized and few-group macroscopic crosssections by employing a whole core model. Moreover, a MGMC model arranged by the already homogenized cross-sections at the assembly level was created to see how biased are the eigenvalue results between the MGMC model, compared to the original heterogeneous MC scenario. This is performed while studying the excess reactivity of a core at the BOC when all control rods are out. On the other hand, a multi-group core PHISICS model was created in parallel. Therefore, the second objective will be to provide a comparison of computed eigenvalues and spatial characteristics of the flux between the deterministic solution given by PHISICS, and the MGMC prototype. In the end, this type of benchmark exercise offers the opportunity of verifying deterministic codes.

Heterogeneous core
At first, a reference model based on the critical MYRRHA version 1.6 at the BoC [2] was performed in OpenMC. Design to operate at a nominal power of 100 MWth, it consists of 108 (18 batches of 6) driving MOX fuel assemblies (FA's) (ranging between 0 and up to 60 MWd/kgHM), 6 control rods (CR's) banks EPJ Web of Conferences 247, 04002 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124704002 and 3 SCRAM rod bundles, 4 in-pile test sections (IPS's), 4 material testing assemblies and 2 outer ring of reflector/shielding sub-assemblies, surrounded by a steel jacket (as depicted in figure 1). What makes this design interesting are the irradiation "thermal" IPS's and the assemblies for fast spectrum material irradiation. Regarding the thermal IPS's (located at the core periphery), they contain irradiation rigs for the production of Mo-99. This target design is based on in-house technology already in use for many years at our BR2 reactor [10]. Analysis has shown [5] that loading this rig at the periphery (four symmetric positions north-east, north-west, south-east and south-west) gives about 300 Ci/g-U235per irradiation cycle, which is enough to be a worthy successor of BR2. The total height of this model was considered to cover from the above core structure and all the way down to the lower nozzle (accounting for a total height of 467.8 cm). Moreover, an excess reactivity calculated with OpenMC at the BoC of 1278 pcm was obtained (e.g. ݇ = 1.01295 ± 11 ‫,݉ܿ‬ computed with a total of 5݁ histories along 200 active cycles with full vacuum boundary conditions). For verification purposes against the same MCNP6.2 [11] model (which gave a ݇ = 1.01312 ± 9 ‫݉ܿ‬ [12]), the excess reactivity relative difference was found to be of 1.28%. Thus, it can be said that the eigenvalue difference between OpenMC and MCNP lies within statistical uncertainty.

Figure 1. MYRRHA 1.6 based-model in OpenMC
Nevertheless, a smaller-scale model was considered for the homogenization purposes and further results of this particular paper. In the end, only a total of 169 hexagonal core locations were modeled in this occasion (compared to 217 ones from the previous case). Radially, the new domain (depicted in figure 2) ended at the first reflector ring while axially, the dimensions were considered from the upper reflector of the fuel pins (i.e. at 44 cm above the origin), and up to the end of the lower plenum, for a total axial length of 94.5 cm. The reason of having a reduced domain was for statistical tallying purposes; as it will be shown later on, performing a global MC-based homogenization is computationally expensive and relative tallying errors are hard to be lowered for large configurations following this methodology. In the end, since the axial geometrical characteristics were greatly shortened compared to the original model and LBE can be considered good at reflecting neutrons, top and bottom reflective boundary conditions were assumed.

Homogeneous core
Homogenized and energy-collapsed macroscopic cross-sections were calculated for the reduced model; radially, at the assembly level (see figure 3) and axially, for a mesh of 19 nodes (for the node at the very top 4 cm were considered, and for the node at the very bottom, 5.5 cm were used instead. Regarding the distance comprised between -45 and 40 cm, a single node was allocated every 5 cm). For the energymesh, a grid of 20 energy groups was deducted from the 33-group ERANOS mesh [13]. The first 8 groups (at low energies) were collapsed, as well as the last 5 high energy groups from the well-known ERANOS mesh, resulting in the grid described below in Table I.  Multi-group cross-sections and scattering matrices (up to P3) were obtained by means of the convolution in space, solid-angle and energy of the respective reaction rates of interest, while being finally weighted with the angular flux. A very good reference describing the OpenMC methodology for obtaining nodal group cross-sections can be found in [7]. Just as a matter of demonstration, equations (1) and (2) describe these procedure for any Legendre-moment and for any group-to-group scattering reaction matrix.
As an example of the resultant parameters computed with OpenMC, figure 4 offers a plot of the scattering cross-section as a sum of Legendre-polynomials both as a function of mu-bar (i.e. average cosine of the scattering angle) and energy-groups for a node located in the fuel from the first ring of the core (group No. 1 represent the fastest group and from there, it descends towards the lowest energy group No. 20. From hereafter, all the multi-group illustrations will have these order, from high to low energies).

Figure 4. Scattering cross-section as a function of angle (up to P3) and energy (node located at the fuel zone)
A global homogenization technique via Monte Carlo would be constrained by its statistical efficiency. Better statistics were seek to be achieved in this new reduced geometrical model, since the tallying of the nodal cross-sections in the original domain would had resulted extremely computational expensive. All the multi-group cross-sections were computed using two-million particles per cycle along 300 cycles (where the first 50 were discarded).
The 20-group mesh allowed to tally reaction rates at energies between 1MeV and 20 MeV in the water zone of the thermal IPS's efficiently. If energy-bins would had been considered within this energy region, bad statistics or even the lack of scoring for some type of reactions could have occurred. On the other hand, the fact of lumping the first 8 low-energy groups from the ERANOS mesh allowed to reduce the statistical uncertainty associated to the computation of the multi-group cross-sections that are far away from the thermal IPS's (e.g. close to the center of the core). To illustrate these, the relative error of the mean value of the total cross-section (defined as ‫ܧܴ‬ = ‫ܵܺ(ܦܶܵ‬ തതതത ) ܺܵ തതതത ⁄ ) for all the 3211 nodes is depicted in figure 5. Node-numbering goes from the central assembly nodes towards the periphery, following an orderly ring configuration (e.g. the last nodes will be the ones located at the reflector ring assemblies). The highest uncertainty happens to be for the lowest energy group (at a maximum value of 25.8%) at EPJ Web of Conferences 247, 04002 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124704002 nodes that are located at the center of the core and up to the third ring. This is contrary to what is observed at peripheral nodes where, due to the presence of the thermal IPS's, such low energy range produces reaction rates that are strong enough to be accurately tallied. This explains the reason behind choosing the aforementioned energy mesh. For other energy intervals out of group No. 20, relative errors lie below the 10% range.

Figure 5. Relative error (%) of the multi-group total cross-section at all nodes
After processing all the correspondent nodal cross-sections and re-building a homogenized core, the MGMC capabilities offered by OpenMC gave an eigenvalue of ݇ = 1.04902 ± 3 ‫݉ܿ‬ (after 5݁ ଼ histories along 250 active cycles). Compared to the heterogeneous reduced case, the difference between eigenvalues is ∆݇ = 489 ‫݉ܿ‬ ± 1 ‫.݉ܿ‬ [9] is a package toolkit design to provide modern analysis for reactor physics research. Originally developed at Idaho National Laboratory, its neutron transport solver known as INSTANT (Intelligent Nodal and Semi-structured Treatment for Advanced Neutron Transport) [14] has been employed during this work as a deterministic core simulator. Fed with the same multi-group macroscopic cross-sections as from the homogenized OpenMC model, the intention is to validate its capabilities in predicting spatial flux distributions by means of an eigenvalue calculation versus the MGMC model. In the end, even if the response of the MGMC model was almost 500 pcm far from the original case, it is still based in the Monte Carlo method and could be used to validate other methodologies that would employ the exact same input model.

PHISICS (Parallel and Highly Innovative Simulation for INL Code System)
The INSTANT solution taken into account in this work follows the discretization of the neutron transport equation based in the Variational Nodal Method [15]. In space, this corresponds to a hybrid finite element method while in angle, it corresponds to a spherical harmonics expansion. Once the discretization takes place and a matricial system is formed, the GMRES (Generalized Minimal Residuals) is employed for its solution by means of an efficient memory usage (although it offers a slower solution compared to other Krylov-subspace solvers). Thus, the INSTANT input model corresponds to a 3D-hexagonal geometry with a Legendre order up to P3. Same as the previous OpenMC models, axial-reflective and radialvacuum boundary conditions were considered.

PHISICS vs. OpenMC MGMC benchmark results
Table II offers both a recompilation of previous Monte Carlo-based ݇ results, along with a comparison to the one computed with PHISICS. It can be seen that between the OpenMC MGMC model and the deterministic one by PHISICS, 171 pcm separate both eigenvalues. Therefore, it can be said that the computation of integral parameters such as the effective multiplication factor with PHISICS is reliable.   Results from figure 6 to figure 9 show a very good agreement between the deterministic solution and the Monte Carlo-based one. For the energy-integrated flux, a maximum of 1% relative difference from the OpenMC MGMC solution is observed for the fast irradiation assembly located at the center of the core. Therefore, the axial distribution of the energy-integrated flux happens to be very close between both methodologies. Regarding the comparison between multi-group fluxes, a maximum relative difference of 15% occurs axially at the center of the fuel assembly located at the first ring and for the lowest energy group. In fact, this region happens to have the largest statistical uncertainty while tallying the homogenized macroscopic cross-sections.

CONCLUSIONS
Utilizing the novel capabilities of the OpenMC code (both in terms of producing global homogenized and multi-group macroscopic cross-sections while employing the angular flux as a weighting function, as well as being able to run MGMC based calculations) were the main novelty of this paper for the purpose of proposing a methodology for verifying deterministic core simulators. This was applied during the study of the neutronics characteristics of a full core based on the MYRRHA-1.6 model (at the BOC and at excess reactivity).
The first conclusion to be derived happens to be at the comparison between the original heterogeneous model and the homogeneous one, both obtained with OpenMC. Almost 500 pcm separate the solutions when an approximate model based on multi-group and assembly-wise homogenized cross-sections were employed. Even if the angular characteristics of the reactions rates and flux (as a weighting function) are being accounted for during the homogenization process, an important difference is observed in the end between the Monte Carlo original and nodal solutions. Perhaps if a homogenization is done at a lower scale (i.e. at the cell-pin level), a direct comparison between heterogeneous and homogenized models would be closer. This discrepancy occurs because reactions rates are not well preserved at the assembly level by the new homogenized model, creating intra-assembly currents (and in the end, the total leakage) not to be physically consistent with the original heterogeneous case. Thus, a correction technique is needed, such as a super homogenization technique (SPH) [16] in order to update the nodal cross-sections and to make the new system to reproduce the physics from the original one. A future work is foreseen to use the INSTANT SPH module to correct the full core model cross-sections.
Other conclusions arise once the same input model is used by both MGMC and deterministic-based methodologies. Even if such nodal cross-sections at this stage could not represent physical results, they could still be used to assess the capabilities of deterministic core simulators in predicting nodal neutronic characteristics. If we assume that that the MGMC offers the true solution (if good statistics are accounted for), then the previous benchmark offered the opportunity to show that the transport solver of PHISICS does indeed a very good job. Therefore, once the cross-sections will be corrected in the future, we can assure that PHISICS will be able to reproduce the MYRRHA model very properly.