Thermal hydraulics model

In the 3DSTH-MSR, the temperature field within the reactor core of the GMC-MSR, which encompassed the graphite assemblies, reflector, and liquid region, was calculated using a 3D convection–diffusion equation: $\begin{array}{l}\frac{\partial \left(\rho C_{p}uT\right)}{\partial x}+\frac{\partial \left(\rho C_{\text{p}}vT\right)}{\partial y}+\frac{\partial \left(\rho C_{\text{p}}wT\right)}{\partial z}\\ =\frac{\partial }{\partial x}\left({\lambda }_{xx}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_{yy}\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left({\lambda }_{zz}\frac{\partial T}{\partial z}\right)+Q,\end{array}$ (1) where ρ, u, v, w, C_p, Q, and T are the density, x-, y-, and z-components of the velocity vector, specific heat capacity at constant pressure, heat source, and temperature, respectively. Thermal conductivity λ is expressed in the tensor form, which is a symmetric matrix. $\lambda =\left[\begin{array}{ccc}{\lambda }_{xx}& 0& 0\\ 0& {\lambda }_{yy}& 0\\ 0& 0& {\lambda }_{zz}\end{array}\right]$ (2) where λxx, λyy and λzz represent the xx, yy, and zz components of tensor λ, respectively. For the solid region, isotropic thermal conductivity was employed, where ${\lambda }_{xx}={\lambda }_{yy}={\lambda }_{zz}$. However, in the liquid region, anisotropic thermal conductivity is utilized. Equation (1) neglects the contributions of viscous dissipation, pressure work, and kinetic energy.

Within the solid region, u, v, and w were each set to zero, leading to the transformation of Eq. (1) into a heat conduction equation with an internal heat source. $\frac{\partial }{\partial x}\left({\lambda }_{xx}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_{yy}\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left({\lambda }_{zz}\frac{\partial T}{\partial z}\right)+Q=0$ (3) To obtain a 1D fuel-salt temperature profile along flow direction z for each channel, the λxx and λyy components of thermal conductivities in the 3DSTH-MSR were assigned infinite values (1×10⁶ W/(m·K)) in the 3D liquid region. Because we set λxx and λyy to infinitely large quantities, the x and y components of temperature gradients ($\frac{\partial T}{\partial x}$ and $\frac{\partial T}{\partial y}$) tended to zero. This led to a nearly uniform distribution of the fuel–salt temperature within each channel along the direction perpendicular to the flow. In addition, the fuel-salt temperature profile of each channel was affected by the z-component of the temperature gradient.

To describe the heat convection along the z direction, velocity components u and v were set to zero, while the z-component of velocity w was set as $w=\frac{m_{\text{flow}}}{\rho (T)\cdot A},$ (4) where m_flow and A denote the mass flow and cross-sectional area of the fuel channel, respectively. $\rho (T)$ denotes the density of fuel salt, which is a function of temperature. The z-component of velocity w was updated after each iteration of the temperature calculation. By modifying the conductivity and velocity in the 3D liquid region, Eq. (1) yields a 1D fuel-salt temperature profile along the flow direction, similar to that of the 1D heat convection equation.

Thermal coupling between the solid and liquid regions was established by adding thermal resistance at the solid–liquid interface, as shown in Fig. 1(a). Thermal resistance causes temperature discontinuities at the interface, as shown in Fig. 1(b). The heat flux at the solid–liquid interface can be written as $q=\frac{\text{Δ}T}{R}\mathrm{,}$ (5) where q, $\text{Δ}T$, and R represent the heat flux, temperature difference at the solid–liquid interface, and specific thermal resistance, respectively. The specific thermal resistance is related to the convective heat transfer coefficient and Nusselt number using the following equation: $h=\frac{Nu\cdot {\lambda }_{zz}}{D}=\frac{1}{R},$ (6) where h, Nu, and D are the convective heat transfer coefficient, Nusselt number, and hydraulic diameter of the fuel channel, respectively. Nu was calculated based on a forced-convection heat-transfer experiment conducted by ORNL. The fitted empirical correlation for the laminar flow has the following form [34]: $Nu=1.89{\left(Re\cdot Pr\cdot \frac{D}{L}\right)}^{0.33}{\left(\frac{{\mu }_l}{{\mu }_s}\right)}^{0.14},$ (7) where Re, Pr, L, μl, and μs represent the Reynolds number, Prandtl number, fuel channel length, dynamic viscosities calculated using the liquid temperature, and dynamic viscosity calculated using the solid surface temperature, respectively.

Fig. 1Full-screen View

The sketch maps of solid–liquid coupling interface (a) and the temperature profile (b)

After imposing thermal resistance at the solid–liquid interface and modifying λxx and λyy for the liquid region, the radial temperature profile of the fuel assembly was similar to that shown in Fig. 1(b). In the 3DSTH-MSR, the adiabatic boundary condition was applied to the outer walls of the solid region, while the temperature boundary condition was set at the inlet of the liquid region.

The pressure losses due to gravity, friction, localized effects, and acceleration effects in the i_th channel are calculated using the following equations: $\text{Δ}{P}_{\text{g}\mathrm{,}i}=\frac{1}{A_i}{\displaystyle \iiint }\rho g\text{d}x\text{d}y\text{d}z\mathrm{,}$ (8) $\text{Δ}{P}_{\text{f}\mathrm{,}i}=\frac{1}{A_i}{\displaystyle \iiint }f\frac{\rho }{2D_i}{\left(\frac{m_{\text{flow}\mathrm{,}i}}{\rho \cdot A_i}\right)}^2\text{d}x\text{d}y\text{d}z\mathrm{,}$ (9) $\text{Δ}{P}_{\text{k}\mathrm{,}i}=\frac{k}{A_i}{\displaystyle \iint }\frac{\rho }{2}{\left(\frac{m_{\text{flow}\mathrm{,}i}}{\rho \cdot A_i}\right)}^2\text{d}x\text{d}y\mathrm{,}$ (10) $\text{Δ}{P}_{\text{a}\mathrm{,}i}=\frac{m_i}{A_i}\left(w_{\text{o}}-w_i\right)\mathrm{,}$ (11) where g, f, k, w_o, and wi represent the acceleration due to gravity, friction factor, coefficient of resistance, and z component of the velocity vector at the channel outlet and inlet, respectively. Subscripts g, f, k, and a represent the gravitational, frictional, local, and acceleration pressure drops, respectively. Total pressure drop $\text{Δ}{P}_{\text{t}\mathrm{,}i}$ for the i_th channel is defined as $\text{Δ}{P}_{t\mathrm{,}i}=\text{Δ}{P}_{g\mathrm{,}i}+\text{Δ}{P}_{f\mathrm{,}i}+\text{Δ}{P}_{k\mathrm{,}i}+\text{Δ}{P}_{a\mathrm{,}i}\mathrm{,}$ (12) The mass flow in different fuel channels was determined by adjusting the mass-flow distribution until a uniform pressure drop was attained across all the fuel channels.

Figure 2 illustrates the computational procedure for the 3DSTH-MSR. First, in the initialization stage, uniform flow rate and temperature distribution were established. These initial conditions were utilized to calculate the z-components of velocity and thermal resistance. Subsequently, Eq. (1) was solved to obtain the temperature distribution. The pressure drops for each channel are updated based on the solution of Eq. (1). Subsequently, the flow distribution was readjusted based on these updated pressure drops, initiating a new iteration of the calculation process. This iterative procedure was repeated until convergence, which ultimately provided the desired solution. The computational procedure of 3DSTH-MSR method, as shown in Fig. 2, is implemented on COMSOL Multiphysics [35]. The Heat Transfer module of COMSOL was used to calculate the temperature distribution in both the solid and liquid phases. The MATLAB interface in COMSOL, was employed to extract the pressure drop values and reset the mass flow distribution for the subsequent iteration of the temperature calculation. The 3DSTH-MSR can perform calculations on complex topologies using either structured or unstructured meshes.

Fig. 2Full-screen View

Flow chart of 3DSTH-MSR

It treats the liquid and solid domains as a collection of geometric objects. In the process of constructing the geometric model, COMSOL automatically detected the interfaces between the solid and liquid domains and established contact pairs to ensure appropriate contact between them. Heat flux, described by Eq. (5), was implemented by applying thermal resistance R at these contact pairs, which allowed coupling of the 3D solid temperature with the liquid temperature. The 3DSTH-MSR required only the solution of a 3D convection–diffusion equation, enabling efficient and automatic thermal coupling between various solid and liquid interfaces. This distinctive characteristic not only yielded substantial reductions in computational time compared with conventional CFD methods, but also eliminated the need for supplementary procedures to couple the 3D solid temperature with the 1D liquid temperature.

Verification of 3DSTH-MSR

The 3DSTH-MSR model developed in this study for GMC-MSR consisted of two primary components. The first component of the model involved the utilization of a single partial differential equation to derive the temperature distributions of both graphite and molten salt. The second component was grounded in the principle of equal pressure drop to ascertain the flow distribution within each flow channel. Regarding the validation of the second part, because it utilized a widely applicable flow-distribution calculation approach, this model has already been employed and verified in a previous study [36]. To verify the temperature field computed using the 3DSTH-MSR model, we employed the RELAP5 code [37].

To verify the program, a circular pipe cooled by a coolant was selected as the simulation object. The pipe had a length of 0.6 m, an inner diameter of 0.02 m, and a wall thickness of 0.01 m. We used FLiBe as the coolant. The thermal conductivity of the pipe wall material was 20 W/(m· K). For the boundary conditions, the inlet mass flow rate of the pipe was set to 0.5 kg/s, and the corresponding inlet temperature was 873.15 K. The total thermal power within the pipe was 23860.0 W. Notably, under these operating conditions, the Dittus–Boelter equation was utilized to calculate the Nusselt number for turbulent flow in both the 3DSTH-MSR and RELAP5 simulations. Due to the thermal power distribution in the GMC-MSR being different from that in traditional reactors, two distinct cases, Cases A and B, were considered. In Case A, 90% of the thermal power was uniformly deposited in the coolant, with the remaining portion being uniformly deposited along the tube wall. In Case B, the opposite configuration was used. For the calculation of the solid region temperature, the RELAP5 code utilized a 1D heat conduction model, disregarding heat conduction along the direction of coolant flow. The conductive heat flux in the axial direction can be neglected because the tube wall was actively cooled.

To accurately verify the 3D temperature field calculated by the 3DSTH-MSR, its results and those obtained using RELAP5 were compared. Figure 3(a) and Fig. 4(a) compare the radial temperature profiles at five positions along the pipe for Cases A and B, respectively. These positions were selected sequentially from the bottom to the top of the pipe. The calculated coolant temperature (in the liquid region) and tube wall temperature (in the solid region) obtained from the 3DSTH-MSR corroborated those calculated using RELAP5. The 3D temperature profile of the pipe obtained using the 3DSTH-MSR for Cases A and B is shown in Fig. 3 (b)) and Fig. 4(b). The verification calculations exhibited satisfactory agreement between the 3DSTH-MSR and RELAP5 code results, confirming the accuracy and reliability of 3DSTH-MSR.

Fig. 3Full-screen View

(Color online) Temperature distribution in Case A. (a) Radial temperature distribution in the pipe. (b) Temperature profile of the pipe (unit: K)

Fig. 4Full-screen View

(Color online) Temperature distribution in Case B. (a) Radial temperature distribution in the pipe. (b) Temperature profile of the pipe (unit: K)

2MW-MSR description

The 2MW-MSR is a small-scale experimental MSR designed by the Shanghai Institute of Applied Physics, Chinese Academy of Sciences. A schematic of the MSR and its dimensional parameters are illustrated in Fig. 5. The designed rated power and mass flow rate of the 2MW-MSR were 2 MW and 59.25 kg/s, respectively. The 2MW-MS was designed with an inlet temperature of 600 ℃ and an outlet temperature of 620 ℃.

Fig. 5Full-screen View

(Color online) Schematic of 2MW-MSR [38]. (a) vertical cross-section; (b) horizontal cross-section; (c) cross-section of the assembly

The reactor core comprises a block of graphite, which can be subdivided into 85 fuel assemblies, as shown in Fig. 5(b). The fuel assembly had a hexagonal cross-section with a circular fuel channel. Its dimensions are shown in Fig. 5(c). The 2MW-MSR utilized graphite as its moderator, with Hastelloy-N alloy as the material for both the reactor vessel and structure. The thermophysical properties of graphite, Hastelloy-N alloy, and fuel salt are listed in Table 1 [36]. The graphite and the reactor vessel in the 2MW-MSR were separated by a 5-mm gap. A portion of the molten salt (the outer salt layer in Fig. 5(a) and (b) flows through this gap, providing a cool environment for the reactor vessel.

The distribution of thermal power within the reactor core of the 2MW-MSR was calculated using the SCALE code [39] in previous research [28]. The results indicated that 65.68% of the thermal power was deposited in the fuel salt region of the assembly, while 3.39% was deposited in the graphite region of the assembly. Additionally, the fractions of thermal power deposited in the fuel salt within the bottom plenum, top plenum, and 5-mm gap were 11.25%, 12.50%, and 4.41%, respectively. The remaining thermal power was released into the structural material and reflector components. Figure 6 depicts the 3D thermal power distribution for the simulation in our study. The impact of the drift of delayed neutron precursors (DNP) and decay heat on the thermal power distribution in the active core was disregarded. Zhang [40] obtained the neutron flux distribution in an MSR for various fuel velocities under steady-state conditions. They observed that the drift of delayed neutron precursors had a negligible impact on the neutron flux. Zhou [41] studied decay heat distribution in a comparable MSR. Their findings revealed that the decay heat released in the active core accounted for 46.7% of the total decay heat. The influence of the drift of the DNPs and delayed heat on the temperature distribution in the MSR was found to be negligible. This is because the delayed neutron contribution accounted for approximately 0.67% of the total neutron population, while the decay heat contributed approximately 7% of the total thermal power.

Fig. 6Full-screen View

(Color online) Power density distribution in the 2MW-MSR (unit: MW/m³)

Results and discussion

The calculation model included all the assemblies, a radial reflector, and flowing salt in the 5-mm gap. However, the plenums and the top and bottom regions of the graphite reflectors were not included in the model. The flow rates in channels except for outer salt layer were calculated using the principle of equal total pressure drop in the 3DSTH-MSR. In addition, the flow rate in the outer salt layer was set to the designed value of 5% of the rated flow.

In the analysis, a free tetrahedral mesh was adopted. A grid independence test was conducted on the calculation model. Table 2 lists the maximum and average temperatures of graphite and fuel salt for different grid numbers. The results indicated that the temperature was not significantly affected by the increase in grid number beyond 0.82 million grids.

The temperature fields of the 2MW-MSR are shown in Fig. 7. Because the energy deposited in graphite was cooled by the flowing fuel salt, the graphite reached a higher temperature than the salt. However, because only a small fraction of the fission power is released in graphite and because of the excellent convective heat transfer coefficient between the fuel salt and graphite, the temperature difference between the graphite and fuel salt is not significant. The maximum temperature reached in the 2MW-MSR was 628.7 ℃, and it was located in the top region of the radial reflector.

Fig. 7Full-screen View

(Color online) Temperature distribution of 2MW-MSR (unit: ℃)

The temperature distribution and isotherms at the horizontal cross-section of the core at half height are presented in Fig. 8. The salt outer layer plays a crucial role in cooling the reflector adjacent to the reactor vessel, causing the reflector near the reactor vessel to have a lower temperature than that near the fuel assemblies. As depicted in Fig. 8(b), the majority of isotherms within the inner region of the reactor core encircle the fuel channels. This observation suggests that most of the thermal energy generated in the graphite of each assembly was dissipated through the surrounding fuel salt within that assembly. Moreover, this pattern indicates insignificant thermal coupling between neighboring assemblies. The fuel assemblies located at the periphery of the core, along with the outer salt layer, play a vital role in dissipating the energy released in the radial reflector. The uniform cooling effect of the outer salt layer results in isothermal lines within the reflector adjacent to the reactor vessel, demonstrating a distinct circular shape and concentric circular distribution. Owing to the combined influence of the radial reflector and the cooling effect of the fuel assemblies, the isotherms near the reflector in the fuel assemblies exhibited noticeable deviations from a circular shape.

Fig. 8Full-screen View

(Color online) Temperature field (a) and isotherm distribution (b) in cut plane z=0.55 m (unit: ℃)

As part of the investigation of heat transfer characteristics of the 2MW-MSR, the conductive heat flux vectors in a 1/4 cross-section at half the height of the core are depicted in Fig. 9. The lengths and directions of the vectors represent the magnitude and direction of heat transfer, respectively. Because the energy deposited in an assembly is primarily dissipated by the fuel salt within the same assembly, the heat flux vectors within an assembly point toward the center of the fuel channel. The conductive heat flux vectors in the reflector exhibit two distinct characteristics. When the radius is less than 0.63 m, the heat flux vectors point toward the center of the reactor. Conversely, when it exceeds 0.63 m, the vectors point outward. This phenomenon indicates that the energy deposited in the reflector within a radius of 0.63 m is primarily removed by the fuel assemblies, while the remaining energy deposited in the reflector is cooled by the outer salt layer.

Fig. 9Full-screen View

(Color online) The conductive heat flux vectors in the 1/4 cross-section (z=0.55 m)

The mass flow distribution in the reactor core of the 2MW-MSR MW is illustrated in Fig. 10. The highest mass flow of 0.73 kg/s was observed at the center of the reactor, while the lowest mass flow of 0.64 kg/s was observed in the fuel channel adjacent to the reflector. The mass flow calculation considered only the gravitational, frictional, and acceleration pressure drops. The local pressure drops at the top and bottom plenums caused by the presence of a distribution plate and changes in the flow section were not considered. Moreover, the density of the fuel salt decreases with increasing temperature; consequently, a greater mass flow was allocated to assemblies with higher energy deposition, and the mass flow distribution closely followed the power density profile.

Fig. 10Full-screen View

(Color online) The mass flow in each channel (unit: kg/s)

Figure 11 depicts the salt outlet temperature in each fuel channel and the outer salt layer. As observed in the figure, the outlet temperature of the 2MW-MSR is relatively uniform, with a maximum deviation of approximately 1 ℃. This is primarily attributed to the higher mass flow distribution in the fuel channels with higher fission energy. In this simulation, a designed mass flow value was specified for the outer salt layer. This layer not only removed 4.41% of the rated power deposited within it, but also effectively dissipated most of the energy in the reflector region. In this design, the highest outlet temperature of the fuel salt was observed in the outer salt layer.

Fig. 11Full-screen View

(Color online) The outlet fuel temperature (unit: ℃)