Results from a multi-physics numerical benchmark for codes dedicated to molten salt fast reactors

Veriﬁcation and validation of multi-physics codes dedicated to fast-spectrum molten salt reactors (MSR) is a very challenging task. Existing benchmarks are meant for single-physics codes, while experimental data for validation are absent. This is concerning, given the importance numerical simulations have in the development of fast MSR designs. Here, we propose the use of a coupled numerical benchmark specifically designed to assess the physics-coupling capabilities of the aforementioned codes. The benchmark focuses on the speciﬁc characteristics of fast MSRs and features a step-by-step approach, where physical phenomena are gradually coupled to easily identify sources of error. We collect and compare the results obtained during the benchmarking campaign of four multi-physics tools developed within the SAMOFAR project. Results show excellent agreement for all the steps of the benchmark. The benchmark generality and the broad spectrum of results provided constitute a useful tool for the testing and development of similar multi-physics codes. (cid:1) 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license


Introduction
Interest in liquid-fuel nuclear reactor research has increased in the last decades (LeBlanc, 2010), especially after the Generation IV International Forum included the Molten Salt Reactor (MSR) in the list of the new generation reactors aiming at delivering a breakthrough in nuclear electricity production in terms of safety, sustainability, and proliferation resistance (Generation International Forum, 2002).The use of a molten salt both as fuel and coolant leads to unique physics phenomena in MSRs: internal heat generation in the coolant; thermal feedback induced by fuel expansion; transport of delayed neutron precursors; and, thus, stronger coupling between neutronics and thermal-hydraulics.These features are absent in traditional solid-fuel reactors, therefore classical codes used in the nuclear community are unsuitable for simulating MSRs behavior ''as they are".Even when some efforts have been made to provide these codes with the capabilities required by MSR modeling (e.g., Zanetti et al., 2015) the code structure (usually developed in the '60s/'70s when strong computational capabilities were not available), along with the often strong modeling assumptions taken (e.g., a 1D modeling) make them unappealing for the study of MSR physics.In addition, the verification and validation campaign that makes these legacy codes suitable for the study of conventional reactors, cannot be considered applicable to the study of MSRs due to the peculiar reactor characteristics.
A comprehensive verification and validation program is necessary in order to increase the confidence in these research codes and bring them close to the level of industrial ones.However, this is a very challenging task, especially when considering tools targeting non-moderated MSR designs, as the Molten Salt Fast Reactor (MSFR) (Allibert et al., 2016;Brovchenko et al., 2019).Validation is in fact not possible, as experimental data are available only as

Contents lists available at ScienceDirect
Annals of Nuclear Energy j o u r n a l h o m e p a g e : w w w .e l s e v i e r .c o m / l o c a t e / a n u c e n e result of the operation of the Molten Salt Reactor Experiment (Haubenreich and Engel, 1970).This was a thermal reactor with graphite moderator and salt channels, characteristics which make it vastly different from the current fast MSR designs.Verification of single-physics codes can be achieved for many applications (see Ghia et al., 1982;Ethier and Steinman, 1994;Botella and Peyret, 1998;Pautz, 2001, for examples in the field of CFD and neutronics), but verification of multi-physics codes remains a difficult task (De Oliveira and Mikityuk, 2018), especially when the coupled problem is solved by iterating different solvers, each one addressing a specific set of equations.
To help overcome the last issues, in this work we propose the use of a numerical benchmark for multi-physics tools dedicated to fast-spectrum molten salt reactors.It was originally developed at LPSC/CNRS-Grenoble (Aufiero, 2015;Laureau, 2015;Aufiero et al., 2018), and for this it is referred to as ''CNRS benchmark" in the following.
The CNRS benchmark is divided into several steps, in which stationary or transient simulations of an MSFR-like system are prescribed.The coupling between the various physics phenomena is gradually introduced, so that it is easy to assess the physics coupling capabilities of the codes, and to highlight possible sources of discrepancy and fix the underlying errors.In this perspective, complex phenomena to model as turbulence or Doppler feedback are also absent.Other results available in literature could be used as benchmark for multi-physics codes dedicated to MSRs.For example, Fiorina et al. (2014) and Brovchenko et al. (2019) reported code-to-code comparisons of MSFR steady and transient calculations performed during the Euratom FP7 EVOL project.However, the complexity of full-reactor calculations hinders the understanding of where possible discrepancies in the results stem from, whether only from specific modeling/data choices (cross section libraries, turbulence models, etc.) or real errors in the code packages.In this perspective, the test case presented here is a more powerful tool for benchmark purposes, especially when the assessed multi-physics code is at early stages of development.
In the CNRS benchmark, the fluid problem assumes the salt to be incompressible, but no prescription is given for the neutronics model, thus making this benchmark suitable to virtually any multi-physics tool developed for fast MSRs.This test case can even be useful to benchmark tools intended for thermal-spectrum MSRs, but no assessment of the capability of such codes to model the peculiar physics phenomena related to the moderator is possible.
The purpose of this work is threefold.First, we present for the first time a complete version of the CNRS benchmark.Then, we collect all results obtained during the benchmarking campaign of four multi-physics tools developed in the framework of the Horizon 2020 Euratom SAMOFAR project (http://samofar.eu/),and we compare them in order to assess the capability of the codes to correctly model the peculiar physics phenomena in a molten salt fast reactor.Finally, and most importantly, the results and the benchmark description itself are meant to be an useful tool for any other institution interested in testing its own multi-physics tool or in developing a new one.In this perspective, the coupled tools assessed here include both diffusion and transport codes, and both finite volume and finite element solvers, thus providing the interested reader with a broad spectrum of results to compare with.
The remainder of the paper is organized as follows.In Section 2, a general description of the benchmark is provided and all the data useful to perform calculations are reported.In Section 3, the benchmark phases and steps are described in detail, together with the definition of the output quantities of interest.Section 4 is devoted to the description of the code packages involved in the benchmark campaign and used to obtain the results collected and analyzed in Section 5. Finally, our conclusions are discussed in Section 6.

General description
The multi-physics benchmark described in this work features a molten salt system whose characteristics (neutron spectrum, strong temperature feedback, salt composition, precursors movement) make it a simple representation of the MSFR.In order to assess the coupling capabilities of the multi-physics tools in a systematic way, facilitating the identification of possible sources of inconsistency in the results, the benchmark is structured into three phases, sub-divided into steps, with Phase 0 as single physics verification, Phase 1 as steady-state coupling, and Phase 2 as time dependent coupling.
The salt flow is considered to be incompressible and laminar, and buoyancy is modeled adopting the Boussinesq approximation.However, no prescription is given regarding the neutronics model.Sources of complexity as turbulence or 3D geometry are explicitly avoided, as these are features other benchmarks deal with.Decay heat is neglected as well.For these reasons, the benchmark is fairly general and suitable to virtually any multi-physics tool dedicated to MSRs, even at early stages of development.

Geometry and boundary conditions
Fig. 1 shows the domain of the problem.It is a 2 m by 2 m cavity filled with molten salt at the initial temperature of 900 K. Pointwise comparison of the observables in the different steps of the benchmark is performed along the cavity centerlines AA 0 and BB 0 .
The domain is treated as a homogeneous, bare reactor.Therefore, the standard vacuum conditions are applied for the neutron flux to each boundary, together with a homogeneous Neumann condition for the delayed neutron precursors.A no-slip boundary condition, with zero fluid velocity, is applied to all walls except the top lid, which moves at U lid .All walls are adiabatic, and salt cooling is simulated with a volumetric heat sink 1 : q 000 ðrÞ ¼ cðT ext À TðrÞÞ; where TðrÞ is the salt temperature at point r, T ext ¼ 900 K, and the volumetric heat transfer coefficient, c, is uniform throughout the entire domain.All steady-state steps of the benchmark are criticality eigenvalue problems, in which the reactor power is normalized to the reference power, P.

Input data
The fuel salt is a LiF-BeF 2 -UF 4 , whose molar composition is reported in Table 1.Its fluid properties are considered constant with temperature and uniform in space.They are reported in Table 2. Composition and properties were optimized to highlight the impact of possible errors in physics coupling strategies.For example, the Schmidt number value is a trade off between two opposite needs: to limit mixing of delayed neutron precursors, thus magnifying the impact of the flow field on their distribution, and to cover the impact of the numerical diffusivity of a specific discretization scheme on the results.
Regarding nuclear data, the JEFF-3.1 library (Koning et al., 2006) at T ref ¼ 900 K is prescribed 2 .Multi-group codes can work with neu-1 Even if the domain is 2D, all quantities are reported with standard 3D units in this work, in order to improve clarity.For the same reason, the heat sink and the heat transfer coefficient are indicated as ''volumetric".For an exact unit match, one should consider the domain to have a depth of 1 m along a pseudo z-axis.However, this is irrelevant from the numerical point of view, since all quantities do not depend on z.
2 For convenience, the same reference temperature is chosen for the Boussinesq approximation.
tronics data generated with Serpent (Leppänen et al., 2015) and condensed into six groups, using the spatial-homogenization and groupcollapsing models described by Leppänen et al. (2016).The complete set of multi-group neutronics data can be found in Appendix A. As the reactor system considered here has a similar spectrum to the one of the full MSFR, the same six energy groups used in previous works to model the MSFR were chosen.This energy-group structure, reported in Table 3, was tested in previous works (Fiorina et al., 2012;Fiorina et al., 2013) against ERANOS (using up to 1968 energy groups) and Serpent, and it proved capable of reproducing with sufficient accuracy the overall MSFR spectrum and relevant neutronics parameters.
In this benchmark, the Doppler effect on microscopic cross sections is negligible; this avoids biases arising from the differences in treatment of this complex effect by each code.Consequently, cross sections and diffusion coefficients are affected only by the saltexpansion feedback and scale with the temperature according to where q fuel is the fuel density.
Finally, delayed neutron precursors are grouped into eight families, whose decay constants and fractions are reported in Appendix A.

Phases and steps of the benchmark
In this section, we describe each phase and step of the benchmark.Each step is schematically described in terms of goal, required input, and the output quantities of interest (the observables).

Phase 0: Single physics benchmark
In this phase, as a preliminary benchmark of the codes, steadystate simulations are carried out without any coupling between the different physics.

Step 0.1: Velocity field
The solution of steady-state incompressible flow is studied.This step is aimed at testing the capability of the codes to get a correct velocity field, which is mandatory to obtain consistent multiphysics coupling.

Observables
Velocity components along centerlines AA 0 and BB 0 .

Step 0.2: Neutronics
The neutronics criticality eigenvalue problem is studied with static and isothermal fuel, given the single physics purpose of this phase.The aim is to verify the neutronics solutions of the codes in simple, static-fuel conditions.Before proceeding to further comparisons, agreement between codes must be verified for the fission source distribution and the estimation of the effective multiplication factor.Minor differences in the effective multiplication factor, k eff , might arise from different neutronics models and approximations adopted and are considered acceptable considering the purpose of the benchmark.The reactor power, P, is set to normalize the neutron flux.Input T ¼ 900 K; and P ¼ 1 GW.

Observables
Fission rate density, R E R f UdE, along AA 0 ; and Reactivity, q.
3.1.3.Step 0.3: Temperature The passive scalar transport capability of the codes is assessed independently from the solution of the fluid flow and the neutron-

Observables
Temperature distribution along centerlines AA 0 and BB 0 .

Phase 1: Steady-state coupling
In this phase, steady-state solutions are found gradually coupling the several physics phenomena characterizing the molten salt fast system.

Step 1.1: Circulating fuel
The criticality eigenvalue problem is solved in the presence of fuel motion, with a fixed velocity field and uniform fuel temperature.This step is aimed at assessing a correct evaluation of the effects of the fluid flow on neutronics.In particular, the consistency of reactivity loss due to fuel motion is verified.

Input
Fixed flow field from Step 0.1; T ¼ 900 K; and P ¼ 1 GW.

Observables
Delayed neutron source, P i k i C i , along AA 0 and BB 0 ; and Reactivity change from Step 0.2, q À q s 0:2 .

Step 1.2: Power coupling
The coupling between neutronics and thermal-hydraulics thermal-hydraulics is investigated in the simple case of a fixed velocity field.The goal is to analyze the flux shape deformation caused by the non-uniform fuel temperature field and the effects of this deformation, in turn, on the temperature distribution.The capability of reproducing correctly this coupling, via the fuel density feedback, is assessed in this step.

Input
Fixed flow field from Step 0.1;

Observables
Temperature distribution along AA 0 and BB 0 ; Reactivity change from Step 1.1, q À q s 1:1 ; and Change of fission rate density along AA 0 and BB 0 with respect to the solution obtained at Step 0.2,

Step 1.3: Buoyancy
The full multi-physics problem is analyzed in the easiest conditions, without an external source of momentum.The fuel flow is driven by buoyancy effects caused by the temperature gradients.The main objective is to test the codes' capability to predict the correct velocity field induced by the fission heat source and the correct reactivity introduction due to the movement of precursors.Discrepancies in the results of this step can be considered mainly related to buoyancy effects or to its modeling, as most of the multi-physics coupling issues for steady-state eigenvalue problems are studied in Step 1.1 and Step 1.2.

Observables
Velocity components along AA 0 and BB 0 ; Temperature distribution along AA 0 and BB 0 ; Delayed neutron source along AA 0 and BB 0 ; and Reactivity change from Step 0.2, q À q s 0:2 .

Step 1.4: Full coupling
This step is thought as representative of realistic reactor simulations in terms of physical phenomena, because it involves the full solution of the multi-physics problem.All the phenomena analyzed separately in the previous steps are introduced here simultaneously: (i) external momentum source (i.e., non-null lid velocity); (ii) buoyancy effects; (iii) delayed neutron precursors motion; and (iv) flux deformation due to an asymmetric temperature distribution. Input P variable in the range ½0; 1 GW with a step of 0.2 GW; and U lid variable in the range ½0; 0:5 m s À1 , with a step of 0.1 m s À1 .

Observables
Reactivity change from Step 0.2, q À q s 0:2 , as a function of P and U lid .

Phase 2: Time dependent coupling
The transient behavior of the reactor is simulated in this last phase, taking into account the full coupling between the physics phenomena.

Step 2.1: Forced convection transient
In order to perform a general assessment, no standard transient is simulated (e.g., reactivity insertion, loss of heat sink), since the comparison would be limited to the specific transient and its characteristic time constant and physical phenomena.The goal is instead to be as general as possible, so we assess the response of the system (in terms of gain and phase shift) to a perturbation in the frequency domain.This approach is limited by the possible influence of the perturbation amplitude, since a linear analysis tool is used on a modeled system characterized by non-linearities.This drawback is overcome using a fixed, small amplitude variation.This type of analysis allows characterizing the outcomes of the multi-physics codes as a function of the excitation frequency, observing in a synthetic, yet quantitative way the different physics at work and their associated dynamics.In this way, in case of discrepancy among the codes, it is quite immediate to understand in which phenomenon the code is failing to reproduce the correct behavior.As a perturbation, we study how an oscillation on the heat transfer coefficient affects the power production.

Input (initial condition)
Steady-state solution from Step 1.4 with U lid ¼ 0:5 m s À1 and P ¼ 1:0 GW; and

Transient description
Starting from the initial condition, the volumetric heat transfer coefficient c is uniformly perturbed according to a sine wave of amplitude 10% and frequency f 2 ½0:0125; 0:025; 0:05; 0:1; 0:2; 0:4; 0:8 Hz.The variation in the cooling of the salt leads to a sinusoidal power trend induced by the negative density feedback coefficient. Observables Power gain and phase-shift as a function of the perturbation frequency.The power gain is defined as where the denominator corresponds to the amplitude of the heat transfer coefficient sine wave (10%), P max is the maximum power, and P avg is the time-averaged power, corresponding to the initial power of 1 GW.

Code packages
In this section, we briefly describe the multi-physics tools participating to the benchmark campaign, and we report the specific choices adopted for the calculations by each partner (e.g., mesh, time step, time discretization scheme).The code packages were developed during the SAMOFAR project at the Centre national de la recherche scientifique-Grenoble (indicated with CNRS from here on), the Politecnico di Milano (PoliMi), the Paul Scherrer Institute (PSI) and the Delft University of Technology (TUD).

CNRS code
A multi-physics (neutronics, thermal-hydraulics, and thermomechanics) tool to perform nuclear reactor studies has been developed at the LPSC/CNRS-Grenoble.This tool uses the C++ finite volume libraries of OpenFOAM (OpenFOAM, 2013) to implement various neutronics models such as diffusion and Simplified P N of first (SP 1 ) and third order (SP 3 ).In addition, the multi-physics tool has been internally coupled to the neutronics code Serpent to allow performing either steady Monte Carlo simulations or transient simulations through a quasi-static approach.
For the present analysis, the SP 1 and SP 3 models were chosen.The thermal-hydraulics model used in all calculations considers an incompressible single-phase flow with the Boussinesq approximation for the mass and momentum balance (Navier-Stokes) and energy conservation solved using OpenFOAM's PIMPLE algorithm (Versteeg and Malalasekera, 2007;Holzmann, 2019).Fig. 2 shows the coupling scheme among the solvers.Within a time step, the neutronics problem is solved first, then the CFD solver is called.A fixed-point iteration process is adopted to resolve the nonlinearities coupling the two problems.More details on this multiphysics tool and the models implemented therein can be found in Blanco et al. (2020).

Specifications for benchmark calculations
In this study, the cavity domain was discretized using a structured, non-uniform mesh (thinner near the walls) containing 200 by 200 cells.An adaptive time step was employed with a Courant number limit of 0.8.In all simulations, the interpolation schemes chosen for the integration of the time derivative, diffusive, and convective terms were Implicit Euler, Gauss linear, and Gauss upwind respectively.

PoliMi code
PoliMi has developed and implemented a multi-physics solver for fast molten salt systems using the finite volume OpenFOAM library.The code implements several thermal-hydraulics and neutronics models (the interested reader is referred to Cervi et al., 2017;Cervi et al., 2018;Cervi et al., 2019;Cervi et al., 2019;Cervi et al., 2019, for more details).For this benchmark, a singlephase, incompressible model was selected for thermal-hydraulics (adopting the Boussinesq approximation to describe buoyancy), and a multi-group diffusion model was chosen for neutronics (Cervi et al., 2019).
The structure and coupling strategy of the multi-physics tool are described in Fig. 3.At the beginning of the time step, the thermal-hydraulics cycle solves for fuel velocity, pressure, and temperature.Picard iterations are performed until convergence of the thermal-hydraulics problem.Then, the neutronics cycle begins, solving for neutron flux and delayed neutron precursors (and decay heat, if necessary).The volumetric power source field is then updated, and the energy equation is solved again.Once the new temperature and density fields of the fuel are calculated, cross sections are updated, and the cycle is repeated with Picard iterations until convergence.Before progressing to the next time step, outer iterations between the thermal-hydraulics and the neutronics cycles can be performed, to resolve the non-linearities between the respective physics.In these cycles, fuel temperature and density are passed from thermalhydraulics to neutronics, in order to evaluate cross sections, while the power density distribution is passed from the neutron-Fig.2. Coupling strategy of the CNRS solver.The neutronics problem is solved first, followed by the thermal-hydraulics one.Fixed-point iterations between the two solvers are performed within a time step to resolve the coupling non-linearities.
ics to the thermal-hydraulics solver, in order to update the temperature.

Specifications for benchmark calculations
All simulations were performed with a 400 by 400 uniform structured mesh, ensuring that results were not affected by further refinement of the grid.In addition, second order schemes were selected to compute the inter-cell fluxes for each physical quantity.The time scheme chosen for transient calculation was the first order implicit Euler scheme, using the time step sizes reported in Table 4.

PSI code
The GeN-Foam multi-physics tool was used by PSI in this study.As an OpenFOAM-based code, it relies on the Finite Volume Method to discretize the system of partial differential equations.Relevant for this benchmark are the solution of the incompressible Navier-Stokes equations and the neutron diffusion equation with delayed neutron precursors transport using a fixed-point iteration scheme to couple the system.All the details on the coupling scheme are described in Fiorina et al. (2015).

Specifications for benchmark calculations 3
The mesh used consisted of a 200 by 200 structured grid gradually refined towards the walls, which was found to give meshindependent results for the problem in the regions and fields of interest.Convective and diffusive terms used upwind biased and centered schemes respectively for spatial discretization and the implicit backward Euler scheme for time discretization.

TU Delft code
The TU Delft multi-physics package consists of two in-house codes: PHANTOM-S N : An S N solver for the multi-group Boltzmann equation coupled with the transport equations of the delayed neutron precursors; DGFlows: A parallel solver for the incompressible Navier-Stokes equations.
Both codes are based on a Discontinuous Galerkin Finite Element Method (DG-FEM) and implicit Backward Differentiation Formulae (BDF) for space and time discretization respectively.
Fig. 4 illustrates the coupling scheme together with the data exchanged between the codes.PHANTOM-S N takes care of correcting cross sections with the density feedback starting from the information on the average temperature on each element received from DGFlows.The neutronics problem is then solved using the velocity field (u) from DGFlows as another input for the precursors equation.Finally, the fission power density (P fiss ) is transferred to the CFD code.In steady-state simulations, the codes are iterated until convergence.In transient simulations, DGFlows is called first, then, after the completion of a time step, PHANTOM-S N is called to handle the neutronics physics.The reader is referred to Tiberga et al. (2019) for a more detailed description of the code package.

Specifications for benchmark calculations
All benchmark simulations were performed on a 50 by 50 uniform structured mesh, adopting a polynomial order P ¼ 2 to approximate the velocity field and P ¼ 1 for all the other quantities.These options proved to ensure mesh independent results.Neutronics simulations were carried out both with order S 2 , qualitatively close to diffusion, thus obtaining results more fairly comparable to those obtained by the other diffusion codes, and S 6 , to get ''full-transport" solutions, able to take into account of the full anisotropy of scattering up to order three.Transient calculations were performed taking a time step equal to 1=200 of the perturbation period and using a BDF2 time discretization scheme.

Results
In this section, we report and compare the results obtained by each partner performing the steps of the benchmark described in Section 3. Whenever an observable is required along the centerlines AA 0 and/or BB 0 , we report a table with its values at 9 equidistant points along the line.The complete 1D profiles (201 equidistant points along the line of interest) can be downloaded from Tiberga et al. (2020).One plot with the complete profile of each observable is reported as well per benchmark step, in order to give a better idea to the reader on how the results (qualitatively) compare along the whole centerline.The plots of some observables in the entire domain for some steps of the benchmark are reported in Appendix B for the sake of completeness, but again they allow for a qualitative comparison of the results only.CNRS performed calculations with two angular discretizations, as described in Section 4.1, indicated with ''CNRS-SP 1 " and ''CNRS-SP 3 ".The same was done by TUD (see Section 4.4), and their results are indicated with ''TUD-S 2 " and ''TUD-S 6 ".
In absence of a reference solution, we report and discuss the discrepancy among the values obtained by each partner.At each point r i , we collect N c values (one per code) of a quantity Q.We define the average value as where N p ¼ 201 is the number of sampling points of the quantity Q.
The average discrepancy is then calculated as ¼ 1  5 reports the velocity components along the centerlines, showing good agreement between the codes.This is not surprising, as this step reproduces a well-known benchmark in the CFD community (Ghia et al., 1982;Botella and Peyret, 1998).The average relative discrepancy is below 0.35% for velocity profiles along AA 0 and rises up to 0.8% for the BB 0 profiles.The profile of u x along BB 0 can be seen in Fig. 5: the agreement is excellent.At the boundaries, the TUD velocity is not exactly null, due to the DG-FEM discretization adopted, which requires a weak imposition of the boundary conditions.It is worth noticing that the second polynomial order chosen by TUD for the velocity field leads to excellent results despite the coarse mesh selected.

Step 0.2: Neutronics
The profile of the fission rate density along AA 0 is reported in Table 6 and in Fig. 6. Results are very similar considering the diffusion codes (PoliMi and PSI) and the CNRS-SP 1 and TUD-S 2 calculations.This is expected, since a diffusion model with P 1 scatter is theoretically equivalent to the SP 1 and S 2 discretizations.Differences increase when considering the higher order transport discretizations.However, the average relative discrepancy is below 0.3%; this is expected in a large homogeneous geometry like the one of this benchmark, where using a more advanced transport model does not lead to great differences in the results.The values obtained by each partner scatter the most at the boundaries, but this is expected due to the difference in the imposition of the vacuum boundary conditions between diffusion and transport models.Finally, the system reactivity is shown in Table 7. Diffusion and CNRS-SP 1 results differ only by 10 pcm, the S 2 model by TUD is only around 70 pcm off the diffusion ones; not surprisingly, the more advanced transport models lead to different results, but only slightly (around 150 pcm of difference from the diffusion ones).

Step 0.3: Temperature
Table 8 reports the temperature profiles along AA 0 and BB 0 .The agreement among the codes is excellent, with an average discrepancy of around 0.1% for both profiles.This is confirmed by the complete plot along BB 0 shown in Fig. 7. Once more, the maximum disagreement between the codes is to be found at the boundaries and can be explained by the different treatment of boundary conditions by the DG-FEM discretization adopted by TUD and the differences in fission density described in Section 5.1.2.Table 8 Step 0.3 -Temperature distribution along centerlines AA 0 and BB 0 .The profiles of the delayed neutron precursors source along the cavity centerlines are shown in Table 9. Very good agreement can be noticed among all codes, with an average discrepancy of 0.35% for the profile along AA 0 and of 0.3% for the one along BB 0 (shown in Fig. 8 as well).Results scatter again at the boundaries, but, taking into account the differences in velocity and fission density analyzed in Sections 5.1.1 and 5.1.2,this is acceptable.
Finally, Table 10 reports the value of reactivity insertion with respect to the static case (Step 0.2).The agreement between the codes is excellent, with a minimal difference between the TUD-S 6 and the others, which is expected, given the more advanced transport model.It can be deduced that all the assessed multi-physics tools can correctly reproduce the effect induced on the effective multiplication factor by the transport of delayed neutron precursors inside the cavity.

Step 1.2: Power coupling
The profiles of the temperature along AA 0 and BB 0 are reported in Table 11.Results are very similar to the ones obtained for Step 0.3 (see Section 5.1.3),with an excellent agreement among all numerical tools.In fact, the average discrepancy is 0.09% for both profiles.Table 11 also reports the values of the change of fission rate density along AA 0 and BB 0 with respect to the solution obtained at Step 0.2.Differences are more considerable here, as can be noticed also in the complete plot shown in Fig. 9. Nevertheless, the average discrepancy is lower than 1.6% for both profiles.PSI and TUD results are those affected by the highest discrepancy (around 2.5% for the former, 2% for the latter).For PSI, this is to be ascribed to the more refined mesh towards the boundary used combined with the vacuum conditions imposed by the diffusion model.For TUD, Table 9 Step 1.1 -Delayed neutron source along AA 0 and BB 0 .Table 10 Step 1.1 -Reactivity change from Step 0.2.Code q À q s0:2 (pcm) this is due to the different approach adopted for correcting cross sections with temperature (element-wise in TUD calculations, point-wise in the other codes).Anyway, the agreement of TUD with the other partners is good in a discrete sense.
Finally, Table 12 reports the value of reactivity insertion with respect to Step 1.1.The agreement between the multi-physics tools is excellent (maximum difference less than 40 pcm) and the same considerations made in Section 5.2.1 hold here.We can conclude that all partners can correctly reproduce the simple ''two-way" coupling between flux shape and temperature distribution, via the fuel-density feedback.

Step 1.3: Buoyancy
Table 13 reports the profiles of the observables along the cavity centerlines AA 0 and BB 0 .Note that the profile of u x along BB 0 is not reported as it is null.In fact, the problem is perfectly symmetric to the vertical centerline.Once more, a very good agreement can be noticed among the codes.Velocity components present an aver-Table 11 Step 1.2 -Profiles along AA 0 and BB 0 of the temperature and the change of fission rate density with respect to the solution obtained at Step 0.2.Table 12 Step 1.2 -Reactivity change from Step 1.1. Code À1122.0 age discrepancy that does not exceed 0.7%.Again, TUD simulations disagree most at the boundaries due to the DG-FEM discretization.Temperature profiles agree excellently, with an average discrepancy lower than 0.08%.This is confirmed by the plots shown in Fig. 10.For what concerns the delayed neutron source, its profiles are characterized by a discrepancy of around 0.5% along AA 0 and 1.2% along BB 0 .For the latter profile, in particular, we notice a good agreement between PSI and CNRS, and between TUD and PoliMi.
We therefore ascribe the small discrepancy in these results to the non-uniform mesh used by both CNRS and PSI, coarser in the middle of the domain.The plot in Fig. 10 corroborates this: a small discrepancy can be noticed only at the center of the cavity.However, considering the differences pointed out in the previous steps as well, the disagreement is considered acceptable.This is confirmed by the results on the reactivity difference with respect to Step 0.2 reported in Table 14.In fact, all values differ less than 50 pcm and are coherent with the ones analyzed in Sec-Table 13 Step 1.3 -Velocity components, temperature distribution, and delayed neutron source along AA 0 and BB 0 .Given the symmetry of the problem, ux is null along BB 0 and so not reported.tion 5.2.2.This proves that all multi-physics codes can correctly reproduce the effects induced by salt buoyancy on velocity, temperature, neutron flux, and delayed neutron precursors fields.

Step 1.4: Full coupling
The reactivity change with respect to Step 0.2 is reported in Table 15, for each reactor power and lid-velocity considered.As expected, higher power levels correspond to more significant insertions of negative reactivity, due to the increased average temperature in the system.The lid-velocity has a dual effect on reactivity, depending on the power level.At low power, higher lidvelocities correspond to higher negative reactivity differences, whereas the effect is reversed from P ¼ 0:6 GW.In fact, the more intense is forced-convection, the more precursors and energy get redistributed in the cavity, and this influences reactivity via two competing mechanisms.Shifting the pick of the precursors concentration further from the cavity center introduces more negative reactivity, with respect to the case U lid ¼ 0, whereas shifting the temperature pick introduces positive reactivity as the temperature gets lower in the region of higher importance.According to the power level, one effect prevails on the other, because at lower power the effect of the temperature redistribution is less significant.
Results are satisfactory.The differences between codes are always limited (the maximum difference varies from 10 to 45 pcm) and can be prescribed to the superposition of the effects already denoted in the previous steps.Tiberga et al. (2020).They were evaluated after the system response reached an asymptotic behavior at all frequencies.
Good agreement is found among the codes.The gain values are characterized by an average discrepancy of around 0.6%.Phaseshift values are more scattered, with a discrepancy of around 2.2%.However, this is considered acceptable.
From a physics point of view, the system dynamic response is as expected: at low frequencies, the power follows the variation of the extracted one (i.e., gain % 1 and minimal phase shift), as precursors can reach an equilibrium during the slow transient; at high frequencies, on the contrary, the contribution of the delayed neutrons to the entire neutron population is filtered, as some precursors families do not ''perceive" the temperature reactivity feedback due to their long half-life; so, the system response amplitude decreases (i.e., gain less than one).Moreover, the phase shift approaches À90°; in fact, when c reaches a maximum, the temperature time-derivative is at its maximum Table 14 Step 1.3 -Reactivity change from Step 0.2.Code q À q s0:2 (pcm) CNRS-SP 1 À1220.5 CNRS-SP 3 À1220.7 PoliMi À1227.0PSI À1219.6 TUD-S 2 À1208.5 TUD-S 6 À1184.4 Table 15 Step 1.4 -Reactivity change from Step 0.2 as a function of P and U lid .(at high frequencies, the time derivative term in the energy equation is dominating over the convective and diffusive ones), and so is the reactivity change and the derivative of the neutrons population, which shifts to a cosine wave.

Conclusions
In this paper, we have presented a coupled neutronics and fluiddynamics benchmark for multi-physics codes targeting fastspectrum molten salt reactors.The problem is based on a simple molten salt reactor system which is representative of the main characteristics of the Molten Salt Fast Reactor (MSFR): strong coupling between thermal-hydraulics and neutronics, fast spectrum, and transport of delayed neutron precursors.
The simplicity of the benchmark and its step-by-step approach, in which physics phenomena are coupled two at a time, make it an excellent tool to test multi-physics tools, because possible sources of discrepancy can be easily detected and fixed (as, in fact, it happened in our experience).This is in contrast with other complex benchmarks currently available in literature.Moreover, the generality of the problem makes this benchmark suitable for virtually any code.Multi-physics tools intended for thermal-spectrum MSRs may benefit from it too, even though the absence of a graphite moderator does not allow for a complete assessment.
We have presented the results collected during the benchmarking campaign involving four multi-physics tools developed in the context of the Horizon 2020 Euratom SAMOFAR project.A CFD-SP N transport tool developed at the Centre national de la recherche scientifique-Grenoble (CNRS), two CFD-diffusion solvers developed at the Politecnico di Milano (PoliMi) and the Paul Scherrer Institute (PSI), and a CFD-S N transport code developed at the Delft University of Technology (TUD).The comparison has shown very good agreement among the partners.Differences vary between a few tenths of a percent and a few percent maximum, depending on the particular benchmark step and quantity compared.They stem mainly from the use of different meshes between participants.Another reason for discrepancy is the neutronics model used.Diffusion solvers have been found to be often in excellent agreement, together with the CNRS-SP 1 and TUD-S 2 transport results, whereas TUD-S 6 transport results tend to have the largest differences.Moreover, the discontinuous Galerkin Finite Element discretization adopted by TUD often introduces slight discrepancies in the solution at the walls, where boundary conditions can be imposed only in a weak sense.The largest differences have been found in the fission density distribution of Step 1.2.They are mainly due to the different approach adopted by TUD in correcting cross sections with density feedback, that is, using the average temperature in each element.The differences are, however, considered acceptable.
In conclusion, the results prove that the recently developed multi-physics codes assessed in this work are able to reproduce accurately the physics phenomena characterizing fast-spectrum MSRs, both at steady-state and during transients.More in general, together with the complete benchmark description itself, they constitute a very useful tool for the future testing and development of any other multi-physics code dedicated to MSRs, as the interested reader is provided with a quite large spectrum of results to compare with.

Disclaimer
The content of this paper does not reflect the official opinion of the European Union.Responsibility for the information and/or views expressed therein lies entirely with the authors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B. Observables fields
This appendix is meant to complete Section 5 with plots of some observables in the entire domain relative to some steps of

Fig. 1 .
Fig. 1.CNRS benchmark: 2 m by 2 m domain.The cavity is insulated, surrounded by vacuum, and filled with molten salt at initial temperature of 900 K. Observables are compared pointwise along the centerlines AA 0 and BB 0 at several benchmark steps.

Fig. 3 .Fig. 4 .
Fig. 3. Coupling strategy of the PoliMi solver.The thermal-hydraulics problem is solved first, and the fuel temperature and density are transferred to the neutronics solver.Then, after solving the neutronics problem, the power density is passed back to the thermal-hydraulics solver.Iterations between the two solvers are performed within a time step until convergence.

5. 3 .
Phase 2: Time dependent coupling 5.3.1.Step 2.1: Forced convection transient Fig. 11 shows the power gain and phase-shift, summarized into Bode diagrams, as a function of the oscillation frequency of the heat transfer coefficient c.Numerical values can be downloaded from

Fig. 11 .
Fig. 11.Step 2.1 -Bode diagrams of power gain and phase-shift as a function of the frequency of the c wave.
Fig. B.12 shows the velocity magnitude field with flow streamlines obtained by each code in three different benchmark steps: Step 0.1 (V); Step 1.3 (B); and Step 1.4 with P ¼ 1 GW and U lid ¼ 0:5 m s À1 (C).

Fig
Fig. B.13. Temperature field with isolines obtained by each code in three different benchmark steps.From left to right: Step 0.3 (V); Step 1.3 (B); and Step 1.4 with P ¼ 1 GW and U lid ¼ 0:5 m s À1 (C).50 K was selected as isoline interval in all figures.Please note that the color bar scale for each step is different, and, for clarity, it has been reported only once per column.

Fig
Fig. B.12. Velocity magnitude field and flow streamlines obtained by each code in three different benchmark steps.From left to right: Step 0.1 (V); Step 1.3 (B); andStep 1.4 with P ¼ 1 GW and U lid ¼ 0:5 m s À1 (C).The streamlines were generated sampling 15 seeds along line BB 0 for case (V), 15 seeds along AA 0 for case (B), and 20 seeds along BB 0 for case (C).Please note that the color bar scale for each step is different, and, for clarity, it has been reported only per column.

Fig
Fig. B.14. Distribution of the first family of delayed neutron precursors (T1=2 ¼ 55:6 s) with isolines obtained by each code in three different benchmark steps.From left to right: Step 1.1 (V); Step 1.3 (B); and Step 1.4 with P ¼ 1 GW and U lid ¼ 0:5 m s À1 (C).5:0 Â 10 16 m À3 was selected as isoline interval in all figures.Please note that the color bar scale for each step is different, and, for clarity, it has been reported only once per column.
Fig. B.13 shows the temperature field obtained by each code in three different benchmark steps: Step 0.3 (V); Step 1.3 (B); and Step 1.4 with P ¼ 1 GW and U lid ¼ 0:5 m s À1 (C).

Table 3
Definition of energy groups.

Table 1
Fuel salt composition.

Table 2
Salt thermodynamic and transport properties.The aim is to compare the temperature distributions obtained by the different codes, fixing the velocity field and the heat source distribution.InputFixed flow field from Step 0.1; Fixed heat source distribution from Step 0.2; and
E R

Table A .
22 reports the fission spectra, the average number of neutrons (prompt and delayed) emitted, and the energy released per fission event.Finally, TableA.23 reports the decay constants and fractions characterizing the 8 families of delayed neutron precursors.

Table A .
18 P1 scattering cross sections.Table A.21 P1-corrected diffusion coefficient and inverse neutron velocities.

Table A .
19P2 scattering cross sections.

Table A .
23Fraction and decay constant for each family of delayed neutron precursors.