Development of a dynamics model for graphite-moderated channel-type molten salt reactor

Development of a dynamics model for graphite-moderated channel-type molten salt reactor 1 2 3 first-author 1 2 1 2 3 1 2 3 corresp guowei@sinap.ac.cn (Wei Guo). guowei@sinap.ac.cn (Wei Guo). guowei@sinap.ac.cn 1 2 1 2 3 corresp Corresponding author: caixz@sinap.ac.cn (Xiang-Zhou Cai) Corresponding author: caixz@sinap.ac.cn (Xiang-Zhou Cai) caixz@sinap.ac.cn Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China. Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China. CAS Innovative Academy in TMSR Energy System, Chinese Academy of Sciences, Shanghai 201800, China CAS Innovative Academy in TMSR Energy System, 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. A molten salt reactor (MSR) is one of the six advanced reactor concepts selected by the Generation IV International Forum because of its advantages of inherent safety, and the promising capabilities of Th-U breeding and transuranics (TRU) transmutation. A dynamics model for the channel-type MSR is developed in this work based on a three-dimensional thermal-hydraulic model (3DTH) and a point reactor model. The 3DTH couples a three-dimensional heat conduction model and a one-dimensional single-phase flow model that can accurately consider the heat conduction between different assemblies. The 3DTH is validated by the RELAP5 code in terms of the temperature and mass flow distribution calculation. A point reactor model considering the drift of delayed neutron precursors (DNP) is adopted in the dynamics model. To verify the dynamics model, three experiments from the Molten Salt Reactor Experiment (MSRE) are simulated. The agreement of the experimental data and simulation results was excellent. With the aid of this model, the unprotected step reactivity addition and unprotected loss of flow of the 2 MWt experimental MSR are modeled, and the reactor power and temperature evolution are analyzed. Molten salt reactor Thermal-hydraulics Point reactor model Thermal coupling 1. Introduction 1 The history of molten salt reactors (MSRs) can be traced back to the 1940s when the USA started a military project, Aircraft Nuclear Propulsion [ 1 1 ]. In the 1950s, the first MSR Aircraft Reactor Experiment (ARE) [ 2 2 ] was built at Oak Ridge National Laboratory (ORNL). The research on MSRs continued in ORNL with the building and operation of the Molten-Salt Reactor Experiment (MSRE) [ 3 3 ]. The successful operation of the MSRE established the basic technologies and led to the conceptual design of a molten-salt breeder reactor (MSBR) [ 4 4 ] at ORNL in the 1970s, which was never built. Since the 1980s, a series of researches on MSRs for Th-U breeding or transuranics (TRU) transmutation have been undertaken in France [ 5 5 , 6 6 ], Japan [ 7 7 ], and Russia [ 8 8 ]. In 2002, the MSR was selected by the Generation IV International Forum as a candidate for the Generation IV Reactor because it meets the criteria of sustainability, economics, safety and reliability, proliferation resistance, and physical protection. In 2011, the Chinese Academy of Sciences (CAS) launched the Thorium-based MSR (TMSR) project with the aim of building a series of MSRs to utilize the thorium resource efficiently [ 9 9 ]. The graphite-moderated channel-type MSR is a typical design of the MSR that adopts a liquid form of fuel. The fuel salt mixing carrier salt and fissile material is used as fuel and coolant for this type of MSR. The reactor core of this type of MSR consists of a large number of graphite stringers that form the fuel channels. The fission reaction occurs when the fuel salt flows through the fuel channel. The graphite-moderated channel-type MSR is considerably different from reactors using solid fuels. First, the fission energy is released predominantly in the fuel salt and removed directly by the flowing fuel salt. The graphite moderator only deposits a small part of the energy produced by gamma and fast neutron radiation. Second, as there is no transverse flow between different fuel channels, the thermal coupling of the fuel salt in this type of MSR is achieved through heat conduction in the graphite moderator. Finally, owing to the liquidity of the fuel salt, a part of the delayed neutron precursors (DNPs) drifts out of the core and emits delayed neutrons (DNs) at the external loop. Therefore, the transient behavior of the MSR requires consideration of the drift effect of the DNPs. Considering the aforementioned characteristics, these dynamics code developed for solid fuel reactors are not suitable for MSR. In recent decades, a series of dynamic codes have been developed for the graphite-moderated channel-type MSR. Křepel et al. developed the DYN3D-MSR [ 10 10 ] code, where the three-dimensional (3D) neutronics and the multichannel thermal hydraulic model were used. Kópházi developed a 3D analysis code DT-MSR [ 11 11 ], where the neutron diffusion model was adopted for neutronics calculations, and the one-dimensional heat convection equation and 3D heat conduction equation were adopted for thermal-hydraulics calculation. Based on the neutron diffusion model and multichannel model, Si et al. and Wei et al. developed the dynamics code TANG-MSR [ 12 12 ] and a two-dimensional dynamics code [ 13 13 ], respectively. Although many studies have been performed, the majority of these codes only adopt the multichannel model for the thermal hydraulics calculation, and some of these codes must assume velocity profiles for the temperature and neutronics calculation. The multichannel model adopts the one-dimensional heat conduction equation to calculate the moderator temperature, and the one-dimensional single-phase flow model for the fuel salt, which neglects the thermal coupling between the assemblies. In the present study, a 3D thermal-hydraulic model (3DTH) that can consider the thermal coupling between different assemblies is developed to analyze the transient features of MSRs. The dynamics model is developed by coupling the 3DTH with a point reactor model. To validate the dynamics model, the protected fuel pump start-up experiment, protected coast-down experiment, and natural circulation experiment of the MSRE are simulated. By applying the dynamics model, the reactor power and temperature evolutions of the 2 MWt experimental MSR (2MW-MSR) [ 14 14 ] are analyzed for several transients. 2. Numerical Model 1 2.1 Thermal hydraulics model 2 In the 3DTH, the temperature of the solid region, which includes the graphite assemblies and reflector, is calculated by a 3D heat conduction equation: <math display="block" id="M1"><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>C</mi><mtext>g</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>T</mi><mtext>g</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mo>∇</mo><msub><mi>λ</mi><mtext>g</mtext></msub><mo>∇</mo><msub><mi>T</mi><mtext>g</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>Q</mi><mtext>g</mtext></msub><mtext> </mtext><mo>,</mo></mrow></math> where ρ g , C g , λ g , T g , and Q g are the density, thermal capacity, thermal conductivity, temperature, and power density, respectively. The exterior surface of the solid region and the inner surface of the fuel channel employ the adiabatic boundary condition and convective heat transfer boundary condition, respectively. Different models are adopted for the fuel salt in different regions. The fuel channels are described individually by the one-dimensional single-phase flow model. The governing equations are the mass, momentum, and energy conservation equations: <math display="block" id="M2"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>f</mi></msub><mi>v</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mn>0</mn><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M3"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mi>v</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mtext>+</mtext><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><msup><mi>v</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mtext>-</mtext><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><mo>∂</mo><msub><mi>p</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><msub><mi>ρ</mi><mtext>f</mtext></msub><mi>g</mi><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M4"><mrow><mi>A</mi><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><msub><mi>C</mi><mtext>f</mtext></msub><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mtext>+</mtext><mi>A</mi><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mi>v</mi><msub><mi>C</mi><mtext>f</mtext></msub><msub><mi>T</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mi>A</mi><msub><mi>Q</mi><mtext>f</mtext></msub><mo>+</mo><msub><mi>Q</mi><mtext>h</mtext></msub><mtext> </mtext><mo>,</mo></mrow></math> where ρ f , v , C f , and T f are the fuel density, velocity, thermal capacity, and temperature, respectively, and A is the cross-section area of the fuel channel. The last two terms in the right of Eq. (3) denote the frictional pressure drop and gravity pressure drop, respectively; Q f represents the volumetric heat released in the fuel salt; and Q h denotes the heat flux at the surface of the fuel channel. The units of Q f and Q h are W/m 3 and W/m, respectively. The temperature of the fuel salt in the plenums is calculated by the lumped parameter model; the equation can be written as <math display="block" id="M5"><mrow><mi>M</mi><msub><mi>C</mi><mtext>f</mtext></msub><mfrac><mrow><mtext>d</mtext><mi>T</mi></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><msub><mi>m</mi><mtext>f</mtext></msub><msub><mi>C</mi><mtext>f</mtext></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mrow><mtext>in</mtext></mrow></msub><mo>−</mo><mi>T</mi><mo stretchy="false">)</mo><mo>+</mo><mi>Q</mi><mtext> ,</mtext></mrow></math> where M , m f , and Q are the mass of the fuel salt in the plenum, mass flow of the reactor, and energy deposited in the plenum, respectively; T in is the inlet temperature of the plenum. For the bottom plenum, T in is the inlet temperature of the reactor. For the top plenum, T in is calculated by the following formula: <math display="block" id="M6"><mrow><msub><mi>T</mi><mrow><mtext>in</mtext></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><msub><mi>m</mi><mi>i</mi></msub><msub><mi>C</mi><mtext>f</mtext></msub><msub><mi>T</mi><mrow><mtext>out</mtext><mo>,</mo><mi>i</mi></mrow></msub></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mi>i</mi></munder><mrow><msub><mi>m</mi><mi>i</mi></msub><msub><mi>C</mi><mtext>f</mtext></msub></mrow></mstyle></mrow></mfrac><mo>,</mo></mrow></math> where m i and T out ,i are the mass flow and outlet fuel salt temperature of channel i , respectively The thermal coupling between each fuel channel and solid region is built with the convective heat transfer boundary condition at the inner surface of the fuel channel. The heat flux between the fuel salt and surface of the fuel channel has the following form: <math display="block" id="M7"><mrow><mi>q</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>N</mi><mi>u</mi><mo>⋅</mo><msub><mi>λ</mi><mtext>f</mtext></msub></mrow><mrow><msub><mi>D</mi><mtext>e</mtext></msub></mrow></mfrac><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>g</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>z</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>T</mi><mtext>f</mtext></msub><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where Nu , λ f , and D e are the Nusselt number, thermal conductivity of fuel salt, and diameter of the fuel channel, respectively. The Nusselt number is calculated by the empirical correlation from the forced-convection heat-transfer experiment in ORNL [ 15 15 ]. The Q h in Eq. (4) is calculated by integrating q ( x,y,z ) along the perimeter of the fuel channel. The line integral of momentum conservation equation (Eq. (3)) along the z direction of the fuel channel can be written as <math display="block" id="M8"><mrow><mfrac><mi>L</mi><mi>A</mi></mfrac><mfrac><mrow><msubsup><mi>m</mi><mi>i</mi><mi>n</mi></msubsup><mo>−</mo><msubsup><mi>m</mi><mi>i</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>a</mtext><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mi>n</mi></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>f</mtext><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>g</mtext><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>k</mtext><mo>,</mo><mi>i</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>,</mo></mrow></math> where L is the length of the fuel channel; m i is the mass flow in channel i , which equals to <math id="M9"><mrow><msub><mi>ρ</mi><mtext>f</mtext></msub><mi>v</mi><mi>A</mi></mrow></math> ; and Δ P a ,i , Δ P i , Δ P f ,i , Δ P g ,i , and Δ P k ,i are the acceleration pressure drop, total pressure drop, frictional pressure drop, gravity pressure drop, and local pressure drop (additional term), respectively. The indexes n and i indicate the time step and channel number, respectively. These pressure drops with indexes ( n – 1) are known quantities, which are calculated by the following equations: <math display="block" id="M10"><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>a</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><msub><mi>A</mi><mi>i</mi></msub><msup><mi>v</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><msub><mi>m</mi><mi>i</mi></msub><mi>v</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mtext> </mtext><mtext> </mtext><mo>=</mo><mfrac><mrow><msub><mi>m</mi><mi>i</mi></msub></mrow><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mi>v</mi><mrow><mtext>out</mtext></mrow></msub><mo>−</mo><msub><mi>v</mi><mrow><mtext>in</mtext></mrow></msub></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M11"><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>f</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mi>f</mi><mfrac><mrow><mi>ρ</mi><msup><mi>v</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><msub><mi>D</mi><mtext>e</mtext></msub></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mi>f</mi><mfrac><mrow><msub><mi>m</mi><mi>i</mi></msub><mi>v</mi></mrow><mrow><mn>2</mn><msub><mi>D</mi><mtext>e</mtext></msub></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mo>=</mo><mfrac><mrow><msub><mi>m</mi><mi>i</mi></msub></mrow><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mi>f</mi><mfrac><mi>v</mi><mrow><mn>2</mn><msub><mi>D</mi><mtext>e</mtext></msub></mrow></mfrac></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mo>,</mo></mrow></math> <math display="block" id="M12"><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>g</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∫</mo><mrow><mi>ρ</mi><mi>g</mi></mrow></mrow></mstyle><mtext>d</mtext><mi>z</mi><mo>,</mo></mrow></math> <math display="block" id="M13"><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>k</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mi>k</mi><mfrac><mrow><mi>ρ</mi><msup><mi>v</mi><mn>2</mn></msup></mrow><mn>2</mn></mfrac><mo>=</mo><mi>k</mi><mfrac><mrow><msub><mi>m</mi><mi>i</mi></msub></mrow><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></mfrac><mfrac><mi>v</mi><mn>2</mn></mfrac><mo>,</mo></mrow></math> The mass flow, velocity, and density in Eqs. (9)–Eq. (12) are known quantities, which are calculated in the former time step ( n – 1). The total pressure drop Δ P a ,i of the different fuel channels should be equal, that is Δ P i = Δ P . For a reactor with N fuel channels, the mass flow in each fuel channel is calculated by [ 16 16 ]: <math display="block" id="M14"><mrow><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mfrac><mi>L</mi><mi>A</mi></mfrac><mfrac><mrow><msubsup><mi>m</mi><mn>1</mn><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>m</mi><mn>1</mn><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mi>Δ</mi><msup><mi>P</mi><mi>n</mi></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>f</mtext><mo>,</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mi mathvariant="normal">g</mi><mo>,</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>k</mtext><mo>,</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>a</mtext><mo>,</mo><mn>1</mn></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>⋮</mo></mtd></mtr><mtr><mtd><mfrac><mi>L</mi><mi>A</mi></mfrac><mfrac><mrow><msubsup><mi>m</mi><mi>N</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>m</mi><mi>N</mi><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mi>Δ</mi><msup><mi>P</mi><mi>n</mi></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>f</mtext><mo>,</mo><mi>N</mi></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>g</mtext><mo>,</mo><mi>N</mi></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>k</mtext><mo>,</mo><mi>N</mi></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>−</mo><msup><mrow><mo>(</mo><mrow><mi>Δ</mi><msub><mi>P</mi><mrow><mtext>a</mtext><mo>,</mo><mi>N</mi></mrow></msub></mrow><mo>)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msup><mtext> </mtext></mtd></mtr><mtr><mtd><mfrac><mrow><mo>∂</mo><msub><mi>M</mi><mrow><mtext>total</mtext></mrow></msub></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>i</mi><mi>N</mi></munderover><mrow><mfrac><mrow><msubsup><mi>m</mi><mi>i</mi><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>−</mo><msubsup><mi>m</mi><mi>i</mi><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac></mrow></mstyle></mtd></mtr></mtable></mrow><mo>,</mo></mrow></math> where M total is the total mass flow in the reactor core, which is a known quantity. Eq. (13) is a system of linear equations with ( N + 1) variables. The unknown variables ( m i and Δ P ) for the time step n are directly calculated by solving the linear equations. The 3DTH is developed based on COMSOL Multiphysics [ 17 17 ] and MATLAB software. COMSOL Multiphysics is a general-purpose platform software for modeling multiphysics applications. The heat conduction and one-dimensional single-phase flow are calculated by the computational fluid dynamics (CFD) module of COMSOL Multiphysics. With the aid of the module "LiveLink for Matlab" provided by COMSOL Multiphysics, MATLAB script was used to build the geometry and thermal coupling between the solid region and fuel salt. 2.2 Neutronics model 2 In this study, a new point reactor model is developed for neutron dynamics. To simulate the drift of DNPs along the primary loop, the reactor core and external loop are divided into several nodes. The schematic representation of the point reactor model is displayed in Fig. 1 Fig. 1 . The point reactor model can be written as <math display="block" id="M15"><mrow><mfrac><mrow><mtext>d</mtext><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><mi>ρ</mi><mo>−</mo><mi>β</mi></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>V</mi><mtext>c</mtext></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover><mrow><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>V</mi><mi>k</mi></msub></mrow></mstyle></mrow><mo>)</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> <math display="block" id="M16"><mrow><msub><mi>V</mi><mi>k</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><msub><mi>V</mi><mi>k</mi></msub><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>V</mi><mi>k</mi></msub><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>L</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>L</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>L</mi></msub><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>k</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>k</mi></msub><mtext> </mtext><mtext> </mtext><mi>k</mi><mo>=</mo><mn>1</mn><mo>,</mo></mrow></math> <math display="block" id="M17"><mrow><msub><mi>V</mi><mi>k</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><msub><mi>V</mi><mi>k</mi></msub><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>V</mi><mi>k</mi></msub><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mrow><mi>k</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>k</mi><mo>−</mo><mn>1</mn><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mrow><mi>k</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>k</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>k</mi></msub><mtext> </mtext><mtext> </mtext><mi>k</mi><mo>=</mo><mn>2</mn><mo>⋅</mo><mo>⋅</mo><mo>⋅</mo><mi>M</mi><mo>,</mo></mrow></math> <math display="block" id="M18"><mrow><msub><mi>V</mi><mi>l</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>V</mi><mi>l</mi></msub><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>M</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mtext>M</mtext><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>M</mi></msub><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>l</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>l</mi></msub><mtext> </mtext><mtext> </mtext><mi>l</mi><mo>=</mo><mn>1</mn><mo>,</mo></mrow></math> <math display="block" id="M19"><mrow><msub><mi>V</mi><mi>l</mi></msub><mfrac><mrow><mtext>d</mtext><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>=</mo><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><msub><mi>V</mi><mi>l</mi></msub><mo>+</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mrow><mi>l</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>l</mi><mo>−</mo><mn>1</mn><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mrow><mi>l</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>−</mo><mfrac><mn>1</mn><mrow><msub><mi>τ</mi><mi>l</mi></msub></mrow></mfrac><msub><mi>C</mi><mrow><mi>l</mi><mo>,</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mi>V</mi><mi>l</mi></msub><mtext> </mtext><mtext> </mtext><mi>l</mi><mo>=</mo><mn>2</mn><mo>⋅</mo><mo>⋅</mo><mo>⋅</mo><mi>L</mi><mo>,</mo></mrow></math> where n is the neutron density; C k,i is the DNP concentration of the i th group in node k of the reactor core; C l,i is the DNP concentration of the i th group in node l of the external loop; ρ is the total reactivity, which contains the compensative reactivity due to the loss of DNs, reactivity feedback caused by the temperature changes, and reactivity input; β is the total DN fraction; V is the fuel salt volume; and τ is the fuel transit time in a node. The subscripts c, k , and l denote the reactor core, node number in the core, and node number in the external loop, respectively. The last two terms in the right of Eqs. (15)–Eq. (18) describe the number of DNPs leaving and entering the node. The point reactor model is a set of nonlinear ordinary differential equations, which is solved by the Mathematics module of COMSOL Multiphysics. 2.3 Coupling scheme 2 The schematic overview of the dynamics model is presented in Fig. 2 Fig. 2 . The 3DTH and point reactor model are coupled by the reactor power and temperature. In the initialization step, the geometry for the 3DTH is built, the power distribution is read, and a steady state is calculated at a constant power. The 3DTH and neutronics are implicitly coupled; the convergence criteria must be satisfied before the calculation of the next time step. The dynamics model developed in this study is based on MATLAB software and COMSOL Multiphysics. MATLAB is used to build the geometry and create the coupling variables. COMSOL Multiphysics is used as a solver. The equations mentioned above are calculated by a specialized module and the Mathematics module of COMSOL Multiphysics. COMSOL Multiphysics provides a segregated solver and coupled solver for multiphysics simulation. The segregated solver divides the multiphysics problem into several equations, where each equation is solved sequentially. As an alternative to the segregated solver, the coupled solver can solve all the physical fields (equations) simultaneously. To reduce the computational time, the coupled solver is adopted in the present study. 3. Validation 1 3.1 Validation of the 3DTH 2 At present, there is no commercial software adopting the same model as the 3DTH. To verify the correctness of the 3DTH, the results from the 3DTH are compared with those calculated by the RELAP5 code [ 18 18 ]. To validate the temperature calculated by the 3DTH, a round pipe with a thin tube wall was chosen as the calculation model. The round pipe was cooled by the molten salt FLiBe. The main parameters are listed in Table 1 Table 1 . Because the axial heat conduction can be neglected, the 3D temperature field of the tube wall calculated by the 3DTH can be accurately validated by calculating the one-dimensional heat conduction in the radial direction. 1 1 Parameter 1 1 Value 1 1 Parameter 1 1 Value 1 1 Pipe length (m) 1 1 0.5 1 1 Thermal conductivity of tube wall (W/(m·K)) 1 1 25.0 1 1 Inner diameter (m) 1 1 0.02 1 1 Total power (W) 1 1 59650.0 1 1 Tube wall thickness (m) 1 1 0.01 1 1 Fraction of power in the coolant (%) 1 1 90.0 1 1 Mass flow (kg/s) 1 1 0.5 1 1 Fraction of power in tube wall (%) 1 1 10.0 1 1 Inlet temperature (K) 1 1 873.15 1 1 Initial temperature (K) 1 1 873.15 The radial temperature distributions in the tube wall at the middle of the pipe ( z = 25.0 cm) calculated by the 3DTH and RELAP5 code are compared in Fig. 3a Fig. 3a . The temperature profiles at different times indicate a reasonably acceptable agreement between the 3DTH and RELAP5 code. The coolant temperature evolution at z = 12.5 cm and 50.0 cm are compared in Fig. 3b Fig. 3b . The results calculated by the 3DTH are consistent with the RELAP5 code. To validate the mass flow distribution calculated by the 3DTH, three identical vertical round pipes, which were cooled by FLiBe, were chosen as the calculation model. The size of these round pipes was identical to the size listed in Table 1 Table 1 . At the beginning ( t = 0 s), the total mass flow was 0.075 kg/s (0.025 kg/s for each pipe) and the energy deposited in each pipe was 3.579 kW. From t = 20 s to t = 30 s, the energy in Pipe 2 linearly increased by 20%, and the energy in Pipe 3 linearly decreased by 20%. From t = 100 s to t = 110 s, the total mass flow linearly increased by 50%. Because the density of FLiBe decreases with a rise in temperature, more mass flow is distributed to the pipe deposited with more energy. Figure 4a Figure 4a indicates the mass flow response in each pipe to the change of energy. Figure 4b Figure 4b represents the mass flow response in each pipe to the change of total mass flow. The mass flow distributions calculated by the 3DTH are consistent with the RELAP5 code. 3.2 Validation of the dynamics model 2 In this work, three MSRE experiments, the protected fuel pump start-up [ 19 19 ], protected fuel pump coast-down [ 19 19 ], and natural circulation experiment [ 20 20 , 21 21 ], were chosen as benchmarks to validate the dynamics model. The protected fuel pump start-up and coast-down experiments were performed in the U235-fueled MSRE core. ORNL DN data [ 22 22 ] was adopted for these simulations. The protected fuel pump start-up and protected coast-down experiments focused on calculating the reactivity, which compensates for the loss of DNs due to the change of velocity of the fuel salt. During these pump-driven transients, a constant power was maintained. In the case of the fuel pump start-up transient, the mass flow changed from zero to the nominal value in 10 s. The DNPs started to drift out of the core and positive reactivity was inserted to compensate for the loss of DNs. In the case of the pump coast-down, the number of DNPs that drifted out of the core reduced, and negative reactivity was induced. The compensative reactivity calculated by the present work was compared with that from the experiment and is shown in Fig. 5 Fig. 5 . The different curves indicate that the reactor core and external loop were divided into a different number of nodes. The C3L6 curve represents the reactor core divided into three nodes and external loop divided into six nodes. As indicated in Fig. 5 Fig. 5 , the shapes and trend of these curves calculated in this study are qualitatively consistent with the measured data. The natural circulation experiment was conducted by ORNL to determine the heat removal characteristics of MSRE by natural-convection flow of the fuel salt. At the beginning of this transient, the reactor power was stabilized at 4.1 kW with limited fuel circulation. The transient was driven by increasing the cooling capacity of the secondary circuit stepwise. During the experiment, reactor power, primary and secondary salt temperature, and primary mass flow rate were measured. The purpose of the benchmark was not to model natural circulation, but to reproduce the evolution of the reactor power. For this benchmark, the reactor inlet temperature and mass flow rate extracted from experimental data were provided as input data. The natural circulation experiment was performed only for the U233-fueled MSRE core. For this simulation, the neutron generation life Λ = 4.0333 <math id="M20"><mo>×</mo></math> 10 -4 s, fuel and graphite temperature feedback coefficient were equal to -9.54 pcm/K and -5.76 pcm/K, respectively, and ORNL DN data [ 22 22 , 23 23 ] was adopted. Figure 6 Figure 6 displays a comparison of this work and the experimental data. The calculated results are extremely consistent with the experiment data. In the 242nd minute, the maximum peak value of 358 kW, calculated by this model, is extremely close to the 354 kW measured by the experiment. 4. 2MW-MSR transient calculation 1 4.1 2MW-MSR description 2 After the validation of the dynamics model, it was applied to analyze several transients of the 2MW-MSR designed by the TMSR center of CAS [ 14 14 ]. A vertical and horizontal view of the 2MW-MSR is presented in Fig. 7a Fig. 7a and Fig. 7b Fig. 7b , respectively. The active core of the 2MW-MSR contains 85 fuel channels. Both the height and diameter of the active core are 1.1 m. The half-pitch of the assembly is 5.5 cm, and the radius of the fuel channel is 2.4 cm, as indicated in Fig. 7c Fig. 7c . In the 2MW-MSR, graphite is chosen as the moderator and reflector, and Hastelloy-N alloy (safety limit is 730 ℃ [ 24 24 ]) is used as the reactor vessel material and other structural material. The inlet and outlet temperatures of the core are 600 ℃and 620 ℃, respectively. The rated mass flow is designed as 59.25 kg/s. In the 2MW-MSR, a small portion of the fuel salt (the outer salt layer in Fig. 7 Fig. 7 ) flowing upward between the reflector and the reactor vessel, is used to cool the structure material. Because of the geometric symmetry of the 2MW-MSR, one-twelfth of the core is chosen for the calculation. Further, only the active core, top and bottom plenums, radial reflector, radial outer salt layer, and reactor vessel are considered in the 3DTH; the top and bottom reflectors are ignored. The cross-sectional view of the calculation model is displayed in Fig. 7d Fig. 7d , and each fuel assembly is numbered based on their radial position. For the thermal-hydraulic calculation, the mass flow in each fuel assembly is calculated based on equal pressure drop, and the mass flow in the outer salt layer is set as 5% of the rated flow. For the 2MW-MSR, owing to the lack of hydraulic experiment, the local pressure drops at the inlet and outlet of each fuel channel were ignored for the present simulation. The power distribution of the 2MW-MSR is calculated by SCALE [ 25 25 ]; the fraction of power deposited in each component is listed in the Table 2 Table 2 . It should be noted that the power distribution calculated by SCALE does not consider the drift of the DNPs and decay heat. Because the fraction of DNP is relatively small and the decay heat accounts for only approximately 7% of the reactor power, the drift of the DNPs and decay heat clearly does not influence the power distribution. For the neutronics calculation, the reactor core and external loop were divided into three and six nodes, respectively. The temperature feedbacks of the following regions: fuel salt of the assembly, graphite of the assembly, top plenum, bottom plenum, and radial reflector, were considered. The temperature coefficients of these regions were -3.820 pcm/K, -4.465 pcm/K, -0.496 pcm/K, -0.447 pcm/K, and -0.129 pcm/K, respectively. Table 3 Table 3 lists the DN data for this simulation. The data mentioned above for the point reactor model was calculated by Monte Carlo N-Particle Transport Code (MCNP) [ 26 26 ] with the continuous energy cross-section library of ENDF/B-VII. 1 1 Component 1 1 Power fraction (%) 1 1 Component 1 1 Power fraction (%) 1 1 Fuel salt of the assembly 1 1 65.68 1 1 Outer salt layer 1 1 4.41 1 1 Graphite of the assembly 1 1 3.39 1 1 Reactor vessel 1 1 0.34 1 1 Top plenum and salt channel 1 1 12.50 1 1 Radial reflector 1 1 1.14 1 1 Bottom plenum and salt channel 1 1 11.25 1 1 Top and bottom reflector 1 1 1.29 1 1 Group 1 1 Decay constant QUOTE λ i (1/s) 1 1 QUOTE β i (×10 -5 ) 1 1 1 1 1 0.0125 1 1 17.07 1 1 2 1 1 0.0319 1 1 122.35 1 1 3 1 1 0.1094 1 1 81.23 1 1 4 1 1 0.3172 1 1 340.97 1 1 5 1 1 1.3538 1 1 91.12 1 1 6 1 1 8.6643 1 1 17.63 4.2 Unprotected step reactivity insertion 2 The transients driven by step reactivity insertion were simulated at nominal power. The transients of 20 pcm, 50 pcm, and 100 pcm reactivity inserted at t = 1.0 s were calculated. The inlet temperature and mass flow were maintained at the nominal value during these transient. The power responses to the reactivity steps are reported in Fig. 8 Fig. 8 . The power increased rapidly at the beginning of these transients. Later, owing to the rapid temperature feedback of the fuel salt, the reactor power began to decrease. Because the graphite moderator (graphite assembly and reflector) has large heat inertia, the negative temperature feedback of the moderator was slow, which caused the slow decrease of the reactor power. The power peaks for 20 pcm, 50 pcm, and 100 pcm reactivity insertions were 1.32, 1.82, and 2.71 times the nominal power, respectively. After the inserted reactivity were compensated by the increases of temperature, the reactor powers tended to be stable and finally decreased to 1.13, 1.34, and 1.68 times the nominal power for 20 pcm, 50 pcm, and 100 pcm reactivity addition, respectively. The mass flow in each fuel channel during the 100 pcm reactivity insertion is displayed in Fig. 9a Fig. 9a . At the initial steady state, the peak of the power density is located at Assembly 1, and the lowest value of the power density is located at Assemblies 7–9. Therefore, Assembly 1 is distributed with the greatest mass flow, and the mass flow distributed to Assemblies 7–9 is the least. Because some neutrons are reflected and moderated by the radial reflector, a small power peak exists in Assemblies 10–11. The values of the mass flow in Assemblies 10–11 are between those in Assemblies 1–2 and those in Assemblies 7–9. With the increases of reactor power, the mass flow in the assembly located in the central region increased, and the mass flow in Assemblies 7–9 decreased. Then, the differences of the mass flow between these assemblies increased. Figure 9b Figure 9b displays the average temperature in different regions of the core and the maximum temperature in Assembly 1 during the 100 pcm reactivity insertion. As the majority of the fission energy is released directly into the salt, the fuel salt temperature rises more quickly than the moderator temperature. The maximum temperature of the fuel salt and moderator in Assembly 1 was 645.3 ℃ and 654.5 ℃, respectively. During this transient, the maximum fuel temperature was located at the top plenum and attained 651.5 ℃, which is below the safety limit of the structural material. The temperature fields at several time points ( t = 0 s, 20 s, 100 s, 500 s, 1000 s, and 2000 s) during the 100 pcm reactivity insertion transient are displayed in Fig. 10 Fig. 10 . At the initial steady state ( t = 0 s, Fig. 10a Fig. 10a ), because energy deposited in the moderator is removed by the fuel salt, the moderator temperature of the assembly is greater than the salt temperature. When the reactor power increases quickly, the fuel salt temperature in the region close to the outlet is greater than the graphite temperature, as displayed in Fig. 10b Fig. 10b and Fig. 10c Fig. 10c . This phenomenon reveals that the direction of the heat flux at the surface of the fuel channel is reversed and part of the power deposited in the fuel salt is removed by the graphite when the reactor power increases rapidly. When the reactor power stabilizes, the temperature distribution is similar to that in the initial state, as indicated in Fig. 10f Fig. 10f . 4.3 Unprotected loss of flow 2 The unprotected loss of flow transient is initiated by the loss of forced circulations in the primary loop due to pump failure, where the mass flow exponentially decreases to 10% of the rated flow in 5 s. The relative power and temperature evolution are displayed in Fig. 11 Fig. 11 . The fraction of DNs and the relative number of DNs in the core are displayed in Fig. 12 Fig. 12 . The reduction of the fuel salt velocity has two effects. First, the loss of DNs in the core is reduced, which introduces positive reactivity. As indicated in Fig. 12 Fig. 12 , at the beginning of this transient, the fraction of DNs and number of DNs in the core are increased, which results in the increase of power at the beginning ( Fig. 11a Fig. 11a ). Secondly, the temperature in the core increases, which introduces negative reactivity. The long-term behavior of reactor power is dominated by the temperature feedback, and the reactor power decreases to a low level. The evolution of the relative number of DNs in the core is similar to that of the reactor power. The average temperature of the graphite and reflector changes slowly, and the total temperature increase is small. In this transient, the maximum temperature is located at Assembly 1 and attains 692.4 ℃, which is less than the safety limit of structural material. At t = 2700 s, both the reactor power and average temperature tend to be stable, as indicated in Fig. 11 Fig. 11 . The 3D temperature distribution at t = 2700 s is presented in Fig. 13 Fig. 13 , which is close to the steady-state temperature distribution for this transient. Compared with the temperature distribution at the initial steady state, the temperature field has clearly changed. The region with the maximum temperature is no longer located in the radial reflector, it moves to the outlet of the fuel channel. 5. Conclusion 1 In this work, a special dynamics model for the channel-type MSR was developed by coupling a 3D thermal-hydraulic model, 3DTH, and a point reactor model. The 3DTH adopts a 3D heat conduction model for the solid region and one-dimensional single-phase flow model for the fuel salt, which can consider the thermal coupling between different assemblies. The point reactor model divides the reactor into several nodes and the drift of the DNP is modeled. The correctness and reliability of the 3DTH is verified by RELAP5 in terms of temperature and mass flow calculations. Based on the simulation of the protected fuel pump start-up and coast-down experiments, and the natural circulation experiment, the accuracy of the dynamics model is validated. The dynamics model is applied to analyze the unprotected step reactivity insertion and unprotected loss of flow for the 2MW-MSR. The results confirm that under both 100 pcm reactivity insertion and the loss of flow transients, the maximum temperature of the reactor is less than the safety limit of the structural material and the reactor remains safe without external control action. References R.C. Briant , C.B. Ellis . Aircraft Nuclear Propulsion Program: Quarterly Progress Report for Period Ending December 10, 1951. ORNL, ORNL-1154, USA (1951) . http://moltensalt.org/references/static/downloads/pdf/ORNL-1154.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-1154.pdf Aircraft Nuclear Propulsion Program: Quarterly Progress Report for Period Ending December 10, 1951. ORNL, ORNL-1154, USA (1951) E.S. BETTIS , R.W. SCHROEDER , G.A. CRISTY , et al . The Aircraft Reactor Experiment—Design and Construction . Nucl. Sci. Eng.  2 , 804 - 825 ( 1957 ). doi: 10.13182/NSE57-A35495 http://doi.org/10.13182/NSE57-A35495 The Aircraft Reactor Experiment—Design and Construction R.C. Robertson . MSRE Design and Operations Report Part I: Description of Reactor Design. ORNL, ORNL-TM-0728, USA (1965) . http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-0728.pdf MSRE Design and Operations Report Part I: Description of Reactor Design. ORNL, ORNL-TM-0728, USA (1965) P.R. Kasten , E.S. Bettis , W.B McDonald , et al . Summary of Molten-Salt Breeder Reactor Design Studies. ORNL, ORNL-TM-1467, USA (1966) . http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-1467.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-1467.pdf Summary of Molten-Salt Breeder Reactor Design Studies. ORNL, ORNL-TM-1467, USA (1966) E. Merle-Lucotte , D. Heuer , C. Le Brun , et al . The TMSR as Actinide Burner and Thorium Breeder . Lucotte , ( 2007 ). http://hal.in2p3.fr/in2p3-00143680/document http://hal.in2p3.fr/in2p3-00143680/document The TMSR as Actinide Burner and Thorium Breeder D. Heuer , E. Merle-Lucotte , M. Allibert , et al . Simulation tools and new developments of the Molten Salt Fast Reactor . RGN Revue Générale Nucléaire . 6 , 95 - 100 ( 2010 ). doi: 10.1051/rgn/20106095 http://doi.org/10.1051/rgn/20106095 Simulation tools and new developments of the Molten Salt Fast Reactor K. Furukawa , K. Arakawa , L.B. Erbay , et al . A road map for the realization of global-scale thorium breeding fuel cycle by single molten-fluoride flow . Energy Conversion & Management . 49 , 1832 - 1848 ( 2008 ). doi: 10.1016/j.enconman.2007.09.027 http://doi.org/10.1016/j.enconman.2007.09.027 A road map for the realization of global-scale thorium breeding fuel cycle by single molten-fluoride flow V. Ignatiev , O. Feynberg , I. Gnidoi , et al . Progress in Development of Li,Be,Na/F Molten Salt Actinide Recycler & Transmuter Concept . In: Proceedings of ICAPP 2007 , Nice, France ( 2007 ). Progress in Development of Li,Be,Na/F Molten Salt Actinide Recycler & Transmuter Concept M.H. Jiang , H.J. Xu , Z.M. Dai . Advanced Fission Energy Program-TMSR Nuclear Energy System . Bulletin of Chinese Academy of Sciences . 27 , 366 - 374 ( 2012 ). (in Chinese) doi: 10.3969/j.issn.1000-3045.2012.03.016 http://doi.org/10.3969/j.issn.1000-3045.2012.03.016 Advanced Fission Energy Program-TMSR Nuclear Energy System 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 J. Kópházi , D. Lathouwers , J.L. Kloosterman . Development of a three-dimensional time-dependent calculation scheme for molten salt reactors and validation of the measurement data of the molten salt reactor experiment . Nucl. Sci. Eng . 163 , 118 - 131 ( 2009 ). doi: 10.13182/NSE163-118 http://doi.org/10.13182/NSE163-118 Development of a three-dimensional time-dependent calculation scheme for molten salt reactors and validation of the measurement data of the molten salt reactor experiment S.Y. Si , Q.C Chen , H. Bei , et al . Modeling and simulation of redesigned thorium molten salt reactor . High Power Laser and Particle Beams . 29 , 71 - 78 ( 2017 ). (in Chinese) doi: 10.11884/HPLPB201729.160260 http://doi.org/10.11884/HPLPB201729.160260 Modeling and simulation of redesigned thorium molten salt reactor Q. Wei , W. Guo , H.L. Wang , et al . Develop and verify coupling program of the neutron physics and thermal hydraulic for MSR . Nucl. Tech . 40 , 100605 ( 2017 ). (in Chinese) doi: 10.11889/j.0253-3219.2017.hjs.40.100605 http://doi.org/10.11889/j.0253-3219.2017.hjs.40.100605 Develop and verify coupling program of the neutron physics and thermal hydraulic for MSR L. He , C.G. Yu , W. Guo , et al . Development of a Coupled Code for Steady-State Analysis of the Graphite-Moderated Channel Type Molten Salt Reactor . Science and Technology of Nuclear Installations , vol. 2018 , 4053254 , ( 2018 ). doi: 10.1155/2018/4053254 http://doi.org/10.1155/2018/4053254 . Development of a Coupled Code for Steady-State Analysis of the Graphite-Moderated Channel Type Molten Salt Reactor J.W. Cooke , B. Cox . Forced-convection heat-transfer measurements with a molten fluoride mixture flowing in a smooth tube. ORNL, ORNL-TM-4079, USA (1973) . http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-4079.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-4079.pdf Forced-convection heat-transfer measurements with a molten fluoride mixture flowing in a smooth tube. ORNL, ORNL-TM-4079, USA (1973) L. Li , Z.J. Zhang . Research on transient flux distribution in parallel channels . Nucl. Power Eng . 31 , 97 - 101 ( 2010 ). (in Chinese) Research on transient flux distribution in parallel channels COMSOL Multiphysics . Introduction To COMSOL Multiphysics. COMSOL Multiphysics , ( 2016 ). https://cdn.comsol.com/documentation/5.3.0.316/IntroductionToCOMSOLMultiphysics.pdf https://cdn.comsol.com/documentation/5.3.0.316/IntroductionToCOMSOLMultiphysics.pdf Introduction To COMSOL Multiphysics. COMSOL Multiphysics C.M. Allison , R.J. Wagner , J.K. Hohorst , M.P. Ferragut , et al . The development of RELAP/SCDAPSIM/MOD4.0 for advanced fluid systems design analysis . The 11th International Topical Meeting on Nuclear Reactor Thermal Hydraulics , Operation and Safety Gyeongju , Korea , October 9-13, 2016 . The development of RELAP/SCDAPSIM/MOD4.0 for advanced fluid systems design analysis B.E. Prince , S.J. Ball , J.R. Engel , et al . Zero- power Physics Experiments on the Molten-salt Reactor Experiment. ORNL, ORNL-4233, USA (1968) http://moltensalt.org/references/static/downloads/pdf/ORNL-4233.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-4233.pdf Zero- power Physics Experiments on the Molten-salt Reactor Experiment. ORNL, ORNL-4233, USA (1968) M.W. Rosenthal , R.B. Briggs , P.R. Kasten . Moltensalt Reactor Program Semiannual Progress Report. ORNL, ORNL-4396, USA (1969) http://moltensalt.org/references/static/downloads/pdf/ORNL-4396.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-4396.pdf Moltensalt Reactor Program Semiannual Progress Report. ORNL, ORNL-4396, USA (1969) M. Delpech , S. Dulla , C. Garzenne , J. Kophazi , J. Krepel , et al . Benchmark of dynamic simulation tools for molten salt reactors . GLOBAL 2003 - Nuclear Science and Technology: Meeting the Global Industrial and R&D Challenges of the 21st Century, Nov 2003 , New Orleans, United States . American Nuclear Society , pp. 2182 - 2187 , 2003 . http://hal.in2p3.fr/in2p3-00020303 http://hal.in2p3.fr/in2p3-00020303 Benchmark of dynamic simulation tools for molten salt reactors S.J. Ball , T.W. Kerlin . Stability Analysis of the Molten-Salt Reactor Experiment. ORNL, ORNL-TM-1070, USA (1965) . http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-1070.pdf http://moltensalt.org/references/static/downloads/pdf/ORNL-TM-1070.pdf Stability Analysis of the Molten-Salt Reactor Experiment. ORNL, ORNL-TM-1070, USA (1965) J. Krepel . Dynamics of Molten Salt Reactors. Ph.D. Thesis , Czech Technical University in Prague , ( 2006 ). Dynamics of Molten Salt Reactors. Ph.D. Thesis D.T. Ingersoll , C.W. Forsberg , L.J. Ott , et al . Status of preconceptual Design of the Advanced High-Temperature Reactor(AHTR). ORNL, ORNL-TM-2004-104 , USA ( 2004 ). doi: 10.2172/861752 http://doi.org/10.2172/861752 Status of preconceptual Design of the Advanced High-Temperature Reactor(AHTR). ORNL, ORNL-TM-2004-104 ORNL . SCALE: A Modular Code System for Performing Standardized Computer Analyses for Licensing Evaluations. ORNL, ORNL-TM-2005, USA (2005) . http://www.wipp.energy.gov/library/CRA/CRA-2014/References/Others/ORNL_2009_SCALE_Version6.pdf http://www.wipp.energy.gov/library/CRA/CRA-2014/References/Others/ORNL_2009_SCALE_Version6.pdf SCALE: A Modular Code System for Performing Standardized Computer Analyses for Licensing Evaluations. ORNL, ORNL-TM-2005, USA (2005) G. Mckinney . MCNP6 User's Manual - Version 1.0. Los Alamos National Laboratory, Report number: LA-CP-13-00634 , ( 2003 ). MCNP6 User's Manual - Version 1.0. Los Alamos National Laboratory, Report number: LA-CP-13-00634 The authors wish to acknowledge Dr. Hai-Ling Wang for providing the computation resources and technical support during the simulation. 10.1007/s41365-018-0541-7 This work was supported by the Chinese TMSR Strategic Pioneer Science and Technology Project (No.XDA02010000), the Frontier Science Key Program of the Chinese Academy of Sciences (No.QYZDY-SSW-JSC016), and the National Natural Science Foundation of China Key Program (No.91326201). 2018-05-31 2018-08-08 2018-08-10 2019-01-07 2019-01-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. 2019 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-018-0541-7 http://dx.doi.org/10.1007/s41365-018-0541-7 . 2019 Nuclear Science and Techniques 01 30