A Simplified Numerical Benchmark for Pool-Type Sodium Fast Reactors

The current work presents a simplified benchmark for a pool-type Sodium Fast Reactor based on a 2D (r-z) geometry for testing tightly-coupled spatial neutron transport, thermal-hydraulics and thermal-mechanics modeling. The new benchmark is motivated by development of the multi-physics OpenFOAM-based GeN-Foam code at the Laboratory for Reactor Physics and Systems Behaviour at the EPFL, Switzerland, and the FAST code system by the Advanced Nuclear Systems group at the Paul Scherrer Institut, Switzerland. Aiming at the improvement of modeling and code-to-code comparison, the benchmark could prove useful for developers of tightly-coupled multi-physics simulation tools for reactor analysis. The benchmark specification and the solutions obtained with GeN-Foam are presented and discussed in the paper.


INTRODUCTION
A simplified numerical benchmark methodology is proposed in the context of the development of tightly-coupled multi-physics simulation tools for pool-type Sodium Fast Reactor (SFR) analysis. The goal of the present work is to: 1) outline and showcase a numerical benchmark methodology for the assessment of SFR multi-physics code performance through selected quantities of interest; 2) provide a tool to enable periodic regression tests of SFR multi-physics codes under development. Such a benchmark is envisioned to improve both the code development and code-to-code comparison processes. On the code development side, tight coupling in multi-physics codes may results in single-physics code updates having unforeseen consequences on the other physics. Thus, a tool for the assessment of code update effects (i.e. regression testing) by a progressive stepby-step testing of the different involved physics is desirable. On the code-to-code comparison side, a consistent numerical benchmark methodology among institutions involved in code development would prove of great importance.
The present numerical benchmark methodology is set in the broader context of the Horizon-2020 ESFR-SMART [1] project, dedicated to the safety assessment and computational tools improvement for the European SFR reactor concept. This work also stems from the internal activities on multi-physics code development of the Laboratory for Reactor Physics and Systems Behaviour at the EPFL, Switzerland, and the Advanced Nuclear System group at the Paul Scherrer Institut, Switzerland, specifically in relation to the OpenFOAM-based GeN-Foam code and the FAST code system [2].
The current preliminary numerical benchmark consists of: 1) a simplified 2-D pool-type SFR model; 2) a methodology to assess the performance of SFR-dedicated multi-physics code. In the present preliminary development stage, the multi-physics code used for the benchmark is the OpenFOAM-based GeN-Foam code, instrumental to demonstrating the methodology. Nonetheless, the rationale of the work is to provide a methodology that is straightforward to implement in any multi-physics platform of interest.
With regard to the present document, an overview of the coupled physics for SFR calculations is presented first. The computational environment of the employed multi-physics GeN-Foam code is then briefly presented. Subsequently, a description of the employed 2-D reactor geometry is provided. Stemming from such considerations, the benchmark methodology is proposed next. The results obtained by the application of the selected code to the benchmark are presented to further clarify the methodology. Conclusion regarding further developments constitute the final section of the present document.

MULTI-PHYSICS COUPLING OVERVIEW
Multi-physics codes for the modelling of nuclear systems rely on adequately coupled single-physics models. The involved physics consist in thermal-hydraulics, neutron transport and thermal-mechanics. A generic, non-code-specific overview of the different physics and their coupling is provided in Fig. 1, followed by a discussion on their interplay. The discussion is presented in terms of idealized single-physics solvers.

Thermal-hydraulics
The individual contributions from a fluid-dynamics solver and an energy solver can be split, as it provides a further layer of understanding of the overall physics coupling.
Fluid-dynamics Given a temperature field T over a computational domain Ω, a fluid-dynamics solver is dedicated to computing the resulting velocity and pressure fields. The temperature field affects the solver via the fluid molecular viscosity µ and fluid density ρ, which are generally temperature dependent quantities.
Energy Given velocity and pressure fields u and p over a computational domain Ω, an energy solver is dedicated to com- puting the resulting temperature distribution in the system. The velocity field affects the solver via advective terms, while the pressure field affects the solver via time derivative terms which can play a significant role in compressible flow scenarios. The source term for the energy equation consists in a volumetric heat source term, that is somehow (i.e. code-dependent) reconstructed from a volumetric power field P vol that can be provided by a neutron transport solver.

Thermal-mechanics
Given a temperature field within the fuel pins and core structures over a computational domain Ω, a thermal-mechanics solver is dedicated to computing the resulting displacement field d. This displacement field is then employed to: 1) provide feedback to a neutron transport solver via the macroscopic cross sections (generically referred to as Σ), as the displacement field can capture fuel expansion phenomena and thus fuel density reduction; 2) displace the computational domain Ω.

Neutron Transport
Given a set of macroscopic cross sections, collectively referred to as Σ, over a computational domain Ω, a neutron transport solver is dedicated to computing the resulting neutron flux Φ. This forms the basis for the computation of the volumetric power field P vol . Macroscopic cross sections are directly affected by the temperature field via the Doppler effect and indirectly by the fuel thermal expansion and subsequent fuel density change. The former effect is captured by parametrising the cross-sections with respect to temperature, while the latter can be captured by a further parametrisation with respect to the displacement field (e.g. axial fuel expansion, i.e. axial fuel displacement).

COMPUTATIONAL ENVIRONMENT
GeN-Foam is an OpenFOAM [3] based multi-physics code developed at the Laboratory for Reactor Physics and Systems Behaviour (LRS), EPFL and at the Paul Scherrer Institut. While a complete code description can be found in [4], [5], [6], the goal of the current section is to provide an overall introduction. This will be instrumental for further discussions regarding the benchmark model.
From the standpoint of thermal-hydraulics, GeN-Foam is capable of both fine and coarse-mesh single-phase modelling. With regard to the coarse-mesh capabilities, these are based on a porous medium approach [7]. The advantage of such an approach is that it allows for a simplified description of complex systems, while the interaction between the fluid phase and the geometrically unresolved sub-scale structures is captured via adequate closure models. In the fluid-dynamics solver, such an interaction is modelled through a volumetric momentum sink term, which relies on correlations for pressure drops. Conversely, a volumetric momentum source term might be prescribed in certain regions to model a pump. In the energy solver, the interaction is modelled through a volumetric heat source/sink term. Heat sources are treated via a sub-scale fuel pin model that relies on a 1-D heat conduction model for the computation of representative fuel and cladding temperatures within each mesh cell. The source term for the 1-D heat conduction model consists in the volumetric fuel power field P vol . Heat sinks can be modelled by prescribing sub-scale structure properties (e.g. surface area per unit volume, a heat transfer coefficient, sink temperature).
From the standpoint of neutron transport, the solver can model: 1) transient scenarios through either a diffusion equation or an SP3 approach [6]; 2) steady-state scenarios through an eigenvalue approach. The macroscopic cross sections are to be provided as an input, and can be parametrised with respect to: 1) axial fuel expansion; 2) radial core expansion; 3) fuel temperature; 4) coolant density; 5) cladding expansion.
From the standpoint of thermal-mechanics, the purpose of the solver is to predict radial core expansion and axial fuel expansion. In particular, the axial fuel and cladding expansion is pre-dicted by a separate 1-D displacement model within the solver. The temperature field for the computation of the displacement fields is provided by the thermal-hydraulics solver.

BENCHMARK MODEL
In this section, the simplified axial-symmetric 2-D ESFR model on which the benchmark is based is presented. The model consists of three different geometries and meshes, each bound to one of the specific physics. The OpenFOAM computational environment allows in fact to take advantage of standardized meshto-mesh mapping functions, which grant a great deal of flexibility in terms of individual mesh and geometry choices.

Thermal-hydraulics
The computational domain over which the fluid-dynamics and energy sub-solvers operate is reported in Fig. 2. All of the regions are modelled as porous media, with the exception of the simplified sodium pools. In fact, one of the advantages provided by a coarse mesh approach is that it enables a uniform system description of both porous and non-porous regions. In the pump, a uniform, downward volumetric momentum source term can be prescribed. In the heat exchanger, diagrid and radial reflector regions, values for the density, heat capacity and volumetric surface areas of the sub-scale structures are prescribed. This allows for the computation of representative sub-scale structure temperatures, which are mapped to the thermal-mechanics geometry for the computation of the displacement field. A heat transfer coefficient and heat sink temperature are prescribed in the heat exchanger for the modelling of a heat sink term. In the inner and outer core regions, the volumetric surface area of the geometrically unresolved fuel pins is prescribed. Coupled with the 1-D heat conduction model for the fuel pins, this provides a volumetric heat source term to the energy solver. Finally, a set of baffles 1 separates the heat exchanger from the radial reflector, the pump from the diagrid and the radial reflector from both the non-fuel regions and the outer core on one side and from the heat exchanger on the other side. A slip boundary condition is applied to all of the domain boundaries for the velocity field, as well as a zero gradient prescription for temperature and pressure.

Neutron Transport
The computational domain over which the neutron transport solver operates is reported in Fig. 3. The lower and upper axial reflectors are modelled separately as they are characterized by different cross sections. Upper and lower gas plena are modelled as well. The neutron flux computed in the inner and outer core regions provides the basis for the computation of the volumetric fuel power distribution. This is FIGURE 2. Computational domain for thermal-hydraulics. The symmetry axis is in the z direction, starting at the left-most system boundary. It is composed of the following regions: 0) liquid sodium, i.e. simplified hot and cold pools; 1) heat exchanger; 2) non-fuel regions, inclusive of axial reflectors and gas plena; 3) radial reflector; 4) inner core; 5) outer core; 6) diagrid; 7) pump; 8) strongback. The baffles that model the inner and outer core barrels are highlighted in black. mapped onto the thermal-hydraulic mesh to provide a heat source term for the 1-D fuel pin model. The neutron flux is prescribed to vanish at the domain boundaries.

Thermal-mechanics
The computational domain over which the thermalmechanics solver operates is reported in Fig. 4. All of the regions which are not relevant in shaping core radial expansion and axial fuel expansion are collectively described by a so-called "soft structure". It is thermo-mechanically described by a null thermal expansion coefficient and low elastic modulus, so to act as a region in which other elements can expand almost freely. Prescribed thermal expansion coefficients in the diagrid and radial reflector regions allow for the computation of the core radial expansion. In fact, from an SFR design perspective, the assemblies are supported by the diagrid, so that a radial expansion of the diagrid directly results in a radial core expan- sion. The inner and outer core regions are further described by a sub-scale model for the computation of the axial fuel expansion, with prescribed thermal expansion coefficients. The resulting displacement field in the domain is then used to displace the mesh. A null displacement field is prescribed at the upper and lower boundaries of the computational domain, while free expansion is allowed at the radial boundary.

BENCHMARK METHODOLOGY
Stemming from what was discussed up to this point, a numerical benchmark methodology is proposed. In the first phase, steady-state results are obtained by the single physics solvers with no coupling over prescribed input fields. Characteristic quantities specific to each physics, which will be introduced in the present section, are thus obtained. In the second phase, steady-state results are obtained by the solvers coupled in different combinations, and the characteristic quantities re-evaluated. This is done to isolate the coupling contributions. The last stage of this phase involves full coupling of all of the physics at steadystate.
In the remainder of the present section, the actual simulation stages are described in terms of the input fields to be provided and the output quantities to be obtained. As a final remark, the term "core" will be extensively used to refer to both the inner and outer core regions, thus exclusive of the axial and radial reflector regions.

Single-physics Simulations I Fluid-dynamics
The fluid-dynamics sub-solver is employed to compute the steady-state velocity field in the system. Fixed fields and parameters: • uniform temperature field, constant in time; • uniform momentum source in the pump (refer to Fig. 2), constant in time. Output: • total mass flow through core (Ṁ); • average sodium velocity (u) in the core; • total pressure drop across the core (∆p, exclusive of hydrostatic component).

II Energy
The energy sub-solver is employed to compute the steady-state temperature distribution in the system. These are inclusive of average core structure temperatures and fuel temperatures. Fixed fields and parameters: • velocity field obtained in stage I, constant in time; • uniform volumetric fuel power field in the core, constant in time. Output: • maximum, minimum and average coolant (T c ) and fuel (T f ) temperatures in the core; • maximum, minimum and average coolant density (ρ) in the core.

III Neutron Transport
The eigenvalue-based neutron transport solver is employed to compute the steady-state volumetric fuel power distribution in the core. Due to the number of feedbacks from different temperature and displacement fields, one-way coupling is tested for different prescribed fields. Thus, the present stage is split in a number of sub-stages, so that in each sub-stage only one particular field is modified with respect the reference field. The reference fields are those for which the employed set of cross sections were evaluated. A total reactor power of 3.6 GW th is prescribed in all sub-stages for flux normalization. The sub-stages are presented hereby in terms of the input fixed fields, constant in time.
• IIIa − Reference spatially uniform temperature field and null displacement field. • IIIb − Spatially uniform fuel temperature field, increased with respect to the reference fuel temperature field. • IIIc − Spatially uniform coolant density field, decreased with respect to the reference coolant density. • IIId − Displacement field with a constant axial gradient in the core, to model a uniform axial core expansion. • IIIe − Displacement field with a constant radial gradient in the whole system (refer to Fig. 3), null gradient in the axial reflector region, to model a uniform radial core expansion. Output: • effective multiplication factor (k e f f ); IV Thermal-mechanics The thermal-mechanics solver is employed to compute the steady-state axial and radial core expansion (i.e. displacement). Fixed fields and parameters: • temperature field obtained in stage II, constant in time.

Coupled Multi-physics Simulations
V Coupled Thermal-hydraulics The interplay between the temperature, velocity and pressure fields is assessed in this stage. Fixed fields and parameters: • uniform volumetric fuel power field in the core, constant in time, same as employed in stage II; • uniform momentum source in the pump, as employed in stage I, constant in time.
The output quantities consists in those outlined in stages I and II.

VI Coupled Thermal-hydraulics: Buoyancy
The interplay between the temperature, velocity and pressure fields is re-assessed in this stage, with particular regard to buoyancy effects which arise from the prescribed volumetric fuel power density. Fixed fields and parameters: • uniform volumetric fuel power field in the core, reduced in magnitude by a factor 20 with respect to stage II.
The volumetric fuel power field is scaled down as natural circulation cannot provide the same cooling power as a forced convection situation. The output quantities consists in those outlined in stages I and II.

VII Coupled Energy and Neutron Transport
The interplay between the temperature and the volumetric fuel power fields is assessed in this stage. It is recalled that cross sections are parametrised, among others, with respect to temperature and coolant density (fuel density reduction, connected to the axial fuel thermal expansion is captured by a parametrisation with respect to the fuel displacement field, not highlighted in this stage). Fixed fields and parameters: • velocity field obtained in Stage V; • pressure field obtained in Stage V; • null displacement field, constant in time.
The output quantities consists in those outlined in stages II and III.

VIII Thermal-hydraulics and Neutron Transport
Stemming from steps V and VII, the effect of the interplay of fluid-dynamics and energy (i.e. thermal-hydraulics) on neutron transport is assessed. For compressible scenarios, coolant density is in fact influence by the pressure field in addition to the temperature field, hence the need to include fluid-dynamics in the coupling with neutronics. Fixed fields and parameters: • null displacement field, constant in time.
The output quantities consists in those outlined in stages I, II and III.

IX Energy, Neutron Transport and Thermalmechanics
The effect of the interplay between the temperature, displacement and volumetric fuel power fields is assessed. Fixed fields and parameters: FIGURE 5. Temperature (surface field, expressed in K) and superficial (or Darcy) velocity (vector field, expressed in m s −1 ) distributions at steady state obtained at Stage X. The superficial velocity is equal to the real velocity field multiplied by the region porosity, i.e. volume of liquid sodium in the region over total region volume, and arises from the coarse-mesh, porous medium treatment of the system. The higher temperature at the core periphery is due to a higher enrichment of the outer core region.
• velocity field obtained in stage VIII; • pressure field obtained in Stage VIII.
The output quantities consists in those outlined in stages I, II and III.

X Fully Coupled Multi-physics Simulation
A fully coupled thermal-hydraulics, neutron transport and thermalmechanics calculations is performed at steady state. Fixed fields and parameters: • uniform pump momentum source, constant in time.
The output quantities consists in those outlined in all of the single-physics stages.

NOMENCLATURE
u Average sodium velocity in the core; M Sodium mass flow through the core; ∆p Pressure drop across the active core length; T c,max , T c,min , T c Maximum, minimum and average coolant temperature in the core; T f ,max , T f ,min , T f Maximum, minimum and average fuel temperature in the core; ρ max , ρ min , ρ Maximum, minimum and average coolant density in the core; k e f f Effective multiplication factor; d a Maximum axial fuel expansion (relative to initial fuel length); d a Maximum radial core expansion (relative to initial core radius).

RESULTS
In the current section, results obtained by the application of the methodology to the GeN-Foam code are reported. A number of assumptions have been used for this preliminary work: • incompressible sodium flow; • Boussinesq approximation for buoyancy effects, temperature dependence of density neglected otherwise; Furthermore, a one-energy-group modelling of neutron transport is employed, as the focus of the work is on the benchmark methodology rather than an unnecessarily accurate modelling of all of the involved physics. For clarity, the temperature and velocity distributions obtained in stage X are reported in Fig. 5.
In Tab. 1, the results obtained by the various simulation stages discussed in the previous sections are presented. Reference parameters for the ESFR reactor concept can be found in [8], [9].
It should be stated that the benchmark model was built on the referenced data. It should also be stated that the current preliminary work will be updated to better represent ESFR concept as new data becomes available in the context of the Horizon-2020 ESFR-SMART project.
With regard to fluid-dynamic characteristic quantities, it should be stated that the reported velocity consists in the real sodium velocity and not the superficial (or Darcy) velocity, on which the employed porous-medium approach is founded. Furthermore, It is recalled that the reported pressure drop is the one across the active core height, and not the entire core height.
With regard to energy characteristic quantities, it should be stated that the density ρ is obtained via the Boussinesq approximation and thus its effect is only reflected in the gravity momentum source term. This aspect, together with the incompressible flow treatment, limits the impact of velocity-energy coupling on one another, as well as on the other physics. This is reflected in the limited difference between the results of stages VII and VIII, IX and X. As a further remark, the sudden increase of the maximum coolant and fuel temperatures after stage V is due to the spatially non-uniform volumetric fuel power field provided by the neutron transport solver. In fact, the outer core region is characterized by a higher enrichment.
With regard to neutron transport, it can be observed from the results of Stage III reported in Tab. 2 that the general expected behaviour of SFRs in terms of reactivity feedbacks is reproduced, namely negative reactivity feedbacks from fuel temperature increase, axial and radial core expansion. The positive reactivity contribution from coolant density reduction is expected as the increase in neutron leakage from the core is compensated by a greater reduction of neutron absorption in the coolant.
With regard to thermal-mechanics, it should be recalled that the radial core expansion is essentially shaped by the diagrid radial expansion. Since the diagrid supports the assembly lattice and is thus at the core inlet, it will be characterized by considerably lower temperatures than those found in the core, and thus by a smaller expansion when compared to the axial fuel expansion in the core.

CONCLUSIONS
A preliminary simplified benchmark for Sodium Fast Reactors (SFRs) was developed in the broader context of the development of multi-physics tools in cooperation with the Horizon-2020 ESFR-SMART project. In light of the different physics that govern the behaviour of SFRs and their interplay, a benchmark methodology for a simplified 2-D axial-symmetric SFR geometry was devised. In the first place, single-physics are tested individually for prescribed input fields in steady-state scenarios and the simulation outcome characterized against a number of selected quantities of interest. Then, the impact of the coupling between the physics is investigated and the quantities of interest reassessed, still in steady-state scenarios. The proposed methodology is suitable for application to multi-physics code regression tests and for code-to-code verification purposes. Nonetheless, further development stages are foreseen within the context of the proposed methodology. In the near term, the benchmark methodology will be applied to further codes such as the FAST system code developed at the Paul Scherrer Institut. On the benchmark development side, some observations regarding future work can be made. In the first place, it should be noted that most of the multi-physics coupling investigation relies on a set of inputs which are progressively provided by the previous simulation stages. However, a set of inputs that do not depend on the previous methodology stages would be desirable for the ultimate purpose of separately resolving all the couplings. Furthermore, future developments will include the extension of the methodology to transient scenarios. The time-stepping algorithm of a multi-physics code and the intra-time step coupling scheme between the physics represent in fact a further layer of complexity that should be resolved for the stated goal of providing a robust benchmark methodology. In addition to these objectives, future work will be also dedicated to extending the methodology to twophase flow scenarios. As a final remark, it should be stressed that due to the numerical nature of the current benchmark, the validation of the employed codes is out of scope.