Simulation of a molten salt fast reactor using the COMSOL multiphysics software

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Simulation of a molten salt fast reactor using the COMSOL multiphysics software

D. H. Daher，

M. Kotb ，

A. M. Khalaf，

Moustafa S. El-Koliel，

Abdelfattah Y. Soliman

Nuclear Science and Techniques

Vol.31, No.12

Article number 115

Published in print 01 Dec 2020

Available online 26 Nov 2020

DOI：10.1007/s41365-020-00833-3

43401

In this study, COMSOL v5.2 Multiphysics software was utilized to perform coupled neutronics and thermal-hydraulics simulations of a molten salt fast reactor, and the SCALE v6.1 code package was utilized to generate the homogenized cross-section data library. The library’s 238 cross-section groups were categorized into nine groups for the simulations in this study. The results of the COMSOL model under no fuel flow conditions were verified using the SCALE v6.1 code results, and the results of the neutronics and thermal-hydraulics simulations were compared to the results of previously published studies. The results indicated that the COMSOL model that includes the cross-section library generated by the SCALE v6.1 code package is suitable for the steady-state analysis and design assessment of molten salt fast reactors. Subsequently, this model was utilized to investigate the neutronics and thermal-hydraulics behaviors of the reactor. Multiple designs were simulated and analyzed in this model, and the results indicated that even if the wall of the core is curved, hotspots occur in the upper and lower portions of the core’s center near the reflectors. A new design was proposed that utilizes a flow rate distribution system, and the simulation results of this design showed that the maximum temperature in the core was approximately 1032 K and no hotspots occurred.

1. Introduction

Molten salt fast reactors (MSFRs) were found to meet the objectives of Generation IV reactors, including sustainability, fuel resource optimization (utilizing thorium instead of uranium enrichment), non-proliferation, safety, and waste management [1]. An initiative named Evaluation and Viability of Liquid Fuel Fast Reactor System (EVOL) was implemented under the European Atomic Energy Community’s Seventh Framework Program [2,3], and during the EVOL project, a benchmark reactor was developed for the pre-conceptual MSFR design. The reactor’s specifications and parameters were studied and analyzed to explore the reactor’s operation and capabilities [4]. This study focuses on the design and performance assessment of MSFRs. First, we validated the model using the benchmark results [3]. Subsequently, we explained the reasons for the difference between our results and the published results with regard to the maximum temperature in the core. Aufiero [5] and Fiorina [6] reported the maximum temperature in the core as 1300 K, whereas Hu [7] and Li [2] reported it as 1200 K and 2100 K, respectively. These results are not valid because the melting point of the nickel alloy utilized for the core structure is approximately 1600 K. Secondly, the optimized geometry of the core of the MSFR design [8,9] was assessed. Finally, a modified design was proposed that utilizes a flow rate distribution system to increase the flow mixing in the core cavity and reduce the temperature peaking factors in the salt and core walls.

2. System description

The MSFR design utilized in this study was a 3000-megawatt thermal reactor with a fuel circuit, intermediate circuit, and transformation system [3]. The fuel circuit (Fig. 1) consists of the salt fuel, which is both the fuel and coolant, key cavity, gas injection system, bubble separators, heat exchangers, tubes, and pumps. Table 1 presents the specifications of the molten salt utilized in the fuel circuit. The salt in the fuel loop is composed of lithium fluoride and actinide fluoride, where actinide fluoride is composed of fissile materials (such as uranium-233[²³³U] and plutonium) and fertile materials (such as thorium). The proportion of actinide fluoride in the fuel salt was fixed at 22.5 mol%, and the total volume of fuel salt in the fuel loop was approximately 18 m³. The fuel salt flows from the bottom to the top of the core cavity. After exiting the core, the fuel salt, which is fed into 16 groups of pumps and heat exchangers located around the core, circulates through the fuel circuit in approximately 4 s [3]. As depicted in Fig. 1, the entire fuel circuit is located inside a reactor vessel, which acts as a second barrier. The three main components of the core are the top and bottom neutron reflectors, the radial fertile blankets (shown in red in Fig. 1), which are part of the radial reflectors, and the fuel reprocessing units. The thick top and bottom reflectors are made of nickel-based alloys and designed to absorb more than 99% of the leaking neutrons. The radial reflectors include a fertile blanket that is approximately 50 cm thick. The blanket is filled with a fertile LiF-ThF₄ salt, initially composed of 22.5 mol% ²³²ThF₄, to increase the MSFR’s breeding ratio. The radial reflectors are enclosed by a 20 cm–thick layer of B₄C, which provides additional neutron protection from the heat exchangers. The fuel circuit includes a salt draining system, which can be used for a planned shutdown or for an emergency shutdown to prevent rapid increases in core temperature. In an emergency, the fuel salt geometry can be passively reconfigured by draining the fuel salt into tanks located under the reactor to produce passive cooling and an acceptable reactivity margin. In this study, two fuel types were utilized at the beginning of the operation, i.e., ²³³U-dependent and transuranic (TRU)-dependent. The initial fuel salt was composed of LiF-ThF₄-(TRU)F₃, which means that the reference MSFR is initiated with a TRU mixture consisting of 87.5% Pu (2.7% ²³⁸Pu, 45.9% ²³⁹Pu, 21.5% ²⁴⁰Pu, 10.7% ²⁴¹Pu, and 6.7% ²⁴²Pu), 6.3% Np, 5.3% Am, and 0.9% Cm. The compositions of the fuels utilized in this study are provided in Table 2, the thermodynamic properties of the fuel and blanket salts are summarized in Table 3, and Fig. 2 shows the benchmark MSFR geometry.

Fig. 1

(Color online) Schematic of the molten salt fast reactor (MSFR) fuel circuit [3]

Fig. 2

(Color online) Benchmark geometry of the MSFR (dimensions in mm) showing the fuel salt (yellow), fertile salt (red), B₄C protection (blue), reflectors, and 20 mm thick walls composed of a Ni-based alloy (green) [3]

Specifications of the molten salt utilized in the fuel circuit [3]

Fuel composition	LiF-[ThF₄-(TRU)F₃]LiF-[ThF₄-²³³UF₄]
Amount of fissile and fertile in molten salt (mole%)	22.5
Volume of fuel salt inside and outside the core (m³)	18
Inlet temperature (K)	923
Outlet temperature (K)	1023
Mean salt temperature (K)	948
Flow rate (m³/s)	4.5
Number of heat exchangers and pump groups	16

Initial composition of the fuels utilized in the ²³³U- and TRU-started MSFRs in the benchmark calculations [3]

²³³U-started MSFR		TRU-started MSFR
Th	²³³U	Th	Actinide
38281 kg	4838 kg	30619 kg	Pu	11079 kg
19.985 mol%	2.515 mol%	16.068 mol%		5.628 mol%
			Np	789 kg
				0.405 mol%
			Am	677 kg
				0.341 mol%
			Cm	116 kg
				0.058 mol%

Fuel salt thermodynamic properties as a function of its temperature (T) in K

Property	Formula	Range (K)
Density, ρ (kg/m³)	4094-0.882 ×(T-1024)	893-1123
Specific heat capacity, C_p (J/K‧kg)	-1111 + 2.78T	867-907
Thermal Conductivity, λ (W/K‧m)	0.928+8.40E-05·T	891-1020
Dynamic viscosity, μ (Pa‧s)	ρ.5.55E-08.exp(3689/T)	898-1119

3. Methodology

The main objectives of this study were to investigate the COMSOL v5.2 Multiphysics (COMSOL) software’s capability of representing and calculating the MSFR’s different parameters by coupling neutronics [10] and thermal-hydraulics [11] models, assess the different MSFR designs, and propose an improved design. COMSOL is a powerful interactive simulation environment used to model and solve various types of scientific and engineering problems; however, to utilize this software in this study, it was necessary to obtain the nuclear data constants, which were calculated using the SCALE v6.1 code package [12]. Additionally, the MSFR geometry was modeled using the Graphically Enhanced Editing Wizard software that is utilized by critical safety specialists who use the SCALE software system developed at the Oak Ridge National Laboratory, which is utilized and accepted worldwide for nuclear safety analyses. The ENDF/B-VII cross-section library was utilized in all the calculations in this study. These codes allow the neutron transport equation to be solved in integral form using multi-group approximation; therefore, it is a valuable tool for the validation process when experimental data are unavailable. Figure 3 illustrates the coupled calculation procedure. The MSFR was modeled utilizing the SCALE v6.1 code package and the T6-DEPL sequence was utilized at a set temperature to calculate the core’s criticality. Then, utilizing different temperatures, the T-DEPL sequence was used to calculate and generate nine groups of cross-sections. From the results, the cross-sections can be described as a function of temperature and used as input data in the neutronics COMSOL model. The statistical error in the neutronics COMSOL model was set to less than 0.001. The MSFR modeled by the COSMOL software was used to validate the neutronics COMSOL model at different temperatures by comparing the criticality with the SCALE v6.1 results. The power density calculated by the neutronics COMSOL model was used as the heat source input in the thermal-hydraulics COMSOL model, which was utilized to produce the temperature and velocity distributions. The temperature and velocity distributions were then inputted into the neutronics COMSOL model to account for the temperature and velocity distribution effects on the different parameters.

Fig. 3

The calculation flowchart

The cross-section library utilized in this study contains microscopic cross-sections from 238 energy groups, which were categorized into nine groups as a benchmark [13]. The 9 groups were classified into three regions, the fast region above 73 keV included 3 groups, the resonance region between 12.9 eV and 73 keV included 5 groups, and the thermal region below 12.9 eV included 1 group. Table 4 presents the energy ranges of the nine groups. Fig. 4 shows the fission cross-sections of ²³³U and TRU. The resonance region of most isotopes is located between 12.9 eV and 73 keV, so it is represented by five groups, while the thermal group (below 12.9 eV) includes only one group.

Fig. 4

(Color online) Neutron fission cross-sections of ²³³U and TRU

Energy ranges of the nine groups [13]

Group number	Neutron energy
	Upper limit	Lower limit
1	20 MeV	1.4 MeV
2	1.4 MeV	573 KeV
3	573 KeV	73 KeV
4	73 KeV	2.29 KeV
5	2.29 KeV	186 eV
6	186 eV	52 eV
7	52 eV	33.25 eV
8	33.25 eV	12.9 eV
9	12.9 eV	1^-5 eV

3.1 Neutronics COMSOL Model

The nine-group diffusion and six-group delayed neutron equations (Equations 1 and 2, respectively) [14] can be represented by the convection and diffusion application model and transport of diluted species model, respectively, as follows:

\begin{array}{l} \frac{1}{ν g} \frac{\partial ϕ^{g}}{\partial t} + \frac{u}{ν g} \cdot \nabla \cdot D^{g} \nabla ϕ^{g} + [\sum_{a}^{g} + \sum_{g' = 1, \neq 9}^{9} \sum_{s}^{g \to g'}] ϕ^{g} \\ = \sum_{g' = 1, \neq 9}^{9} \sum_{s}^{g' \to g} ϕ^{g'} + \frac{1}{k} (1 - β_{eff}) χ_{p}^{9} \sum_{g' = 1}^{9} ν \sum_{f}^{g'} ϕ^{g'} \\ + \sum_{i = 1}^{6} χ_{i}^{g} λ_{i} C_{i}, \end{array}

(1)

\frac{\partial C_{i}}{\partial t} + u . \nabla C_{i} + λ_{i} C_{i} = \frac{1}{k} β_{i} \sum_{g = 1}^{9} ν Σ_{f}^{g} ϕ^{g} .

(2)

In both equations, the left- and right-hand sides define the loss and production mechanisms, respectively. The subscripts a, s, and f denote absorption, scattering, and fission, respectively. The superscript g is the number of neutron energy groups, i is the precursor, $v^{g}$ is the neutron speed of group g (cm/s), φ and is the neutron flux (n/cm²‧s¹), D is the diffusion constant (cm), ∑ is the macroscopic removal cross-section (cm^-1), χ is the fraction of delayed neutrons, $ν$ is the average number of neutrons produced by fission induced by a neutron, λ is the decay constant (s^-1), β is the fraction of delayed neutrons produced, $C_{i}$ is the precursor concentration (cm^-3), and $u$ is the fuel’s velocity vector (cm/s).

When the stationary diffusion model includes the time-dependent coefficient partial differential equation [15], it is represented as follows:

e_{a} \frac{\partial^{2} u}{\partial t^{2}} + d_{a} \frac{\partial u}{\partial t} + \nabla \cdot (- c \nabla u - α u + γ) + β \cdot \nabla u + a u = f,

(3)

where u is the dependent variable, e_a is the mass coefficient, d_a is the damping/mass coefficient, c is the diffusion coefficient, α is the conservative flux convection coefficient, γ is the conservative flux source term, β is the convection coefficient, a is the absorption coefficient, and f is the source term.

The coefficient form of the partial differential equation is used to represent Equation 1, which must be solved. First, the equation must be solved as an eigenvalue problem by assuming that power is only produced via neutron fission reactions and k_eff is equal to 1, which produces Eq. (4) as follows:

λ^{2} e_{a} u^{2} + λ d_{a} u + \nabla \cdot (- c \nabla u - α u + γ) + β \cdot \nabla u + a u = f

(4)

The diffusion equation’s final form is that of Eq. (5) and the eigenvalue solution is utilized as the initial guess to the stationary problem. The convergence rate of stationary problems is strongly dependent on the initial guess; therefore, the solution converges faster if the initial guess is close to the solution.

λ ϕ^{g} - \nabla \cdot D^{g} \nabla ϕ^{g} + [Σ_{a}^{g} + \sum_{g^{'} = 1, \neq g}^{9} Σ_{s}^{g \to g^{'}}] ϕ^{g} = \sum_{g^{'} = 1, \neq g}^{9} Σ_{s}^{g^{'} \to g} ϕ^{g^{'}} + \sum_{g^{'} = 1}^{9} ν Σ_{f}^{g^{'}} ϕ^{g^{'}}

(5)

The second step is to solve Eqs. (1) and (2) utilizing the initial guess generated by the eigenvalue solution under the condition that the reactor’s power is equal to 3 GW utilizing the global equations in the module.

The same symmetry and vacuum boundary condition can be represented by Eqs. (6) and (7), respectively, as follows:

- n \cdot (- D \nabla ϕ) = 0

(6)

- n \cdot (- D \nabla ϕ) = - \frac{1}{2.1312} ϕ,

(7)

where D is the diffusion coefficient.

The symmetry and outer boundary conditions of the fuel and fertile circuits can be represented by Eqs. (8) and (9), respectively, and utilized to solve the precursor's concentration equation (Equation 2).

- n \cdot (- \nabla C_{i}) = 0

(8)

C_{i} = 0

(9)

The COMSOL software utilizes the finite element method to introduce test functions that are defined through a computational mesh. For each computational cell or mesh element, several test functions are locally defined. Additionally, as part of the finite element method, shape functions are defined that are used to represent the candidate solutions.

TRITON software utilizes the T-DEPL module to generate the macroscopic neutron cross-section, the energy released per fission, delayed neutron fraction, decay constants, and delayed neutron production of the precursor groups for the different regions (i.e., fuel, blanket, and absorber and structure materials). The nuclear constants required to solve the neutron multi-group diffusion equation in the COMSOL model were generated utilizing the ENDIF/B-VII library in the T-DEPL module. In this study, the 238-energy group library was adopted and the 238-group macroscopic cross-sections were categorized into nine energy groups. The group constants generated at different temperatures must be interpreted as a function of temperature. The criticality calculated by the T6-DEPL module was used to validate the neutronics COMSOL model at different temperatures for both fuel types and the nine energy groups were chosen from the benchmark.

3.2 Thermal-hydraulic COMSOL model

Turbulent flow (k-ε model) and heat transfer (fluids model) models were utilized in the COMSOL software to study the thermal-hydraulics of MSFRs, which are based on the fundamental laws of the mass, momentum, and energy conservation equations used in the computational fluid dynamics method. The laminar and Euler models were initially incorporated and then extended to account for turbulence effects using the k-ε model. The propensity for an isothermal flow to become turbulent is measured by the Reynolds number as follows:

R_{e} = \frac{ρ U L}{μ},

(10)

where μ is the dynamic viscosity (Pa‧s), ρ is the density (kg/m³), U is the velocity(m/s), and L is the scale of the flow’s length (m).

Assuming that the fluid is incompressible and Newtonian, the Navier-Stokes equations [15,16] become Eqs. (11) and 12 as follows:

ρ \frac{\partial u}{\partial t} + ρ (u \cdot \nabla) u = \nabla \cdot - p I + μ (\nabla u + {(\nabla u)}^{T}) + F,

(11)

ρ \nabla \cdot u = 0,

(12)

where ρ is the density (kg/m³), u is the velocity vector (m/s), P is the pressure (Pa), and F is the volume force vector (N/m³).

The k-ε model is one of the most widely used turbulence models for industrial applications. The T-DEPL module utilized in this study includes the standard k–ε model, which introduces two additional transport equations (Eqs. (7) and (8)) and two dependent variables, turbulent kinetic energy, k, and turbulent dissipation rate, ε. The turbulent viscosity is modeled as follows:

μ_{T} = ρ C_{μ} \frac{k^{2}}{ε},

(13)

where C_μ is a model constant.

The transport equations for k and e (Eqs. (14) and (15), respectively) are as follows:

ρ \frac{\partial k}{\partial t} + ρ u \cdot \nabla k = \nabla \cdot ((μ + \frac{μ_{T}}{σ_{k}}) \nabla k) + P_{k} - ρ ε,

(14)

ρ \frac{\partial ε}{\partial t} + ρ u \cdot \nabla ε = \nabla \cdot ((μ + \frac{μ_{T}}{σ_{ε}}) \nabla ε) + C_{ε 1} \frac{ε}{k} P_{k} - C_{ε 2} ρ \frac{ε^{2}}{k},

(15)

where the production term is calculated by Eq. (16) as follows:

P_{k} = μ_{T} (\nabla u : (\nabla u + {(\nabla u)}^{T}) - \frac{2}{3} {(\nabla \cdot u)}^{2}) - \frac{2}{3} ρ k \nabla . u .

(16)

The heat transfer in the solid and liquid interfaces is obtained by solving Eqs. (17) and (18) as follows:

ρ C_{p} (\frac{\partial T}{\partial t} + u_{trans} \cdot \nabla T) + \nabla \cdot (q + q_{r}) = - α_{p} T ∶ \frac{\partial S}{\partial t} + Q

(17)

ρ C_{p} (\frac{\partial T}{\partial t} + u \cdot \nabla T) + \nabla \cdot (q + q_{r}) = α_{p} T (\frac{\partial p}{\partial t} + u \cdot \nabla p) + τ ∶ \nabla u + Q,

(18)

where C_p is the specific heat capacity at constant stress (J/kg‧K), T is the absolute temperature (K), u_trans is the velocity vector of translational motion (m/s), α is the coefficient of thermal expansion (1/K), S is the second Piola-Kirchhoff stress tensor (Pa), τ is the viscous stress tensor (Pa), Q contains additional heat sources (W/m³), and q and q_r are the heat flux by conduction and radiation (W/m²), respectively.

4. Results

The MSFR geometry was modeled using COMSOL software and the SCALE v6.1 code package (Fig. 5); however, due to its symmetry, only one-sixteenth of the core was modeled. To validate the neutronics model created by the COMSOL software, the results from the neutronics COMSOL model were compared with the SCALE v6.1 results. In the neutronics COMSOL model, the number of meshes is equal to 17615. In the thermal-hydraulics COMSOL model, three different numbers of meshes in the core were utilized to study the effect of mesh size on the results due to discrepancies in the previously published study results [2,5–7]. The number of meshes in the core utilized was 8665, 25804, and 217105 in Cases I, II, and III, respectively, which is shown in Fig. 6; however, to reduce the computation time, a small number of meshes were utilized in the solid regions.

Fig. 5

(Color online) 3-D models of the fuel circuit (a) COMSOL v5.2 Multiphysics and (b) SCALE v6.1. The red, blue, yellow, and green areas represent the fuel, fertile blanket, absorber, and structure, respectively.

Fig. 6

Number of meshes in the core utilized in the thermal-hydraulics COMSOL model: (a) Case I = 865, (b) Case II = 25804, and (c) Case III = 217105

Table 5 shows the effective multiplication factors, k_eff, calculated utilizing the SCALE v6.1 and COMSOL software under the condition of no flow at different temperatures for the ²³³U- and TRU-started MSFRs.

k_eff calculated by SCALE v6.1 and COMSOL

Temperature (°K)	²³³U-started MSFR		TRU-started MSFR
	SCALE v6.1	COMSOL	SCALE v6.1	COMSOL
923	0.9904 ± 0.0011	0.9871754±0.00049	1.0254 ± 0.0010	1.02695989±0.00041
948	0.9868 ± 0.0013	0.9848586±0.00032	1.0241 ± 0.0012	1.02579067±0.00063
1000	0.9847 ± 0.0012	0.9810806±0.00064	1.0237 ± 0.0015	1.02343186±2.6e-06

For the ²³³U-started MSFR, the differences between the SCALE v6.1 and COMSOL results were approximately 322, 194, and 361 pcm at 923, 948, and 1000 K, respectively. These differences are attributed to the different calculation methods utilized by the software and because the nuclear data were categorized into nine groups in the COMSOL model.

For the TRU-started MSFR, the differences between the SCALE v6.1 and COMSOL results were approximately 156, 169, and 27 pcm at 923, 948, and 1000 K, respectively. These differences are attributed to the same reasons as the differences in the ²³³U-started MSFR results. The results were in good agreement, which validated the neutronics COMSOL model; therefore, it was utilized in this study to represent the MSFR.

The neutronics COMSOL model calculates the local power densities, neutron flux, precursor concentration, and reactivity. In TRU- and ²³³U-started MSFRs, the maximum power density is located at the core’s center, as shown in Fig. 7. The power density at the core’s center in a TRU-started MSFR is greater than in a ²³³U-started MSFR, which implies that the fast neutron flux is greater in a TRU-started MSFR compared to that in a ²³³U-started MSFR. Additionally, the total fission energy release from a fissile element in TRU is greater than that in ²³³U, except for ²³⁹Pu and ²³⁸Pu, as shown in Fig. 8.

Fig. 7

(Color online) Power densities (W/m³) of (a) ²³³U-started and (b) TRU-started MSFR

Fig. 8

Total energy released from fission

The power densities generated by the neutronics COMSOL model were utilized as the heat source input values in the thermal-hydraulics COMSOL model, which was used to predict the temperature and velocity distributions. The temperature and velocity distributions in the fuel circuits of the ²³³U- and TRU-started MSFRs calculated using the thermal-hydraulics COMSOL model at a fuel flow rate of 4.5 m³/s for Cases I, II, and III are presented in Figs. 9, 10, and 11, respectively. In all cases, the maximum temperature was located near the core wall where the flow became vortex, as can be seen in the figures. In the ²³³U-started MSFR, the maximum temperatures were 1460.37, 1612.37, and 2223.95 K for Cases I, II, and III (which had 8655, 25804, and 217105 numbers of meshes), respectively. It is evident from the results that the maximum temperature increased with an increase in the number of meshes in the core. The same conclusion was determined for Cases I, II, and III for the TRU-started MSFR, where the maximum temperatures were 1448.22, 1612.45, and 2224.20 K, respectively. The difference in the number of meshes explains the different maximum temperature results. The maximum temperature was published by Aufiero [5] and Fiorina [6] of 1300 K and Hu’s result [7] of 1200 K. These results are relatively consistent with the Case I results in this study. Our results, however, were closer to the results of Li’s study [2], which determined a maximum temperature of 2100 K. Therefore, the discrepancy was determined to be caused by the difference in the number of meshes utilized in the model, as revealed by Cases I, II, and III, because the heat transfer analysis depends strongly on the flow field, particularly when vortex flow is present. The velocity at the center of the vortex is approximately zero and the heat transfer from convection is minimal; therefore, only heat transfer from conduction occurs. Consequently, when the element size is decreased, the heat transfer from conduction decreases, and the heat is then generated from fission. This flow is more accurately described when using a smaller element size (greater number of meshes) [17]; therefore, to obtain accurate results, the element size must be very small. As such, these results were unacceptable; therefore, the reference core geometry must be optimized to prevent vortex regions from occurring, which cause high temperatures.

Fig. 9

(Color online) Temperature and velocity distributions for Case I (number of mesh = 8665). (a) ²³³U-started MSFR temperature distribution, (b) TRU-started MSFR temperature distribution, (c) ²³³U-started MSFR velocity distribution, (d) TRU-started MSFR velocity distribution.

Fig. 10

(Color online) Temperature and velocity distributions for Case II (number of mesh = 25804). (a) ²³³U-started MSFR temperature distribution, (b) TRU-started MSFR temperature distribution, (c) ²³³U-started MSFR velocity distribution, (d) TRU-started MSFR velocity distribution.

Fig. 11

(Color online) Temperature and velocity distributions for Case III (number of mesh = 217105). (a) ²³³U-started MSFR temperature distribution, (b) TRU-started MSFR temperature distribution, (c) ²³³U-started MSFR velocity distribution, (d) TRU-started MSFR velocity distribution.

As expected, the outlet temperature was approximately 100 K higher than the inlet temperature, which corresponds with the benchmark and the temperature predicted utilizing Equation 11. The maximum temperature in the core, however, was approximately 1300.95 K higher than the inlet temperature calculated utilizing Equation 11 as follows:

d T = \frac{P}{ρ \times C_{p} \times flow rate},

(19)

where P is the reactor power (3 GW), ρ is the fuel salt density (4124.87 kg/m³ at T=973 K), C_p is the specific heat of the fuel salt (1593.94 J/kg‧K at T=973 K), and the flow rate was 4.5 m³/s.

The temperature and velocity distributions determined by the thermal-hydraulics COMSOL model were used as input data in the neutronics COMSOL model to study the effect of fuel flow on temperature. The macroscopic cross-sections are highly temperature-dependent and the emission location of the delayed neutrons is affected by the fluid's motion due to their half-life time. In contrast, prompt neutrons have a very small half-life. Table 6 compares the k_eff values calculated via the COMSOL software for the ²³³U- and TRU-started MSFRs utilizing a flow of 4.5 m³/s for the different numbers of meshes in the core to those of the EVOL benchmark [3]. For the lowest mesh number (8665), the ²³³U- and TRU-started COMSOL results were similar to the EVOL results. The ²³³U-started differences were 946, 1137, and 1807 pcm when the number of meshes was 8665, 25804, and 217105, respectively. For the TRU-started results, the differences were 55, 194, and 642 pcm when the number of meshes was 8665, 25804, and 217105, respectively. As previously discussed in Section 3, the temperature increases where the k_eff values decrease with the largest number of meshes because the cross-sections are highly temperature-dependent. Therefore, to obtain more accurate results, the number of meshes must be increased. Additionally, the reactor’s reactivity decreases based on the fuel flow, because large amounts of delayed neutrons are emitted from the core. The maximum losses were approximately 1614 and 1030 pcm for the ²³³U- and TRU-started MSFRs, respectively, which are in agreement with the EVOL results for a low number of meshes.

Comparison of k_eff values for the COMSOL models and EVOL benchmark

²³³U-started MSFR			TRU-started MSFR
EVOL	COMSOL		EVOL	COMSOL
	No. of mesh	k_eff		No. of mesh	k_eff
0.98301	8665	0.97355	1.01955	8665	1.01900
	25804	0.97164		25804	1.01761
	217105	0.96494		217105	1.01313

The prompt neutron distributions are presented in Figs. 12, 13, 14, 15, 16, 17, 18, 19, 20.

Fig. 12

(Color online) Neutron flux for the first energy group. (a) ²³³U-started MSFR; (b) TRU-started MSFR.

Fig. 13

(Color online) Neutron flux for the second energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 14

(Color online) Neutron flux for the third energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 15

(Color online) Neutron flux for the fourth energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 16

(Color online) Neutron flux for the fifth energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 17

(Color online) Neutron flux for the sixth energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 18

(Color online) Neutron flux for the seventh energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 19

(Color online) Neutron flux for the eighth energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

Fig. 20

(Color online) Neutron flux for the ninth energy group. (a) ²³³U-started MSFR, (b) TRU-started MSFR.

For the ²³³U-started MSFR, the maximum fluxes of the first, second, third, fourth, and fifth neutron energy groups occurred at the core’s center with values of 7.98383×10¹⁴, 1.05851×10¹⁵, 2.06449×10¹⁵, 5.69202×10¹⁵, and 1.05053×10¹⁵ neutrons/(cm²‧s), respectively. The maximum fluxes of the sixth, seventh, and eighth neutron energy groups occurred at the top and bottom of the core with values of 9.06652×10¹³, 1.4110×10¹³, and 2.78732×10¹³ neutrons/(cm²‧s), respectively. The maximum flux of the ninth neutron energy group occurred at the fertile material with a value of 1.35349×10¹³ neutrons/(cm²‧s).

For the TRU-started MSFR, the maximum fluxes of the first, second, third, fourth, and fifth neutron energy groups occurred at the core’s center with values of 8.63153×10¹⁴, 1.11209×10¹⁵, 2.13570×10¹⁵, 5.57955×10¹⁵, and 7.85002×10¹⁴ neutrons/(cm²‧s), respectively. The maximum fluxes of the sixth, seventh, and eighth neutron energy groups occurred at the top and bottom of the core with values of 6.45237×10¹³, 9.79504×10¹³, and 1.90963×10¹³ neutrons/(cm²‧s), respectively. The maximum flux of the ninth neutron energy group occurred at the fertile material with a value of 1.14321×10¹³ neutrons/(cm²‧s).

From these results, it is evident that for the ninth neutron energy group, which has neutron energy in the thermal energy range of 12.9 to 10^-5 eV, the maximum neutron flux occurs inside the fertile blanket. The flux in the fertile blanket was more thermal than in the core region because the fast fission neutrons generated in the core thermalized as they moved towards the outer regions, becoming maximum in the fertile blanket, and a high thermal flux was observed in the fertile blanket and structure’s material.

Figure 21 shows the total neutron fluxes in the core, blanket, absorber, and structure regions of the ²³³U- and TRU-started MSFRs. For the ²³³U-started MSFR, the average flux in the core was 4.03697×10¹⁵ neutrons/(cm²‧s). The maximum flux in the core was 1.07494×10¹⁶ neutrons/(cm²‧s), which steadily decreased towards the outer region of the core until it reached 1.67329×10¹⁵ neutrons/(cm²‧s) at the fuel core/structural material boundary. The total flux in the blanket was 1.24907×10¹⁵ neutrons/(cm²‧s) and decreased to 2.90720×10¹³ neutrons/(cm²‧s) at the fertile blanket/second structural material boundary.

Fig. 21

(Color online) (a)²³³U-started MSFR, (b) TRU-started MSFR, (c) Total neutron flux for all energy groups

The same behavior was observed for the TRU-started MSFRs. The average flux in the core was 3.62538×10¹⁵ neutrons/(cm²‧s) and the maximum flux in the core was 1.04980×10¹⁶ neutrons/(cm²‧s), which steadily decreased towards the outer region of the core until it reached 1.65252×10¹⁵ neutrons/(cm²‧s) at the fuel core/structural material boundary. The total flux in the blanket was 1.25863×10¹⁵ neutrons/(cm²‧s) and decreased to 2.99966×10¹³ neutrons/(cm²‧s) at the fertile blanket/second structural material boundary.

²³³U-started MSFR: In the core region, it is evident that the largest contribution in total flux comes from the fourth neutron energy group (neutron energy between 2.29 and 73 KeV), which contributes approximately 53% of the total flux. The third neutron energy group has the second-largest contribution of 19%, followed by the fifth, first, and second neutron energy groups with contributions of 10%, 7.3%, and 9.5%, respectively. The contributions of the sixth, seventh, eighth, and ninth groups are less than 0.8%.

In the blanket region, the largest flux contribution in total flux comes from the fourth neutron energy group, contributing 42% to 53% of the total flux. The fifth neutron energy group contributes the second-largest amount ranging from 20% to 31%, followed by the sixth, third, and ninth neutron energy groups with contributions up to 23%, 15%, and 11%, respectively. The maximum contribution of the remaining groups was 11%.

TRU-started MSFR: In the core region, it is evident that the largest contribution in total flux comes from the fourth neutron energy group, which contributes approximately 53% of the total flux. The third neutron energy group has the second-largest contribution of 20%, followed by the second, first, and fifth neutron energy groups with contributions of 11%, 8%, and 7%, respectively. The contributions of the sixth, seventh, eighth, and ninth groups are less than 0.25%.

In the blanket region, the largest contribution comes from the fourth neutron energy group, contributing 40% to 53% of the total flux. The fifth neutron energy group contributes the second-largest amount ranging from 18% to 31%, followed by the third, ninth, and sixth neutron energy groups with contributions up to 14%, 8.8%, and 8%, respectively. The maximum contribution of the remaining groups was 11%.

The presence of light elements in the salt, such as lithium and fluorine, produced a softer neutron spectrum in the MSFR core compared with other solid-fuel fast reactors, which is the reason the highest neutron flux values were from the epithermal energy groups.

The precursor concentration distributions for the six groups of the ²³³U-started MSFR are presented in Fig. 22. The maximum concentrations of the first through sixth groups were 6.71018×10¹⁰, 5.76196×10¹⁰, 2.08157×10¹⁰, 1.53499×10¹⁰, 1.87656×10⁹, and 4.7346×10⁸ cm^-3, respectively. The maximum precursor concentrations of the first through fourth groups occurred near the core’s wall and the maximum precursor concentrations of the fifth and sixth groups occurred near the top of the core’s center, which was attributed to the fuel’s flow becoming vortex near the wall (Figs. 22c and d). This phenomenon was more apparent in precursor groups with small decay constants. When the salt's circulation time was 4 s, the half-life of the precursors of delayed neutrons with a small decay constant was approximately 55.46 s. For delayed neutron precursors with a short half-life of 0.22 s, this effect was almost nonexistent. For all groups, the maximum precursor concentration decreased with the decay constant.

Fig. 22

(Color online) Precursor concentrations (precursors/cm³) for the ²³³U-started MSFR. (a) First group, (b) Second group, (c) Third group, (d) Fourth group, (e) Fifth group, (f) Sixth group

The precursor concentration distributions for the six groups of the TRU-started MSFR are presented in Fig. 23. The maximum concentrations of the first through sixth groups were 4.08564×10¹⁰, 1.24218×10¹¹, 2.77181×10¹⁰, 2.25090×10¹⁰, 2.88048×10⁹, and 5.4025×10⁸ cm^-3, respectively. The maximum precursor concentrations of the first through fifth groups occurred near the core’s wall and the maximum precursor concentration of the sixth group occurred near the top of the core’s center, which was attributed to the fuel’s flow becoming vortex flow near the wall. The half-life values of the six precursor groups generated via SCALE v6.1 are presented in Table 7 for the ²³³U- and TRU-started MSFRs. These results show that the delayed neutron’s half-life for the TRU-started MSFR is greater than that of the ²³³U-started MSFR; therefore, its effect differs based on the fuel flow. Notably, the prompt neutrons were not affected by the fuel-salt flow because they were created via fission within 10⁻¹³ s after the fission event [18]; however, the distribution of precursors, and consequently the delayed neutrons, were affected by the fuel flow. This is a unique feature of the MSFR, wherein some of the precursors exit the core, leading to a decrease in the reactor’s criticality.

Fig. 23

(Color online) Precursor concentrations (precursors/cm³) for the TRU-started MSFR. (a) First group, (b) Second group, (c) Third group, (d) Fourth group, (e) Fifth group, (f) Sixth group.

Decay constants and Beta fractions of the six precursor groups

Precursor group	TRU-started MSFR		²³³U-started MSFR
	Decay constant,λ (1/s)	Fraction(β)	Decay constant,λ (1/s)	Fraction(β)
1	9.65E-03	7.51E-05	1.25E-02	2.91E-04
2	2.33E-02	5.23E-04	3.59E-02	6.66E-04
3	1.00E-01	4.17E-04	1.38E-01	7.66E-04
4	2.48E-01	7.30E-04	3.18E-01	1.14E-03
5	9.48E-01	2.46E-04	1.22E+00	3.27E-04
6	2.45E+00	7.46E-05	3.15E+00	1.37E-04

Based on these results, it was concluded that this model is suitable for studying MSFRs; however, the model’s design must utilize a large number of meshes and prevent the occurrence of vortex regions. Therefore, a thermal-hydraulic model was utilized [8,9]. The model’s geometry, which is presented in Fig. 24a, included curved radial core walls (in contact with the fertile blanket) to prevent the vortex region that occurred in the benchmark geometry from forming.

Fig. 24

(a) Geometry of the MSFR and (b) computational fluid dynamic mesh utilized in the model

To obtain a balance between calculation efficiency and accuracy, the maximum bulk cell size was set at 3 cm. The mesh utilized in the model is illustrated in Fig. 24b.

The fuel salt temperature and velocity distributions for the model are presented in Fig. 25. The maximum velocities occurred near the lower and upper reflectors due to the curvature of the inlet and outlet, respectively. The minimum velocities occurred at the top and bottom of the core’s center and near the wall, which are also where the hotspots occurred with temperatures of approximately 1225, 1105, and 1154 K, respectively.

Fig. 25

(Color online) (a) Temperature distribution and (b) velocity distribution

From these results, it was concluded that this design does not prevent the occurrence of hotspots. Consequently, two modified designs, Design I and II, were studied in which the upper and lower walls were curved. These designs are presented in Fig. 26.

Fig. 26

Geometries of the two modified MSFRs: Design (a) I and (b) II

The temperature distributions for the Design I and II models are presented in Fig. 27. It is evident that hotspots still occur on the top and bottom of the core’s center and the wall between the core and blanket in both designs. The highest temperatures in Design I were approximately 1102 and 1148 K and in Design II were approximately 1091 and 1157 K, which occurred on the top of the core’s center and the wall between the core and blanket, respectively, in both designs.

Fig. 27

(Color online) Temperature distributions for Design (a) I and (b) II

Based on these results, hotspots still occur even when both the upper and lower walls of the core are curved. Therefore, a new design was determined that can distribute the flow proportionately to the energy generated. Fig. 28 shows the proposed design, which includes 12 channels with different cross-sectional areas at the inlet.

Fig. 28

Geometry of the proposed MSFR

To perform the grid sensitivity analysis, we used different element sizes in the mesh. Table 8 shows the maximum temperature results for each element size and the differences between the maximum temperatures do not exceed 4.2 degrees. The fuel salt temperature and velocity distributions for the proposed model are presented in Fig. 29. The maximum and minimum velocities occurred near the curvatures of the 12 channels near the wall between the core and blanket, respectively. The highest temperature was approximately 1032 K, while the temperature out of the core is equal to 1023. So that the highest temperature is eight degrees higher than the outlet temperature. Based on these results, the proposed design prevents hotspots from occurring in the core.

Fig. 29

(Color online) (a) Temperature and (b) velocity distributions of the proposed MSFR design

Maximum temperatures utilizing different element sizes

Element size(cm)	Maximum Temperature(K)
3	1031.840
4	1027.650
5	1029.223
6	1027.744
7	1028.949
8	1028.956

5. Conclusion

The steady-state of the MSFR was investigated based on the neutronics benchmark of the MSFR developed by the EVOL project. We used the SCALE v6.1 code package to generate the nuclear constants required to perform this study and validate the neutronics COMSOL model. We performed a coupled neutronics and thermal-hydraulics analysis of the MSFR core utilizing COMSOL v5.2 Multiphysics software and the use of turbulence models was determined more suitable for this simulation. The main objective of this study was to validate the reactor core model created by COMSOL. First, the model was validated under the condition of no flow by comparing the k_eff results with the SCALE v6.1 results, which were found to be in good agreement. Second, the model was validated using EVOL data and the k_eff results, which determined that the shape of the fluxes in the model and reference were in good agreement. The discrepancies in the temperature distributions in previously published studies were caused by the difference in the number of meshes utilized in the core. The heat transfer and computational fluid dynamics analysis results depend strongly on fluid flow, especially when vortex flow occurs. Vortex flow is more accurately described using a smaller element size (increasing the number of meshes); therefore, the reference core geometry must be modified to prevent the formation of vortex flow, which causes hotspots. The maximum total flux values in the ²³³U- and TRU-started MSFRs were 1.07494×10¹⁶ and 1.04980×10¹⁶ neutrons/(cm²‧s), respectively, which occurred at the core’s center in both MSFRs. The largest total flux contribution comes from the neutron energy groups with energy ranging between 2.29 and 73 KeV (the epithermal energy groups). The presence of light elements in the salt, such as lithium and fluorine, produced a softer neutron spectrum in the MSFR core compared to other solid-fuel fast reactors, which is the reason the highest neutron flux contributions were from the epithermal energy groups. The maximum thermal flux occurred in the fertile region and the maximum precursor concentrations of the six groups occurred near the wall of the core and at the top of the core’s center. This is attributed to the fuel’s flow becoming vortex flow near the wall and this phenomenon was more prominent in precursor groups that have a small decay constant. The computational fluid dynamics and heat transfer results showed that the maximum temperatures occur near the walls, which is caused by the fuel’s flow becoming vortex near the wall. From the results, it was determined the original reference geometry needed to be optimized to prevent the formation of vortex flow regions, which leads to the occurrence of hotspots. Finally, it was determined that the coupled neutronics and thermal-hydraulics model created by the COMSOL software was suitable to model MSFRs; however, the modification of the model’s geometry to utilize a large number of meshes and prevent the occurrence of vortex regions was required. By studying multiple designs, it was concluded that the optimum design must have a flow distribution system to prevent the occurrence of hotspots and the final design proposed in this study was able to prevent the occurrence of hotspots in the core.

References

[1]

D. Buckthorp.

Introduction to Generation IV nuclear reactors

. Structural Materials for Generation IV Nuclear Reactors. Elsevier, 2017. p. 1-22. https://doi.org/10.1016/B978-0-08-100906-2.00001-X.

Baidu Scholar

Google Scholar

[2]

R. Li, S. Wang, A. Rineiski, et al.

Transient analyses for a molten salt fast reactor with optimized core geometry

. Nucl. Eng. Design 292:164-176 (2015). https://doi.org/10.1016/j.nucengdes.2015.06.011

Baidu Scholar

Google Scholar

[3]

M. Brovchenko, J.-L. Kloosterman, L. Luzzi, et al.

Neutronic benchmark of the molten salt fast reactor in the frame of the EVOL and MARS collaborative projects

. EPJ Nucl. Sci. Technol. 5:2 (2019). https://doi.org/10.1051/epjn/2018052

Baidu Scholar

Google Scholar

[4]

E. Merle-Lucotte, M. Allibert, M. Brovchenko, et al.

Introduction to the Physics of Thorium Molten Salt Fast Reactor (MSFR) Concepts

. Thorium Energy for the World, Cham: Springer International Publishing; 2016. p. 223-231. https://doi.org/10.1007/978-3-319-26542-1_34

Baidu Scholar

Google Scholar

[5]

M. Aufiero, A. Cammi, O. Geoffroy, et al.

Development of an OpenFOAM model for the Molten Salt Fast Reactor transient analysis

. Chem.Engineering Sci. 111:390-401 (2014). https://doi.org/10.1016/j.ces.2014.03.003

Baidu Scholar

Google Scholar

[6]

C. Fiorina, D. Lathouwers, M. Aufiero, et al.

Modelling and analysis of the MSFR transient behaviour

. Ann. Nucl. Energy 64: 485-498 (2014). https://doi.org/10.1016/j.anucene.2013.08.003

Baidu Scholar

Google Scholar

[7]

T. Hu, L. Cao, H. Wu, et al.

Coupled neutronics and thermal-hydraulics simulation of molten salt reactors based on OpenMC/TANSY

. Ann. Nucl. Energy 109:260-276 (2017). https://doi.org/10.1016/j.anucene.2017.05.002

Baidu Scholar

Google Scholar

[8]

A. Laureau, M. Aufiero, P.R. Rubiolo, et al. Coupled neutronics and thermal-hydraulics transient calculations based on a fission matrix approach: Application to the Molten Salt Fast Reactor. 2015.

[9]

A. Laureau, P.R. Rubiolo, D. Heuer, et al.

Coupled neutronics and thermal-hydraulics numerical simulations of a Molten Fast Salt Reactor (MFSR)

. Joint International Conference on Supercomputing in Nuclear Applications + Monte Carlo 2014. https://doi.org/10.1051/snamc/201402307

Baidu Scholar

Google Scholar

[10]

D. Chandler, J. Freels, G. Maldonado, et al. COMSOL-based Nuclear Reactor Kinetics Studies at the HFIR. 2011.

[11]

B.R. Betzler, Powers JJ, Worrall A.

Molten salt reactor neutronics and fuel cycle modeling and simulation with SCALE

. Ann. Nucl. Energy 101:489-503 (2017). https://doi.org/10.1016/j.anucene.2016.11.040

Baidu Scholar

Google Scholar

[12]

S.M. Bowman.

SCALE 6: Comprehensive Nuclear Safety Analysis Code System

. Nucl. Technol. 174:126-148 (2011). https://doi.org/10.13182/NT10-163

Baidu Scholar

Google Scholar

[13]

M. Brovchenko, E. Merle. Optimization of the pre-conceptual design of the MSFR - EVOL project. 2013;2.

[14]

M.A. Schultz.

Dynamics of Nuclear Reactors

.Nucl. Technol.12:415 (1971). https://doi.org/10.13182/NT71-A30993

Baidu Scholar

Google Scholar

[15]

COMSOL. Introduction to COMSOL Multiphysics 5.3. 2014;

[16]

E.L. Paul, V.A. Atiemo-Obeng, S.M. Krest. Handbook of Industrial Mixing: Science and Practice - Wiley Online Library. John Wiley Sons. 2004.

[17]

S. Mulligan. Experimental and numerical analysis of three-dimensional free-surface turbulent vortex flows with strong circulation. 2015.

[18]

US Department of Energy.

DOE FUNDAMENTALS HANDBOOK NUCLEAR PHYSICS 2, 1993

. https://doi.org/DOE-HDBK-1019/2-93

Baidu Scholar

Google Scholar