1. Introduction
In 2000, the United States Department of Energy (DOE) identified the need for advanced nuclear technologies and convened the Generation IV International Forum (GIF) to expand the mission of nuclear energy and further the development of future nuclear systems. The next generation advanced reactors seek to achieve (1) optimal resource utilization and waste minimization; (2) high proliferation resistance; (3) safety and reliability through passive safety systems that eliminate the likelihood of fission product release during accident scenarios [
1].
Liquid-fueled molten-salt reactors (LFMSRs) have the potential to reduce the cost of electricity generation with an overall cycle efficiency that is higher than conventional reactors [
2]. There is also no need to manufacture fuel pellets and rods since online refueling ensures that the reactor can continue operating. These objectives are realized by LFMSR systems since different salt compositions and fuel forms provide flexibility in temperature margins, neutron spectrum, and power density [
1]. Other advantages to LFMSR systems include high-temperature process heat, minimal high-level nuclear waste, negative void and thermal reactivity feedback, and low-pressure operation.
Researchers need advanced modeling approaches to obtain results with a reasonable computational efficiency while maintaining accuracy owing to the unique characteristic of the design such as liquefied fuel, complex core geometry, fast spectrum, and moving precursors. The key is to identify and model these multi-physics phenomena that impose the potential challenges to modeling and simulation of the reactors. In this study, GeN-Foam [
3] is adopted as the main modeling and simulation tool for LFMSR. GeN-Foam is based on the Open-Source Field Operation and Manipulation (OpenFOAM) C++ library [
4]. The LFMSR can be modeled using the GeN-Foam code to avoid a detailed Reynolds-Averaged Navier-Stokes (RANS) model of components such as heat exchangers and pumps. The flow and heat transfer in these components can be modeled using a porous k-epsilon approach, thereby saving a lot of computational expense and time.
Most of the recent research efforts have been focused on studying the LFMSR, developing and verifying multi-physics codes for the LFMSR simulation [
5]. The Safety Assessment of the molten-salt fast reactor (SAMOFAR) project coordinates research efforts on the molten-salt fast reactor (MSFR) concept in Europe. The MSFR is a fast-spectrum breeder reactor with a large negative power coefficient. The fast spectrum allows the reduction in the reprocessing requirements and enables a better reactor breeding ratio [
6]. The reference MSFR concept is a 3000 MW
th reactor with three different circuits: the fuel circuit, the intermediate circuit, and the power conversion system.
Figure 1 illustrates a basic layout of the MSFR. The components of the primary fuel circuit are the fuel salt, the core cavity, the heat exchangers, the pumps, the inlet and outlet pipes, and the salt-bubble separators. The initial composition of the MSFR is 77.5 mol% of LiF and a 22.5% mol mixture of AcF
3 and AcF
4 (Ac indicates actinides, mostly
232Th or
233U). The total volume of fuel salt is about 18 m
3 in the fuel loop, and the average core salt temperature is approximately 725 °C. The fuel salt, upon exiting the core cavity from the top, flows into 16 different loops, each containing a pump and heat exchangers located around the core. The fuel salt circulates the loop in about 3–4 s [
7].
The MSFR design has some potential advantages as discussed but also poses new challenges in terms of reactor design and safety. The presence of several phenomena inherent to molten salts that do not exist in other reactor types requires the development of better modeling tools and new methodologies to perform the safety evaluation. One such research project is the Euratom Evaluation and Viability of Liquid Fuel Fast Reactor System (EVOL) project. In their recent work, a cylindrical core with the diameter equal to the height, surrounded by steel reflectors axially and a fertile blanket, a boron carbide layer, and a reflector radially is studied. The fuel salt used in their study is a fluoride-based salt and is a mixture of lithium fluoride, thorium fluoride, and actinide fluoride salts.
Though there are optimized core geometry designs available, rigorous studies have not yet been performed in this regard. The optimized core geometry available in the literature [
8,
9] is only based on one factor, i.e., hot spots inside the core. These hot spots occur mainly due to the formation of recirculation and stagnation zones which increase the salt temperature considerably. This undesirable phenomenon needs to be avoided or minimized as it leads to higher temperature gradients inside the core. There are also other factors that need to be considered when optimizing the design. Extensive studies need to be performed to optimize the capital cost and performance of major systems. These analyses must consider the multi-physics design impacts on plant economic performance such as pumping requirements, sizing of heat exchangers, sizing of system pressure boundaries, and coolant inventory requirements.
The objective of this work is to develop a framework for multi-physics modeling, design optimization, and uncertainty propagation of fast-spectrum liquid-fuel molten-salt reactors. GeN-Foam is capable of performing coupled fine-/coarse-mesh thermal-hydraulics calculations, multi-group neutron diffusion calculations, and thermo-mechanics calculations using both structured and unstructured meshes. This multi-fidelity approach, together with DAKOTA, will address the need for a new structure for advanced reactor analysis.
This article next presents the materials and methods used in its associated GeN-Foam simulations, optimization, and uncertainty analysis which are loosely based on a work by Brovchenko et al. [
10]. The results are divided into five main sections: primary loop optimization, steady-state results, transient results, uncertainty propagation, and sensitivity analysis. Core design, heat exchanger placement, and pump position are included in the primary loop optimization results. The uncertainty propagation includes scalar analysis, functional analysis, and three-dimensional analysis.
3. Results
The results of this work have been structured into five sections according to separate analyses. The following sections describe each of these analyses, with some being divided into several sub-sections.
3.1. Primary Loop Optimization
Three components were optimized for the proposed LFMSR design: reactor core dimensions, heat exchanger size, and pump position. The optimization of each of these components is described in the following sections.
3.1.1. Sizing of Primary Heat Exchanger
The Printed Circuit Heat Exchanger (PCHE) has been chosen for the proposed LFMSR system because a conventional shell and tube heat exchanger may not be sufficient to achieve the required heat transfer rate in the proposed design. Printed circuit heat exchangers use chemical etching and a diffusion bonding process to create a heat exchanger core with no joints, welds, or points of failure [
23]. Because of this, the reactor core dimensions become the main set of constraints for PCHE size. Several design iterations were performed by building 3D CAD models in SOLIDWORKS for both single-loop and sixteen-loop models to tentatively pick the PCHE dimensions. The assumptions in the proposed PCHE design include the neglection of heat conduction along axial direction, heat transfer from heat exchanger walls, and counter-flow heat transfer.
Table 5 lists the proposed PCHE design parameters. Thermophysical properties of the salt on the primary and secondary side are known, and the temperatures on both sides of the PCHE are fixed. The secondary inlet temperature is fixed as well, and the secondary outlet temperature is varied to determine the total number of channels needed. This is to verify that the proposed dimensions of the PCHE system are well within the design limits imposed by the number of channels. Reactor power was held constant during steady-state calculations at 2400 MW; for the 1/16 slice modeled, 150 MW was applied to the computational domain.
For a fixed secondary outlet temperature (870 K), the overall heat transfer coefficient and the pressure losses are calculated with GeN-Foam and included in the proposed LFMSR design.
Table 6 shows the remaining proposed PCHE design parameters. A constant heat transfer coefficient has been assumed for the heat exchanger region in all simulations.
The size of the primary heat exchanger depends on four main parameters: reactor power, reactor core dimensions, and thermophysical properties of primary and secondary fluid [
24]. The reactor power, core dimensions, and thermophysical properties are fixed. Heat exchanger design for a liquid-fueled reactor must also consider how to limit the residence time of fuel salt in inactive regions, the possibility of low salt inventory, lower pressure drops, and resistance to high temperature (973–1073 K) [
25].
3.1.2. Pump Position
Pump position in the reactor system plays a significant role as it dictates the flow characteristics in the reactor system. The steady-state velocity magnitude distribution is shown in
Figure 4. The velocity contours are similar for both the pump positions, but the velocity magnitudes are quite different. The average velocity and maximum velocity (
Table 7) of the salt is higher in the reactor system when the pump is placed after the HEX. Achieving higher velocities inside the reactor system for the same pump momentum source is ideal as it minimizes the risk of formation of recirculation zones leading to better thermal efficiency and performance.
When the pump is placed after the heat exchanger, it copes with the pressure lost during the heat transfer to the secondary side, and if the pressure lost is not too high, then the pump can be placed before the heat exchanger too. The governing metric for the pump position will be flow distribution, average velocity, and maximum velocity of the fuel salt in the system.
3.2. Steady-State Results of the Optimized Design
The steady-state simulation of the 3D LFMSR model was performed in two steps. The core region is 4.1 m in height and has a radius of 3.0 m. The first step is a steady-state case with only thermal hydraulics active. A second steady-state simulation is launched with neutronics and energy equations active and starting from the results of the first step. The spatial velocity and pressure distribution for the 3D reactor system is shown in
Figure 5a,b. All steady-state results shown were gathered after the neutronics and energy solvers converged.
The pressure profile is higher at the bottom because of the salt static pressure and the pump momentum head. The spatial-temperature distribution is shown in
Figure 6a along with the density profile in
Figure 6b. Linear dependence of density on temperature can be observed. As expected, the temperature is higher near the core outlet and the upper wall and lower at the core inlet.
The precursor distributions for delayed neutron groups (families) 1 and 8 are shown in
Figure 7. These results are available because GeN-Foam solves the diffusion equation for the entire primary loop domain. Since the Group 1 precursors have the longest half-life (55.6 s), they do not decay immediately and are convected by the fuel salt from the core into the heat exchanger, then into the pump, and then back into the core again where they accumulate near the outlet in the upper region.
3.3. Unprotected Loss of Flow (ULOF) Transient Accident Simulation
The ULOF accident scenario can occur in an LFMSR when there is a sudden plant outage and the back-up generators fail. This scenario is modeled as a null transient since GeN-Foam does not currently support transient input. To simulate the following transient in GeN-Foam, the pump momentum source is immediately set to zero after the reactor system reaches steady-state conditions as described in the previous section. GeN-Foam iteratively couples the neutronics and thermal-hydraulics calculations. Unlike the steady-state simulations, the transient models do not have a constant applied power. The power distribution is calculated by the neutronics solver and is passed on to the thermal-hydraulics solver which uses the power distribution to calculate the temperature distribution.
The two parameters of interest during the ULOF transient are velocity and temperature. It is crucial that the maximum temperature does not approach the boiling temperature (1673 K) of the fuel salt and that the system attains a non-zero natural circulation velocity. Upon initiating the transient, salt velocity decreases rapidly in the first 25 s, after which it stabilizes due to the establishment of buoyancy-driven flow.
A significant increase in temperature occurs in this transient which decreases the density of the fuel. Fuel-salt temperature is stabilized after the ULOF transient by negative temperature feedback effects. Parameterized cross-sections for 700 K and 1100 K were prepared for these simulations according to the GeN-Foam requirements. The spatial distribution of the velocity and the temperature at the end of the simulated transient (after 300 s) are shown in
Figure 8a,b, respectively. The overall maximum velocity decreases as expected, and its contour changes considerably from the steady-state solution in
Figure 5. During the ULOF transient, the momentum source slowly changes from the pump region to the core region. The predicted velocity magnitudes shown in
Figure 8 seem large for natural circulation and may be due to the large core region that is contributing energy and therefore fluid momentum.
3.4. Uncertainty Propagation and Correlation Analysis
Uncertainty propagation and correlation analyses were performed on the quantities of interest: scalar analysis, one-dimensional functional analysis, and three-dimensional functional analysis. The quantities of interest are represented by the following function:
where the sample number is
and the node number is
Here,
is equal to 100, and M is equal to 511, 568. The timestep variable is represented by
where
is the final timestep. A Python script was written to automate the preprocessing, simulation job, and postprocessing of each case for DAKOTA.
Spatial averages and ensemble averages were applied to the results of uncertainty propagation. Spatial average is defined in Equation (1) as the statistical average over the functional phase-space (space and time). The standard deviation over this functional phase-space is calculated by Equation (2).
A correlation analysis using Spearman rank correlation coefficient (SRCC) is conducted to nonparametrically measure the strength and direction of association between the two ranked variables, i.e., the input and output data. The SRCCs capture the monotonic correlation between any input and output of interest. The input data considered for this study are the parameters through which uncertainties are propagated, and they are the delayed neutron precursor (DNP), delayed precursor decay constant (), liquid fuel density (), liquid fuel dynamic viscosity (), heat transfer coefficient (HTC), pump momentum (PM), specific heat (), and thermal expansion (). The output parameters of interest for this study are temperature and velocity. The correlation analysis study is conducted at the three different levels also investigated for the uncertainty propagation (scalar, one-dimensional, and three-dimensional) to understand how the correlations change with time and space. At the scalar level, the correlations are computed between the input data and the ensemble averaged system average/maximum of the output quantities of interest is computed at the final timestep. The one-dimensional functional correlation analysis focuses on how the correlations between the input and output data evolves with time. This functional study is important as it gives more information about the significance of the parameters before and after the transient.
3.4.1. One-Dimensional Scalar Analysis Results
The two quantities of interest in this scalar analysis are velocity and temperature.
Figure 9 and
Figure 10 show the Empirical Cumulative Distribution Function (ECDF) plots for maximum temperature and average velocity. These plots were generated using the entire data set of the three-dimensional domain. It is clear from
Figure 9 that most of the data for the maximum temperature lie below 1000 K. The temperatures of the system are well below the boiling temperature of the salt, thereby not imposing any design safety risks. The respective standard deviation is 10.73 K. Similarly,
Figure 10 shows that 99% of the data for average velocity lie between 0.4 and 0.5 m/s. These velocities are attained for a natural circulation-driven system that can manage the reactor system during accident scenarios by keeping the salt flowing.
In the scalar sensitivity correlation analysis study, the SRCCs are computed for the ensemble average of system maximum and system average temperature and velocity. The rank correlations for maximum temperature and average velocity are shown in
Figure 11 and
Figure 12, respectively. Kinetic parameters do not have much of an impact on the temperature or velocity. Specific heat, thermal expansion, and heat transfer coefficient have the most important correlations at the end of the transient.
3.4.2. One-Dimensional Functional Analysis Results
The functional sensitivity analysis presented next is like the scalar analysis except the SRCCs are computed at every timestep. The statistical behaviors of temperature and velocity over time are shown in
Figure 13 and
Figure 14, respectively. In
Figure 13, the relative values (dashed line) occur at the same numerical value as the absolute values (solid line).
The time-dependent rank correlations for maximum temperature and average velocity are shown in
Figure 15 and
Figure 16. The thermal expansion coefficient correlation gains more importance as the transient begins. At 150 s, when the pump is turned off completely, the temperature of the system changes, and
determines the magnitude by which the density change occurs as shown in Equation (3).
3.4.3. Three-Dimensional Analysis Results
The evolution of quantities of interest with respect to time is studied as well but not included in this work for brevity. Since not all details are captured through the ECDF plots, this study is extended to every single nodal point in the spatial domain of the reactor system but limited to the last timestep.
To perform a 3D analysis and understand the reactor system behavior spatially, it is important to calculate the ensemble average/statistical mean of the quantities of interest at every nodal point m over different cases at the final timestep (
). The two quantities of interest in this 3D analysis are velocity and temperature. The spatial distribution of the statistical means temperature and velocity in the reactor system is shown in
Figure 17a and
Figure 18a. The standard deviations of temperature and velocity are shown in
Figure 17b and
Figure 18b, respectively.
Recirculation zones do not seem to appear during the transient, and the maximum temperature in the reactor system is well below the boiling temperature of the fuel. The gradients in velocity magnitude near the HEX and pump are due to the sharp porosity changes in those regions. The optimized reactor core design appears to be a good option for minimizing recirculation zones for both steady-state and transient simulation.
The SRCCs are computed in the 3D sensitivity study for the ensemble average of temperature and velocity at every nodal point in space at the end of the transient. The advantage of a 3D SRCC study is that it helps in understanding input parameter correlations with the output quantity of interest at every point in space.
Figure 19,
Figure 20 and
Figure 21 illustrate the 3D results of specific heat capacity, thermal expansion coefficient, and heat transfer coefficient sensitivity with temperature and velocity, respectively.
4. Discussion
The primary loop components have lower significance for the neutron kinetics calculation, and hence, the neutron flux is smaller by several orders of magnitude compared to the core region. The neutron flux is concentrated in the energy groups 2, 3, and 4, with the maximum in energy group 3. The flux concentration depends on the salt composition and the energy boundaries used. Most prompt fission neutrons of U-235, U-238, and Pu-239 isotopes appear in the 0.7 MeV energy range, while the average prompt neutron energy is around 2.2 MeV. Also, for the molten-salt mixture in the energy range given by three energy groups, fewer absorptions occur as a result of a bigger neutron population.
The group 1 precursors can live long enough to move through the whole system many times, and their fraction is low (β1 = 6.186 × 10−5). Their concentration does not depend anymore on where the precursors are generated, and it becomes relatively homogeneous compared to, e.g., the group 8 precursor concentration. The concentration of group 1 is higher at the core’s upper region because of the vicinity of the generation area and the direction of flow. Group 8 will mostly decay in the place where it is generated, and as a result, its concentration is proportional to the neutron flux. GeN-Foam lacks a mathematical steady-state modeling feature, and all the results obtained for the steady simulation are with a time-dependent solver and a system that is not perturbed.
In the reactor system after the pump is turned off, the difference in density is what drives the flow of fuel salt and helps the system attain natural circulation, thereby dictating the temperature distribution. When the flow is driven by pump momentum, the effect of is suppressed as pump momentum is the dominant factor that influences the velocity and temperature of the fuel salt. To understand the significance of specific heat capacity correlations with input parameters, it is important to understand the differences between single side () and total heat transfer ().
The stable nature of these transient results may be due to the simplified geometry and the fact that symmetry boundaries were used on the cell faces that would connect adjacent loops in the actual design. The symmetry boundary condition enforces a zero gradient, so a periodic boundary or a model of two or more adjacent loops may show oscillatory behavior.
The heat gained or lost in the fuel salt depends on the inherent characteristics of the salt such as specific heat; as a result, it would affect the outlet temperature of the primary fuel salt exiting the heat exchanger and in general the overall temperature of the salt. After the transient is initiated, specific heat ceases to be a dominating factor and its effect on temperature is suppressed as soon as density and thermal expansion take over as they start dictating the temperature of the fuel salt in the reactor system. Since specific heat still has an impact on temperature, it is just not very pronounced; this impacts the velocity as density and temperature difference drive the flow upon attaining natural circulation. This suggests specific heat capacity correlation with velocity increases after initiating the transient.
Parameter dependence on temperature decreases upon initiating the transient. As soon as the transient is initiated, it becomes a density-driven flow, and specific heat capacity and the thermal expansion coefficient start dictating the temperature and velocity of the system, whereas the velocity correlation with the heat transfer coefficient increases after initiating the transient. The two behaviors are interlinked: the amount of cooling (single-side heat transfer) achieved in the primary fuel salt depends on the velocity of the system, and the heat transferred from the primary fuel salt to the secondary coolant depends on the heat transfer coefficient. The velocity of the fuel salt determines the residence time of the salt in the heat exchanger, and as a result, the heat transfer coefficient’s correlation with the velocity of the system increases.
Specific heat, thermal expansion, and heat transfer coefficients have the biggest impact on temperature and velocity of the system and only these parameters are taken into consideration in this study. The minimum and maximum temperatures inside the core occur along the centerline and these temperatures are impacted by the specific heat of the molten salt. As a result, there is a strong correlation between temperature and specific heat observed along the centerline of the core (
Figure 17). Since this is a natural circulation-driven flow, the temperature differences inside the core drive the fuel-salt flow and specific heat contributes towards these temperature differences within the core; as a result, a strong correlation between specific heat and velocity of the fuel salt can be observed.
The thermal expansion coefficient has a similar effect on the temperature of the salt as specific heat and is strongly correlated near the centerline of the core. As soon as the pump is deactivated, the system becomes density-driven and the thermal expansion coefficient determines the magnitude by which the density change occurs, thereby impacting the fuel-salt velocity. The amount of cooling (single-side heat transfer) achieved in the primary fuel salt depends on the velocity of the system, and the heat transferred from the primary fuel salt to the secondary coolant depends on the heat transfer coefficient. The velocity of the fuel salt determines the residence time of the salt in the heat exchanger, and as a result, a strong correlation between heat transfer coefficient and velocity occurs.
As the hot fuel salt enters the HEX, the temperature of the fuel salt depends primarily on the heat transfer coefficient input parameter, but as the fuel salt moves through the core, other input parameters such as the thermal expansion coefficient, specific heat, and density come into play, and the uncertainties propagated through these parameters start to add up, thereby leading to higher uncertainties just before the entrance to the HEX. With regards to velocity, the maximum uncertainty is observed near the curved wall of the reactor where the velocities are very close to zero, whereas in other regions, the standard deviation is relatively lower by many orders of magnitude.
As heat is transferred to and from the system, the temperature of the fuel salt changes, and this temperature change depends on two factors: the number of molecules in the system and the heat capacity of the system. Since the volume of the system is fixed, the fuel-salt temperature strongly depends on the specific heat capacity of the salt. The fuel-salt temperature shows the expected negative correlation with specific heat capacity. Upon initiating the ULOF transient, as the temperature of the system increases, the thermal expansion of fluid causes an increase in the density difference between the cold and hot sections which leads to the increase in driving head, thereby increasing fuel-salt velocity. Since it is now a density-driven system, it impacts the maximum temperature of the core as shown in
Figure 18. The SRCC of maximum temperature and average velocity are very similar to the average temperature and maximum velocity, and for brevity, they are not discussed in detail.
5. Conclusions
Through the scalar analysis, 99% of the data for the maximum and average temperature were below 1000 K and 875 K, and the respective standard deviations were 10.73 K and 5.025 K, thereby not imposing any design safety risks. Since a transient scenario was evaluated, it is important to ensure that reasonable temperature and velocity was maintained during the entire duration of the transient. Therefore, this core design seems feasible during the transient evaluated.
One-dimensional functional analysis results revealed that the standard deviations of both temperature and velocity increased right after the transient was initiated. This is exactly when the reactor system was most unstable. As the system attained a natural circulation regime, the standard deviations decreased. The rate of convergence of output parameters can also play a role in the large standard deviations observed as they could vary depending on the input perturbations.
The relative standard deviation of results in 3D were higher near the inlet of the heat exchanger for temperature predictions. This is the result of an accumulation of uncertainties contributed by the input design parameters such as the thermal expansion coefficient and specific heat that take effect as the fuel salt moves through the core.