PBMR-400 BENCHMARK SOLUTION OF EXERCISE 1 AND 2 USING THE MOOSE BASED APPLICATIONS: MAMMOTH, PRONGHORN

High temperature gas cooled reactors (HTGR) are a candidate for timely Gen-IV reactor technology deployment because of high technology readiness and walk-away safety. Among HTGRs, pebble bed reactors (PBRs) have attractive features such as low excess reactivity and online refueling. Pebble bed reactors pose unique challenges to analysts and reactor designers such as continuous burnup distribution depending on pebble motion and recirculation, radiative heat transfer across a variety of gas-filled gaps, and long design basis transients such as pressurized and depressurized loss of forced circulation. Modeling and simulation is essential for both the PBR’s safety case and design process. In order to verify and validate the new generation codes the Nuclear Energy Agency (NEA) Data bank provide a set of benchmarks data together with solutions calculated by the participants using the state of the art codes of that time. An important milestone to test the new PBR simulation codes is the OECD NEA PBMR-400 benchmark which includes thermal hydraulic and neutron kinetic standalone exercises as well as coupled exercises and transients scenarios. In this work, the reactor multiphysics code MAMMOTH and the thermal hydraulics code Pronghorn, both developed by the Idaho National Laboratory (INL) within the multiphysics object-oriented simulation environment (MOOSE), have been used to solve Phase 1 exercises 1 and 2 of the PBMR-400 benchmark. The steady state results are in agreement with the other participants’ solutions demonstrating the adequacy of MAMMOTH and Pronghorn for simulating PBRs.


INTRODUCTION
Software quality assurance requires new codes to undergo rigorous benchmarking to demonstrate their capabilities in simulating problems of interest for the design and assessment of pebble bed reactors (PBRs). Verification and Validation (V&V) for PBRs is currently ongoing for the MAM-MOTH [1] reactor multiphysics code and the Pronghorn [2][3][4] thermal-hydraulics code, that are developed at INL under the multiphysics object oriented simulation environment (MOOSE) [5]. MAMMOTH is built on the Rattlesnake radiation multiphysics code [6,7]. To date, MAMMOTH has been validated against the HTR-10 benchmark [8]; while Pronghorn has been validated against the SANA benchmark [3]. In this work, MAMMOTH and Pronghorn are used for solving the PBMR-400 steady-state benchmark exercises [9]. As PBMR-400 is a numerical benchmark, success is measured by comparing the obtained results with results submitted by the other benchmark participants. Participants from 15 different countries contributed 14 different solutions for the Neutron-Kinetic (NK) exercise 1 using ZIRKUS [10], PARCS [11], NEM [12], VSOP [13], PEBBED [14], and 9 independent solutions using THERMIX [13], TINTE [15], WIMS [16], and GCR [17] for the Thermal-Hydraulic (TH) exercise 2. In section 2 a brief description of the used codes and the models are presented, in section 3 the analysis of the results and the comparison with the other participants' solutions are reported.

CODES AND MODELS
As MOOSE based applications, both MAMMOTH and Pronghorn are developed following modern standards and benefit from the plethora of existing multiphysics models and algorithms in the MOOSE ecosystem such as unstructured higher-order meshes, massive parallelization, dimension agnosticism, etc.

The MOOSE based applications MAMMOTH and Pronghorn
MAMMOTH is a MOOSE-based reactor multiphysics application. It is designed as a single entry point for multiphysics analysis of advanced nuclear reactors and allows seamless coupling to the MOOSE ecosystem [18][19][20]. MAMMOTH uses capabilities implemented in the Rattlesnake application for solving the multigroup radiation transport equation. Rattlesnake provides various discretization schemes, including the self-adjoint angular flux (SAAF) and least squares (LS) formulation with continuous finite element method (FEM), and first-order formulation with discontinuous FEM. The two directional variables are discretized using the discrete ordinates method (SN), the spherical harmonics expansion method (PN) or using the diffusion approximation. Rattlesnake solves steady-state, transient and eigenvalue problems with arbitrary order of scattering anisotropy. Pronghorn is a coarse-mesh thermal hydraulics code based on MOOSE. Pronghorn solves the compressible Euler equations (mass, momentum, and energy conservation) in 1, 2, and 3D augmented by closure relations to capture the average interaction effects of a mixture of fluid and solid phases. Pronghorn's treatment of thermal hydraulics is suitable for applications involving complicated fluid-solid structures that can be treated in an average sense. The current closure relationships available in Pronghorn are specific to beds of spherical particles; within the nuclear engineering community, the target application area is thermal-fluid analysis of Pebble Bed Reactors (PBRs) with either gas or liquid coolant. However, porous media methods can be applied to many diverse areas of nuclear engineering, such as prediction of heat transfer in quenched corium heaps or thermal hydraulics of complex heat exchangers.

The PBMR-400 Neutron-Kinetic and Thermal-Hydraulic Models
The benchmark design has been derived from the 400 MW Pebble Bed Modular Reactor (PBMR) demonstration power plant. The PBMR-400 is a modular thermal reactor moderated with graphite and cooled with helium. The 9.6 % uranium dioxide fuel is encapsulated by different shells of pyrolytic graphite and ceramic material called TRISO particle. Almost 15,000 TRISO particles are contained in every one of the almost 452,000 pebbles. The pebbles are contained in an annular region, they enter the core at the top and leave it through the defueling chute at the bottom. At the bottom outlet, pebbles are either recirculated or discarded based on their burnup; a pebble makes an average of six passes through the core. During normal operation 192.7 kg/s of helium at an inlet temperature of 773 K flow trough the pebbles from the top to the bottom reaching roughly 1173 K at the outlet.
Several simplifications have been applied to the benchmark problem geometry to reduce the number of approximations that every participant have to introduce in their models, reducing the discrepancy of the solutions due to the different codes capabilities. A prime example is the inlet of the Helium into the pebble bed. Helium usually enters the upper plenum well above the pebble bed, but due to shortcomings in the geometry representation in V.S.O.P. Helium enters just below the top of the pebble bed. In addition, the central reflector is not cooled in the benchmark model, while the actual PBMR-400 design has an engineered coolant flow from through the central reflector.
As shown in Fig. 1 (left) two axysimmetric 2D models have been defined preserving the geometry of the most relevant core components respectively for the TH and the NK. The phase 1 exercise 1 NK model, Fig. 1 (center), includes the pebble bed, a small portion of reflector above and below the active region, the top cavity and the full radial reflector from the center to the core barrel including the gas gap between side reflector and barrel. All the regions far from the core, where flux solutions may be problematic have been excluded and replaced by void boundary conditions. In total only four material regions exist: the pebble bed (red), the reflector (yellow), the control rod (green), and the void areas (cyan). In order to take into account temperature and spectrum effects the two group cross section library has been generated subdividing the material regions into spectral regions selected in the VSOP model resulting in 190 different materials, 110 of which are used in the Pebble Bed. Since special treatment is required to correctly simulate the neutron streaming effect in the void regions, directional diffusion coefficients for the top cavity and the gas gap are provided in Ref. [9]. For the current study an average of the two directional diffusion coefficient has been used. A Preconditioned Jacobian-Free Newton-Krylov (PJFNK) method initialized by two free power iterations has been used to solve the steady state eigenvalue problem.
As shown in Fig. 1 (right) the TH model for the phase 1 exercise 2 includes more components than the NK model. The arrows indicate the flow path. Helium traverses different porous components including the inlet lower collector, the outlet collector (20% porosity) and the pebble bed (39% porosity). The fluid flow domain is limited to just these components, while the solid conduction equations are solved on the whole domain depicted in Fig. 1 (right) to take into account the heat conducted from the core to the outside boundary trough the reflector, the barrel, and Reactor Pressure Vessel (RPV). Special radiation conduction components have been used to take into account the heat exchange between the two annular cavities around the core barrel. The top, bottom, center, and centerline of the axysimmetric model are considered adiabatic, where the normal heat flux is set to zero. On the right boundary a thermal resistance model taking into account the radiation and the conduction trough the stagnant (benchmark assumption) air gap between the RPV and the Reactor Cavity Cooling System (RCCS) at 293.15 K has been imposed. Finally 192.7 kg/s of inlet mass flow rate, 773 K of inlet helium temperature, and 9 MPa of outlet pressure have been imposed for the fluid domain. The geometry has been discretized using 7904 regular quadrilateral elements paying particular attention to the regions with step changes in porosity. A fully implicit time advancement scheme with an adaptive time step and a PJFNK method have been used to advance the system to the stable steady state solution, i.e. the steady-state solution is obtained via a pseudo-transient solve. The adaptive time step algorithm changes the next time step proportionally to the inverse of the number of non-linear iterations needed to reach convergence for the current time step. Steady-state is assumed if the change of the solution with time becomes negligible.

MAMMOTH Neutron Kinetic Results Comparison
The two-dimensional flux and power density profiles are computed by Rattlesnake and shown in Fig. 2; Rattlesnake's results match physical intuition for PBRs. The fast flux is concentrated in the top part of the pebble bed region where the fuel is cooler, and the thermal flux reaches its maximum in the central graphite reflector close to the pebble bed. The power density mimics the thermal flux profile and reaches its maximum at the top left of the pebble bed, it starts to decrease toward the center of the pebble bed annulus and increases again close to the side reflector thanks to its moderating effect.
In Fig. 3, we demonstrate that global parameters including  The steady-state helium pressure, temperature and velocity streamlines are shown in Fig. 10 (leftcenter). The helium temperature increases once it enters into the pebble bed where the power is generated. The spatial temperature distribution shows a temperature peak of 1, 300 K close to the bottom center of the core. This value is not visible in the axial and radial profiles because they show radial and axial averages of the temperature, respectively, and only a small volume reaches temperatures close to the peak temperature. The streamlines show how most of the helium flows directly into the bed and just a small amount flows into the top cavity. Helium originating from different radial positions in the core mixes in the outlet plenum to produce the average helium outlet temperature. In Fig. 10 (right) the pebble surface temperature and the graphite structures temperature are shown. In normal operation, most of the heat generated in the pebble bed is removed by the helium flow but the central reflector graphite has no active cooling (this is a benchmark simplification in reality active cooling is provided); therefore its temperature is at least 500 K higher than the rest of the graphite in the reactor.
In order to cross-validate the code with the other participants' solutions, an extensive comparison of the results has been performed. In Fig. 11 we show how well Pronghorn matches the other participants' solutions for the following parameters: inlet Helium temperature, outlet Helium temperature, inlet Helium pressure, outlet Helium pressure, total pressure drop, pebble bed pressure drop, and average pebble bed Helium temperature. The values computed by Pronghorn for all  Fig. 12, 13, 14 we show that the Pronghorn solution is situated within the range of the other participants' solutions especially for the axial pressure drop profile. Although the helium temperature is within the range of the solutions, it is situated at the lower end of the range of the submitted results. It shows a similar behavior as the helium average temperature in Fig. 11. The radial (axially averaged) profile of relevant properties has been compared in Fig. 15, 16. In this case the pebble surface and helium temperature profile is slightly different from most of the other participant results, showing a minimum not at the outer circumference of the bed but in the region immediately to the left of it. This behavior is explained by Ref. [21,22]; it is caused by neglecting the effect of flow braiding in the effective thermal conductivity correlation used for the pebble bed. As demonstrated in [21], all the codes used for the phase 1 exercise 2 (DIREKT, THERMIX-KONVEK, TINTE) except for KAERI-MARS-GCR use a "dispersive" helium conductivity instead of the default KTA rule. The dispersive helium conductivity is significantly larger than the thermal conductivity of Helium leading to a flatter temperature profile close to the outer wall. For this reason, a new calculation applying the same code modification has been performed. As shown in Fig. 15, 16 the new calculations labelled "INL-PRONGHORN-CORR" is consistent with the axial and average values and shows the same radial profile as the other participants with the exception of a small negative constant bias. This conductivity should (in theory) compensate for the difference between the real (winding) flow past the pebbles and the idealized (smoothed) flow corresponding to the porous medium approximation. The default correlations used for this study are the official KTA rules as explained in [2]. Using the conductivity correction, the global results are acceptable, but further investigations are needed to understand the validity and the nature of the correction that was causing the temperature profile discrepancy. Initially computed radial Helium temperature profiles showed some discrepancy with the legacy THERMIXlike codes because a modified thermal conductivity that accounted for flow braiding was used in the THERMIX-like codes. A second calculation implementing the same modification has been performed and the new results show a much better agreement with the other participants' solutions demonstrating the code simulation capabilities and flexibility for this type of design. The future research efforts will focus first on developing a coupled TH-NK model validating it with the other participants results for the phase 1 exercise 3 and then use the same model to perform the phase 2 transient scenario simulations.

ACKNOWLEDGEMENTS
This manuscript has been authored by Battelle Energy Alliance, LLC under Contract No. DE-AC07-05ID14517 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.