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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

Long He，

Cheng-Gang Yu，

Rui-Min Ji，

Wei Guo，

Ye Dai，

Xiang-Zhou Cai

Nuclear Science and Techniques

Vol.30, No.1

Article number 18

Published in print 01 Jan 2019

Available online 07 Jan 2019

DOI：10.1007/s41365-018-0541-7

47000

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 reactorThermal-hydraulicsPoint reactor modelThermal coupling

1. Introduction

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]. In the 1950s, the first MSR Aircraft Reactor Experiment (ARE) [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]. The successful operation of the MSRE established the basic technologies and led to the conceptual design of a molten-salt breeder reactor (MSBR) [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,6], Japan [7], and Russia [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].

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] 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], 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] and a two-dimensional dynamics code [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] are analyzed for several transients.

2. Numerical Model

2.1 Thermal hydraulics model

In the 3DTH, the temperature of the solid region, which includes the graphite assemblies and reflector, is calculated by a 3D heat conduction equation:

ρ_{g} C_{g} \frac{\partial T_{g}}{\partial t} = \nabla λ_{g} \nabla T_{g} (x, y, z) + Q_{g},

(1)

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:

\frac{\partial ρ_{f}}{\partial t} + \frac{\partial (ρ_{f} v)}{\partial z} = 0,

(2)

\frac{\partial (ρ_{f} v)}{\partial t} + \frac{\partial (ρ_{f} v^{2})}{\partial z} = - \frac{\partial p}{\partial z} - \frac{\partial p_{f}}{\partial z} - ρ_{f} g,

(3)

A \frac{\partial (ρ_{f} C_{f} T_{f})}{\partial t} + A \frac{\partial (ρ_{f} v C_{f} T_{f})}{\partial z} = A Q_{f} + Q_{h},

(4)

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³ 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

M C_{f} \frac{d T}{d t} = m_{f} C_{f} (T_{in} - T) + Q,

(5)

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:

T_{in} = \frac{\sum_{i} m_{i} C_{f} T_{out, i}}{\sum_{i} m_{i} C_{f}},

(6)

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:

q (x, y, z) = \frac{N u \cdot λ_{f}}{D_{e}} (T_{g} (x, y, z) - T_{f} (z)),

(7)

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]. 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

\frac{L}{A} \frac{m_{i}^{n} - m_{i}^{n - 1}}{Δ t} + {(Δ P_{a, i})}^{n - 1} = {(Δ P_{i})}^{n} - {(Δ P_{f, i})}^{n - 1} - {(Δ P_{g, i})}^{n - 1} - {(Δ P_{k, i})}^{n - 1},

(8)

where L is the length of the fuel channel; m_i is the mass flow in channel i, which equals to $ρ_{f} v A$ ; 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:

Δ P_{a, i} = \frac{1}{A_{i}} \int \frac{\partial (ρ A_{i} v^{2})}{\partial z} d z = \frac{1}{A_{i}} \int \frac{\partial (m_{i} v)}{\partial z} d z = \frac{m_{i}}{A_{i}} [v_{out} - v_{in}],

(9)

Δ P_{f, i} = \int f \frac{ρ v^{2}}{2 D_{e}} d z = \frac{1}{A_{i}} \int f \frac{m_{i} v}{2 D_{e}} d z = \frac{m_{i}}{A_{i}} \int f \frac{v}{2 D_{e}} d z,

(10)

Δ P_{g, i} = \int ρ g d z,

(11)

Δ P_{k, i} = k \frac{ρ v^{2}}{2} = k \frac{m_{i}}{A_{i}} \frac{v}{2},

(12)

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]:

{\begin{cases} \frac{L}{A} \frac{m_{1}^{n + 1} - m_{1}^{n}}{Δ t} = Δ P^{n} - {(Δ P_{f, 1})}^{n - 1} - {(Δ P_{g, 1})}^{n - 1} - {(Δ P_{k, 1})}^{n - 1} - {(Δ P_{a, 1})}^{n - 1} \\ ⋮ \\ \frac{L}{A} \frac{m_{N}^{n + 1} - m_{N}^{n}}{Δ t} = Δ P^{n} - {(Δ P_{f, N})}^{n - 1} - {(Δ P_{g, N})}^{n - 1} - {(Δ P_{k, N})}^{n - 1} - {(Δ P_{a, N})}^{n - 1} \\ \frac{\partial M_{total}}{d t} = \sum_{i}^{N} \frac{m_{i}^{n + 1} - m_{i}^{n}}{Δ t} \end{cases},

(13)

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] 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

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. The point reactor model can be written as

Fig. 1

Schematic of point reactor model.

\frac{d n (t)}{d t} = \frac{ρ - β}{Λ} n + \frac{1}{V_{c}} \sum_{i = 1}^{N} (\sum_{k = 1}^{M} λ_{i} C_{k, i} V_{k}),

(14)

V_{k} \frac{d C_{k, i} (t)}{d t} = \frac{β_{i}}{Λ} n V_{k} - λ_{i} C_{k, i} V_{k} + \frac{1}{τ_{L}} C_{L, i} (t) V_{L} - \frac{1}{τ_{k}} C_{k, i} (t) V_{k} k = 1,

(15)

V_{k} \frac{d C_{k, i} (t)}{d t} = \frac{β_{i}}{Λ} n V_{k} - λ_{i} C_{k, i} V_{k} + \frac{1}{τ_{k - 1}} C_{k - 1, i} (t) V_{k - 1} - \frac{1}{τ_{k}} C_{k, i} (t) V_{k} k = 2 \cdot \cdot \cdot M,

(16)

V_{l} \frac{d C_{l, i} (t)}{d t} = - λ_{i} C_{l, i} V_{l} + \frac{1}{τ_{M}} C_{M, i} (t) V_{M} - \frac{1}{τ_{l}} C_{l, i} (t) V_{l} l = 1,

(17)

V_{l} \frac{d C_{l, i} (t)}{d t} = - λ_{i} C_{l, i} V_{l} + \frac{1}{τ_{l - 1}} C_{l - 1, i} (t) V_{l - 1} - \frac{1}{τ_{l}} C_{l, i} (t) V_{l} l = 2 \cdot \cdot \cdot L,

(18)

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

The schematic overview of the dynamics model is presented in 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.

Fig. 2

Flowchart of dynamics model.

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

3.1 Validation of the 3DTH

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].

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. 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.

Parameters used for the calculation model.

Parameter	Value	Parameter	Value
Pipe length (m)	0.5	Thermal conductivity of tube wall (W/(m·K))	25.0
Inner diameter (m)	0.02	Total power (W)	59650.0
Tube wall thickness (m)	0.01	Fraction of power in the coolant (%)	90.0
Mass flow (kg/s)	0.5	Fraction of power in tube wall (%)	10.0
Inlet temperature (K)	873.15	Initial temperature (K)	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. 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. The results calculated by the 3DTH are consistent with the RELAP5 code.

Fig. 3

(a) Temperature profiles of tube wall at z = 25.0 cm and (b) Temperature of fuel salt at z = 12.5 cm and 50.0 cm.

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. 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 indicates the mass flow response in each pipe to the change of energy. 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.

Fig. 4

Mass flow history for each pipe.

3.2 Validation of the dynamics model

In this work, three MSRE experiments, the protected fuel pump start-up [19], protected fuel pump coast-down [19], and natural circulation experiment [20, 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] 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. 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, the shapes and trend of these curves calculated in this study are qualitatively consistent with the measured data.

Fig. 5

(Color online) Compensative reactivity inserted during (a) Pump start-up and (b) Coast-down.

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 $\times$ 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, 23] was adopted. 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.

Fig. 6

Power history during 360 minutes of the natural circulation experiment.

4. 2MW-MSR transient calculation

4.1 2MW-MSR description

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]. A vertical and horizontal view of the 2MW-MSR is presented in Fig. 7a and 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. In the 2MW-MSR, graphite is chosen as the moderator and reflector, and Hastelloy-N alloy (safety limit is 730 ℃ [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) flowing upward between the reflector and the reactor vessel, is used to cool the structure material.

Fig. 7

(Color online) Sketch of the experimental reactor and calculation model: (a) Vertical view, (b) Horizontal view, (c) Fuel assembly, (d) Cross-sectional view of the calculation model.

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, 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]; the fraction of power deposited in each component is listed in the 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 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] with the continuous energy cross-section library of ENDF/B-VII.

Power Distribution of 2MW-MSR.

Component	Power fraction (%)	Component	Power fraction (%)
Fuel salt of the assembly	65.68	Outer salt layer	4.41
Graphite of the assembly	3.39	Reactor vessel	0.34
Top plenum and salt channel	12.50	Radial reflector	1.14
Bottom plenum and salt channel	11.25	Top and bottom reflector	1.29

Delayed neutron data for the 2MW-MSR.

Group	Decay constant QUOTE λi(1/s)	QUOTE βi(×10^-5)
1	0.0125	17.07
2	0.0319	122.35
3	0.1094	81.23
4	0.3172	340.97
5	1.3538	91.12
6	8.6643	17.63

4.2 Unprotected step reactivity insertion

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. 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.

Fig. 8

Relative power after 20 pcm, 50 pcm, and 100 pcm reactivity step.

The mass flow in each fuel channel during the 100 pcm reactivity insertion is displayed in 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 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.

Fig. 9

(Color online) (a) Mass flow distribution (number indicates the Assembly number) and (b) Temperature evolutions (Tf_av: average temperature of fuel salt in assembly, Tg_av: average temperature of graphite in assembly, Tr_av: average temperature of radial reflector, Tb_av: average temperature of fuel salt in bottom plenum, Tt_av: average temperature of fuel salt in top plenum, Tf_max: maximum temperature of fuel salt in Aassembly 1, and Tg_max: maximum temperature of graphite in Assembly 1).

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. At the initial steady state (t = 0 s, 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 and 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. 10

(Color online) Temperature distribution of 2MW-MSR at (a) 0 s, (b) 20 s, (c) 100 s, (d) 500 s, (e) 1000 s, and (f) 2000 s (Temperature unit: ℃).

4.3 Unprotected loss of flow

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. The fraction of DNs and the relative number of DNs in the core are displayed in Fig. 12.

Fig. 11

(Color online) (a) Relative power and (b) Temperature evolution during loss of flow transient.

Fig. 12

Fraction of DNs and relative number of DNs in the core.

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, 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). 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. The 3D temperature distribution at t = 2700 s is presented in 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.

Fig. 13

(Color online) Temperature distribution of the 2MW-MSR at 2700 s (Temperature unit: ℃).

5. Conclusion

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