Multiphysics Modeling of Precursors in Molten Salt Fast Reactors using PROTEUS and Nek5000

The goal of this work was to calculate the impact of the delayed neutron precursor drift in fast spectrum Molten Salt Reactors (MSRs) using coupled solutions from the neutronics code PROTEUS and the computational fluid dynamics code Nek5000. Specifically, using a multiphysics approach to solve the effective delayed neutron fraction (βeff) or delayed neutron precursor distribution for reactors with flowing fuel salts would provide valuable information for transient simulations and safety assessments. Given the multiple options for the flux solution and geometric resolution/fidelity in PROTEUS, two approaches were developed and applied to various test cases: PROTEUS-NODAL/Nek5000 and PROTEUS-SN/Nek5000. For the former, the precursors are tracked in the built-in precursor drift model in PROTEUS-NODAL, whereas in the latter, Nek5000 directly tracks the precursors. Both approaches were used to solve a single test channel problem and showed excellent agreement in the calculated βeff. Separately, a 3D hourglass-shaped core was modeled using the PROTEUS-SN/Nek5000 approach. This problem was designed to demonstrate the capability of the discrete ordinates (SN) solver and Nek5000 to model complex core designs with axially varying geometries and the ability for Nek5000 to track the precursors and calculate the resulting βeff. In addition, the Nek5000 calculations revealed the presence of recirculation zones in the hourglass design, which could lead to significant temperatures in the fuel salt and surrounding materials. These first coupled solutions show why these approaches may be necessary for not only predicting the precursor drift effect in fast MSRs but also for reactor design and performance assessments.


INTRODUCTION
In recent years, there has been a resurgence in interest in molten salt reactors (MSR) as indicated by many reactor vendor startups and established companies pursuing the development and deployment of these nontraditional reactor technologies. To support this interest, the Nuclear Energy Advanced Modeling and Simulation (NEAMS) program within the U.S. Department of Energy's Office of Nuclear Energy (DOE-NE) has been developing and applying its computational tools that can be used to perform both design and high-fidelity modeling and simulation work for these MSRs. One of several problems that the NEAMS program is looking to address is using multiphysics approaches to model the spatial distribution of delayed neutron precursors (DNPs) in fast spectrum MSRs. Specifically, this is important for designs that involve flowing fuel salt that leaves the core region because the loss of these precursors from the core reduces the effective delayed neutron fraction, βeff, of the system and can potentially impact safety performance.
To tackle this complex problem, two multiphysics approaches were developed to track the movement of these DNPs, each with a thermal-hydraulics component and a neutronics component. For the thermalhydraulics component of both approaches, the computational fluid dynamics (CFD) code Nek5000 [1] was used. CFD was deemed necessary to account for the complex fuel salt cross flow inside a fast spectrum MSR core (or tank) in which flow is not restricted to channels. In terms of the reactor physics component, the first multiphysics approach uses the NODAL flux solver [2] within the neutronics code PROTEUS [2,3,4], while the second approach uses the second order discrete ordinates or SN flux solver [3] within PROTEUS. These two multiphysics approaches, (1) PROTEUS-NODAL/Nek5000 and (2) PROTEUS-SN/Nek5000, were developed simultaneously to enable a larger range of modeling suitability for different fast MSR designs and to provide both medium and high fidelity solutions given the user's preference or computational resources. PROTEUS also includes a third flux solver using the Method of Characteristics or MOC [4], but coupling it to Nek5000 was not pursued since the SN solver provided similarly high fidelity solutions but also has the capability to model axially non-uniform geometries. The PROTEUS-NODAL option already has the complex precursor drift model built-in and has a light computational burden that may be preferred for design calculations, so it was a good choice to couple with Nek5000, with the understanding that some mesh interpolation would be required. However, it is limited to diffusion/SP3 and simpler regular geometries (no finite element mesh) so PROTEUS-SN was also chosen for coupling to offer an option with transport and as-built 3D geometries. Both approaches were applied to a single channel test problem and the PROTEUS-SN/Nek5000 approach was applied to a 3D problem consisting of an hourglass-shaped fast MSR.
For the single channel test problem, 33-group cross sections were generated for an infinite homogeneous mixture of fuel using the cross section generation code MC 2 -3 [5]. For the hourglass problem, a two-step cross section generation procedure was adopted which involved generating spatially-and spectrallydependent fine group cross sections based on a RZ (pseudo 3D) representation of the core and then using the local spectra to collapse the cross sections into the 33 groups. ENDF/B-VII.0 cross sections were used. Due to the complexity of the task of coupling these codes, for this initial study, a few approximations were made to test the feasibility of these approaches on a few sample problems. No cross sections were updated with new temperatures calculated from CFD; this will be done in future work. The impact of delayed neutrons on the flux and power distributions was neglected; this impact was tested on the single channel test problem and it was found to be indeed negligible (a maximum of 0.16% difference in local power).  Fig. 1 shows a schematic for the coupled solution approach for PROTEUS-NODAL/Nek5000. Initially, PROTEUS-NODAL is used to solve the steady-state power distribution of the fast MSR core and passes this information via script to Nek5000. Nek5000 then calculates the velocity field of the fuel salt for the entire core and passes it back to PROTEUS-NODAL. PROTEUS-NODAL then uses this velocity field as input in its precursor drift model [6] and neutronics solver to iteratively calculate the precursor distribution and initial precursor source distribution until convergence. After convergence, the updated power distribution can then be sent back to Nek5000 to confirm convergence. Fig. 2 shows the PROTEUS-SN/Nek5000 approach, which uses PROTEUS-SN to solve for the initial power distribution and initial precursor source distribution under steady-state conditions. This information is then passed to Nek5000, which performs not only the heat transfer and fluid dynamics calculations but also the precursor tracking. For this approach, the precursor distribution calculated by Nek5000 is the final answer because feeding this information back to PROTEUS-SN to recalculate the flux (rigorous solution) would require significant code development, and it was previously stated that the impact of the precursors on the flux solution is assumed to be minimal.

SINGLE CHANNEL TEST
A simple test case was developed to identify the requirements for information transfer between PROTEUS and Nek5000. For this purpose, the geometry of the test case was intentionally made simple. It consists of only a long square tube with length of 128 cm and width of 4 cm. The square tube contains only fuel salt composed of UCl3 and NaCl with 20% enriched 235 U. The density of the fuel salt is assumed to be constant. The thermal power is selected to be 33.25 kW with the intention to facilitate the Nek5000 calculation (at this power level, the flow is in the laminar regime). A single precursor family model is defined for simplicity; extension to 6 or more families is simple for the NODAL approach but requires additional calculations to generate the precursor source term for the SN approach. It was assumed that the delayed neutron fraction is 650 pcm, the half-life is 5 seconds, and the time that the fuel salt spends outside the core is 10 seconds.
Both the PROTEUS-NODAL/Nek5000 and PROTEUS-SN/Nek5000 approaches were applied to this test case. In addition, an OpenMC [7] model of the test case with stationary fuel was also developed for verification purposes. The calculated eigenvalues and effective delayed neutron fractions with and without flow, i.e., precursor drift, for this test case are shown in Table I. The precursor concentrations corresponding to the solutions in Table I are shown in Fig. 3. These results will be explained in more detail in the following sub-sections, but the main observation is that the two multiphysics approaches showed excellent agreement in their calculated βeff values (130 and 140 pcm), which is reasonable given the flow rate. PROTEUS-SN/Nek5000** see note** 140 * A βeff calculation was not automically generated in PROTEUS-SN ** Since there is currently no feedback from Nek5000 to PROTEUS-SN in this approach, the k-eff with flow was not re-calculated (N1) (N2) (S1) (S2) Figure 3. Precursor distributions corresponding to cases N1, N2, S1, and S2 from Table I.

PROTEUS-NODAL + Nek5000 Result
For this test case, the PROTEUS-NODAL calculated eigenvalues for Step 1 and 2 of the approach as shown in Fig. 1 are listed in Table I as Cases N1 and N2, respectively. The relatively large difference between the OpenMC k-eff and PROTEUS-NODAL (N1) result is largely due to the diffusion approximation in NODAL. In contrast, the no-flow PROTEUS-SN (S1) and OpenMC solutions have a difference of only 59 pcm. However, with a more realistic and larger reactor core, this difference between the N1 and OpenMC results would be reduced. The ~600 pcm lower k-eff in the Case N2 result compared to that from Case N1 indicates the impacts of using NODAL's precursor drift model and the velocity profile from Nek5000. Under the assumptions of this test, most of the precursors decay outside of the core. This is also illustrated in Fig. 3b, which shows that when the precursor drift model was enabled, the precursor distribution concentrates near the outlet of the core region instead of the axial center, as was the case for the stationary fuel as shown in Fig. 3a. The βeff also decreases significantly from 650 to 130 pcm for the same reason. The large difference in the precursor distribution between Cases N1 and N2 for the PROTEUS-NODAL calculations has little impact on the power distribution for the test problem. The maximum relative difference in power distribution is about 0.16% and occurs near the outlet of the core. This shows that the power distribution is mainly dominated by the prompt neutron flux, as expected. The small fraction of the delayed neutron flux contributes insignificantly to the total power distribution. Thus, the iteration stops at this point or may not need to be performed at all as an approximation.
In Nek5000, the fuel density, viscosity, thermal conductivity, and specific heat were all considered constant for this analysis and the Reynolds number was maintained at a low value to keep the flow laminar. All transverse surfaces used no-slip wall boundary conditions for the fluid velocity and insulated boundaries for the temperature and precursor family. The upper and lower surfaces used a recycling boundary condition for the velocity and inlet/outlet conditions for the temperature and neutron precursor. Due to the assumption of constant fuel density and viscosity, the velocity is fully decoupled from the temperature solution. For the PROTEUS-NODAL/Nek5000 coupled solution, the velocity field predicted by Nek5000, using the initial power distribution from PROTEUS-NODAL, was provided back to PROTEUS-NODAL to predict the βeff according to the algorithm shown in Fig. 1.

PROTEUS-SN + Nek5000 Result
Unlike PROTEUS-NODAL, there is currently no precursor drift model implemented in PROTEUS-SN. Thus, the precursor distribution was solved in Nek5000 based on the precursor source distribution calculated by PROTEUS-SN, as shown in Fig. 3c, as well as the power distribution. The precursor source distribution is currently not a direct output from PROTEUS-SN calculations so some post-processing was required. The vertex-wise multi-group neutron flux from PROTEUS-SN outputs were used to calculate the vertex-wise precursor source. It needs the isotopic fission cross section from the cross section dataset ISOTXS and the delayed neutron information from the DLAYXS file. In addition, a Python script was developed to extract the ߥΣ ߶ nu-fission rates from the .h5 file generated by PROTEUS-SN along with the coordinates of the vertices. The nu-fission rate is multiplied by the delayed neutron fraction for the single precursor family assumed for the test problem to yield the precursor source distribution. All these data are written to ASCII file that is passed to Nek5000.
In order to facilitate the data transfer from PROTEUS to Nek5000, a python script was developed to interpret the data structure of a Nek5000 binary restart file. By using the developed script, the ASCII data from PROTEUS can be translated into the Nek5000 native restart file format. This allows Nek5000 to use its own scalable spectral interpolation tool to interpolate heat and precursor generation rates onto its own mesh, allowing both PROTEUS and Nek5000 to use their own optimal meshes. This makes the data transfer process, as well as the class of problems that can be simulated significantly more flexible and robust.

Test Problem Summary
This simple exercise was completed using both coupling approaches that obtained similar answers for the effective delayed neutron fraction and helped shed light on the data exchange requirements. As a result of this exercise, two external scripts were developed to facilitate the coupling between PROTEUS and Nek5000 that can be used for realistic cores. This single channel test did not include cross flow, which is a major phenomenon in fast MSRs without channels and necessitates the use of CFD codes like Nek5000. The 3D core example in the proceeding section will include the cross flow effect among other details that were approximated here.

3D HOURGLASS CORE
An hourglass shaped core developed to evaluate the MSR coupling capability between PROTEUS-SN and Nek5000, as the overall design shape is similar in concept to some of those being explored by private industry. It also represents a good test case, as cross flow will be significant due to the complex geometry, and the axial variation of the geometry requires PROTEUS-SN to model the core as-is without axial extrusion of a 2D plane. The core shape is represented by a cylinder of 2.2 m diameter and 3 m height with circular sections removed from the tops and sides. This is shown in Figure 4a. The sections were removed from the tops and sides in order to reduce recirculation regions.
The simulated test case was for a reactor of 10 MW total power using a chloride-based fuel salt of density 3430 kg/m 3 with an inlet temperature of 750°C. The chosen half-life for the delayed neutron precursor family was 5 seconds. These conditions are consistent with the fuel properties used in the single channel test problem. In the hourglass core design, liquid fuel is pumped in radially at the bottom of the core. It then flows up through the center and exits radially at the top of the core. This results in a complex flow pattern. Flow through the core was simulated at a low Reynolds number of 1000. This was done in order to keep the thermal hydraulics analysis relatively simple and more effort focused on developing the coupled PROTEUS-SN/Nek5000 capability.

PROTEUS-SN + Nek5000 Result
PROTEUS-SN was used calculate the power distribution, as shown in Fig. 4b, which was passed along with the initial precursor source distribution to Nek5000. A hexahedral mesh was generated to facilitate the mesh-wise data transfer from PROTEUS-SN to Nek5000. For the Nek5000 simulation, a cylindrical mesh was generated using preNek, the native mesh generation software included with the Nek5000 distribution. The cylindrical mesh was then distorted into the hourglass shape and smoothed at runtime using Nek5000's built-in mesh manipulation capability. The previously described python script was used to translate the data from PROTEUS-SN to the Nek5000 native restart file format. This allowed the power and precursor generation rates to be interpolated from the PROTEUS domain onto the Nek5000 domain. These were then used as the source terms for the energy conservation and precursor family transport equations.
Results for the velocity profile in the core are shown in Fig. 5. Blue regions in the velocity magnitude represent where the flow is moving slowly. The regions where the core has been removed, resulting in the unique shape, are observed to be where the flow would move slowest. This is significant in liquid fueled MSRs, as the fuel flowing through the areas with lowest flow rates have the longest residence time in the core and, because the fuel is also the coolant, these areas will become the hottest. The temperature and precursor concentration profiles are shown in Fig. 6. It can be seen by comparing to the velocity profile that the regions of highest temperature coincide with the regions of the lowest velocity, and despite core design attempts to minimize recirculation zones, areas of very high temperature associated with recirculation zones are present. These areas arise because the heat cannot be convected away by the flow and must rely on conduction for cooling. In particular, the two most significant areas are at the bottom in the center of the core and on the outside wall just below the outlet. In these regions, the temperature of the fuel salt is predicted to rise to over 1000 °C, while the regions with flow remain practically at the inlet temperature. This underscores the importance of accurately predicting recirculation zones. A similar trend is shown in the precursor concentration profile, shown in Fig. 6b. However, the decay term included as part of the transport equation results in significantly different overall distributions. The precursor concentration does not build up as significantly in the recirculation zones. Using the presented conditions, the βeff has been calculated for this coupled simulation as 260 pcm.

CONCLUSIONS
To solve the effective delayed neutron fraction, βeff, or precursor distribution resulting from precursor drift in fast spectrum MSRs, two high-fidelity coupled solutions involving PROTEUS-NODAL, PROTEUS-SN, and Nek5000 were developed and applied to various test cases. Both approaches were used to solve a single test channel problem and showed excellent agreement in the calculated βeff. Separately, a 3D hourglass-shaped core was modeled using the PROTEUS-SN/Nek5000 approach and its βeff was calculated. These calculations also revealed the presence of recirculation zones in the core that could lead to significant temperatures in the fuel salt and surrounding materials.