In the neutronics calculations we employ the spatial neutron diffusion code (Griffin) and the 0-D PK solver (Squirrel) to model the effect of the flowing fuel and the redistribution of DNPs during steady-state and transient calculations. The solutions of Squirrel, which is a modified PK solver for flowing-fuel MSRs, are verified against the solutions of Griffin’s diffusion solver. Squirrel can account for the time-dependent spatial distribution of the DNPs in the reactor by considering a proper weighting function.

Pronghorn is a coarse-mesh finite-volume thermal-hydraulics solver implemented in the MOOSE framework. It supports the development of nuclear reactors at the engineering scale (Lindsay et al., 2021).

In this work, Pronghorn solves for the fluid flow distribution in the reactor, the DNP scalar field, and the fluid and solid enthalpy fields. The equations solved for thermal-hydraulics modeling in Pronghorn are the single-phase weakly compressible Navier-Stokes and energy conservation in porous media. The porous medium formulation assumes that the liquid fuel and the internal reactor structures occupy the same homogenized volume. It is primarily used to model the reactor core as a homogeneous system. The porosity γ is defined as the ratio between fluid volume and total volume.

In this context, the conservation of mass is given by γ r ∂ ρ t , r ∂ t + ∇ ⋅ ρ t , r U t , r = 0 , (19) with the porosity γ ( r ) , the fluid density ρ ( t , r ) , and the superficial velocity U ( t , r ) .

The conservation of momentum is given by ∂ ρ t , r U t , r ∂ t + ∇ ⋅ ρ t , r 1 γ r U t , r ⊗ U t , r = − γ r ∇ p t , r + ∇ ⋅ μ + ρ ν t ∇ U t , r + ∇ U t , r T + γ r ρ t , r g − W r ρ t , r U t , r + γ r F t , r , (20) where μ is the dynamic viscosity, ν t is the turbulent viscosity using the 0-D turbulent mixing length model (Lindsay et al., 2023), a time-dependent forcing term F ( t , r ) that represents the pumping force, the gravitational vector g , the pressure p ( t , r ) , and the pressure drop coefficient W ( r ) using the Darcy and Forchheimer friction models.

The conservation of energy for the fluid phase is given by ∂ γ r ρ t , r H f t , r ∂ t + ∇ ⋅ ρ t , r H f t , r U t , r = ∇ ⋅ γ r κ f + ρ α t ∇ T t , r + h V T f t , r − T s t , r + γ r q ‴ f t , r , (21) and the conservation of energy for the solid phase is given by ∂ 1 − γ r ρ s t , r H s t , r ∂ t = ∇ ⋅ 1 − γ r κ s ∇ T s t , r + h V T s t , r − T f t , r + 1 − γ r q ‴ s t , r , (22) where H f is the specific enthalpy of the fluid, H s is the specific enthalpy of the solid, T f is the fluid temperature, T s is the solid temperature, ρ s is the solid density, κ f is the thermal conductivity of the fluid, α t is the turbulent heat diffusivity, κ s is the solid thermal conductivity tensor, h V is the fluid-porosity-averaged volumetric heat exchange coefficient between the liquid and solid phases, q ̇ l ‴ is the heat source being deposited directly in liquid fuel (e.g., fission heat source), and q ̇ s ‴ is the heat source in the solid (e.g., residual power production in structures). For an assembly modeled as a porous medium, the fission source is defined as the total power of the assembly divided by its total volume.

The convective heat transfer between the porous core region and the solid core barrel is also included in the model through a computed heat transfer coefficient. The value of the porosity in the core region is 0.222. The pressure in the closed loop is fixed by fixing one pressure point at the reactor outlet. To avoid the appearance of a checkerboard pattern in the pressure field caused by the separation of velocity and pressure, Pronghorn employs the collocated formulation of the Rhie-Chow interpolation method. In this method the velocity flux at the faces is computed using the difference between the direct and cell-interpolated pressure gradients, suppressing the formation of strong pressure gradients at the cell centers, which enables the discontinuities in the porosity and body forces like drag or a pump force to be simulated without resulting in velocity oscillations. A more detailed description of the implementation can be found in Lindsay et al. (2021), and a more general discussion of the topic can be found in Moukalled et al. (2016). Pronghorn also calculates the DNP production, decay, and transport for Griffin Equation 2 and for Squirrel Equation 16.

The MSRE was a thermal spectrum reactor in which liquid fuel salt flowed through graphite moderator channels. The main reactor parameters considered in this work are listed in Table 1. A diagram of the reactor layout, including the core, the primary circuit, the primary heat exchanger, and the secondary system, is shown in Figure 1 (Kerlin et al., 1971a). The MSRE lattice was made of vertical graphite stringers with a 5.08 by 5.08 cm cross section, while the fuel salt flowed through a rectangular channel of 3.05 by 1.016 cm with round corners of radius 0.508 cm on the sides of the stringers (Kedl, 1970). The MSRE core configuration is shown in Figure 2. The left and top-right plots show the MSRE reactor assembly, while the bottom-right plot shows the MSRE 4-halves fuel salt channels with graphite stringers (Jaradat and Ortensi, 2023).

The MSRE’s liquid fuel salt was a FLiBe base salt bearing U and Zr, composed of LiF-BeF₂-ZrF₄-UF₄. The main fissile isotope was ²³⁵U, which was later replaced with ²³³U. The atomic fractions of the fuel salt isotopes are provided in Table 2 (Jaradat, 2021). The thermophysical properties of the fuel salt and solid moderator are provided in Table 3. This data was collected from several design and analysis reports of the MSRE (Gabbard, 1970; Compere et al., 1975). All the material properties considered in this model are temperature independent except for the fuel salt density ρ , which can be represented as a function of the fuel salt temperature T as follows: ρ T = 2263.0 − 0.4798 × T − 923.0 . (23)

The MSRE model was developed using a 2-D axisymmetric domain in R-Z coordinates for both the spatial-dynamics neutronics and thermal-hydraulics calculations based on a model of the MSRE previously developed by Idaho National Laboratory (Schunert et al., 2023; Jaradat and Ortensi, 2023). The spatial-dynamics neutronics model in Griffin uses the multigroup diffusion approximation of the linearized Boltzmann transport equation in a homogeneous medium. It accounts for the effect of the advection- and diffusion-driven drift of the DNPs and their decay in the outer loop by obtaining the DNP distributions from thermal-hydraulic calculations performed in Pronghorn using the porous media approximation.

A cut-through view of the reactor core vessel is provided on the left panel of Figure 3 (Robertson, 1965); it shows the graphite stringers with the fuel channels, the lower and upper plenum, the downcomer, and the reactor core vessel. On the right-hand side of Figure 3, the developed 2-D R-Z model for neutronics and thermal-hydraulics calculations is presented. This model includes the graphite-moderated core, lower and upper plenum, riser, fuel salt pump, primary heat exchanger, outer loop pipe, and the downcomer.

The MSRE reactor geometry was simplified so the simulation and analysis could be efficient. The control rods and the outside regions of the active core, including the reactor vessel, the insulator, the gap between the reactor vessel and insulator, and the thermal shield, were not developed as part of the core configuration. The secondary loop was also not considered in the model. The temperature of the secondary side of the heat exchanger was set to a constant. The azimuthal symmetry of the MSRE model ignores some of the details of the real reactor system, but it preserves the important quantities of the core. The modeled core barrel has a radius of 0.70485 m, while the graphite-moderated core has a radius of 0.6914 m and a height of 1.6637 m. The lower plenum is 0.12954 m high, and the upper plenum is 0.21336 m high. The total salt volume in the reactor is 1.64 m 3 .

The fuel’s total circulation time in the whole system is 25.2 s; the fuel transitions through the core, lower plenum, and upper plenum in 17 s. In the pump block, the fuel is accelerated by a uniformly distributed momentum source in the momentum equation with a total value equal to 1.15 × 1 0 8 m/s². The volumes of the in-core and out-of-core regions reflect the volumes and circulation times reported in Refs. Ball and Kerlin (1965), Robertson (1965), Kerlin et al. (1971b). The fuel circulates upward through the core and is heated by the released fission power. At full power, the fuel temperature rises 23 K from the inlet to the outlet of the core. The fuel leaves the upper plenum through the riser, and heat ( q ‴ ) is removed in the heat exchanger region according to the following relation: q ‴ = α T − T H X , (24) where α is the volumetric heat transfer coefficient with a value of 2.22 × 1 0 7 W/m³K, and T H X is the secondary side heat exchanger temperature, which is equal to 819.15 K. These values were selected to ensure the models matched the reported steady-state full-power fuel salt inlet temperature of 905 K.

A full 3D core model was developed for the MSRE using the Monte Carlo code OpenMC (Romano et al., 2015) to generate multigroup microscopic cross sections by accounting for the neutron spectral changes in the reactor core. This cross-section-generation model was adopted from Refs. Jaradat (2021), Jaradat and Sik Yang (2023). Additionally, the multigroup energy structure and the number of groups used to generate cross sections were selected based on a previous study (Jaradat, 2021; Jaradat and Sik Yang, 2023) that was performed to optimize the energy group structure and the number of groups. In this work, a 16-energy-group-structure was used to generate multigroup cross sections at several fuel salt and moderator temperatures using ENDF/B-VII.1 neutron data (Chadwick et al., 2011).

The fuel salt and graphite volume fractions in the reactor core region are 22.29% and 77.71%, respectively. This ratio is considered in obtaining cross-section homogenization and in the porous flow model implemented in Pronghorn. Additionally, the model considers the heat deposited in the solid graphite due to neutron thermalization and gamma heating. The considered value is 5.467 % of the total power produced in the reactor, which was calculated from previous OpenMC simulations presented in Jaradat (2021). The delayed neutron fraction and the decay constant for six delayed neutron groups of the ²³⁵U fuel and ²³³U fuel were also obtained from the same simulation, as reported in Table 4.

The multiphysics model of the MSRE was developed using the MOOSE MultiApps system, in which the coupling between the thermal-hydraulics solver and neutronics solvers is designed to allow for flexibility in exchanging problem variables using different neutronics models coupled to the same thermal model. Figure 4 shows the coupling scheme and the transferred parameters between the main application (Pronghorn) and the sub-applications (Griffin and Squirrel).

Pronghorn is the main application, and it solves the Navier-Stokes equations, the DNP field for Griffin, and the reduced DNP field for Squirrel. Griffin provides the power density and fission source distributions to Pronghorn. The power density distribution is used to obtain fuel and moderator temperatures, density, pressure, and velocity fields. The fuel and moderator temperatures and the fuel density are used to update the cross sections for Griffin calculations. The fission source distribution is used in the calculation of the DNPs distribution to determine the delayed neutron source and reconstruct the total fission neutron source in the core region, considering their decay in the outer loop.

Squirrel relies on the initial neutron flux, power density, fission source density, and the kinetics parameters calculated by the steady-state Griffin simulation. Squirrel scales the power density and neutron flux during the transient calculations with a power amplitude function, where the steady-state power distribution or shape function is used throughout the entire transient simulation.

Since the neutronics solvers follow two different approaches, the key aspect of the coupling is to maintain the flexibility to swap each solver in each transient while maintaining the same Pronghorn simulation. This approach eliminates differences in the performance of the diffusion based neutronics solvers and the PK solver originating from thermal-hydraulics solutions.

In this work, we performed validation tests encompassing steady-state, pump start-up, pump coast-down, and natural circulation tests. The details of the coupling need to be different for each simulation. In this paper, the following three cases are evaluated:

1. Steady state: The eigenvalue problem is solved by Griffin, considering the steady-state diffusion equation with DNP concentration distributions based on the thermal field and DNP concentration distributions provided by Pronghorn, assuming a constant power level. Pronghorn solves a time-dependent problem until a steady-state solution is achieved at each Picard iteration.

In this sub-section, we provide the steady-state solutions of the MSRE model obtained from coupled Griffin and Pronghorn calculations for isothermal temperature coefficients and effective delayed neutron fractions of stationary and circulating fuels, along with steady-state distributions of the power, temperature, and DNP concentrations, which Squirrel needs to initiate coupled transient calculations.

Table 5 provides the calculated values of the isothermal temperature coefficients of the ²³⁵235U and ²³³U fuel salts. Due to the temperature dependency assumed for density in Equation 23, the temperature coefficients involve both the density and Doppler feedback. The isothermal temperature coefficient is calculated by globally varying the temperature in the reactor between 850 and 1000 K using Equation 5 and considering stationary fuel. The calculated values are also compared to values measured during operations of the MSRE, as reported in Prince et al. (1968).

Figure 5 shows the reactor field variables obtained from coupled steady-state Griffin-Pronghorn calculations for power density, fuel salt temperature, and fuel salt velocity. The majority of the power is generated at the center of the active core region, while the outer core does not get heated. There are two apparent peaks in power density in the upper and lower plenum regions due to the massive increase in fuel salt volume and the absence of the graphite moderator in these regions.

The middle plot of Figure 5 shows the fuel salt temperature distribution. The temperature distribution is skewed toward the top of the core, as expected, with an approximate maximum fuel temperature of 950 K at the top of the graphite block, close to the core’s center. In the upper plenum, colder salt from the outer regions of the core is drawn toward the riser located above the core’s center. There, the cooler salt mixes with the hotter salt from the core’s center, significantly reducing the fuel salt temperature.

The right plot of Figure 5 shows the fuel salt superficial velocity with streamlines indicating the fuel salt path. The salt enters the lower plenum from the downcomer and flows radially inward into the graphite block region, then flows upward without any lateral movement. In the upper plenum, the salt flows again radially inward toward the riser. In the riser, the flow rate is the highest since the cross-sectional area is the smallest. In the pump, the salt is forced to move through the elbow into the downcomer and back to the lower plenum again, with a total circulation time of 25.2 s.

For the flowing-fuel case, the DNPs are advected in the reactor core and the primary loop, resulting in the redistribution of the DNPs based on their half-lives, which leads to reactivity losses as DNPs decay in regions of lower importance. Outside of the active core, the importance of the DNPs is almost zero. Figure 6 depicts the distribution of the six groups of DNPs within the fuel salt for the reactor operating with ²³³U fuel salt. The longest-lived DNP group (Group 1) decays outside the active core region and peaks toward the top region of the core. The undecayed DNPs flow back into the core through the lower plenum, resulting in almost homogeneous mixing in the core. Groups 2 and 3 of the DNPs decay mostly in the outer region of the core, resulting in larger delayed neutron losses. For the shortest-lived group (Group 6), the DNPs decay at almost the same location as their initial position and peak toward the core center, resulting in smaller delayed neutron losses.

Table 6 provides the effective delayed neutron fraction β of stationary and flowing fuels, and the loss in effective delayed neutron fraction β loss for both ²³⁵U and ²³³U fuel salts. The reported MSRE values of the ²³⁵U fuel salt are from Prince et al. (1968); reactivity loss measurements were not reported for the ²³³U fuel salt. The effective delayed neutron fraction of the ²³⁵U fuel salt is 654 pcm, which is more than twice that of the ²³³U fuel salt, which is 298 pcm. The reactivity losses are almost 220 pcm for ²³⁵U fuel salt and 110 pcm for ²³³U fuel salt, with the relative loss in reactivity being similar for both fuels at around 35%. Squirrel-calculated reactivity losses due to fuel flow agree very well with Griffin values, with a difference of 3 pcm in both fuel cases. Also, both codes’ values agree well with the measured reactivity losses of 212 pcm for ²³⁵U fuel salt. For the stationary fuel case, Squirrel relies on the effective delayed neutron fraction calculated by Griffin.

In this section we discuss the results of the validation test against the experimental data of the pump start-up and coast-down tests. These tests were performed with ²³⁵U fuel salt at a low power level (zero power or cold conditions). A detailed description of the experimental setup of this transient can be found in Prince et al. (1968). The mass flow rate of the fuel salt was adjusted through the primary salt pump, which resulted in reactivity perturbations. These reactivity changes were compensated for by adjusting the control rod position to match positive or negative reactivity insertion and by maintaining the reactor at a constant low power level (around 10 W) throughout the test. Therefore, the temperature reactivity feedback could be neglected in the simulations. The reported position of the control rods represents the reactivity change due to DNP redistribution and decay with changes in the fuel salt flow rate.

During the pump start-up test, as the fuel salt mass flow rate started increasing, the DNPs flowed outside the active core region and decayed in regions of lower importance and in the outer loop, resulting in increased reactivity losses. In the pump start-up simulation, the pump force in the thermal hydraulics model was adjusted to achieve a fuel mass flow rate change from 0% to 100% in approximately 8 s, while in the pump coast-down simulation the pump force was adjusted to achieve a fuel mass flow rate change from 100% to 0% in approximately 20 s. The slower coast-down of the fuel flow can be explained by the hydraulic inertia that the system had at maximum flow rate. The measured and calculated mass flow rates for both test cases are displayed in Figure 7, along with the normalized pump force values used in Pronghorn.

The estimated uncertainty for the measured reactivity considered in this work originates from the following parameters:

• Control rod position indicator: The uncertainty in the control rod position measurement during the zero power tests was estimated to be ± 0.01 inches.

Figures 8, 9 display, on the left-hand side, the experimental values of the reactivity along with simulation results for pump start-up and coast-down, respectively. During the pump start-up test, the redistribution and decay of the DNPs led to a maximum reactivity loss of 284 pcm after 13 s. As the undecayed DNPs reentered the active core region, the measured reactivity loss was reduced again and started oscillating at approximately 212 pcm until it reached steady-state conditions. During the pump coast-down test, the reactivity loss decreased as the fuel salt flow rate decreased due to the fact that the DNPs decay rate in the active core region increases. Finally, the reactivity loss reaches a zero value at zero-flow conditions, as the system returns to its initial state.

The simulation results of both codes show very good agreement for both pump tests. For the pump start-up test, both codes are in good agreement with the measured experimental values and captured the initial increase in reactivity loss and the final oscillation in reactivity loss, although they failed to capture the measured peak reactivity loss of 284 pcm by 34 pcm. The reason for this discrepancy is not fully understood, but it could be related to one or more of the following:

• Channel effect: The fuel channels are not resolved in the homogenized core model investigated here. The flow pattern in the fuel channels of the MSRE might cause a larger reactivity effect that the porous media approximation can not capture.

For the pump coast-down test, both codes are in good agreement with the measured values for the entire period of the test. Additionally the difference in the simulation results obtained with Squirrel and Griffin is shown on the right axis of Figures 8, 9. The maximum difference between the two codes is 3.5 pcm for the pump start-up and 2.5 pcm for the coast-down. Since DNP advection has a minimal effect on the fundamental mode, good agreement between the codes is expected, as the assumption posed in Equation 9 is mostly fulfilled.

These transients increase our confidence that the circulation time and fuel volume are used properly, and they demonstrate that the neutronics solvers correctly evaluated the effect the shifted spatial distribution of the DNPs had on the reactivity.

In Figure 10, the density of the longest-lived DNP group (Group 1) is shown during the pump start-up test at different times up to 50 s. Initially, the DNP shape matches the fission source density in the core since there is no fuel flow. As the pump starts, the DNPs move upward into the riser (at approximately 10 s) and then down through the downcomer. At 20 s, the DNPs reenter the core from the bottom. This return of the DNPs into the core is responsible for the reduction in reactivity loss observed in Figure 8. The DNPs that are present at 30 s are advected out of the core, and at this time the returning DNPs have not yet reached the upper region of the core, leaving a region of lower DNP density in the middle of the reactor. At 50 s, the maximum DNP density is inside the core. Throughout the process, the DNPs are being mixed in the core.

In Figure 11, the density of the longest-lived DNP group (Group 1) is shown during the pump coast-down test at different times up to 50 s. Initially, the DNPs are well mixed throughout the core. As the fuel flow is reduced (pump force reduced), more DNPs decay in the active core region, resulting in a decrease in reactivity losses. At 20 s, significantly fewer DNPs reenter the core from the bottom, clearly showing that fewer DNPs are advected through the outer core. At 50 s, the DNP shape matches the fission source density in the core. During the coast-down transient, the DNP density in the core increases steadily with minor oscillations.

In this test, the MSRE’s natural circulation was examined. The fuel pump was turned off and the fuel salt circulated through the system due to natural circulation. The main test characteristics can be summarized as follows:

• The reactor was operated with ²³³U fuel salt with the primary pump turned off.

The reduction in the fuel salt temperature led to a positive reactivity insertion (due to the negative temperature coefficient); therefore, the power of the reactor increased as the temperature decreased. This increase in power and the resulting heating of the fuel increased the fuel temperature. The hotter fuel caused the power to stabilize at a new level. This process of reducing the fuel temperature was repeated several times, and multiple power plateaus were reached. The maximum power was 351 kW, and the final power was 324 kW. During the transient, the reactor power developed fully unprotected, with the control rod remaining in its initial position, and the reactor was fully driven by temperature feedback during the test. The effect of the DNP movement in the circulating fuel was minimal due to the very low fuel-flow rate (Rosenthal et al., 1969).

During the simulation, the inlet temperature was prescribed by adjusting the secondary temperature of the heat exchanger. The model used is not able to reproduce the reported fuel flow with natural circulation, most likely because it lacks the correct height difference between the core and the heat exchanger. Therefore, the mass flow rate was adjusted by changing the pump force in Pronghorn. The transient was induced by the prescribed change in inlet temperature and the adjustment of the pump force. With colder fuel entering the core, the power increases due to the positive temperature feedback. The increasing power heats up the core. The reactor reaches a steady power when the inlet temperature is kept constant, which happened multiple times during the transient.

Figure 12 shows the inlet mass flow rate and the inlet and outlet fuel temperatures alongside the measured inlet and outlet fuel temperatures reported in Rosenthal et al. (1969). Additionally, Figure 13 shows the reactor power evolution as measured and calculated by Griffin and Squirrel. The simulation results and the experimental data match reasonably well for the entire transient. Both the Griffin and Squirrel solutions match very well throughout the transient, with the maximum observed difference in total reactor power being 7.5 kW. The good agreement between the codes can be explained by the transient’s slow evolution and the minimal deformation of the fundamental mode. The largest difference in power measurements for the simulation and experimental data is seen at the beginning of the transient. After that, the power is well matched, resulting in good agreement between the simulated and experimental outlet temperatures. The assumption of thermal equilibrium conditions (steady state) at the beginning of the simulation may explain the temperature differences in the data for the beginning of the transient.

Figure 14 shows the spatial power density in the reactor core at different times during the natural circulation test, as calculated by Griffin. Initially, the power density is low at 19.0 kW/m³; during the transient the power density increases up to 1.5 MW/m³ at 240 min. After that, it decreases slightly. At all times, the power density is highest in the center region of the core, and there is no measurable shift in the shape of the power density. Figure 15 shows the fuel salt temperature at different times during the natural circulation test. At the beginning of the transient, the inlet temperature increased up to 910 K, and it reached 960 K at the end of the transient. The difference in the simulated power obtained with Squirrel and Griffin is shown on the right axis of Figure 15. The maximum difference between the codes is 7.5 kW.