A multi-physics solver for liquid-fueled fast systems based on the discontinuous Galerkin FEM discretization

Abstract Performing accurate numerical simulations of molten salt reactors is challenging, especially in case of fast-spectrum designs, due to the unique physics phenomena characterizing these systems. The limitations of codes traditionally used in the nuclear community often require the development of novel high-fidelity multi-physics tools to advance the design of these innovative reactors. In this work, we present the most recent code developed at Delft University of Technology for multi-physics simulations of liquid-fueled fast reactors. The coupling is realized between an incompressible RANS model and an S N neutron transport solver. The models are implemented in two in-house codes, based on the discontinuous Galerkin Finite Element discretization, which guarantees high-quality of the solution. We report and discuss the results of preliminary simulations of the Molten Salt Fast Reactor at steady-state and during a Total Loss of Power transient. Results prove our code has capabilities for steady-state and transient analysis of non-moderated liquid-fueled reactors.


Introduction
Interest in molten salt reactors (MSR) has increased over the past two decades (LeBlanc, 2010), after the inclusion of this technology among the six Generation IV nuclear reactors in 2002 (Generation IV International Forum, 2002). Many MSR concepts are currently being investigated worldwide by companies and research institutes (Dolan, 2017), encouraged by the success of the Molten Salt Reactor Experiment (MSRE), which operated at Oak Ridge National Lab in the late sixties demonstrating the safety and the feasibility of the technology (Haubenreich and Engel, 1970;MacPherson, 1985).
Numerical simulations play an essential role in developing these innovative reactor designs, given the lack of substantial operational experience and know-how compared to the field of light-water reactors. However, simulating an MSR is a challenging task, as the use of a liquid salt both as fuel and coolant leads to unique physics phenomena which strongly (and non-linearly) couple neutronics and thermal-hydraulics : distributed, internal heat generation in the coolant; transport of delayed neutron precursors; feedback induced by fuel density variations. Moreover, in non-moderated designs, as the Molten Salt Advanced Reactor Transmuter (MOSART) (Ignatiev et al., 2014) or the Molten Salt Fast Reactor (MSFR) (Allibert et al., 2016), the traditional core shape with repeated structures (e.g., fuel pins) is lost.
Some efforts have been made to equip classical codes used in the nuclear community to model traditional solid-fueled reactors with the features necessary to address MSR simulations (e.g., Zanetti et al., 2015). However, legacy codes were usually developed decades ago, when strong computational capabilities were unavailable, and are often characterized by modeling assumptions (e.g., 1D flow model, point kinetics) that limit their range of usability. For these reasons, accurate simulations of molten salt reactor systems (especially fast-spectrum designs) need the development of new dedicated tools. We give hereafter a (non-exhaustive) overview of the previous work in the field before introducing our contribution.
The first decade of the century saw a succession of preliminary studies and modeling approaches based on simplifying assumptions (see Cammi et al., 2011, for a complete overview). First investigations on the coupling between neutron dynamics and fuel motion were performed by Lapenta et al. (2001) and Dulla et al. (2004), using simple neutronics models and considering only a 1D imposed velocity field with no temperature feedback. A prescribed salt velocity was considered also in the multi-group diffusion code Cinsf1D (Lecarpentier and Carpentier, 2003), used for the analysis of the AMSTER system (Vergnes and Lecarpentier, 2002). The temperature feedback on cross sections was taken into account by solving the energy equation in the fluid and graphite structures coupled through empirical heat transfer correlations. graphite moderator) model implemented in the same simulation environment (COMSOL Multiphysics®). It was used to get a deeper insight into the steady-state and transient characteristic of the Molten Salt Breeder Reactor (Robertson, 1971). At TU Delft, Van der Linden (2012) coupled the multi-group diffusion code DALTON with a full RANS/kÀ ε thermal hydraulic model. In the framework of the Euratom FP7 EVOL project (EVOL, 2010(EVOL, -2013, the two aforementioned multi-physics codes were used to conduct an extensive investigation of the behavior of the MSFR, both at steady-state and during accidental transient scenarios (Fiorina et al., 2014). A 2D axisymmetric geometry was adopted, given the cylindrical shape of the MSFR core at that time. Aufiero et al. (2014b) developed the first OpenFOAM® multi-physics tool for MSRs, coupling a diffusion model with the RANS/kÀ ε equations. Being capable of dealing with full-core 3D analysis, the tool was employed to investigate for the first time an asymmetric loss of flow accident inside the MSFR. A single-group diffusion model was adopted to relieve the computational burden though. A similar tool is the OpenMC/TANSY code system developed by Hu et al. (2017) for the analysis of the MOSART. More innovative from the neutronics point of view is the MSFR multi-physics tool of Laureau et al. (2017). The Transient Fission Matrix method was in fact chosen to reproduce transient Monte Carlo-like calculations with a reduced computational cost. Fiorina et al. (2015) presented GeN-Foam, a multi-physics code coupling a RANS/k À ε model, extended to coarse-mesh applications thanks to the adoption of a porous-medium approach, and a multi-group diffusion neutronics solver. The code is also equipped with a displacement-based thermal-mechanics solver and a finite-difference model for the temperature field in the fuel, since it targets homogeneous liquid-fueled reactors as well as light-water or liquid-metal reactors. This tool was recently used by Altahhan et al. (2019) to perform preliminary design investigations of the primary loop of a chloride-based molten salt fast reactor. All previously mentioned multi-physics models consider the fuel salt to be an incompressible fluid. However, Aufiero et al. (2017) pointed out that taking into account the fuel compressibility is essential to correctly model the fuel-expansion feedback during fast, super-prompt critical transients. This was confirmed by their simulations, performed on a simplified MSFR geometry, using a novel multi-physics tool coupling a compressible flow model, available in OpenFOAM®, with the Monte Carlo code Serpent (Lepp€ anen et al., 2015). Further investigations have been recently performed by Cervi et al. (2019a), using a multi-physics tool consisting of a multi-group diffusion model coupled to a two-phase compressible RANS solver (both implemented in Open-FOAM®), able therefore to model the effects of non-uniform distributions of gas bubbles (helium) in the fluid mixture. The neutronics model was later improved adopting a Simplified P 3 approximation of the neutron transport equation (Cervi et al., 2019b).
In this paper, our goal is to present the most recent advances in the field of multi-physics simulations of (non-moderated) MSRs made at TU Delft. In the framework of the H2020 SAMOFAR project (http://samofa r.eu/), we developed a novel multi-physics tool realized by coupling an incompressible RANS model with an S N neutron transport solver, implemented in two in-house codes. In both of them, the discontinuous Galerkin Finite Element method (DG-FEM) is adopted to discretize in space the model equations. The method was chosen as it combines the advantages of local conservation, as in finite volumes, with the highorder discretization and high geometric flexibility (thus easily handling complex shapes as the MSFR core) of finite elements, guaranteeing high accuracy of simulations.
The MSFR can be considered a highly homogeneous system when compared to solid-fueled reactors, so diffusion models provide very good results in standard scenarios (as those considered in this work). Nevertheless, the simulation of reactor conditions characterized by significant material discontinuities will have to be carried out at a later development stage. For example, it will be necessary to consider draining scenarios of the reactor or the blanket, or malfunctions of the helium bubbling system (Delpech et al., 2009) leading to large and non-uniform concentration of bubbles in the fuel salt. In this situations, a full-transport model is necessary to describe more accurately the neutronics behavior of the reactor (Duderstadt and Hamilton, 1976). In this perspective, the combination of S N with the DG-FEM discretization leads to a state-of-the-art deterministic method to solve the Boltzmann equation (Reed and Hill, 1973;Wareing et al., 2001) which guarantees high-quality of the solution even in presence of large discontinuities in material properties. Our tool constitutes therefore a useful contribution to the set of available multi-physics codes dedicated to the analysis of molten salt fast reactors.
The paper is outlined as follows. In Section 2, we describe in detail the thermal-hydraulics and the neutronics models and the coupling between them. Then, as an example of the capabilities of the coupled code, which has been benchmarked against other similar numerical tools in a recent work of some of the authors , we perform preliminary simulations of the MSFR behavior at steady-state and during the first 30s of an unprotected Total Loss of Power (TLOP) accident. For this reason, we first describe the MSFR design and the modeling choices adopted in Section 3. Then, we report and analyze the results of our simulations in Section 4, drawing some conclusions in Section 5.

Description of the multi-physics tool
In this section, we describe in details the models constituting our multi-physics code. The coupling is realized between two in-house Fortran codes: DGFlows for CFD, and PHANTOM-S N for neutron transport.

Fluid dynamics: DGFlows code
DGFlows is designed to solve the set of Navier-Stokes equations in the low-Mach number limit (Hennink et al., 2020). Turbulent flows are handled through Reynolds Averaged Navier-Stokes (RANS) models. The RANS model equations read (omitting dependencies for clarity) (Patankar, 1980) ∂ρh ∂t Here, ρ is the fluid density, m the mass flux and u ¼ m= ρ the velocity, p a pseudo-pressure, and h the specific enthalpy. The shear stress tensor, τ, is defined as where I is the identity tensor, while μ and μ t are the molecular and turbulent (or eddy) dynamic viscosities. κ in Equation (1a) is the molecular thermal conductivity, while κ t is the turbulent one. The latter is computed as κ t ¼ μ t c p =Pr t , where c p is the salt specific heat capacity and Pr t is the turbulent Prandtl number. Finally, S h is the energy source term (e.g., fission heat), F buoy is the buoyancy force, and F m includes all other possible momentum sources.
In this work, as we did not address any super-prompt critical transient, we safely considered the salt density to be constant, thus reducing System (1) to the set of incompressible RANS, and we adopted the Boussinesq approximation to model buoyancy: where T ref is the Boussinesq reference temperature, g is the gravity acceleration, β buoy is the fluid thermal expansion coefficient, and ρ ref is the salt density computed at the reference temperature. System (1) has to be closed with a turbulence model and proper initial and boundary conditions. In this work we adopted the standard k À ε model (Launder and Spalding, 1974), as customary in the context of MSRs multi-physics simulations (e.g., Cammi et al., 2011;Fiorina et al., 2014;Aufiero et al., 2014b;Laureau et al., 2017). The model was modified to take into account the turbulence production by buoyancy (Van der Linden, 2012;Van Maele and Merci, 2006) and cast into a logarithmic form to ensure the positivity of turbulence quantities (Ilinca et al., 1998): Here, K ¼ lnðkÞ and E ¼ lnðεÞ, μ t ¼ ρC μ e 2K À E . Turbulence is produced by both shear stress and buoyancy: Pr t β buoy g⋅rT: The latter term can actually be negative in case of thermal stratification. All remaining closure coefficients are reported in Table 1.
All boundary and initial conditions chosen for the simulations reported in this study are described in Section 3.2.6.

Spatial and temporal discretization and solution of the linear systems
The discontinuous Galerkin FEM is chosen to discretize in space the model equations. We briefly describe here only the key ingredients of the discretization, and refer to Hennink et al. (2020) (and the references cited therein) for more details.
Instead of primitive quantities, as usual, we solve for the set of conservative variables, to fully exploit the local conservation of DG. All unknowns are represented in the Galerkin space using a hierarchical set of modal basis functions up to order P (for vector quantities) and P À 1 (for all other scalars). Diffusive terms are discretized with the Symmetric Interior Penalty (SIP) method, which is consistent and stable and leads to optimal spatial convergence rates and compact stencil size (Shahbazi et al., 2007). The Lax-Friedrichs numerical flux is chosen for the convective terms (Cliffe et al., 2010).
The general transport equation for quantity ζ is discretized implicitly in time with the backward differentiation formula (BDF) of order 2 (with constant time step Δt): where γ 0 ¼ 3=2, γ 1 ¼ À 2, and γ 2 ¼ 1=2, and the superscript n þ 1 indicates the new time step. Matrix A collects all contributions deriving from the discretization of diffusive and convective terms.
The discretized RANS/k À ε system is solved starting from the energy equation. Then, the coupled momentum-continuity equations are solved in a segregated way using a second-order time accurate pressure correction scheme (Van Kan, 1986); the splitting is done at algebraic level, as explained by Shahbazi et al. (2007), to avoid the imposition of pressure boundary conditions. Finally, the equations for the turbulence quantities are solved in sequence. All non-linearities (e.g., in the convective terms, or in the k À ε equations) are solved by initializing the interested quantity, ζ * , with a second-order extrapolation from the previous time steps, that is, ζ * ¼ 2ζ n À ζ nÀ 1 , or by simply choosing ζ * ¼ ζ nþ1 , when the latter is known. Moreover, if necessary, inner iterations between the momentum and the pressure-Poisson equation and outer iterations between the momentum-continuity block and the k À ε block can be performed.
DGFlows has capabilities for parallel computation. We use METIS to partition the mesh (Karypis and Kumar, 1998) and the MPI-based software library PETSc (Balay et al., 2018) to assemble and solve all linear systems. The pressure-Poisson equation is solved with the conjugate gradient method and an additive Schwarz preconditioner, with one block per process and an incomplete Cholesky factorization on each block; all other linear systems are solved with the GMRES method and a block Jacobi preconditioner, with one block per process and an incomplete LU factorization on each block. All linear solvers are initialized with a second-order time extrapolated solution from the previous time steps, in order to accelerate convergence.

Neutronics: PHANTOM-S N code
PHANTOM-S N is a solver for the multi-group Boltzmann neutron transport equation. This in-house code, already capable of calculating criticality and time eigenvalues (K� oph� azi and Lathouwers, 2012), of performing both regular and generalized perturbation analysis (Perk� o, 2015, Ch. 3) and goal-oriented mesh refinement (Lathouwers, 2011a, b), was extended to solve time-dependent problems and, in particular, to handle the transport of delayed neutron precursors.
The time-dependent multi-group model equations read (omitting Table 1 Values of the closure parameters of the k À ε turbulence model (Launder and Spalding, 1974;Van Maele and Merci, 2006 (Lewis and Miller, 1993) 1 Here, ϕ and Φ are the angular and scalar fluxes, Σ t and Σ f are the total and fission macroscopic cross sections, while Σ s ðΩ' ⋅ΩÞ is the scattering differential cross section. ν indicates the average number of neutrons emitted per fission event, χ is the fission spectrum (superscript p and d indicate the prompt and delayed spectra), Ω is the neutron travel direction, while v is the neutron speed. Subscripts g and g' span all energy groups. C i , β i , and λ i are the concentration, fraction, and decay constant of the i th -family of delayed neutron precursors, and β ¼ P i β i . In Equation (7b), ν and ν t are the molecular and turbulent kinematic viscosities, while Sc and Sc t are the molecular and turbulent Schmidt numbers. Standard vacuum (i.e., zero incoming flux) or reflective boundary conditions typically close Equation (7a).
In case of steady-state calculations, System (7) is cast into a criticality eigenvalue problem, by removing the time derivative terms and dividing ν by the multiplication factor k eff . The resulting flux is then normalized to the desired reactor power level.
When simulating long-term transient scenarios in nuclear reactors, it is paramount to take decay heat into account. In PHANTOM-S N , we implemented transport equations to model the evolution of the fission products responsible for the decay heat. Following Aufiero et al. (2014b), the approach is very similar to what is done for the delayed neutron precursors. The decay heat fission products are divided into "families", each one characterized by a decay constant λ dh j , and a fraction β dh j . The transport equation for family j reads where E f is the average energy released by each fission event. This equation models the balance of d j , the concentration of decay heat precursors in family j multiplied by the average energy released by that decay family. As the delayed neutron precursors, also the decay heat fission products are drifted by the flow in liquid-fueled systems, hence the presence of the convective and diffusive terms in (8).
When decay heat is modeled, the power density term to be passed to Equation (1a) becomes

Spatial and temporal discretizations and solution of the linear systems
After the expansion of the scattering term in a series of spherical harmonics, Equation (7a) is discretized in angle by the method of discrete ordinates, selecting directions and weights based on the levelsymmetric quadrature set (Lewis and Miller, 1993). The DG-FEM is adopted for spatial discretization, upwinding the flux values in all face integrals of the weak formulation. This is a well established approach in the transport community (e.g., Reed and Hill, 1973;Wareing et al., 2001). Moreover, all the details of it have been reported by K� oph� azi and Lathouwers (2012), so they are not repeated here. Equations (7b) and (8) are discretized in space with the DG-FEM, adopting the same choices as in DGFlows (see Section 2.1.1). Finally, in transient calculations, all time derivatives are discretized using the implicit BDF2 scheme (Equation (6)).
The coupled discrete System (7) can be cast into the following linear form: where A t ϕϕ , A f ϕϕ , and A s ϕϕ represent respectively the discretization of the transport, fission, and scattering operators in Equation (7a); A ϕC represents the delayed neutron precursors decay in Equation (7a), A CC ðu * Þ the precursors time-derivative, convection, diffusion, and decay, and A f Cϕ is the fission source term in Equation (7b). The right hand side vector collects the known terms deriving from the discretization of the time derivatives.
In steady-state calculations, System (10) is cast into the following eigenvalue problem: which is solved with a standard power method approach (Golub and Van Loan, 2013). Linear Systems (10) and (11) are solved by a flexible Krylov method (e.g., GCR), with a physics-based preconditioning approach, which depends on the type of problem. Let r be the residual and z the preconditioned residual vectors at iteration l of the outer Krylov method. For steady-state calculations, System (11) is preconditioned by a standard block Jacobi method: This is a good enough preconditioner, taking into account that the contribution of the delayed precursors to the neutron transport equation is small, at steady-state. Moreover, System (11) is typically solved with large tolerances, in order to quickly move to the next iteration of the power method.
In transient calculations, on the other hand, System (10) must be solved with sufficient accuracy, but a standard block Gauss-Seidel method does not perform well, as the contribution of the precursors to Equation (7a) is essential. For this reason, we adopted an "improved" version of a block Gauss-Seidel preconditioner. As explained by Pautz and Birkhofer (2003), if the fuel were solid, one could easily diagonalize the time-discretized System (7) by inverting the left hand side of the precursors equation (thanks to the absence of convective and diffusive terms) and plugging C nþ1 i into Equation (7a). This would lead to an equation for ϕ g only with a modified fission spectrum, where , and an additional term to the right hand side: At preconditioner step, we modify the neutron transport equation in a similar fashion, as if the fuel were solid. In fact, we take a modified equation for the precursors' preconditioned residual: in which Ã CC is an approximation of A CC ðu * Þ, without convection and diffusion terms. Hence, we easily invert it and plug z lþ1 C into the equation for z lþ1 ϕ . Therefore, the full preconditioner is This is an improved Gauss-Seidel method, as it takes into account the contribution of the precursors in the flux equation (even if just approximately), thus leading to faster convergence of the outer Krylov method. This preconditioner is the more effective the smaller Δt, as lim Δt→0 A CC ðu * Þ ¼Ã CC (and, of course, the lower is the salt flow rate). To conclude, it is worth to underline that the flexibility of the outer Krylov method is fundamental for convergence, given the variability of the described preconditioners at each iteration.
Independently from the chosen preconditioner, the multi-group flux subsystem is solved with Gauss-Seidel iterations. Each group fixedsource problem is in turn solved as proposed by Warsa et al. (2004), that is, by applying a Krylov method preconditioned with a directional sweep procedure that inverts the group transport operator. Then, as the molten salt is a highly scattering medium, we apply a Diffusion Synthetic Acceleration preconditioner as well. Our DSA scheme is partially consistent and based on an SIP-DG discretization of the diffusion-like equation (Wang and Ragusa, 2010). A more detailed description of the solution technique for the multi-group problem can be found in K� oph� azi and Lathouwers (2012). Precursors subsystems (as well as the decay heat precursors systems) are solved with GMRES and ILU preconditioner implemented in PETSc. Finally, as in DGFlows, all solvers are initialized with a second-order time extrapolated solution in transient calculations, to speed-up convergence.

Coupling strategy and cross sections treatment
Fig. 1 displays the structure of the multi-physics tool and the data exchanged between the codes. The fission power density is transferred to DGFlows as it appears on the right hand side of the energy equation. Velocity and turbulent viscosity fields influence the precursors distribution, so they follow the inverse route. Moreover, DGFlows computes an average temperature on each element (T E ) and exports it to the routines in PHANTOM-S N devoted to the computation of cross sections. Based on these temperatures, cross sections are interpolated starting from a set of libraries at prescribed temperatures, generated by Monte Carlo or deterministic codes. In case only a single library is prescribed (as for the calculations reported in this work), cross sections are corrected by assuming a linear dependence on the fuel density and a logarithmic dependence on the temperature to model Doppler feedback (Waltar et al., 2012;Aufiero et al., 2014b): where T 0 is the library reference temperature, ρ ref is the reference density at which macroscopic cross sections are evaluated, 1 and α r is a coefficient.
In transient simulations, a loose-coupling strategy is adopted, where DGFlows is run first and, after completion of a time-step, PHANTOM-S N is called to handle the neutronics physics. In order to have a nonlinearconsistent coupling scheme (Ragusa and Mahadevan, 2009), thus preserving a global second order time accuracy, the power density passed to DGFlows at the beginning of a new time step is extrapolated using a second order scheme from the values at the previous time steps. In steady-state simulations, the codes are iterated until convergence.

Mesh generation and manipulation
Both codes can handle structured and unstructured meshes with tetrahedral or hexahedral elements (triangles or quadrangles in 2D), generated with the open source tool Gmsh (Geuzaine and Remacle, 2009). Starting from the same "master" mesh, each code independently can perform a local, hierarchical refinement. However, in practice the neutronics mesh corresponds to the master mesh, while DGFlows refines it (even more than once) in regions where the flow gradients have to be resolved but the neutron importance is low (e.g., close to walls, or in portions of the fuel circuit far from the core cavity). Fig. 2 illustrates the hierarchic refinement on a simple triangular mesh. 2 The hierarchy of the refinement makes the exchange of data between the two meshes easy, exploiting Galerkin projection routines. Let E R be an element on the refined mesh and E C the corresponding mother on the coarse mesh; let also ζ C and ζ R be the vectors of coefficients of the DG-FEM expansion for quantity ζ on E C and E R , respectively. The coarse-torefined mapping is performed as follows: For each element of the refined mesh, 1. if E R is equal to E C (i.e., there is no actual refinement), then ζ R ¼ ζ C ; otherwise: 2. compute the quadrature points (ξ R ) and weights (w R ) on E R and map them onto E C ; 3. compute the FEM mass matrix on E R (M R ), and solve the following linear system: where the right hand side vector is evaluated by Here, θ R and θ C indicate the basis functions on E R and E C , while ξ R→C the mapped quadrature point; index q spans all the quadrature points, and indices i and k span the degrees of freedom on E R and E C respectively.
This mapping is exact. The refined-to-coarse mapping is performed as follows: For each element of the coarse mesh, do 1. if E C is equal to E R , then ζ C ¼ ζ R ; otherwise: 2. compute the FEM mass matrix on E C (M C ) 3. initialize ζ C ¼ 0. Then, for each child element on the refined mesh, (a) compute the quadrature points and weights on E R and map them onto E C ; Fig. 1. Computational scheme of the multi-physics tool. DGFlows, the CFD solver, and PHANTOM-S N , the neutron transport code, exchange data at each iteration to resolve the coupling between the various physical phenomena characterizing a molten salt fast reactor.
(b) solve the linear system with the right hand side vector evaluated by Here, indices i and k span the degrees of freedom on E C and E R respectively, while q spans all quadrature points on E R ; This mapping inevitably leads to loss of information. This has no significant impact when transferring the eddy viscosity field to the neutronics mesh. However, the loss of accuracy on the velocity field is detrimental for the precursors solution, because the velocity field no longer satisfies the divergence-free constraint in a discrete sense on the neutronics mesh. For this reason, in PHANTOM-S N the velocity field is not mapped to the neutronics mesh when discretizing the convection integrals in the weak form of Equations (7b) and (8). Instead, each integral on an element of the neutronics mesh is computed as the sum of the integrals on the corresponding children elements of the CFD mesh. The velocity is computed at quadrature points on the latter, so that no loss of accuracy is introduced. The same technique is also used to compute the average temperature on each element of the neutronics mesh, without loss of information.

Verification and validation
Activities to validate and verify the single-physics codes were carried out in previous works (K� oph� azi and Lathouwers, 2012;Hennink et al., 2020).
Recently, the coupled-code took part in a benchmarking campaign for multi-physics simulation tools dedicated to liquid-fuel fast reactors, whose results were published in Tiberga et al. (2020). The system under investigation consisted of a molten salt reactor whose characteristics make it a simple representation of the MSFR. Several steady-state or transient problems were simulated, in which physics phenomena were progressively coupled to easily identify potential sources of error. We compared our results with those obtained by other partners of the SAMOFAR project. The very good agreement proved our code can correctly reproduce and model the unique physics phenomena characterizing fast-spectrum MSRs both at steady-state and during transients.

The MSFR: design and modeling approach
In this section, we describe the MSFR design under investigation and the modeling choices for the simulations reported in Section 4.

Design description
The current MSFR concept (Allibert et al., 2016;G� erardin et al., 2017) is a non-moderated reactor operating in the thorium fuel cycle. The molten salt mixture acts both as fuel and thermal carrier, and it is characterized by a strong negative temperature feedback coefficient (Brovchenko et al., 2013;Heuer et al., 2014) that increases the reactor safety and makes the presence of control rods in the core cavity unnecessary. Fig. 3 schematically displays a cross section of the fuel circuit design. The molten salt enters the core cavity at the bottom, rises, and then flows out towards sixteen identical loops called sectors. Each sector is equipped with a pump and a heat exchanger, where the fuel salt releases heat to an intermediate circuit. Both pump and heat exchanger designs are not fully specified at the moment, though a preliminary analysis of a Printed Circuit Heat exchanger (PCHE) was carried out during the SAMOFAR project (Di Ronco et al., 2019). Dispersion of helium bubbles in the salt mixture is foreseen to remove gaseous and metallic fission products and as potential reactivity control mechanism (Delpech et al., 2009).
The toroidal core shape is the result of optimization studies carried out during the EVOL project and aimed at homogenizing the fuel temperature in the core cavity thus eliminating large hot spot regions Fig. 2. Illustration of the hierarchic refinement of a triangular master mesh. Element E 1 is refined twice. Galerkin projection procedures are used to map quantities between the two meshes. Fig. 3. Schematic view of the MSFR fuel circuit. (Rouch et al., 2014). The central torus is surrounded by a blanket, containing a fertile mixture of lithium and thorium fluorides, which improves the reactor's breeding capabilities. Nickel-based alloy reflectors are present at the top and bottom of the core to improve neutron economy. The gap between the blanket and the external sectors is filled with a boron-carbide shield . Table 2 reports the MSFR design parameters relevant for this study. Fig. 4 shows the MSFR geometry used in this work. Given the symmetry of the problem, only a single recirculation loop was simulated, in order to reduce the computational burden. Moreover, heat transfer in reflectors and blanket was not modeled. However, these regions are included in our neutron transport calculations.

Geometry and mesh
Two sets of meshes were generated, to study the variation of the main results upon refinement: Set "M1" The same master mesh is used for both neutronics and CFD calculations. It has 68001 tetrahedra, of which 60468 within the boundaries of the fuel salt domain. The two meshes are shown in Fig. 5; Set "M2" A master mesh having 28574 tetrahedra (23032 in the fuel salt domain) is used for neutronics. This mesh is finer than M1 in the core, but much coarser in the external sector, where neutron importance is low, thus optimizing neutronics calculations from a computational point of view. For this reason, the master mesh is refined for CFD calculations: once, uniformly in the outer-core region (i.e., legs, pump, and heat exchanger), and then once more close to the wall boundaries, to better resolve the boundary layers, resulting into 189464 elements. The two meshes are shown in Fig. 6;

Pump modeling
The MSFR pump design has not been defined yet. For this reason, it was modeled as a simple momentum source: where H pump ¼ 0:1m is the height of the portion of the outlet leg identified with the pump (see Fig. 4), Δp pump is a parameter suitably tuned to obtain the desired volumetric flow rate of m 3 s À 1 , and b e z is the unit vector indicating the direction of the z-axis.

Heat exchanger modeling
The heat exchanger was modeled through a porous medium approach. In fact, a detailed CFD simulation of the component would be unfeasible, from a computational point of view; moreover, its design has not been finalized yet. However, based on preliminary studies during SAMOFAR (Di Ronco et al., 2019), we considered only a portion of the sector "box" to be active (Fig. 4), 0:69m long. Resistance to flow was reproduced by the following force: where K loss is a loss coefficient suitably tuned to obtain the desired pressure drop of 4:5bar. Salt cooling was modeled via a volumetric heat sink term: where γ is a volumetric heat transfer coefficient suitably tuned in order to achieve the target minimum salt temperature in the system of 650 ∘ C. This coefficient does not depend on the salt velocity, because the small size of the channels of the PCHE limits the salt flow to laminar conditions, thus leading to a constant Nusselt number (Di Ronco et al., 2019).

Salt physical properties
Two fuel salt compositions are currently being investigated for the MSFR (Allibert et al., 2016;G� erardin et al., 2017). In this work, as reported in Table 2, we opted for the TRU-based composition, due to the availability of values and equations of state for the physical properties. These are reported in Table 3. Moreover, we considered a turbulent Prandtl number of 0.85, and, in absence of detailed data on the diffusion of species in the salt mixture, we assumed a turbulent Schmidt number of 0.85 (Aufiero et al., 2014b).

Table 2
Main MSFR design parameters considered in the present study (Allibert et al., 2016;G� erardin et al., 2017 Fig. 4. MSFR geometry used for simulations: only 1/16th of the core was modeled, with the associated sector.   6. M2 mesh set. The neutronics mesh (left) has 28574 tetrahedra, while the CFD mesh (right) has 189464 elements. The latter is derived from the former by refining it once uniformly in the outer-core region, then again close to all wall boundaries.

Neutronics data
A set of six-groups condensed neutronics data was used, evaluated at temperature T 0 ¼ 900K with Serpent (Lepp€ anen et al., 2015), selecting the JEFF-3.1.1 nuclear data library (Santamarina et al., 2009). The energy group structure is reported in Table 4. Despite the few groups, this structure proved to be sufficient to reproduce the overall MSFR spectrum and the relevant neutronics parameters when compared to Monte Carlo calculations (Fiorina et al., 2012). The coefficients α r in Equation (17) were evaluated by logarithmic interpolation of Serpent cross sections between 900K and 1200K.
Eight families of delayed neutron precursors were modeled. Their fractions and decay constants, taken from the JEFF-3.1.1 library, are reported in Table 5. Decay heat was modeled through Equation (8), taking into account three families of precursors. As shown by Aufiero et al. (2014b), the superposition of three exponentially decaying terms is an acceptable approximation for transients limited to a few minutes, as the one described in Section 4.2. The considered decay heat precursors fractions and time constants are reported in Table 6.

Boundary conditions
Boundary conditions for CFD calculations included symmetry conditions at the wedge sides and standard logarithmic wall-functions plus adiabatic conditions at all walls. Neglecting the heat transfer with reflectors, blanket, and external environment is an acceptable approximation and is conservative (Rouch et al., 2014).
For neutronics, reflective boundary conditions were imposed at the wedge sides and vacuum conditions everywhere else. Surfaces delimiting the gap between blanket and heat exchanger can be reasonably assumed to be facing void, given the presence of the boron-carbide layer (a strong neutron absorber). The convective and diffusive terms in Equations (7b) and (8) were included only in the mesh portion corresponding to the fuel salt domain. Hence, homogeneous Neumann and no-inflow conditions were imposed at all walls.

Results of preliminary MSFR simulations
In this section, as an example of the capabilities of our code-system, we report and discuss the results of preliminary simulations aimed at investigating the steady-state and transient behavior of the MSFR.
All simulations were run adopting a polynomial approximation of order P ¼ 2 for the mass flux and P ¼ 1 for all other quantities. Order N ¼ 4 was chosen for the discrete-ordinates discretization, and the scattering anisotropy was taken into account up to first order. To investigate the effect of mesh refinement, we report the results obtained on both meshes described in Section 3.2.1. Transient calculations were run choosing a coupling time step size Δt ¼ 0:025s on the M1 mesh and of Δt ¼ 0:01s on the finer one; to keep the CFL number in CFD calculations sufficiently small and avoid numerical instabilities, sub-steps were performed in DGFlows: 10 on the M1 mesh, and 20 on the M2 mesh.

Steady-state solution
As explained in Section 2.3, the steady-state solution was sought by iterating PHANTOM-S N and DGFlows until convergence, and adjusting the pump head, the friction factor and the heat transfer coefficient in the heat exchanger to obtain the desired design specifications reported in Table 2. The final values of these parameters are reported in Table 7. Table 8 reports some design parameters of interest derived from the steady-state solution (comparing them on meshes M1 and M2). Table 9 reports the history per coupling iteration of the same parameters evaluated on M2. As the initial step is a k eff calculation with static and isothermal fuel conditions, the non-linearities related to thermal feedback and precursors motion are resolved mostly in the second iteration. Four iterations in total were sufficient to reach convergence. Fig. 7 shows the power density and temperature fields obtained on meshes M1 and M2. The power density (Fig. 7a) has the expected shape: almost a cosine in the axial direction and a Bessel J 0 function in the Table 5 Fraction and decay constant for each family of delayed neutron precursors.

Parameter (unit) Value on M1 Value on M2
Δppump (Pa) 5.513 � 10 5 5.518 � 10 5 K loss (À ) 2.811 � 10 2 2.811 � 10 2 γ (Wm À 1 K À 1 ) 1.995 � 10 7 1.995 � 10 7 Effective multiplication factor k eff -1.00998 1.00994 Table 9 History per coupling iteration of the design parameters of interest evaluated on mesh M2. Four iterations were sufficient to reach convergence.  2.479 � 10 À 2 5 5.531 � 10 À 3 6 7.485 � 10 À 4 radial one, deformed by the non-uniform temperature field and the presence of reflectors and blanket. This is confirmed by the axial and radial power profiles shown in Fig. 8. As expected, the power quickly drops in the legs and heat exchanger but it does not reach zero, due to the distribution of the fission products which decay and release heat throughout the domain. Contrary to solid-fueled reactors, these large power gradients are not a concern in the MSFR, thanks to a homogeneous fuel irradiation ensured by recirculation and the absence of safety criteria limiting the power shape to preserve the integrity of the solid fuel pins. Actually, these gradients are beneficial, as they lead to lower radiation-induced damage on the core vessel. From Fig. 7a, no relevant differences can be observed between the two meshes. The k eff reported in Table 8 indicates the reactor is quite far from criticality, so the concentration of fissile material should be adjusted. However, a criticality search was not in the scope of the present study. The difference in reactivity between the two meshes is less than 4pcm. The salt temperature increases in the core cavity, but not uniformly. From Fig. 7b, one can see that the salt is the hottest at the center of the core cavity, at the boundary with the upper reflector, where the fluid is almost stagnant. The upper reflector is in general subject to high temperatures and strong thermal gradients. On the contrary, there is no hot spot at the interface with the blanket, at core inlet, confirming the effectiveness of previous core-shape optimization studies (Rouch et al., 2014). The maximum salt temperature is around 70K higher than the average temperature of the salt in the outlet leg. In the heat exchanger, as reported in Table 8, the temperature drops less than the reference 100K (Allibert et al., 2016), due to a peak in the turbulent diffusivity, as explained in the following. The average salt temperature is around 700 ∘ C, as expected.
From Fig. 7b and Table 8, only minor discrepancies can be noticed in the results on the two meshes. Considering the solution on M2 as reference, the error on ΔT hx and T avg calculated on M1 is below 0:1%. The better boundary layer resolution on M2 leads to a T max 3K higher than on M1, but the difference is only 0:3% in relative terms.
Having proved the mesh-convergence of the results, we continue the analysis considering only the simulations on M2, which led to a better flow resolution. Fig. 9 shows the velocity field in the domain. The salt Reynolds number is around 4.6 � 10 5 in the middle of the core cavity and 1:75 � 10 5 in the pump region. However, the salt velocity increases considerably (almost threefold) towards the outer leg, due to the significant reduction of cross section. The contribution of natural circulation to the total mass flow rate is negligible at nominal conditions, due to the large pressure drop in the heat exchanger. At the core inlet, boundary layer detachment can be noticed, but this does not lead to local hot spots. Another detachment is present at the entrance of the heat exchanger "box", due to the sudden expansion. However, the flow in this region is of less interest, given the simple and unrealistic geometry used to represent this portion of fuel circuit.
The large recirculation at the entrance of the core leads to a peak in the turbulent kinetic energy and diffusivity, as shown in Fig. 10. The temperature distribution suppresses the salt turbulence while the salt rises through the core, but then turbulence increases again when the salt accelerates towards the outlet leg. The highest values of k and ν t are found in the heat exchanger and inner leg regions, but, once more, these are of less interest. A more realistic design of the heat exchanger "box" would eliminate the large vortex at the entrance and thus the spike in Fig. 7. Power density (left) and temperature fields (right) obtained at steady-state with the two meshes (mid-plane cuts). Only small differences in the fields can be noticed between the two meshes, confirming the good convergence of the results. turbulence production. This would also lead to a drop in the salt temperature across the heat exchanger closer to the reference 100K.
To conclude the steady-state analysis, Fig. 11 shows the distribution of the first delayed neutron precursors family on the left side and the one of the last family on the right. Both distributions are largely distorted by the flow field. However, the short-lived precursors manage to decay mostly within the core cavity, contrary to the long-lived ones, which can be found in high concentration in the heat exchanger region. This reduces the margin to prompt criticality (Aufiero et al., 2014a).

Transient: unprotected Total Loss of Power accident
As an example of transient scenario, to show the potential of our code package, we chose one of the most (potentially) severe accidents that can happen in the MSFR: the unprotected Total Loss of Power (TLOP). This accident could occur in case of total blackout; the absence of electricity would stop the pumps in the fuel, intermediate, and energy conversion circuits; moreover, the removal of the fission power from the core would become impossible.
In this work, the pump stop was modeled through an exponential decay of F pump in Equation (22), to take into account the pump inertia. In absence of precise design specifications, a time constant of 5s was chosen. The intermediate and energy conversion circuits would continue, in reality, to remove some power from the fuel circuit, thanks to their thermal inertia; passive decay heat removal systems are foreseen as well and their design is under study; moreover, a non-negligible fraction of power would be exchanged with the external environment by passive heat transfer mechanisms, as conduction or radiative transfer. However, throughout this accident scenario, we conservatively assumed that no heat is removed from the heat exchanger (achieved by imposing γ ¼ 0 in Equation (24)) and that the adiabatic boundary conditions described in Section 3.2.6 continue to hold. Considering one recirculation loop only is reasonable, as this accident affects all 16 loops in the same way. The initial condition is represented by the steady-state solution described in the previous section. Since the reactor is not critical at steady-state with the current salt composition, the fission operator was scaled by a factor 1=k eff . Fig. 12 shows the time trends of various characteristic quantities during the first 30s of the accident. Up to 1s, the reactor power (Fig. 12a) does not vary significantly, thanks to the precursors hold-back, but then it quickly drops to less than one-third after 4s; at 5s, the descent slows down, and the power reaches almost a plateau around 7s; from roughly Fig. 10. Eddy viscosity (left) and turbulent kinetic energy (right) fields obtained at steady-state on the M2 mesh (mid-plane cuts). Turbulence is generated at the entrance of the core, due to the presence of a vortex, and at the core outlet. Salt thermal stratification suppresses turbulence in the core cavity. 8s onward, finally, the power decreases monotonically. At t ¼ 30s, the total power is reduced to 208:5MW, 40% of which is constituted by decay heat (see Fig. 12b).
The plateau in the power trend around 7s is due to a peculiar phenomenon of liquid-fueled reactors. In the first seconds of the transient, before the power drops, a certain amount of delayed neutron precursors exits the core. The long lived portion does not decay in the sector but only after it re-enters the core. This injects positive reactivity and slows down the power decay provoked by the increase in T avg . The time scale at which this happens depends on the recirculation time. Here, it is longer than the nominal 4s (Table 2) due to the reduction in the mass flow rate, as can be seen in Fig. 12e. The flow rate decay is exponential, with a time constant of around 10s, as expected. In fact, the momentum forces in the loop are roughly proportional to the square of the flow rate (see, e.g., Equation (23)). Fig. 12c and 12d show the evolution of the salt average and maximum temperature, respectively. T avg can only increase in time, due to the absence of any heat exchange with the external environment. A quick rise of 50K in the first 2:5s is followed by a milder temperature increase, due to the drop in the total power. The maximum salt temperature has a more interesting trend. A first slight increase, due to the increased power-to-flow ratio, is followed by a drop of 30K, due to the abrupt power drop and a first partial homogenization of the temperature in the core cavity, determined by the still quite high turbulent diffusivity. Then, T max rises again, due to the plateau in the power around 7s. The same pattern is repeated in the range 8 À 10s, after the power starts Fig. 11. Distribution of the first (T 1=2 ¼ 55:6s, on the left) and the last (T 1=2 ¼ 0:195s, on the right) family of delayed neutron precursors, obtained at steady-state with the M2 mesh (mid-plane cuts). The short-lived precursors decay mostly within the core cavity, contrary to the long-lived ones. to decrease again. After 20s, both T avg and T max rise at a similar rate (1:68Ks À 1 for the former and 1:63Ks À 1 for the latter), almost constant, due to the almost constant residual power in the system. With these heating-up rates, the grace time before the salt reaches the critical temperature of 1200 ∘ C, at which the reactor incurs structural damage, lies between 235 s and 250 s (evaluated considering T max and T avg , respectively). This is half the grace time of 480s estimated by Brovchenko et al. (2013). Considering the temperature increase rate as constant is too conservative in this case, so longer simulations are necessary to better estimate this important safety parameter. In fact, the grace time strongly influences the design of safety devices such as the freeze-plugs (Tiberga et al., 2019): passive valves of frozen fuel salt that should melt fast enough to allow for the salt draining in large tanks underneath the core (see Fig. 3), to avoid structural damage.
No relevant differences can be noticed between the two meshes and the two time-step sizes, proving the space-time convergence of the results. The differences in the T max trend are to be prescribed to the higher boundary layer resolution of M2, as explained in the previous section, and are always limited to less than 0:3%.
Finally, Fig. 13 shows snapshots of the temperature distribution in the fuel circuit at different times. As time passes, the salt temperature rises and homogenizes, due to the reduction of the power density and the mass flow rate (the difference T max À T avg is reduced to around 25K after 30s). At t ¼ 2:5s one can see that the maximum temperature in the topleft corner of the core cavity is lower than at t ¼ 0, and that the temperature field is smoother.

Conclusions
In this paper, we have presented a novel multi-physics tool for molten salt fast reactors developed at Delft University of Technology.
The code couples an incompressible RANS/k À ε model with a solver for the S N multi-group neutron transport equation and the delayed neutron precursors equations. Decay heat is modeled as well through a set of balance equations that take the drift of the fission products into account. The models are implemented in two in-house codes, based on the discontinuous Galerkin Finite Element method for spatial discretization of the model equations. Second-order implicit schemes are adopted for time discretization. This guarantees high accuracy of simulations. The codes are iterated until convergence in steady-state simulations, while a loose-coupling strategy is adopted during transients. The quantities exchanged between the codes are properly time extrapolated, thus preserving its global second order time accuracy.
To show the capabilities of our multi-physics tool, we have reported the results of preliminary simulations aimed at studying the Molten Salt Fast Reactor (MSFR) at steady-state and during an unprotected Total Loss of Power (TLOP) accident, taking into account the most recent updates to the design brought in the framework of the H2020 SAMOFAR project.
The reactor steady-state conditions revealed that previous studies on the optimization of the core toroidal shape were effective in eliminating hot spots in the core cavity. However, a large recirculation region is still present at inlet, and this induces localized pressure drops and makes the flow less predictable in the core cavity. To avoid this, further improvements of the core design are necessary. Moreover, particular care has to be taken in designing the upper wall, which separates the core from the upper reflector. In fact, it is subject to high temperatures and large thermal gradients, which enhance corrosion effects and induce large mechanical stresses. Finally, the reference fuel composition should be adjusted, reducing the uranium enrichment, thus eliminating the reactivity excess of almost 1000pcm found at steady-state.
During a TLOP accident, the reactor power quickly drops to less than one-third due to the strong negative temperature feedback coefficient, but after 30s it is still at a few hundred MW, due to decay heat. Our code proved able to reproduce the temporarily injection of positive reactivity induced by the re-entrance in core of delayed neutron precursors, which is a peculiar feature of molten salt systems. No conclusions could be retrieved though in terms of grace time before the salt reaches the maximum allowed temperature. Longer simulations are required to analyze this aspect further. This, together with the assessment of the MSFR behavior during a broader spectrum of accidental scenarios is the subject of a future publication, which will also compare our results with those obtained by other partners of the SAMOFAR project. Fig. 13. Evolution of salt temperature distribution in the fuel circuit during the TLOP accident (on the M2 mesh).