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

Simulation of a molten salt fast reactor using the COMSOL multiphysics software 1 2 first-author https://orcid.org/0000-0002-9525-4008 1 corresp Corresponding author, mahmoudkotb@azhar.edu.eg Corresponding author, mahmoudkotb@azhar.edu.eg mahmoudkotb@azhar.edu.eg https://orcid.org/0000-0001-7261-1343 1 2 3 https://orcid.org/0000-0002-4831-9251 Department of Physics, Faculty of Science, Al-Azhar University, Cairo, Egypt Department of Physics, Faculty of Science, Al-Azhar University, Cairo, Egypt Department of Reactors, Nuclear Research Center, Atomic Energy Authority, Cairo, Egypt Department of Reactors, Nuclear Research Center, Atomic Energy Authority, Cairo, Egypt Nuclear Engineering Department, King Abdulaziz University, Saudi Arabia Nuclear Engineering Department, King Abdulaziz University, Saudi Arabia 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 1 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 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 2 , 3 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 4 ]. This study focuses on the design and performance assessment of MSFRs. First, we validated the model using the benchmark results [ 3 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 5 ] and Fiorina [ 6 6 ] reported the maximum temperature in the core as 1300 K, whereas Hu [ 7 7 ] and Li [ 2 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. S econdly, the optimized geometry of the core of the MSFR design [ 8 8 , 9 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 1 The MSFR design utilized in this study was a 3000-megawatt thermal reactor with a fuel circuit, intermediate circuit, and transformation system [ 3 3 ]. The fuel circuit ( Fig. 1 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 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[ 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 3 . 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 3 ]. As depicted in Fig. 1 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 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 4 salt, initially composed of 22.5 mol% 232 ThF 4 , to increase the MSFR’s breeding ratio. The radial reflectors are enclosed by a 20 cm–thick layer of B 4 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., 233 U-dependent and transuranic (TRU)-dependent. The initial fuel salt was composed of LiF-ThF 4 -(TRU)F 3 , which means that the reference MSFR is initiated with a TRU mixture consisting of 87.5% Pu (2.7% 238 Pu, 45.9% 239 Pu, 21.5% 240 Pu, 10.7% 241 Pu, and 6.7% 242 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 Table 2 , the thermodynamic properties of the fuel and blanket salts are summarized in Table 3 Table 3 , and Fig. 2 Fig. 2 shows the benchmark MSFR geometry. 1 1 Fuel composition 1 1 LiF-[ThF 4 -(TRU)F 3 ]LiF-[ThF 4 - 233 UF 4 ] 1 1 Amount of fissile and fertile in molten salt (mole%) 1 1 22.5 1 1 Volume of fuel salt inside and outside the core (m 3 ) 1 1 18 1 1 Inlet temperature (K) 1 1 923 1 1 Outlet temperature (K) 1 1 1023 1 1 Mean salt temperature (K) 1 1 948 1 1 Flow rate (m 3 /s) 1 1 4.5 1 1 Number of heat exchangers and pump groups 1 1 16 2 1 233 U-started MSFR 3 1 TRU-started MSFR 1 1 Th 1 1 233 U 1 1 Th 2 1 Actinide 1 1 38281 kg 1 1 4838 kg 1 1 30619 kg 1 1 Pu 1 1 11079 kg 1 1 19.985 mol% 1 1 2.515 mol% 1 1 16.068 mol% 1 1 1 1 5.628 mol% 1 1 1 1 1 1 1 1 Np 1 1 789 kg 1 1 1 1 1 1 1 1 1 1 0.405 mol% 1 1 1 1 1 1 1 1 Am 1 1 677 kg 1 1 1 1 1 1 1 1 1 1 0.341 mol% 1 1 1 1 1 1 1 1 Cm 1 1 116 kg 1 1 1 1 1 1 1 1 1 1 0.058 mol% 1 1 Property 1 1 Formula 1 1 Range (K) 1 1 Density, ρ (kg/m 3 ) 1 1 4094-0.882 ×(T-1024) 1 1 893-1123 1 1 Specific heat capacity, C p (J/K‧kg) 1 1 -1111 + 2.78 T 1 1 867-907 1 1 Thermal Conductivity, λ (W/K‧m) 1 1 0.928+8.40E-05· T 1 1 891-1020 1 1 Dynamic viscosity, μ (Pa‧s) 1 1 ρ.5.55E-08.exp(3689/ T ) 1 1 898-1119 3. Methodology 1 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 10 ] and thermal-hydraulics [ 11 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 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 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. 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 13 ]. The 9 groups were classified into thr ee 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 Table 4 presents the energy ranges of the nine groups. Fig. 4 Fig. 4 shows the fission cross-sections of 233 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. 1 1 Group number 2 1 Neutron energy 1 1 1 1 Upper limit 1 1 Lower limit 1 1 1 1 1 20 MeV 1 1 1.4 MeV 1 1 2 1 1 1.4 MeV 1 1 573 KeV 1 1 3 1 1 573 KeV 1 1 73 KeV 1 1 4 1 1 73 KeV 1 1 2.29 KeV 1 1 5 1 1 2.29 KeV 1 1 186 eV 1 1 6 1 1 186 eV 1 1 52 eV 1 1 7 1 1 52 eV 1 1 33.25 eV 1 1 8 1 1 33.25 eV 1 1 12.9 eV 1 1 9 1 1 12.9 eV 1 1 1 -5 eV 3.1 Neutronics COMSOL Model 2 The nine-group diffusion and six-group delayed neutron equations (Equations 1 and 2, respectively) [ 14 14 ] can be represented by the convection and diffusion application model and transport of diluted species model, respectively, as follows: <math display="block" id="M1"><mtable columnalign="left"><mtr><mtd><mfrac><mn>1</mn><mrow><mi>ν</mi><mi>g</mi></mrow></mfrac><mfrac><mrow><mo>∂</mo><msup><mi>ϕ</mi><mi>g</mi></msup></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mi>u</mi><mrow><mi>ν</mi><mi>g</mi></mrow></mfrac><mo>⋅</mo><mo>∇</mo><mo>⋅</mo><msup><mi>D</mi><mi>g</mi></msup><mo>∇</mo><msup><mi>ϕ</mi><mi>g</mi></msup><mo>+</mo><mrow><mo>[</mo><mrow><msubsup><mo>∑</mo><mtext>a</mtext><mi>g</mi></msubsup><mo>+</mo><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn><mo>,</mo><mo>≠</mo><mn>9</mn></mrow><mn>9</mn></msubsup><msubsup><mo>∑</mo><mi>s</mi><mrow><mi>g</mi><mo>→</mo><mi>g</mi><mo>'</mo></mrow></msubsup></mrow><mo>]</mo></mrow><msup><mi>ϕ</mi><mi>g</mi></msup></mtd></mtr><mtr><mtd><mo>=</mo><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn><mo>,</mo><mo>≠</mo><mn>9</mn></mrow><mn>9</mn></msubsup><msubsup><mo>∑</mo><mi>s</mi><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup><msup><mi>ϕ</mi><mrow><mi>g</mi><mo>'</mo></mrow></msup><mo>+</mo><mfrac><mn>1</mn><mi>k</mi></mfrac><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msub><mi>β</mi><mrow><mtext>eff</mtext></mrow></msub><mo stretchy="false">)</mo><msubsup><mi>χ</mi><mi>p</mi><mn>9</mn></msubsup><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn></mrow><mn>9</mn></msubsup><mtext> </mtext><mi>ν</mi><msubsup><mo>∑</mo><mtext>f</mtext><mrow><mi>g</mi><mo>'</mo></mrow></msubsup><mtext> </mtext><msup><mi>ϕ</mi><mrow><mi>g</mi><mo>'</mo></mrow></msup></mtd></mtr><mtr><mtd><mo>+</mo><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mn>6</mn></msubsup><msubsup><mi>χ</mi><mi>i</mi><mi>g</mi></msubsup><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo>,</mo></mtd></mtr></mtable></math> <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>C</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>u</mi><mo>.</mo><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub><mo>+</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo>=</mo><mfrac><mn>1</mn><mi>k</mi></mfrac><msub><mi>β</mi><mi>i</mi></msub><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mn>9</mn></munderover><mi>ν</mi><msubsup><mtext>Σ</mtext><mtext>f</mtext><mi>g</mi></msubsup><mo> </mo><msup><mi>ϕ</mi><mi>g</mi></msup><mo> </mo><mo> </mo><mo>.</mo><mo> </mo></mrow></math> 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, <math id="M3"><mrow><mo> </mo><msup><mi>v</mi><mi>g</mi></msup><mo> </mo></mrow></math> is the neutron speed of group g (cm/s), φ and is the neutron flux (n/cm 2 ‧s 1 ), D is the diffusion constant (cm), ∑ is the macroscopic removal cross-section (cm -1 ), χ is the fraction of delayed neutrons, <math id="M4"><mi>ν</mi></math> 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, <math id="M5"><mrow><msub><mi>C</mi><mi>i</mi></msub></mrow></math> is the precursor concentration (cm -3 ), and <math id="M6"><mi>u</mi></math> is the fuel’s velocity vector (cm/s). When the stationary diffusion model includes the time-dependent coefficient partial differential equation [ 15 15 ], it is represented as follows: <math display="block" id="M7"><mrow><msub><mi>e</mi><mtext>a</mtext></msub><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>u</mi></mrow><mrow><mo>∂</mo><msup><mi>t</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><msub><mi>d</mi><mtext>a</mtext></msub><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mo>−</mo><mi>c</mi><mo>∇</mo><mi>u</mi><mo>−</mo><mi>α</mi><mi>u</mi><mo>+</mo><mi>γ</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>β</mi><mo>⋅</mo><mo>∇</mo><mi>u</mi><mo>+</mo><mi>a</mi><mi>u</mi><mo>=</mo><mi>f</mi><mtext>,</mtext></mrow></math> 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: <math display="block" id="M8"><mrow><msup><mi>λ</mi><mn>2</mn></msup><msub><mi>e</mi><mtext>a</mtext></msub><msup><mi>u</mi><mn>2</mn></msup><mo>+</mo><mi>λ</mi><msub><mi>d</mi><mtext>a</mtext></msub><mi>u</mi><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mo>−</mo><mi>c</mi><mo>∇</mo><mi>u</mi><mo>−</mo><mi>α</mi><mi>u</mi><mo>+</mo><mi>γ</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>β</mi><mo>⋅</mo><mo>∇</mo><mi>u</mi><mo>+</mo><mi>a</mi><mi>u</mi><mo>=</mo><mi>f</mi></mrow></math> 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. <math display="block" id="M9"><mrow><mi>λ</mi><msup><mi>ϕ</mi><mi>g</mi></msup><mo>−</mo><mo>∇</mo><mo>⋅</mo><msup><mi>D</mi><mi>g</mi></msup><mo>∇</mo><msup><mi>ϕ</mi><mi>g</mi></msup><mo>+</mo><mrow><mo>[</mo><mrow><msubsup><mtext>Σ</mtext><mtext>a</mtext><mi>g</mi></msubsup><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn><mo>,</mo><mo>≠</mo><mi>g</mi></mrow><mn>9</mn></munderover><msubsup><mtext>Σ</mtext><mtext>s</mtext><mrow><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>′</mo></msup></mrow></msubsup></mrow><mo>]</mo></mrow><msup><mi>ϕ</mi><mi>g</mi></msup><mo>=</mo><mo> </mo><mo> </mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn><mo>,</mo><mo>≠</mo><mi>g</mi></mrow><mn>9</mn></munderover><msubsup><mtext>Σ</mtext><mtext>s</mtext><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>→</mo><mi>g</mi></mrow></msubsup><mo> </mo><msup><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msup><mo>+</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mn>9</mn></munderover><mi>ν</mi><msubsup><mtext>Σ</mtext><mtext>f</mtext><msup><mi>g</mi><mo>′</mo></msup></msubsup><mo> </mo><msup><mi>ϕ</mi><msup><mi>g</mi><mo>′</mo></msup></msup></mrow></math> 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: <math display="block" id="M10"><mrow><mo>−</mo><mi>n</mi><mo>⋅</mo><mrow><mo>(</mo><mrow><mo>−</mo><mi>D</mi><mo>∇</mo><mi>ϕ</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn></mrow></math> <math display="block" id="M11"><mrow><mo>−</mo><mi>n</mi><mo>⋅</mo><mrow><mo>(</mo><mrow><mo>−</mo><mi>D</mi><mo>∇</mo><mi>ϕ</mi></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mrow><mn>2.1312</mn></mrow></mfrac><mi>ϕ</mi><mo> </mo><mtext>,</mtext><mo> </mo></mrow></math> 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). <math display="block" id="M12"><mrow><mo>−</mo><mi>n</mi><mo>⋅</mo><mrow><mo>(</mo><mrow><mo>−</mo><mo>∇</mo><msub><mi>C</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo> </mo></mrow></math> <math display="block" id="M13"><mrow><msub><mi>C</mi><mi>i</mi></msub><mo>=</mo><mn>0</mn><mo> </mo></mrow></math> 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 2 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: <math display="block" id="M14"><mrow><msub><mi>R</mi><mi>e</mi></msub><mo>=</mo><mfrac><mrow><mi>ρ</mi><mi>U</mi><mi>L</mi></mrow><mi>μ</mi></mfrac><mo> </mo><mtext>,</mtext><mo> </mo></mrow></math> where μ is the dynamic viscosity (Pa‧s), ρ is the density (kg/m 3 ), 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 15 , 16 16 ] become Eqs. (11) and 12 as follows: <math display="block" id="M15"><mrow><mi>ρ</mi><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>ρ</mi><mrow><mo>(</mo><mrow><mi>u</mi><mo>⋅</mo><mo>∇</mo></mrow><mo>)</mo></mrow><mi>u</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mo>−</mo><mi>p</mi><mtext>I</mtext><mo>+</mo><mi>μ</mi><mrow><mo>(</mo><mrow><mo>∇</mo><mi>u</mi><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mi>u</mi></mrow><mo>)</mo></mrow></mrow><mi>T</mi></msup></mrow><mo>)</mo></mrow><mo>+</mo><mi>F</mi><mtext>  </mtext><mo>,</mo></mrow></math> <math display="block" id="M16"><mrow><mi>ρ</mi><mo>∇</mo><mo>⋅</mo><mi>u</mi><mo>=</mo><mn>0</mn><mtext>,</mtext></mrow></math> where ρ is the density (kg/m 3 ), u is the velocity vector (m/s), P is the pressure (Pa), and F is the volume force vector (N/m 3 ). 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: <math display="block" id="M17"><mrow><msub><mi>μ</mi><mtext>T</mtext></msub><mo>=</mo><mi>ρ</mi><msub><mi>C</mi><mi>μ</mi></msub><mfrac><mrow><msup><mi>k</mi><mn>2</mn></msup></mrow><mi>ε</mi></mfrac><mo>,</mo></mrow></math> where C μ is a model constant. The transport equations for k and e (Eqs. (14) and (15), respectively) are as follows: <math display="block" id="M18"><mrow><mi>ρ</mi><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>ρ</mi><mi>u</mi><mo>⋅</mo><mo>∇</mo><mi>k</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mtext>T</mtext></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>k</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>P</mi><mi>k</mi></msub><mo>−</mo><mi>ρ</mi><mi>ε</mi><mo>,</mo></mrow></math> <math display="block" id="M19"><mrow><mi>ρ</mi><mfrac><mrow><mo>∂</mo><mi>ε</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>ρ</mi><mi>u</mi><mo>⋅</mo><mo>∇</mo><mi>ε</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mtext>T</mtext></msub></mrow><mrow><msub><mi>σ</mi><mi>ε</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>ε</mi></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>C</mi><mrow><mi>ε</mi><mn>1</mn></mrow></msub><mfrac><mi>ε</mi><mi>k</mi></mfrac><msub><mi>P</mi><mi>k</mi></msub><mo>−</mo><msub><mi>C</mi><mrow><mi>ε</mi><mn>2</mn></mrow></msub><mi>ρ</mi><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mi>k</mi></mfrac><mo>,</mo></mrow></math> where the production term is calculated by Eq. (16) as follows: <math display="block" id="M20"><mrow><msub><mi>P</mi><mi>k</mi></msub><mo>=</mo><msub><mi>μ</mi><mtext>T</mtext></msub><mrow><mo>(</mo><mrow><mo>∇</mo><mi>u</mi><mo>:</mo><mrow><mo>(</mo><mrow><mo>∇</mo><mi>u</mi><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mi>u</mi></mrow><mo>)</mo></mrow></mrow><mi>T</mi></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><msup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mo>⋅</mo><mi>u</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mi>ρ</mi><mi>k</mi><mo>∇</mo><mo>.</mo><mi>u</mi><mo> </mo><mo>.</mo></mrow></math> The heat transfer in the solid and liquid interfaces is obtained by solving Eqs. (17) and (18) as follows: <math display="block" id="M21"><mrow><mi>ρ</mi><msub><mi>C</mi><mtext>p</mtext></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><msub><mi>u</mi><mrow><mtext>trans</mtext></mrow></msub><mo>⋅</mo><mo>∇</mo><mi>T</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>q</mi><mo>+</mo><msub><mi>q</mi><mtext>r</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><msub><mi>α</mi><mtext>p</mtext></msub><mi>T</mi><mo>∶</mo><mfrac><mrow><mo>∂</mo><mi>S</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>Q</mi></mrow></math> <math display="block" id="M22"><mrow><mi>ρ</mi><msub><mi>C</mi><mtext>p</mtext></msub><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>u</mi><mo>⋅</mo><mo>∇</mo><mi>T</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>q</mi><mo>+</mo><msub><mi>q</mi><mtext>r</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>α</mi><mtext>p</mtext></msub><mi>T</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>u</mi><mo>⋅</mo><mo>∇</mo><mi>p</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>τ</mi><mo>∶</mo><mo>∇</mo><mi>u</mi><mo>+</mo><mi>Q</mi><mo> </mo><mtext>,</mtext></mrow></math> 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 3 ), and q and q r are the heat flux by conduction and radiation (W/m 2 ), respectively. 4. Results 1 The MSFR geometry was modeled using COMSOL software and the SCALE v6.1 code package ( Fig. 5 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 2 , 5 5 – 7 7 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 Fig. 6 ; however, to reduce the computation time, a small number of meshes were utilized in the solid regions. Table 5 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 233 U- and TRU-started MSFRs. 1 1 Temperature (°K) 2 1 233 U-started MSFR 2 1 TRU-started MSFR 1 1 1 1 SCALE v6.1 1 1 COMSOL 1 1 SCALE v6.1 1 1 COMSOL 1 1 923 1 1 0.9904 ± 0.0011 1 1 0.9871754±0.00049 1 1 1.0254 ± 0.0010 1 1 1.02695989±0.00041 1 1 948 1 1 0.9868 ± 0.0013 1 1 0.9848586±0.00032 1 1 1.0241 ± 0.0012 1 1 1.02579067±0.00063 1 1 1000 1 1 0.9847 ± 0.0012 1 1 0.9810806±0.00064 1 1 1.0237 ± 0.0015 1 1 1.02343186±2.6e-06 For the 233 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 233 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 233 U-started MSFRs, the maximum power density is located at the core’s center, as shown in Fig. 7 Fig. 7 . The power density at the core’s center in a TRU-started MSFR is greater than in a 233 U-started MSFR, which implies that the fast neutron flux is greater in a TRU-started MSFR compared to that in a 233 U-started MSFR. Additionally, the total fission energy release from a fissile element in TRU is greater than that in 233 U, except for 239 Pu and 238 Pu, as shown in Fig. 8 Fig. 8 . 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 233 U- and TRU-started MSFRs calculated using the thermal-hydraulics COMSOL model at a fuel flow rate of 4.5 m 3 /s for Cases I, II, and III are presented in Figs. 9 Figs. 9 , 10 10 , and 11 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 233 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 5 ] and Fiorina [ 6 6 ] of 1300 K and Hu’s result [ 7 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 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 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. 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: <math display="block" id="M23"><mrow><mi>d</mi><mi>T</mi><mo>=</mo><mfrac><mi>P</mi><mrow><mi>ρ</mi><mo>×</mo><msub><mi>C</mi><mtext>p</mtext></msub><mo>×</mo><mtext>flow rate</mtext></mrow></mfrac><mo> </mo><mtext>,</mtext><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mo> </mo><mrow><mo/><mrow><mn/></mrow><mo/></mrow></mrow></math> where P is the reactor power (3 GW), ρ is the fuel salt density (4124.87 kg/m 3 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 3 /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 Table 6 compares the k eff values calculated via the COMSOL software for the 233 U- and TRU-started MSFRs utilizing a flow of 4.5 m 3 /s for the different numbers of meshes in the core to those of the EVOL benchmark [ 3 3 ]. For the lowest mesh number (8665), the 233 U- and TRU-started COMSOL results were similar to the EVOL results. The 233 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 233 U- and TRU-started MSFRs, respectively, which are in agreement with the EVOL results for a low number of meshes. 3 1 233 U-started MSFR 3 1 TRU-started MSFR 1 1 EVOL 2 1 COMSOL 1 1 EVOL 2 1 COMSOL 1 1 1 1 No. of mesh 1 1 k eff 1 1 1 1 No. of mesh 1 1 k eff 1 1 0.98301 1 1 8665 1 1 0.97355 1 1 1.01955 1 1 8665 1 1 1.01900 1 1 1 1 25804 1 1 0.97164 1 1 1 1 25804 1 1 1.01761 1 1 1 1 217105 1 1 0.96494 1 1 1 1 217105 1 1 1.01313 The prompt neutron distributions are presented in Figs. 12 Figs. 12 , 13 13 , 14 14 , 15 15 , 16 16 , 17 17 , 18 18 , 19 19 , 20 20 . For the 233 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 14 , 1.05851×10 15 , 2.06449×10 15 , 5.69202×10 15 , and 1.05053×10 15 neutrons/(cm 2 ‧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 13 , 1.4110×10 13 , and 2.78732×10 13 neutrons/(cm 2 ‧s), respectively. The maximum flux of the ninth neutron energy group occurred at the fertile material with a value of 1.35349×10 13 neutrons/(cm 2 ‧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 14 , 1.11209×10 15 , 2.13570×10 15 , 5.57955×10 15 , and 7.85002×10 14 neutrons/(cm 2 ‧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 13 , 9.79504×10 13 , and 1.90963×10 13 neutrons/(cm 2 ‧s), respectively. The maximum flux of the ninth neutron energy group occurred at the fertile material with a value of 1.14321×10 13 neutrons/(cm 2 ‧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 Figure 21 shows the total neutron fluxes in the core, blanket, absorber, and structure regions of the 233 U- and TRU-started MSFRs. For the 233 U-started MSFR, the average flux in the core was 4.03697×10 15 neutrons/(cm 2 ‧s). The maximum flux in the core was 1.07494×10 16 neutrons/(cm 2 ‧s), which steadily decreased towards the outer region of the core until it reached 1.67329×10 15 neutrons/(cm 2 ‧s) at the fuel core/structural material boundary. The total flux in the blanket was 1.24907×10 15 neutrons/(cm 2 ‧s) and decreased to 2.90720×10 13 neutrons/(cm 2 ‧s) at the fertile blanket/second structural material boundary. The same behavior was observed for the TRU-started MSFRs. The average flux in the core was 3.62538×10 15 neutrons/(cm 2 ‧s) and the maximum flux in the core was 1.04980×10 16 neutrons/(cm 2 ‧s), which steadily decreased towards the outer region of the core until it reached 1.65252×10 15 neutrons/(cm 2 ‧s) at the fuel core/structural material boundary. The total flux in the blanket was 1.25863×10 15 neutrons/(cm 2 ‧s) and decreased to 2.99966×10 13 neutrons/(cm 2 ‧s) at the fertile blanket/second structural material boundary. 233 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 233 U-started MSFR are presented in Fig. 22 Fig. 22 . The maximum concentrations of the first through sixth groups were 6.71018×10 10 , 5.76196×10 10 , 2.08157×10 10 , 1.53499×10 10 , 1.87656×10 9 , and 4.7346×10 8 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 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. The precursor concentration distributions for the six groups of the TRU-started MSFR are presented in Fig. 23 Fig. 23 . The maximum concentrations of the first through sixth groups were 4.08564×10 10 , 1.24218×10 11 , 2.77181×10 10 , 2.25090×10 10 , 2.88048×10 9 , and 5.4025×10 8 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 Table 7 for the 233 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 233 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 −13 s after the fission event [ 18 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. 1 1 Precursor group 2 1 TRU-started MSFR 2 1 233 U-started MSFR 1 1 1 1 Decay constant,λ (1/s) 1 1 Fraction(β) 1 1 Decay constant,λ (1/s) 1 1 Fraction(β) 1 1 1 1 1 9.65E-03 1 1 7.51E-05 1 1 1.25E-02 1 1 2.91E-04 1 1 2 1 1 2.33E-02 1 1 5.23E-04 1 1 3.59E-02 1 1 6.66E-04 1 1 3 1 1 1.00E-01 1 1 4.17E-04 1 1 1.38E-01 1 1 7.66E-04 1 1 4 1 1 2.48E-01 1 1 7.30E-04 1 1 3.18E-01 1 1 1.14E-03 1 1 5 1 1 9.48E-01 1 1 2.46E-04 1 1 1.22E+00 1 1 3.27E-04 1 1 6 1 1 2.45E+00 1 1 7.46E-05 1 1 3.15E+00 1 1 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 8 , 9 9 ]. The model’s geometry, which is presented in Fig. 24a 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. 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 Fig. 24b . The fuel salt temperature and velocity distributions for the model are presented in Fig. 25 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. 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 . The temperature distributions for the Design I and II models are presented in Fig. 27 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. 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 Fig. 28 shows the proposed design, which includes 12 channels with different cross-sectional areas at the inlet. To perform the grid sensitivity analysis, we used different element sizes in the mesh. Table 8 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 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. 1 1 Element size(cm) 1 1 Maximum Temperature(K) 1 1 3 1 1 1031.840 1 1 4 1 1 1027.650 1 1 5 1 1 1029.223 1 1 6 1 1 1027.744 1 1 7 1 1 1028.949 1 1 8 1 1 1028.956 5. Conclusion 1 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 233 U- and TRU-started MSFRs were 1.07494×10 16 and 1.04980×10 16 neutrons/(cm 2 ‧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 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 https://doi.org/10.1016/B978-0-08-100906-2.00001-X . Introduction to Generation IV nuclear reactors 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 https://doi.org/10.1016/j.nucengdes.2015.06.011 Transient analyses for a molten salt fast reactor with optimized core geometry 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 https://doi.org/10.1051/epjn/2018052 Neutronic benchmark of the molten salt fast reactor in the frame of the EVOL and MARS collaborative projects 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 https://doi.org/10.1007/978-3-319-26542-1_34 Introduction to the Physics of Thorium Molten Salt Fast Reactor (MSFR) Concepts 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 https://doi.org/10.1016/j.ces.2014.03.003 Development of an OpenFOAM model for the Molten Salt Fast Reactor transient analysis 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 https://doi.org/10.1016/j.anucene.2013.08.003 Modelling and analysis of the MSFR transient behaviour 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 https://doi.org/10.1016/j.anucene.2017.05.002 Coupled neutronics and thermal-hydraulics simulation of molten salt reactors based on OpenMC/TANSY 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 . Coupled neutronics and thermal-hydraulics transient calculations based on a fission matrix approach: Application to the Molten Salt Fast Reactor 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 https://doi.org/10.1051/snamc/201402307 Coupled neutronics and thermal-hydraulics numerical simulations of a Molten Fast Salt Reactor (MFSR) D. Chandler , J. Freels , G. Maldonado , et al . COMSOL-based Nuclear Reactor Kinetics Studies at the HFIR . 2011 . COMSOL-based Nuclear Reactor Kinetics Studies at the HFIR 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 https://doi.org/10.1016/j.anucene.2016.11.040 Molten salt reactor neutronics and fuel cycle modeling and simulation with SCALE S.M. Bowman . SCALE 6: Comprehensive Nuclear Safety Analysis Code System . Nucl. Technol . 174 : 126 - 148 ( 2011 ). https://doi.org/10.13182/NT10-163 https://doi.org/10.13182/NT10-163 SCALE 6: Comprehensive Nuclear Safety Analysis Code System M. Brovchenko , E. Merle . Optimization of the pre-conceptual design of the MSFR - EVOL project . 2013 ; 2 . Optimization of the pre-conceptual design of the MSFR - EVOL project M.A. Schultz . Dynamics of Nuclear Reactors . Nucl. Technol . 12 : 415 ( 1971 ). https://doi.org/10.13182/NT71-A30993 https://doi.org/10.13182/NT71-A30993 Dynamics of Nuclear Reactors COMSOL . Introduction to COMSOL Multiphysics 5.3 . 2014 ; Introduction to COMSOL Multiphysics 5.3 E.L. Paul , V.A. Atiemo-Obeng , S.M. Krest . Handbook of Industrial Mixing: Science and Practice - Wiley Online Library . John Wiley Sons . 2004 . Handbook of Industrial Mixing: Science and Practice - Wiley Online Library S. Mulligan . Experimental and numerical analysis of three-dimensional free-surface turbulent vortex flows with strong circulation . 2015 . Experimental and numerical analysis of three-dimensional free-surface turbulent vortex flows with strong circulation US Department of Energy . DOE FUNDAMENTALS HANDBOOK NUCLEAR PHYSICS 2, 1993 . https://doi.org/DOE-HDBK-1019/2-93 https://doi.org/DOE-HDBK-1019/2-93 DOE FUNDAMENTALS HANDBOOK NUCLEAR PHYSICS 2, 1993 10.1007/s41365-020-00833-3 2020-07-26 2020-10-21 2020-10-22 2020-11-26 2020-12-01 2021-04-10 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2020 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-020-00833-3 http://dx.doi.org/10.1007/s41365-020-00833-3 . 2020 Nuclear Science and Techniques 12 31