Flow field effect of delayed neutron precursors in liquid-fueled molten salt reactors

Flow field effect of delayed neutron precursors in liquid-fueled molten salt reactors 1 2 corresp first-author *mscheng@sinap.ac.cn *mscheng@sinap.ac.cn mscheng@sinap.ac.cn https://orcid.org/0000-0002-0723-266X 1 2 https://orcid.org/0000-0002-5570-0384 1 2 1 2 1 2 Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China In molten salt reactors (MSRs), the liquid fuel salt circulates through the primary loop and a part of the delayed neutron precursors (DNPs) decays outside the reactor core. To model and analyze the flow field effect of DNPs in channel-type liquid-fueled MSRs, a three-dimensional space-time dynamics code, named ThorCORE3D, that couples neutronics, core thermal-hydraulics, and a molten salt loop system was developed and validated with the Molten Salt Reactor Experiment (MSRE) benchmarks. The effects of external loop recirculation time, fuel flow rate, and core flow field distribution on the delayed neutron fraction loss of MSRE at steady-state were modeled and simulated using the ThorCORE3D code. Then, the flow field effect of the DNPs on the system responses of the MSRE in the reactivity insertion transient under different initial conditions was analyzed systematically for the channel-type liquid-fueled MSRs. The results indicate that the flow field condition has a significant effect on the steady-state delayed neutron fractions and will further affect the transient power and temperature responses of the reactor system. The analysis results for the effect of the DNP flow field can provide important references for the design optimization and safety analysis of liquid-fueled MSRs. Molten Salt Reactor Delayed Neutron Precursor Nodal expand method Coupled neutronics and thermal-hydraulics Introduction 1 Liquid-fueled molten salt reactors (MSRs), first designed, constructed, and operated by the Oak Ridge National Laboratory (ORNL) in the United States during the 1940s, were selected by the Generation IV Nuclear Energy Systems International Forum as one of six Generation IV reactor systems [ 1 1 ]. Unlike traditional solid fuel nuclear reactors, liquid-fueled MSRs use high-temperature molten salt as the fuel and coolant. High-temperature fuel salt circulates in the primary loop, and a part of the delayed neutron precursors (DNPs) decays in the primary loop outside the core, which can cause delayed neutron loss. Therefore, the neutroni cs behavior of the liquid-fueled MSR is tightly coupled with the thermo-hydraulics of the system because the reactivity contribution of the DNP is strongly influenced by the fluid flow conditions. As the delayed neutron fraction (DNF) can significantly affect the dynamic characteristics of nuclear reactors, DNP effects are very important in reactor transient and safety analyses. Considering the unique flow characteristics of liquid fuel salt in liquid-fueled MSRs, the DNP flow effects have become a focus in contemporary liquid-fueled MSR research and development. Owing to the introduction of liquid fuel, the delayed neutron precursors circulate in the primary loop with the carrier molten salt. The motion of the DNPs is one of the main challenges. Meanwhile, the fission heat directly released in the liquid fuel. Therefore, the thermal fluid with internal heat source is another key challenge. As a result, the physical model and corresponding codes for traditional solid-fuel nuclear reactors are not suitable for liquid-fueled MSRs, and it is necessary to re-derive the mathematical model and develop a new dynamic code for liquid-fueled MSRs. Currently, many researchers have developed new models and codes for liquid-fueled MSR and conducted flow effect analyses of the DNP. Lapenta et al. [ 2 2 ] deduced a new point kinetics equation for liquid-fueled MSRs by considering a 1D imposed velocity field without temperature feedback. Shi and Li et al. [ 3 3 , 4 4 ] successively extended the RELAP5/MOD4.0 code, added 0D and 1D DNP models with the primary loop, validated with the Molten Salt Reactor Experiment (MSRE) benchmarks, and simulated the different transient scenarios of the Molten Salt Breeder Reactor (MSBR) including the load demand change, primary flow, secondary flow, and reactivity transients. Diniz et al. [ 5 5 ] derived analytical solutions and a numerical method to analyze the influence of recirculation time through the external loop on the effective DNFs. Zhu et al. [ 6 6 ] added the DNP flow model to the Monte Carlo N-Particle Transport (MCNP) code and calculated the DNF of a small molten salt reactor under different steady-state flow conditions. Furthermore, Aufiero et al. [ 7 7 ] extended the SERPENT-2 Monte Carlo code to analyze the effective DNF at steady-state of the Molten Salt Fast Reactor (MSFR). Krepel et al. [ 8 8 ] developed a core dynamics code, DYN3D-MSR, based on the Forschungszentrum Rossendorf (FZR) in-house light water reactor code, DYN3D, which allows transient simulation by 3D neutronics and parallel channel thermal-hydraulics, and simulated the steady-state DNF loss and several transient scenarios of MSRE and MSBR. Zhuang et al. [ 9 9 ] developed a nodal expand method (NEM) code, MOREL, and performed dynamic analyses on MSRE and the Thorium Molten Salt Reactor (TMSR) under perturbations of fuel pump start-up and coast-down, and by overheating and overcooling the inlet fuel temperature [ 10 10 ]. Zhang et al. [ 11 11 ] developed a coupled code, COUPLE, by solving the multi-group neutron diffusion equations and the incompressible Navier-Stokes equations with the standard k- ε model under 2D cylindrical geometry. This code was applied to the conceptual analysis of MSFR and the Molten Salt Actinide Recycler & Transmuter (MOSART) [ 12 12 , 13 13 ]. Cammi et al. [ 14 14 - 16 16 14-16 ] developed a multi-physics simulation code based on the commercial software COMSOL Multiphysics, which was used to model and analyze the steady-state and transient characteristics of MSBR, MSRE, and MSFR. Cervi et al. [ 17 17 ] implemented an SP3 neutron transport solver and a neutron diffusion solver in OpenFOAM to analyze MSFR. A framework, named GeN-ROM [ 18 18 ], was developed using OpenFOAM and employs a proper orthogonal decomposition-aided reduced-basis technique (POD-RB). This framework was tested using a 2D multi-physics MSFR model with steady-state and transient scenarios. Tiberga et al. [ 19 19 ] from the Delft University of Technology developed a code coupling an incompressible Reynolds-averaged Navier-Stokes (RANS) model and an SN neutron transport solver and performed preliminary steady-state and transient simulations of MSFR. An open-source MSR simulation code, Moltres, was developed based on the Multi-physics Object-Oriented Simulation Environment (MOOSE) framework. In Moltres, the neutron diffusion, temperature, and precursor equations are fully coupled and solved with implicit time stepping. Moltres was utilized to study the steady-state and transient behavior of MSFR [ 20 20 ]. The PROTEUS-NODAL code was modified for MSR and coupled with the SAM system code under the MOOSE framework [ 21 21 ]. This coupled code system was verified and validated using the MSFR benchmark problem and MSRE experimental data. With these developed codes, the unique phenomena of fuel flow in liquid-fueled MSRs have been well studied, especially in no-moderator fast-spectrum MSR concepts, such as MSFR [ 7 7 , 16 16 , 12 12 , 17 17 ]. However, on the basis of the models and findings of the researchers mentioned above, the entire primary and secondary loop models need to be further considered and extended in channel-type liquid-fueled MSRs. To more systematically analyze the flow field effect of the DNP in channel-type liquid-fueled MSRs, in this study, the ThorCORE3D code with coupled 3D core dynamics and 1D loop dynamics models of channel-type liquid-fueled MSRs was developed and validated with MSRE benchmarks. The DNF of MSRE under different flow conditions, including external loop recirculation time, fuel salt flow rate, and core flow field distribution, were modeled and simulated using the ThorCORE3D code. Subsequently, the flow field effects of DNPs in the reactivity insertion transient of MSRE under different flow field conditions were modeled and analyzed in detail. Theory and method 1 Neutronics model 2 In liquid-fueled MSRs, liquid fuel salt circulates through the primary loop. The flow of the fuel salt has little effect on the neutron flux because the velocity is much lower than the neutron velocity [ 22 22 ]; however, the DNPs circulate with the liquid fuel salt and partially decay in the primary loop outside the reactor, which results in a reduction in the effective DNF and loss of reactivity. The multi-group diffusion theory is used to deduce the neutronics model of the liquid-fueled MSRs, which consists of 3D space-time multi-group neutron diffusion equations and DNP balance equations [ 23 23 ], given in Eq. (1) and (2), for node n. <math display="block" id="M1"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mo>∂</mo><msubsup><mi>Φ</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><msubsup><mi>v</mi><mi>g</mi><mi>n</mi></msubsup><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><mo>∇</mo><mo>⋅</mo><msubsup><mi>D</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>∇</mo><msubsup><mi>Φ</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mtext>t</mtext><mn>,</mn><mi>g</mi></mrow><mi>n</mi></msubsup><mrow><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mi>Φ</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>=</mo><mfrac><mrow><msub><mi>χ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msup><mi>β</mi><mi>n</mi></msup></mrow><mo>)</mo></mrow><msubsup><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mi>Φ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>      </mtext><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>Σ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>→</mo><mi>g</mi></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mi>Φ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>g</mi></mrow></msub><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>          </mtext><mi>g</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>G</mi></mtd></mtr></mtable></math> <math display="block" id="M2"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mo>∂</mo><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><msup><mi>u</mi><mi>n</mi></msup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>       </mtext><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>D</mi><mrow><mtext>c</mtext><mn>,</mn><mi>i</mi></mrow></msub><mo>+</mo><msub><mi>D</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow><mo>∇</mo><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>     </mtext><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>β</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mi>Φ</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>,</mn><mi>j</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>M</mi></mtd></mtr></mtable></math> where Φ is the scalar flux; v is the neutron velocity; D g is the neutron diffusion coefficient; ∑ t , ∑ f , and <math id="M3"><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msub><mrow/><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>→</mo><mi>g</mi></mrow></msub></mrow></mstyle></mrow></math> represent the macroscopic total, macroscopic fission, and macroscopic scattering cross-sections, respectively; k eff is the effective multiplication factor; C is the DNP concentration; λ is the DNP decay constant; β is the total fraction of delayed neutrons; β j is the fraction of delayed neutrons in each group; χ p and χ d, j represent the energy spectrum of prompt and delayed neutrons, respectively; u is the fuel salt velocity; D c, i and D t are the DNP molecular and turbulent diffusion coefficients. The time-dependent terms of the neutron multi-group diffusion equations are discretized using the backward Euler method together with an exponential transformation technique. These equations are then solved using NEM [ 24 24 ]. As the fuel flows through different separated channels in the core, the 3D DNP balance equations are simplified to a 1D form. Additionally, as molecular and turbulent diffusion are much smaller than convection transport in channel-type liquid-fueled MSRs, they can be neglected. These simplified 1D DNP balance equations are as follows: <math display="block" id="M4"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><mo>∂</mo><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msup><mi>u</mi><mi>n</mi></msup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mo>−</mo><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>     </mtext><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>β</mi><mtext>j</mtext><mi>n</mi></msubsup><msubsup><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mi>Φ</mi><mi>g</mi><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>,</mn><mi>j</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>M</mi></mtd></mtr></mtable></math> Applying the backward Euler method for temporal discretization and first-order upwind scheme for spatial discretization [ 25 25 ]: <math display="block" id="M5"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mrow><mo>∂</mo><msup><mi>C</mi><mi>n</mi></msup></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>≈</mo><mfrac><mrow><msubsup><mi>C</mi><mi>t</mi><mi>n</mi></msubsup><mo>−</mo><msubsup><mi>C</mi><mrow><mi>t</mi><mo>−</mo><mi>Δ</mi><mi>t</mi></mrow><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M6"><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo stretchy="false">(</mo><msup><mi>u</mi><mi>n</mi></msup><msup><mi>C</mi><mi>n</mi></msup><mo stretchy="false">)</mo><mo>≈</mo><mfrac><mrow><msup><mi>u</mi><mi>n</mi></msup><msup><mi>C</mi><mi>n</mi></msup><mo>−</mo><msup><mi>u</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><msup><mi>C</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup></mrow><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac></mrow></math> Then, Eq. (3) can be written as follows: <math display="block" id="M7"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mi>u</mi><mi>n</mi></msup></mrow><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><mo>+</mo><mi>λ</mi></mrow><mo>)</mo></mrow><msup><mi>C</mi><mi>n</mi></msup><mo>−</mo><mfrac><mrow><msup><mi>u</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup></mrow><mrow><mi>Δ</mi><mi>z</mi></mrow></mfrac><msup><mi>C</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>=</mo><msup><mi>Q</mi><mi>n</mi></msup><mo>+</mo><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><msubsup><mi>C</mi><mrow><mi>t</mi><mo>−</mo><mi>Δ</mi><mi>t</mi></mrow><mi>n</mi></msubsup></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M8"><mrow><msup><mi>Q</mi><mi>n</mi></msup><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>β</mi><mi>g</mi><mi>n</mi></msubsup><mi>ν</mi><msubsup><mi>Σ</mi><mrow><mtext>f</mtext><mn>,</mn><mi>g</mi></mrow><mi>n</mi></msubsup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msubsup><mtext>Φ</mtext><mi>g</mi><mi>n</mi></msubsup><mn>,</mn><mtext> core</mtext></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mtext>       </mtext><mn>0</mn><mo>,</mo><mtext>external loop</mtext></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> The DNP of the different channels is homogenously mixed at the core outlet. At each time step, Eq. (6) is solved for each channel in the core and external loop using the Gauss-Seidel iteration method. The same space discretization mesh is adopted by the 3D neutron diffusion equations and 1D DNP balance equations, and the node average DNP concentration is transferred to the 3D neutron diffusion equations in the same node. To calculate the reactor kinetics parameters and DNF loss, adjoint equations including adjoint DNP balance equations are proposed for liquid-fueled MSRs, as follows: <math display="block" id="M9"><mtable columnalign="left"><mtr><mtd><mo>−</mo><mo>∇</mo><mo>⋅</mo><msubsup><mi>D</mi><mi>g</mi><mi>n</mi></msubsup><mo>∇</mo><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mi>Σ</mi><mrow><mtext>t</mtext><mn>,</mn><mi>g</mi></mrow><mi>n</mi></msubsup><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msup><mi>β</mi><mi>n</mi></msup></mrow><mo>)</mo></mrow><msubsup><mrow><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mi>g</mi><mi>n</mi></msubsup></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msub><mi>χ</mi><mrow><mi>p</mi><mn>,</mn><msup><mi>g</mi><mo>'</mo></msup></mrow></msub><msubsup><mtext>Φ</mtext><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>'</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>Σ</mi><mrow><mi>g</mi><mo>→</mo><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><msubsup><mtext>Φ</mtext><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mi>g</mi><mi>n</mi></msubsup><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msubsup><mi>β</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mtext>        </mtext><mi>g</mi><mo>=</mo><mn>1,2,...,</mn><mi>G</mi></mtd></mtr></mtable></math> <math display="block" id="M10"><mtable><mtr><mtd><mo>−</mo><mfrac><mrow><mtext>d</mtext><mrow><mo>(</mo><mrow><msup><mi>u</mi><mi>n</mi></msup><msubsup><mi>C</mi><mi>j</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow><mrow><mtext>d</mtext><mi>z</mi></mrow></mfrac><mo>=</mo><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi></mrow></msub><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mi>j</mi><mo>=</mo><mn>1,2,</mn><mo>…</mo><mn>,</mn><mi>M</mi></mtd></mtr></mtable></math> where <math id="M11"><mrow><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></math> and <math id="M12"><mrow><msubsup><mi>C</mi><mi>j</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></math> are the adjoint flux and adjoint DNP concentration, respectively. The adjoint equations are in the same form as the forward equations, and the same solver can be used to solve these equations, with several coefficients replaced. Then, the effective delayed neutron fractions β eff and effective DNP decay constants λ eff can be calculated as follows: <math display="block" id="M13"><mrow><msub><mi>β</mi><mrow><mtext>eff</mtext><mn>,</mn><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mo>〈</mo><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo>|</mo><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi></mrow></msub><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo>〉</mo></mrow></mstyle></mrow><mi>F</mi></mfrac></mrow></math> <math display="block" id="M14"><mrow><msub><mi>λ</mi><mrow><mtext>eff</mtext><mn>,</mn><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mo>〈</mo><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo>|</mo><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi></mrow></msub><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo>〉</mo></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></msubsup><mrow><mo>〈</mo><msubsup><mtext>Φ</mtext><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo>|</mo><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi></mrow></msub><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup><mo>〉</mo></mrow></mstyle></mrow></mfrac></mrow></math> <math display="block" id="M15"><mtable><mtr><mtd><mi>F</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mrow><mo>〈</mo><mrow><msubsup><mi>Φ</mi><mi>g</mi><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mo>|</mo><msub><mi>χ</mi><mrow><mi>d</mi><mn>,</mn><mi>j</mi></mrow></msub><msubsup><mi>λ</mi><mi>j</mi><mi>n</mi></msubsup><msubsup><mi>C</mi><mi>j</mi><mi>n</mi></msubsup></mrow><mo>〉</mo></mrow><mo>+</mo></mtd></mtr><mtr><mtd><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi></mrow><mo>)</mo></mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>'</mo><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover><mrow><mrow><mo>〈</mo><mrow><msubsup><mi>Φ</mi><msup><mi>g</mi><mo>′</mo></msup><mrow><mo>†</mo><mi>n</mi></mrow></msubsup><mrow><mo>|</mo><mrow><msub><mi>χ</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub></mrow></mrow><msubsup><mrow><mrow><mo>(</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup><msubsup><mi>Φ</mi><mrow><msup><mi>g</mi><mo>'</mo></msup></mrow><mi>n</mi></msubsup></mrow><mo>〉</mo></mrow></mrow></mstyle></mtd></mtr></mtable></math> Homogenized few-group parameters 2 The homogenized few-group parameters required for neutronics calculation include few-group macroscopic cross-sections, neutron energy spectrum, delayed neutron parameters, etc. Considering the thermal feedback effect, the lattice code DRAGON [ 26 26 ] is used to generate few-group parameters at different fuel and graphite temperatures. During calculation, temperature-dependent few-group parameters are obtained using the bilinear interpolation method [ 27 27 ]. Core thermal-hydraulics model 2 A parallel multi-channel core thermal-hydraulic model is implemented in the ThorCORE3D code to account for thermal feedback and fuel flow effects. In this model, the fuel channel and surrounding graphite components are approximated as cylindrical tubes, as shown in Fig. 1 Fig. 1 . Owing to the very high boiling point of molten salt, a single-phase flow model is adopted in the core thermal-hydraulics model. The 1D fluid momentum, energy, and mass balance equations can be written as follows: <math display="block" id="M16"><mrow><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>u</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><msup><mi>u</mi><mn>2</mn></msup><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>p</mi><mrow><mtext>fric</mtext></mrow></msub></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mi>g</mi><mi>ρ</mi><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mn>0</mn></mrow></math> <math display="block" id="M17"><mrow><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>H</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>u</mi><mi>H</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><msub><mi>Q</mi><mtext>f</mtext></msub><mo>+</mo><msub><mi>Q</mi><mtext>w</mtext></msub></mrow></math> <math display="block" id="M18"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>u</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mn>0</mn></mrow></math> where ρ is the density of the fuel salt; u is the fuel salt velocity; p is the pressure; H is the specific enthalpy of the fuel salt; g is the acceleration of gravity; p fric is the frictional pressure drop; Q w is the heat transfer between the fluid and graphite components; Q f is the volumetric heat source in the fuel salt, which can be calculated as follows: <math display="block" id="M19"><mrow><msub><mi>Q</mi><mtext>f</mtext></msub><mo>=</mo><msub><mi>α</mi><mtext>f</mtext></msub><mo>⋅</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>g</mi><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msub><mi>Σ</mi><mrow><mtext>kf</mtext><mn>,</mn><mi>g</mi></mrow></msub><msub><mi>Φ</mi><mi>g</mi></msub></mrow></math> where ∑ kf, g is the macroscopic kappa-fission cross-section of group g , and α f is the proportion of fission energy released in the fuel salt. The graphite components are approximated as hollow circular tubes with an inner radius of r 0 and an outer radius of R . Similar to the heat transfer model of fuel rods in solid fuel nuclear reactors, the heat conduction equation of the graphite components in cylindrical geometry is defined as follows: <math display="block" id="M20"><mrow><msub><mi>ρ</mi><mi>g</mi></msub><msub><mi>c</mi><mrow><mi>p</mi><mn>,</mn><mi>g</mi></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mi>g</mi></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><msub><mi>λ</mi><mi>g</mi></msub></mrow><mi>r</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mi>r</mi><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mi>g</mi></msub></mrow><mrow><mo>∂</mo><mi>r</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>Q</mi><mi>g</mi></msub></mrow></math> with the following boundary conditions: <math display="block" id="M21"><mrow><msub><mrow><mrow><mrow><msub><mi>λ</mi><mi>g</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>T</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>r</mi></mrow></mfrac></mrow><mo>|</mo></mrow></mrow><mrow><mi>r</mi><mo>=</mo><msub><mi>r</mi><mn>0</mn></msub></mrow></msub><mo>=</mo><mi>h</mi><mrow><mo>(</mo><mrow><msub><mi>T</mi><mi>g</mi></msub><mo stretchy="false">(</mo><msub><mi>r</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>−</mo><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow></math> <math display="block" id="M22"><mrow><msub><mrow><mrow><mrow><msub><mi>λ</mi><mi>g</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>T</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>r</mi></mrow></mfrac></mrow><mo>|</mo></mrow></mrow><mrow><mi>r</mi><mo>=</mo><mi>R</mi></mrow></msub><mo>=</mo><mn>0</mn></mrow></math> where T g denotes the graphite temperature, ρ g denotes the graphite density, c p,g denotes the graphite specific heat capacity, λ g denotes the graphite thermal conductivity, r denotes the distance from the central line, h denotes the convective heat transfer coefficient, and Q g denotes the volumetric heat source in the graphite. The momentum and mass balance equations are solved using the backward Euler and finite difference methods. The energy balance equation is solved using the method of characteristics (MOC) [ 28 28 ], and the fuel temperature is used as the boundary condition for graphite heat conduction. It is assumed that the fuel salt flows separately in the lower and upper plenums, meaning that the fuel salt is distributed to each channel at the core inlet and mixed at the core outlet. The numerical models in the lower and upper plenums are the same as those in the channels. For the heat conduction equation, the circular graphite component is divided into a series of concentric circular tubes of equal areas. The relationship between the average temperatures of neighboring tubes can be obtained analytically. Steady-state and transient equations can be transformed into the same form, which can then be solved using the tridiagonal matrix algorithm (TDMA). Molten salt loop model 2 Generally, liquid-fueled MSRs consist of two molten salt loops, that is, the fuel salt (primary) loop and the secondary coolant salt loop. Therefore, it is important to consider the molten-salt loop model in the reactor system response. With the fuel salt in the liquid-fueled MSRs circulating in the primary loop, the DNP balance equations and the corresponding numerical methods in the primary loop are described in Sect. 2.1. In the loop model, the pipes and heat exchangers are approximately 1D and the heat loss in the pipes is ignored. As shown in Fig. 2 Fig. 2 , the fuel salt flows in the direction opposite to that of the cooling salt in the heat exchanger. The energy balance equation in the fuel and coolant salts is defined as follows: <math display="block" id="M23"><mrow><mi>ρ</mi><msub><mi>c</mi><mrow><mtext>p</mtext><mn>,</mn><mtext>f</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mi>ρ</mi><mi>u</mi><msub><mi>c</mi><mrow><mtext>p</mtext><mn>,</mn><mtext>f</mtext></mrow></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><msub><mi>Q</mi><mtext>f</mtext></msub><mo>+</mo><msub><mi>Q</mi><mtext>w</mtext></msub></mrow></math> where c p,f denotes the fluid specific heat capacity. Q f =0 is always satisfied because the decay heat is ignored. Q w is the heat transfer between the fluid and steel wall; in the pipes, Q w =0, while in the heat exchanger, Q w was calculated as follows: <math display="block" id="M24"><mrow><msub><mi>T</mi><mtext>f</mtext></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>T</mi><mtext>w</mtext></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>φ</mi><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo></mrow><mi>h</mi></mfrac><mo>=</mo><mo>−</mo><mfrac><mrow><msub><mi>Q</mi><mtext>w</mtext></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mi>A</mi></mrow><mrow><mi>π</mi><mi>D</mi><mi>h</mi></mrow></mfrac></mrow></math> where φ ( z ) denotes the heat flux through the wall, A denotes the flow area, and D is the hydraulic diameter. Using the fuel salt temperature and velocity from the core outlet as the inlet boundary conditions, the finite difference method is used to solve Eq. (20) and (21) at steady-state, and the MOC method is used to solve the equations in the transient state. Next, the fuel salt temperature in the loop outlet is transferred as the inlet boundary condition of the core thermal-hydraulics model. Coupling scheme 2 The calculation flow scheme of ThorCORE3D is shown in Fig. 3 Fig. 3 . Before steady-state or transient calculations, the few-group parameters must be generated at different operating temperatures. At the beginning of the calculation, a uniform power distribution is assumed, and the first thermal-hydraulic calculation of the entire system is performed. Subsequently, the few-group parameters are interpolated and updated according to the calculated thermal fluid parameters. The neutronics calculations, including the DNP balance equations, are conducted using the updated parameters, including cross-sections and flow velocity, and a new power distribution is obtained. The iteration is carried out until the parameters (e.g., k eff and temperatures) converge. The adjoint calculation is performed after the above iteration. The transient calculation can begin after a steady-state calculation. As in the steady-state calculation, the neutronics module and thermal-hydraulics module are also calculated iteratively until convergence, and then advanced to the next time step until the given time is reached. Code validation 2 This section describes a preliminary validation of the ThorCORE3D code that was conducted based on MSRE benchmarks. Description and modeling of the MSRE 3 The MSRE is a molten salt experimental reactor designed, built, and operated by ORNL in the 1960s and is currently the only available source of experimental data for molten salt reactors with liquid fuel. The design power of the MSRE was 10 MWth, and the actual operating power was 8 MWth. The U-235 and U-233 fuels were used in the early and later stages, respectively. The main design parameters and material properties are shown in Table 1 Table 1 . 1 1   1 1 Parameters 1 1 Value 1 1 Core 1 1 Reactor power 1 1 8 MW 1 1   1 1 Mass flow rate 1 1 187.5 kg/s 1 1   1 1 Fuel inlet temperature 1 1 908.15 K 1 1   1 1 Fuel transit time in external primary loop 1 1 10 s 1 1 Fuel Salt 1 1 Molar composition of U-235 fuel salt 1 1 LiF-BeF 2 -ZrF 4 -UF 4 1 1   1 1   1 1 (65-29.2-5-0.8) 1 1   1 1   1 1 35% U-235, 99.99% Li 1 1   1 1 Molar composition of U-233 fuel salt 1 1 LiF-BeF 2 -ZrF 4 -UF 4 1 1   1 1   1 1 (64.50-30.17-5.19-0.138) 1 1   1 1   1 1 91.5 % U-233, 99.99% Li 1 1   1 1 Density 1 1 2322×(1-2.12×10 -4 ×( T (K)-922))kg/m 3 1 1   1 1 Specific heat capacity 1 1 1.968 kJ/(kg·K) 1 1   1 1 Thermal Conductivity 1 1 5.56 W/(m·K) 1 1   1 1 Dynamic Viscosity 1 1 0.009 Pa·s 1 1 Graphite 1 1 Density 1 1 1860 kg/m 3 1 1   1 1 Thermal Conductivity 1 1 58.84 W/(m·K) 1 1 Control rod 1 1 Composition 1 1 30%(wt) Al 2 O 3 , 70%(wt) Gd 2 O 3 1 1   1 1 Density 1 1 5873 kg/m 3 A ThorCORE3D model of the MSRE system is shown in Fig. 4 Fig. 4 , and the two loops of the MSRE were modeled. The MSRE core model is shown in Fig. 5 Fig. 5 and 6 6 . The upper and lower plenums were equivalent to a cylinder with a height of 17.15 cm, and the height of the active core was 166.37 cm. The reactor was divided into 20 nodes along the axial direction, as shown in Fig. 5 Fig. 5 . Each fuel channel and the surrounding graphite were radially divided into a square segment, and each control rod channel and sample channel were equivalent to four channels, with a total of 1372 channels in the radial direction, as shown in Fig. 6 Fig. 6 . The lattice code DRAGON with the ENDF/B-VII.0 data library was used to generate the two-group parameters required by the ThorCORE3D code. Two-group parameters were generated at different fuel and graphite temperatures, and the expansion effects were also considered. The DNP parameters were obtained from the ORNL report, as listed in Table 2 Table 2 . 1 1   1 1   1 1 Total 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 U-235 [ 29 29 ] 1 1 Delayed neutron fraction β (10 -5 ) 1 1 640.5 1 1 21.1 1 1 140.2 1 1 125.4 1 1 252.8 1 1 74.0 1 1 27 1 1   1 1 Decay constant λ (s -1 ) 1 1   1 1 0.0124 1 1 0.0305 1 1 0.111 1 1 0.301 1 1 1.14 1 1 3.01 1 1 U-233 [ 32 32 ] 1 1 Delayed neutron fraction β (10 -5 ) 1 1 271.2 1 1 22.55 1 1 80.28 1 1 67.48 1 1 76.86 1 1 14.94 1 1 9.10 1 1   1 1 Decay constant λ (s -1 ) 1 1   1 1 0.0126 1 1 0.0337 1 1 0.139 1 1 0.325 1 1 1.13 1 1 2.5 Loss of DNF at steady-state operation 3 The effective DNF under the flow and static conditions of the fuel salt were calculated using Eq. (10), and the difference between the two states corresponds to the DNF loss caused by the molten fuel salt flow effect. As shown in Tables 3 Tables 3 and 4 4 , a comparison of the ThorCORE3D numerical results, MSRE experimental values, ORNL numerical results, and other code numerical results with the loads of U-235 fuel and U-233 fuel are listed. The DNF losses calculated by the ThorCORE3D code agree well with the MSRE experimental values and are within the range of the values calculated by the other codes. The numerical results of the ThorCORE3D code using the DNP parameters from the ENDF/B-VII.0 library are close to the results obtained using parameters from the ORNL report. Hence, the DNP parameters from the ORNL report were directly applied in subsequent simulations. 1 1 U-235 salt 1 1 Total 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 MSRE 1 1 212 1 1 ORNL 1 1 222 1 1 EDF 1 1 207.6 1 1 12.8 1 1 65.6 1 1 55.3 1 1 71.4 1 1 2.5 1 1 0 1 1 FZK 1 1 258.3 1 1 16.4 1 1 83 1 1 68.3 1 1 85.8 1 1 4.5 1 1 0.3 1 1 FZR 1 1 223 1 1 10.2 1 1 74.6 1 1 60.5 1 1 75.1 1 1 2.6 1 1 0 1 1 MOREL 1 1 246.3 1 1 13.6 1 1 88.9 1 1 67.7 1 1 73.6 1 1 2.5 1 1 0 1 1 ThorCORE3D_ENDF 1 1 228.5 1 1 13.7 1 1 85.8 1 1 60.6 1 1 64.6 1 1 3.4 1 1 0.3 1 1 ThorCORE3D_ORNL 1 1 224.5 1 1 13.3 1 1 84.2 1 1 59.6 1 1 63.6 1 1 3.5 1 1 0.3 1 1 U-233 salt 1 1 Total 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 MSRE 1 1 100.5 1 1 ORNL 1 1 100.5 1 1 12.39 1 1 43.88 1 1 28.2 1 1 15.6 1 1 0.35 1 1 0.05 1 1 EDF 1 1 107.8 1 1 12.4 1 1 44.1 1 1 30.8 1 1 20 1 1 0.5 1 1 0 1 1 FZK 1 1 125 1 1 14.5 1 1 51.7 1 1 35.1 1 1 22.7 1 1 0.8 1 1 0.1 1 1 FZR 1 1 120.4 1 1 14.3 1 1 50.2 1 1 33.4 1 1 21.1 1 1 1.4 1 1 0 1 1 MOREL 1 1 102.2 1 1 13.5 1 1 38.2 1 1 24.0 1 1 26.1 1 1 0.7 1 1 0.0 1 1 ThorCORE3D 1 1 110.3 1 1 14.2 1 1 47.9 1 1 29.3 1 1 18.0 1 1 0.7 1 1 0.1 Protected pump start-up and coast-down experiments 3 The protected pump start-up and coast-down experiments were part of the MSRE zero-power experiments in ORNL, and U-235 fuel was used during the experiment. During the pump start-up transient, the fuel pump of the primary loop was started at time zero, and the inlet flow rate started to change at 2 s and increased to the nominal value within 8 s. During the pump coast-down transient, the fuel flow rate decreased from the nominal values to zero in 20 s, as shown in Fig. 7 Fig. 7 . As the molten-fuel salt flowed out of the core, the distribution of DNP concentration changed, causing a loss of reactivity. A power of 1 kWth was maintained in the protected pump start-up and coast-down experiments using a control rod, and the change in the position of the control rod with time was recorded. With the core power low and basically constant, the temperature feedback was ignored when the protected transients were simulated using the ThorCORE3D code. By searching for the position of the control rod to maintain the core power, the change in the control rod position was directly calculated, as shown in Fig. 8 Fig. 8 . Then, the rod position was converted into reactivity using the control rod value and compared with the experimental values and calculated values of other codes [ 8 8 , 9 9 ], as shown in Fig. 9 Fig. 9 and 10 10 . As these figures show, the numerical results of the ThorCORE3D code are in good agreement with the experimental values. Natural circulation experiment 3 During the natural circulation experiment, MSRE was operated using the U-233 fuel. The primary molten salt pump was kept closed and no control rod movement was considered. The fuel molten salt flow was driven only by natural circulation. At the beginning of the natural circulation transient, the reactor power was maintained at 4.1 kWth with a limited fuel mass flow rate. Subsequently, the cooling capacity of the secondary loop increased gradually; thus, the core inlet temperature decreased, and the core power increased owing to the reactivity negative feedback effect of temperature [ 30 30 ]. During the simulation, changes in the inlet temperature and fuel salt mass flow rate were used as inputs, as shown in Fig. 11 Fig. 11 . The variation in the core power was calculated and is shown in Fig. 12 Fig. 12 . There are a few different behaviors in the first 2000 s. This is because the experiment did not start from a steady state [ 33 33 ], whereas the Thor-CORE3D code did. However, the numerical results at other times are in good agreement with the experimental values. Reactivity insertion experiment 3 At ORNL, a series of experiments on MSRE with the U-233 fuel were performed to test the power response to reactivity insertion. In the 8 MWth case, a reactivity of 1.39× 10 -4 was inserted in 1 s and the reactor power increased rapidly. Subsequently, the core temperature increased and caused a reactivity negative feedback effect, resulting in a slower increase in power and finally dropping to the initial level [ 34 34 ]. To fully consider the variation in the core inlet temperature with time, the MSRE molten salt loops were modeled with the ThorCORE3D code, and the temperature fields of the primary and secondary loops were solved. The numerical results were compared with the experimental values of MSRE, as shown in Fig. 13 Fig. 13 . The calculated results are in good agreement with the experimental values. Results and discussion 1 Flow field effect on DNP at steady-state 2 In liquid-fueled MSRs, DNPs circulate with the fuel salt and partially decay in the primary loop outside the core. The DNF loss is directly affected by the flow field. In this section, the DNF loss of the MSRE with U-235 fuel under different flow conditions, including the external loop recirculation time, fuel salt flow rate, and core flow field distribution, is modeled and simulated using the ThorCORE3D code. Effect of the recirculation time 3 In this scenario, the fuel molten salt flow rate and flow field distribution remain unchanged, and the loss of DNF on recirculation time from 0 to 20 s was simulated. As shown in Fig. 14 Fig. 14 , the DNF loss of each DNP group increases with recirculation time and gradually reaches saturation. The numerical results for recirculation times of 0 s, 5 s,10 s, and 20 s are listed in Table 5 Table 5 , and the total loss increases from 162.2 pcm at 0 s to 239.8 pcm at 20 s. To obtain the DNF loss ratio of each DNP group and evaluate the effect of recirculation time on different DNP groups, the DNF loss value was normalized using the corresponding effective DNF of the non-flowing state, and the result is shown in Fig. 15 Fig. 15 . The losses of each DNP group are affected by the recirculation time differently. The first DNP group has a higher loss percentage, whereas the sixth group reaches saturation faster, which is because the sixth DNP group has a higher decay constant, most of which decay directly in the core, and an increase in the outer loop length has little effect on it. 1 1 Recirculation time(s) 1 1 Total(pcm) 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 0 1 1 162.2 1 1 10.1 1 1 61.9 1 1 41.4 1 1 46.0 1 1 2.6 1 1 0.2 1 1 5 1 1 203.1 1 1 12.5 1 1 77.3 1 1 53.2 1 1 57.0 1 1 3.0 1 1 0.3 1 1 10 1 1 221.3 1 1 13.9 1 1 86.3 1 1 58.7 1 1 59.2 1 1 3.0 1 1 0.3 1 1 20 1 1 239.8 1 1 15.7 1 1 97.5 1 1 63.7 1 1 59.8 1 1 2.9 1 1 0.2 Notably, at the limit situation of zero recirculation time, there should be no DNPs decay outside the core, but the results show that the DNF loss is not zero. This is because the flow of molten salt leads to the homogenization of the DNP in the core. The density distributions of the six DNP groups with zero recirculation time and static fuel are presented in Fig. 16 Fig. 16 and 17 17 . More delayed neutrons are released at the core edge as the molten salt flowed, where the neutron importance is lower, leading to reactivity loss. Effect of fuel flow rate 3 The DNF loss on the fuel salt mass flow rate from zero to twice the nominal mass flow rate (187.5 kg/s) was simulated with a recirculation time of 10 s. As Fig. 18 Fig. 18 shows, the DNF loss of each DNP group increases with the fuel flow rate and gradually reaches saturation. When the flow rate is zero, there is no DNF loss in any DNP group, as in solid-fuel nuclear reactors. The results for 50 %, 100 %, and 200 % of the nominal flow rate are listed in Table 6 Table 6 . The total loss reaches 286.6 pcm in the 200 % mass flow state. Then, the DNF loss is normalized using the corresponding effective DNF of the non-flowing state. As Fig. 19 Fig. 19 shows, the DNF loss of each DNP group is affected differently by the fuel salt flow. The first DNP group has a higher percentage loss and reached saturation faster, but the last three DNP groups still do not reach saturation in the 200 % mass flow state. This is because the decay constant of the first DNP group is small and the decay is slower. When the fuel salt flow increases, the effect of the DNP returning to the core in the first DNP group is clearer, reaching the saturation state faster. 1 1 Flow rate fraction(%) 1 1 Total(pcm) 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 50 1 1 160.9 1 1 13.4 1 1 79.0 1 1 40.0 1 1 27.3 1 1 1.1 1 1 0.1 1 1 100 1 1 221.3 1 1 13.9 1 1 86.3 1 1 58.7 1 1 59.2 1 1 3.0 1 1 0.3 1 1 200 1 1 286.6 1 1 14.1 1 1 90.3 1 1 71.3 1 1 102.4 1 1 7.9 1 1 0.7 Effect of the core flow field distribution 3 The mass flow distribution of the fuel salt in the radial direction of the core was linearly fitted as follows: <math display="block" id="M25"><mrow><msub><mi>G</mi><mtext>m</mtext></msub><mo>=</mo><mi>k</mi><msub><mi>r</mi><mtext>c</mtext></msub><mo>+</mo><mn>0.149</mn></mrow></math> where G m (kg/s) denotes the mass flow rate of a single channel, r c (cm) denotes the distance from the channel to the core central line, and k (kg/(cm·s)) denotes the slope. As Fig. 20 Fig. 20 shows, the upper limit of the absolute value of k is 2.04× 10 -3 kg/(cm·s) because the flow rate is zero in the center or edge channel at this time. k < 0 indicates that the flow rate in the central area of the reactor is relatively large, and k =0 indicates a uniform mass flow rate for each channel. DNF loss was calculated using the flow distribution represented by the shaded area in Fig. 20 Fig. 20 . As Fig. 21 Fig. 21 shows, DNF loss decreases approximately linearly as k increases. Owing to the higher neutron flux density in the central core area, the DNP concentration is also higher than that at the edge. When k increases, the flow rate in the central area of the core decreases, causing a decrease in the amount of DNP flowing out of the core, thereby decreasing DNF loss. The results of the three special cases are listed in Table 7 Table 7 , and the DNF losses differ by at most 68.2 pcm under different flow field distributions. 1 1 k (10 -3 kg/(cm·s)) 1 1 Total(pcm) 1 1 Group 1 1 1 Group 2 1 1 Group 3 1 1 Group 4 1 1 Group 5 1 1 Group 6 1 1 -2.04 1 1 248.7 1 1 14.8 1 1 92.2 1 1 65.0 1 1 72.2 1 1 4.1 1 1 0.3 1 1 0 1 1 210.9 1 1 13.5 1 1 84.0 1 1 55.9 1 1 54.5 1 1 2.7 1 1 0.2 1 1 2.04 1 1 180.5 1 1 12.2 1 1 75.8 1 1 47.6 1 1 42.6 1 1 2.0 1 1 0.1 To accurately determine the influence of the flow field distribution, the DNF loss was normalized using the DNF loss under uniform flow conditions. As Fig. 22 Fig. 22 shows, the DNP loss of the first DNP group is approximately 10% affected by the flow field distribution, whereas it reaches 60% for the sixth DNP group. This is because when the fuel flow rate is near the nominal flow rate, the DNF loss of the first DNP group that changes with the flow rate is not obvious, while the DNF loss of the sixth DNP group continues to rise. The change in the flow field distribution is reflected in the change in the flow rate of each channel; therefore, the sixth DNP group is more significantly affected by the flow distribution. Flow field effect in reactivity insertion transient 2 In this section, the MSRE system responses with U-235 fuel under different initial conditions are calculated for the reactivity insertion transient. The transient starts with a nominal power of 8 MWth, and the reactivity of 26 pcm is inserted in 1 s so that the reactor power increases rapidly; however, the power increase caused by the reactivity jump is slowed down by the reactivity negative feedback caused by the temperature increase. First, the influence of different external loop recirculation times of the fuel salt on the system response in reactivity insertion transient was simulated and analyzed. Three cases with recirculation times of 5 s, 10 s, and 20 s were selected for this study, whereas the other parameters were the same. These power responses are shown in Fig. 23 Fig. 23 . The peak value of the power curve at the recirculation time of 20 s is approximately 20% higher than that of the normal case of 10 s, whereas the power of the 5 s case is lower because the DNPs decay in the external loop increase when the recirculation time increases, resulting in a lower effective delayed neutron fraction of the reactor core, as discussed in Sect. 3.1.2. DNF is an important parameter in reactor control. When introducing the same reactivity, the peak value in the process of the power response of the reactor with a lower effective delayed neutron fraction is higher. In Fig. 24 Fig. 24 , the responses of the average fuel salt and graphite temperatures are shown. At the beginning of the transient, the average fuel salt temperature increases owing to the increase in the power. When the hot molten salt flows through the external loop and re-enters the core, the average molten salt temperature further increases and reaches its highest value at approximately 40 s. In the 20 s recirculation time case, the initial temperature increases faster because of the faster power increase; however, owing to the longer external loop, the hot molten salt re-enters the core later, and the average molten salt temperature of the 20 s recirculation time case is overtaken by the other two cases. The curve fluctuates slightly because of the inflow and outflow of molten salt in the primary circuit. Finally, the average molten salt temperature and graphite temperature of the three cases are stable at the same level. Next, the influence of the fuel salt flow rate on the system response in the transient state was simulated and analyzed. Keeping the other parameters constant, three cases with mass flow rates of 93.75 kg/s, 187.5 kg/s, and 375 kg/s were selected. As shown in Fig. 25 Fig. 25 , the peak of the power curve is higher when the fuel salt flow is larger; this higher peak is caused by two factors. The first reason is that when the flow rate increases, the effective DNF decreases, and thus, the system response is more sensitive to the reactivity jump. Another reason is that, when the fuel flow rate is high, the fuel temperature increases slowly, as shown in Fig. 26 Fig. 26 ; thus, the reactivity negative feedback introduced is smaller. When the fuel flow rate is low, the heat could not be removed quickly, and the temperature of the molten salt increases extremely rapidly, which results in a rapid decrease in the core power, as shown by the red line in Fig. 25 Fig. 25 and 26 26 . Finally, the influence of the core fuel salt flow field distribution on the system response in reactivity insertion transient was calculated. Three fuel flow distribution cases, with a flow rate higher in the central area, a uniform flow rate in the whole core, and a flow rate higher in the peripheral area, were selected for analysis, corresponding to k =-2.04, 0, 2.04× 10 -3 kg/(cm·s) in Sect. 3.1.3. As shown in Fig. 27 Fig. 27 , the power peaks occur at approximately 20 s in all cases, and the peak value decreases with an increase in k . This is because when k is smaller, the DNP loss is larger, and the effective DNF of the reactor is lower, which leads to a higher power response, and the peak value of the molten salt temperature is also higher, as shown in Fig. 28 Fig. 28 . Finally, the temperature and power of the three transient conditions became stable basically at the same level. Conclusion 1 The effective DNF is an important parameter in the safety analysis of liquid-fueled MSRs, and the DNF loss is closely related to the state of the molten salt flow field in the reactor core and the external primary loop. In this study, to model and analyze the flow field effect of the delayed neutron precursor in channel-type liquid-fueled MSRs more systematically, a 3D dynamics code, ThorCORE3D, was developed and validated with MSRE benchmarks. Subsequently, the effects of external loop recirculation time, fuel flow rate, and core flow field distribution on DNF loss under steady-state conditions were simulated using the ThorCORE3D code. The results show that the external loop recirculation time, fuel salt flow rate, and distribution of the fuel salt flow field significantly affect the steady-state effective DNF loss of the liquid-fueled MSRs, and the effect on DNF of different DNP groups differ. Both the DNP out-of-core decay and the change in DNP core distribution could cause a loss of DNF in the core. Ultimately, the reactivity insertion transients of the MSRE under different external loop recirculation times, fuel flow rates, and core molten salt flow field distributions were simulated using the ThorCORE3D code, and the power and temperature responses were observed and analyzed in detail. The results show that the reactor system has different and instantaneous responses owing to different steady-state DNFs. Generally, the peak values of the reactor power and average fuel and graphite temperatures reach higher values when the initial DNF is smaller. Owing to the reactivity negative feedback effect, the power and temperature are eventually stabilized at the same level in different external loop recirculation times and different flow field distribution cases. When the fuel flow rates differ, the final powers and temperatures significantly differ owing to the combined effect of different DNF losses and temperature feedback. References D.L. Zhang , L.M. Liu , M.H. Liu et al , Review of conceptual design and fundamental research of molten salt reactors in China . Int. J. Energ. Res . 42 ( 5 ), 1834 ( 2018 ). doi: 10.1002/er.3979 http://doi.org/10.1002/er.3979 Review of conceptual design and fundamental research of molten salt reactors in China G. Lapenta , F. Mattioda , P. Ravetto , Point kinetic model for fluid fuel systems . Ann. Nucl. Energy . 28 , 1759 - 1772 ( 2001 ). doi: 10.1016/S0306-4549(01)00012-3 http://doi.org/10.1016/S0306-4549(01)00012-3 Point kinetic model for fluid fuel systems C.B. Shi , M.S. Cheng , G.M. Liu. Extending and verification of RELAP5 code for liquid fueled molten salt reactor . Nucl. Power Eng . 37 , 16 - 20 ( 2016 ). doi: 10.13832/j.jnpe.2016.03.0016 http://doi.org/10.13832/j.jnpe.2016.03.0016 .( in Chinese ) G.M. Liu. Extending and verification of RELAP5 code for liquid fueled molten salt reactor R. Li , M.S. Cheng , Z.M. Dai. Improvement and validation of the delayed neutron precursor transport model in RELAP5 code for liquid fuel molten salt reactor . Nucl. Tech . 44 ( 6 ): 060603 ( 2021 ). doi: 10.11889/j.0253-3219.2021.hjs.44.060603 http://doi.org/10.11889/j.0253-3219.2021.hjs.44.060603 ( in Chinese ) Z.M. Dai. Improvement and validation of the delayed neutron precursor transport model in RELAP5 code for liquid fuel molten salt reactor R.C. Diniz , F.S. de S. da Rosa , A. da C. Gonçalves , Calculation of delayed neutron precursors’ transit time in the external loop during a flow velocity transient in a Molten Salt Reactors . Ann. Nucl. Energy . 165 , 108640 ( 2022 ). doi: 10.1016/j.anucene.2021.108640 http://doi.org/10.1016/j.anucene.2021.108640 Calculation of delayed neutron precursors’ transit time in the external loop during a flow velocity transient in a Molten Salt Reactors G.F. Zhu , R. Yan , H.H. Pengu et al , Application of Monte Carlo method to calculate the effective delayed neutron fraction in molten salt reactor . Nucl. Sci. Tech . 30 , 34 ( 2019 ). doi: 10.1007/s41365-019-0557-7 http://doi.org/10.1007/s41365-019-0557-7 Application of Monte Carlo method to calculate the effective delayed neutron fraction in molten salt reactor M. Aufiero , M. Brovchenko , A. Cammi et al , Calculating the effective delayed neutron fraction in the Molten Salt Fast Reactor: Analytical, deterministic and Monte Carlo approaches . Ann. Nucl. Energy . 65 , 78 - 90 ( 2014 ). doi: 10.1016/j.anucene.2013.10.015 http://doi.org/10.1016/j.anucene.2013.10.015 Calculating the effective delayed neutron fraction in the Molten Salt Fast Reactor: Analytical, deterministic and Monte Carlo approaches J. Křepel , U. Rohde , U. Grundmann et al , DYN3D-MSR spatial dynamics code for molten salt reactors . Ann. Nucl. Energy . 34 , 449 - 462 ( 2007 ). doi: 10.1016/j.anucene.2006.12.011 http://doi.org/10.1016/j.anucene.2006.12.011 DYN3D-MSR spatial dynamics code for molten salt reactors K. Zhuang , L.Z. Cao , Y.Q. Zheng et al , Studies on the molten salt reactor: code development and neutronics analysis of MSRE-type design . J. Nucl. Sci. Technol . 52 , 251 - 263 ( 2015 ). doi: 10.1080/00223131.2014.944240 http://doi.org/10.1080/00223131.2014.944240 Studies on the molten salt reactor: code development and neutronics analysis of MSRE-type design L.Z. Cao , K. Zhuang , Y.Q. Zheng et al , Transient analysis for liquid-fuel molten salt reactor based on MOREL2.0 code . Int. J. Energ. Res . 42 , 261 - 275 ( 2018 ). doi: 10.1002/er.3828 http://doi.org/10.1002/er.3828 Transient analysis for liquid-fuel molten salt reactor based on MOREL2.0 code D.L. Zhang , Z.G. Zhai , X.N. Chen et al , COUPLE, A coupled neutronics and thermal- hydraulics code for transient analyses of Molten Salt Reactors . Transactions 108 ( 1 ) 921 - 922 ( 2013 ). COUPLE, A coupled neutronics and thermal- hydraulics code for transient analyses of Molten Salt Reactors D.L. Zhang , Z.G. Zhai , A. Rineiski et al , COUPLE, A time-dependent coupled neutronics and thermal-hydraulics code, and its application to MSFR . Proceedings of the 2014 22nd International Conference on Nuclear Engineering . https://doi.org/10.1115/ICONE22-30609 COUPLE, A time-dependent coupled neutronics and thermal-hydraulics code, and its application to MSFR D.L. Zhang , L.M. Liu , M.H. Liu et al , Neutronics/thermal-hydraulics coupling analysis for the liquid-fuel MOSART concept . Energy Procedia 127 , 343 - 351 ( 2017 ). doi: 10.1016/j.egypro.2017.08.075 http://doi.org/10.1016/j.egypro.2017.08.075 Neutronics/thermal-hydraulics coupling analysis for the liquid-fuel MOSART concept A. Cammi , V. Di Marcello , L. Luzzi et al ., A multi-physics modelling approach to the dynamics of Molten Salt Reactors . Ann. Nucl. Energy 38 , 1356 ( 2011 ). doi: 10.1016/j.anucene.2011.01.037 http://doi.org/10.1016/j.anucene.2011.01.037 A multi-physics modelling approach to the dynamics of Molten Salt Reactors M. Zanetti , A. Cammi , C. Fiorina et al ., A geometric multiscale modelling approach to the analysis of MSR plant dynamics . Prog. Nucl. Energ . 83 , 82 ( 2015 ). doi: 10.1016/j.pnucene.2015.02.014 http://doi.org/10.1016/j.pnucene.2015.02.014 A geometric multiscale modelling approach to the analysis of MSR plant dynamics C. Fiorina , D. Lathouwers , M. Aufiero et al ., Modelling and analysis of the MSFR transient behaviour, Ann . Nucl. Energy . 64 , 485 ( 2014 ). doi: 10.1016/j.anucene.2013.08.003 http://doi.org/10.1016/j.anucene.2013.08.003 Modelling and analysis of the MSFR transient behaviour, Ann E. Cervi , S. Lorenzi , L. Luzzi et al ., Multiphysics analysis of the MSFR helium bubbling system: A comparison between neutron diffusion, SP3 neutron transport and Monte Carlo approaches . Ann. Nucl. Energy . 132 , 227 ( 2019 ). doi: 10.1016/j.anucene.2019.04.029 http://doi.org/10.1016/j.anucene.2019.04.029 Multiphysics analysis of the MSFR helium bubbling system: A comparison between neutron diffusion, SP3 neutron transport and Monte Carlo approaches P. German , M. Tano , C. Fiorina et al ., GeN-ROM—An OpenFOAM®-based multiphysics reduced-order modeling framework for the analysis of Molten Salt Reactors . Prog. Nucl. Energ . 146 , 104148 ( 2022 ). doi: 10.1016/j.pnucene.2022.104148 http://doi.org/10.1016/j.pnucene.2022.104148 GeN-ROM—An OpenFOAM®-based multiphysics reduced-order modeling framework for the analysis of Molten Salt Reactors M. Tiberga , D. Lathouwers , and J.L. Kloosterman , A multi-physics solver for liquid-fueled fast systems based on the discontinuous Galerkin FEM discretization . Prog. Nucl. Energ . 127 , 103427 ( 2020 ). doi: 10.1016/j.pnucene.2020.103427 http://doi.org/10.1016/j.pnucene.2020.103427 A multi-physics solver for liquid-fueled fast systems based on the discontinuous Galerkin FEM discretization S. M. Park and M. Munk , Verification of moltres for multiphysics simulations of fast-spectrum molten salt reactors . Ann. Nucl. Energy 173 , 109111 ( 2022 ). doi: 10.1016/j.anucene.2022.109111 http://doi.org/10.1016/j.anucene.2022.109111 Verification of moltres for multiphysics simulations of fast-spectrum molten salt reactors G. Yang , M.K. Jaradat , W. Sik Yang et al ., Development of coupled PROTEUS-NODAL and SAM code system for multiphysics analysis of molten salt reactors . Ann. Nucl. Energy 165 , 108889 ( 2022 ). doi: 10.1016/j.anucene.2021.108889 http://doi.org/10.1016/j.anucene.2021.108889 Development of coupled PROTEUS-NODAL and SAM code system for multiphysics analysis of molten salt reactors D.L. Zhang , S.Z. Qiu , G.H. Su et al. , Development of a steady state analysis code for a molten salt reactor . Ann. Nucl. Energy 36 , 590 - 603 ( 2009 ). doi: 10.1016/j.anucene.2009.01.004 http://doi.org/10.1016/j.anucene.2009.01.004 Development of a steady state analysis code for a molten salt reactor X.D. Zuo , M.S. Cheng , Z.M. Dai , Development and validation of a three-dimensional dynamics code for liquid-fueled molten salt reactors . Nucl. Tech . 45 ( 3 ), 030603 ( 2022 ). doi: 10.11889/j.0253-3219.2022.hjs.45.030603 http://doi.org/10.11889/j.0253-3219.2022.hjs.45.030603 ( in Chinese ) Development and validation of a three-dimensional dynamics code for liquid-fueled molten salt reactors M.S. Cheng , M. Lin , X.D. Zuo et al ., Development and validation of a three-dimensional hexagonal nodal time-spatial kinetics code based on exponential transform . Nucl. Tech . 41 ( 6 ), 060604 ( 2018 ). doi: 10.11889/j.0253-3219.2018.hjs.41.060604 http://doi.org/10.11889/j.0253-3219.2018.hjs.41.060604 . ( in Chinese ) Development and validation of a three-dimensional hexagonal nodal time-spatial kinetics code based on exponential transform D.L. Zhang , S.Z. Qiu , C.L. Liu et al ., Nuclear calculation and program development for Molten Salt Reactor . At. Energy Sci. Technol . 42 ( 12 ), 1103 - 1108 ( 2008 ). doi: 10.7538/yzk.2008.42.12.1103 http://doi.org/10.7538/yzk.2008.42.12.1103 . ( in Chinese ) Nuclear calculation and program development for Molten Salt Reactor G. Marleau , A. Hébert , R. Roy , A user guide for dragon version 4 , Institute of Genius Nuclear, Department of Genius Mechanical, School Polytechnic of Montreal ( 2011 ). A user guide for dragon Bilinear interpolation , https://en.wikipedia.org/wiki/Bilinear_interpolation https://en.wikipedia.org/wiki/Bilinear_interpolation ; 2022 . U. Grundmann , U. Rohde , S. Mittag et al. , DYN3D version 3.2-code for calculation of transients in light water reactors (LWR) with hexagonal or quadratic fuel elements-description of models and methods, Forschungszentrum Rossendorf ( 2005 ). DYN3D version 3.2-code for calculation of transients in light water reactors (LWR) with hexagonal or quadratic fuel elements-description of models and methods, Forschungszentrum Rossendorf B.E. Prince , S.J. Ball , J.R. Engel et al. , Zero-power physics experiments on the molten-salt reactor experiment . Oak Ridge Natl Lab ( 1968 ). doi: 10.2172/4558029 http://doi.org/10.2172/4558029 Zero-power physics experiments on the molten-salt reactor experiment M.W. Rosenthal , R.B. Briggs , P.R. Kasten , Molten-Salt reactor program semiannual progress report for period ending February 28, 1969 . Oak Ridge Natl Lab ( 1969 ). doi: 10.2172/4780471 http://doi.org/10.2172/4780471 Molten-Salt reactor program semiannual progress report for period ending R.B. Briggs , Molten-salt reactor program semiannual progress report for period ending July 31, 1964 , Oak Ridge Natl Lab ( 1964 ). doi: 10.2172/4676587 http://doi.org/10.2172/4676587 Molten-salt reactor program semiannual progress report for period ending R.C. StefYy , P.J. Wood , Theoretical dynamic anatysis of the MSRE with 233-U fuel , Oak Ridge Natl Lab ( 1969 ). doi: 10.2172/4771215 http://doi.org/10.2172/4771215 Theoretical dynamic anatysis of the MSRE with 233-U fuel M. Delpech , S. Dulla , C. Garzenne et al , Benchmark of dynamic simulation tools for molten salt reactors , in: GLOBAL 2003-Nuclear Science and Technology: Meeting the Global Industrial and R&D Challenges of the 21st Century , American Nuclear Society . pp. 2182 - 2187 ( 2003 ). Benchmark of dynamic simulation tools for molten salt reactors R.C. Steffy , Experimental dynamic analysis of the MSRE with 233-U fuel , Oak Ridge Natl Lab ( 1969 ). doi: 10.2172/4132458 http://doi.org/10.2172/4132458 Experimental dynamic analysis of the MSRE with 233-U fuel All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xian-Di Zuo, Mao-Song Cheng, Yu-Qing Dai, and Kai-Cheng Yu. The first draft of the manuscript was written by Xian-Di Zuo and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01084-0 This work was supported by Strategic Pilot Science and Technology Project of Chinese Academy of Sciences (No. XD02001005). 2022-04-10 2022-05-30 2022-06-26 2022-08-12 2022-09-26 © 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. 2022 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-022-01084-0 http://dx.doi.org/10.1007/s41365-022-01084-0 . Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2022 Nuclear Science and Techniques 08 33