Due to their reduced core sizes, Small Modular Reactors (SMRs) typically exhibit larger heterogeneities and gradients in the neutron flux distribution with respect to traditional large Light Water Reactors (LWRs) (Jose' & Reyes, 2017). This can be partly countered with aggressive profiling of the fuel loading; however the significance of the core boundary fuel assemblies (FAs) that exhibit the steepest internal flux gradients is an inescapable fact as the boundary FAs make up a larger proportion of the reactor core in small cores when compared to the large ones. On one hand, the state-of-the-art LWR analysis tools are usually based on an assembly-level (nodal) coupling of the different physics and may face challenges in resolving local intra-assembly details of the coupled problem. On the other hand, advanced low-order methods with rod-level resolved coupled physics are already available (although beyond the state-of-the-art) and trade the fast runtimes of the traditional tools with an increased accuracy in the prediction of rod-level safety parameters. High-fidelity methods for reactor analysis are becoming more and more attractive due to their potential to cope with the more strict requirements of the nuclear safety standards and with the goals of the reactor designers’ and operators’ to reduce safety margins for improved commercial performances. The use of these advanced simulation methodologies at high resolution is possible nowadays due to availability of HPC architectures at a relatively low cost. As a part of this trend, novel high-fidelity methods have been recently developed that couple continuous energy Monte Carlo neutronics tools with sub-channel thermal hydraulics solvers (Sanchez-Espinoza, et al., 2021). These beyond the state-of-the-art solvers have the drawback of the significant runtime, but due to their fewer approximations on the core physics when compared with the low-order tools, they are used to provide reference solutions for those cases where little or no experimental data are available. In fact, one of the challenges in the licensing of new SMRs concepts is the lack of operational data from the concepts that could be used for validating the day-to-day reactor analysis tools used in the deterministic safety analyses.

In this paper, the simulation of a Rod Ejection Accident (REA) in a SMR core is presented and discussed. This work was performed within the framework of the EU H2020 McSAFER Project (Sanchez & al., 2021) whose main objective is to advance the safety research for SMRs by combining safety-relevant thermal-hydraulic experiments with advanced numerical simulations at different spatial resolution levels.

The core design considered in this study is the KSMR core developed at the Karlsruhe Institute of Technology (KIT), based on the SMART concept (Lee & Zee, 2003) and optimized to produce 330 MWth (Alzaben, et al., 2019). The selected transient is assumed to initiate at Hot Zero Power condition, resulting in a rapid reactivity insertion that may challenge the fuel-cladding integrity. In this study two different solutions of the selected transient are presented. Firstly, the transient has been analyzed by means of the standard nodal diffusion two-steps approach using the PARCS code and this will provide the “state-of-the-art-solution” against which any new solver will be compared to in order to demonstrate the new capability regarding the prediction capability of the safety parameters. In a second step, the transient has been modeled with the high-fidelity coupled code Serpent2/SCF allowing performing pin-by-pin transient simulations by coupling neutronics and thermal-hydraulics at the sub-channel level.

The concept of boron-free PWRs that has been proposed several decades ago (Alex, et al., 1986) is nowadays regaining popularity, especially within the framework of water cooled SMRs, due to the availability of advanced designs for burnable absorbers (BA) capable of guarantee enough small excess reactivity without compromising the neutronic and safety performances (Yahya & Kim, 2017) (Kim, et al., 2016) (Merwe & Hah, 2018).

The KSMR core is the results of an innovative boron-free core design developed at KIT that fits into a generic SMART-like SMR cooled by light water and that incorporates well established LWRs technologies with regard to FAs design and type of materials (Alzaben, et al., 2019). It is optimized to produce 330 MWth and loaded with 57 fuel assemblies (FAs) (Figure 1). The soluble-boron-free operation concept has been adopted in order to have a higher moderator temperature coefficient and also to exclude the possibility of a boron dilution accident. The basic FA design is based on a 17 × 17 rods pattern containing 24 guide tubes and a central instrumentation tube. Each FA has a different number of BA rods (either 20 or 24) depending on their position inside the core (Figure 1). These BA are made of Al₂O₃ mixed with B₄C and allow the reduction of the HFP excess reactivity at the beginning of life as well as the core peaking power. The fuel pins are cladded with Zr-4 and their enrichment is lower than 5%. In order for the radial and axial peaking factors to be as smooth as possible, six different types of FAs have been conceived with different radial and axial variation of enrichment and BAs (Figure 1). The central FA is the one with the lowest enrichment and the FAs at the periphery do not contain any BAs but they have dummy rods filled with Al₂O₃. The final optimized core design configuration together with the axial enrichment distribution is shown in Figure 1.

As far as the control rods (CRs), they were designed to provide fast reactivity insertion as well as reactivity adjustment as required for power manoeuvring and transient compensation. Figure 2 shows the CRs configuration that satisfies these goals. The final core design (Figure 1) has 53 rodded FAs divided in two banks, namely, for regulation and shutdown. All the additional details of the KSMR core design can be found in (Alzaben, et al., 2019).

The transient scenario considered in this study consists in a Rod Ejection Accident (REA) (Szilard, et al., 2011), (Abdo, 2014), (Khalil, 2011), (Kim, et al., 2011). This a typical design basis event where a CR assembly is rapidly completely ejected out of the core due to a failure of the Control Rod Mechanism (CRDM) housing structure (OECD/NEA, 2010). In safety analysis, the case studied for this type of scenario is the one corresponding to the ejection of the rod with the highest reactivity worth. In particular, in our analysis this is investigated at the beginning of cycle (BOC) and at Hot Zero Power (HZP) conditions. The initial core configuration for the REA transient is shown in Figure 3 and the initial conditions are summarized in Table 1. In this arrangement, all the boron carbide (B₄C) control rods are fully extracted whereas the remaining ones are fully inserted. This configuration is not in a critical state and as a consequence in the solutions presented in this paper the fission source is accordingly normalized. The reason for choosing this type of condition is because in this way the most reactive CR to be ejected is worth approximately 1.5$, thus enabling the possibility of studying a super-prompt critical transient with an asymmetrical power distribution.

In this study the solution of the REA transient at the nodal level has been performed by means of the PARCS code (Downar, 2017) coupled with the Subchanflow code (U. Imke, U.; V.H. Sanchez and al, 2012).

PARCS is a three-dimensional core solver able to simulate static and dynamic reactor conditions in both Cartesian and Hexagonal geometries using either the multi-group diffusion approximation or the SP3 simplified transport formulation. In particular, for the PARCS simulations carried out in this work, the hybrid ANM/NEM nodal scheme has been used, which is based on the finite difference method in the x and y directions and on the well-established nodal expansion method in the z direction (Lee, et al., 2004).

Subchanflow (SCF) is a four equations (mass, energy, axial and lateral momentum) subchannel thermal-hydraulic code developed at KIT. It can be used for both static and transient simulations and is widely adopted and validated for LWRs applications.

In order to generate the PARCS solution of the REA transient, the first step consists in the transport calculation at the lattice level in order to produce the multi-group homogenized cross sections. In our approach this task has been performed using the Serpent2 Monte Carlo code (J. Leppänen et al., 2015).

Detailed Serpent 2 2D FA models with explicit pin-wise geometrical description have been created taking into consideration the variation of the material compositions as a function of the core height. Subsequently, spatial homogenization and energy-group condensation have been carried-out for each model using reflective boundary conditions (BCs) and a two-group energy structure with the typical cut-off at 0.625 eV. For each model, the homogenized cross-sections have been generated for branch variations originating from the variation of the moderator temperature, the fuel temperature and the CRs position (Table 2). As far as the statistics in the Monte Carlo simulations, for each branch a total of 10⁹ active particle histories divided into 2000 active cycles are used with 500 inactive cycles in order to make sure that the fission source is converged. The main reason for such a large number of active neutron histories is the need to reduce the statistical noise associated with the generation of the group-wise form functions that are subsequently used in PARCS for the pin powers reconstruction.

The group constants for the top and bottom axial reflectors have been produced according to the same two group-structure used for the FAs energy condensation through a Serpent2 3D FA model using reflective BCs in the radial direction and black BCs in the axial direction. In the spatial homogenization, for the top axial reflector the average core outlet temperature is considered while for the bottom part the coolant conditions are set to the core inlet temperature. For the radial reflector, a Serpent2 3D full core model has been considered, in which all regions are set to the average core temperatures for fuel and coolant. As shown in Figure 4, each node in which the radial reflector is divided is adjacent to a FA. Each homogenized region has the same dimension of a FA where the core baffle is homogenized with the coolant around. As a consequence, nine reflector group constants data are generated for each unique position in a quarter of the core (Figure 4). Therefore, with this methodology it is possible to get the detailed spectrum used in the reflector energy-group condensation process.

After having generated the 2-group XSs with the approach described above, these are subsequently converted into the PARCS Macroscopic XS (PMAXS) format via the GenPMAXS code (Ward, et al., 2016) in order to be used for the solution of the REA transient.

The coupling between PARCS and SCF has been implemented at KIT through the ICoCo interface (salome-platform, 2007). The main advantage of this coupling approach is that it does not interfere with the syntax of the codes, i.e., inputs have to be made for PARCS and SCF as usual. A MED mesh, as available in the SALOME platform (salome-platform, 2007) is generated using a MED pre-processor (Garcia, 2019a), for the codes to store the variables fields, and, via ICoCo routines, both codes communicate using get and set functions for fields and time step definition.

The core neutronic model consists radially in one node per FA or reflector element (it may contain sub-nodes) and has an axial discretization of 27 elements: 20 for the active height, 3 for the bottom reflector and 4 for the upper reflector. In SCF only the active height has been modelled with the same radial and axial discretization as the neutronic model. Each FA has been modelled as a single channel with a representative rod considering cross flow between neighbour nodes.

The transient calculation for the REA scenario was performed with an explicit time scheme between neutronics and TH with a time step of 0.0005 s without relaxation factors applied during the feedback fields’ exchange. Given that the core is not initially in a critical configuration, the source has been normalized at the beginning of the transient calculations.

As an alternative to the standard two-steps method usually adopted in reactor safety analysis, an increasing effort is currently being observed in the nuclear community to develop high accurate multi-physics approaches based both on high-fidelity Monte Carlo and deterministic codes coupled with subchannel thermal hydraulic codes. This trend is also driven by large and cheap HPC-clusters for massive parallel computing. In this context and within the framework of the EU H2020 McSAFE Project (Sanchez-Espinoza, et al., 2021), several attempts have been undertaken to develop Monte Carlo-based high-fidelity coupling schemes capable to simulate pin-by-pin full cores’ behaviour both at steady state but also with burnup and time dependent capabilities. In particular, a new beyond state-of-the-art master-slave coupled code was developed based on Serpent2 (as master) and SCF (as slave) (Ferraro et al., 2020a; Ferraro et al., 2020b; Ferraro et al., 2020c). This scheme is strongly based on the use of the advanced capabilities available in the two codes, including the possibility to handle variable temperatures and density profiles in a convenient way for coupled calculations through superimposed meshes (namely, Multiphysics Interface Files (IFC)) and the Serpent2 capability to model time dependent geometry transformations (such as CR movements) as well as prompt and delayed neutrons. The fields’ exchange (TH parameters and fission power) is carried out through the values from this IFC, where proper mapping files are included to consider the FAs’ position. In the developed Serpent2/SCF coupling, on one side a number of so called open door functions are offered by Serpent2 which allow managing the typical issues of a coupled scheme (i.e., initialization, iteration, termination, deallocation, etc.) and on the other side SCF is used as a shared library, where ad–hoc high–level functions are used (Garcia, 2019b). In order to connect these two aspects, specific routines were implemented to properly consider their interaction. Thanks to this type of implementation, both Serpent2 and SCF are independently maintainable and the tool can be easily updated with the next releases of the codes.

As far as for the transient capabilities of the Serpent2/SCF coupling, the standard approach for time-dependent Monte Carlo is used, in which a complex number of different issues need to be taken into consideration. In general, since the position and the energy of any neutron at each interaction point are known, generally these simulations are performed as an external source with a proper time steps discretization and a population control between each time step (Leppänen, 2013). One of the most important challenges is also the large difference in the time scales between the prompt neutrons and the delayed neutrons precursors which need also to be considered. In order to cope with this problem, the approach adopted in Serpent2 (as well as in other MC codes) is to track the precursors’ producing the delayed neutrons instead of directly sampling the delayed neutron from fission (Valtavirta, 2016). As a consequence, two different external sources have to be taken into consideration, namely, the one that can be identified with the so called live neutrons and a second one representative of the precursors. In order to properly generate these two sources a prior criticality calculation need to be performed. Therefore, Serpent2 time-dependent calculations are based on a two-steps approach, as shown in Figure 5.

Since SCF has already time-dependent capabilities by its own (Imke et al., 2012; .; Sanchez and al, 2012), from the TH point of view of the actual coupling approach, the main issue is to properly update the power fields at each time bin. From the neutronic side, in Serpent2 a time-dependent calculation is carried out by tallying the fission power at each time bin, a population control for live neutrons and precursors is applied and the IFC TH fields are updated using the SCF’s results. In order to generate the two initial source distributions, a steady-state criticality calculation is performed in which the positions of live neutrons and delayed neutron precursors are recorded for the converged steady state critical configuration into dedicated files. In the second step the actual transient calculation is performed based on the values stored in the previous step. Moreover, the converged power from the criticality coupled calculation in the first step is used for an initial SCF run in order to consider the right TH fields for the first time step.

In this study the modeling of a reactivity insertion accident for the KSMR boron-free core has been presented. This type of investigation is of relevance as part of the deterministic safety analysis in support to the licensing process for most of the future SMRs. The scenario considered, representative of a super-prompt critical design basis accident, has been simulated at two different fidelity levels, namely, the classical two-steps diffusion based nodal approach and a more accurate high-fidelity approach based on the Serpent2/SCF coupling tool. The exercise represents the first verification of the nodal PARCS/SCF solution against a Serpent2/SCF high-fidelity solution. Based on the comparison of key parameters predicted by both solution approaches it can be concluded that the nodal solution is not conservative. It means that the high-fidelity solution predicts a higher local power peak, higher cladding temperature and lower MDNBR compared to the one of the nodal approach. Both solutions predict similar evolution of the total power, reactivity, cladding temperature but the timing is not the same. The high quality of the PARCS/SCF results with the pin power reconstruction (PPR) approach is achieved through the use of detailed Serpent2 FAs and core models for the fuel’s and reflector’s group constants generation including discontinuity factors and form functions for the PPR-approach. It is worth to highlight that the Serpent2/SCF simulation allows identifying the rod with the highest power and the subchannel with the highest coolant temperature and the local respective safety parameters, which cannot be obtained by a nodal solution. From a safety analysis point of view, the investigations allow to conclude that the integrity of the fuel rod cladding is not challenged during the REA (fuel and cladding peak temperatures are below the melting points). Further investigations including a more realistic description of the irradiated fuel by incorporating a fuel performance code in the coupled system will allow us to quantify the impact of burnup on the final results for end of cycle (EOC) conditions.