The nuclear thermal propulsion (NTP) reactor has become a potential leading power option for NASA’s crewed landing on Mars due to its much shorter transit time in a long-distance space mission (Hibberd and HEIN, 2021). NTP reactor designs are grouped into high-enriched uranium (HEU) and low-enriched uranium (LEU) designs according to the different enrichment of the nuclear fuel used. A few decades of the cold war between the United States and the Soviet Union had witnessed the amount of mature and effective numerical and experimental work (Belair et al., 2013; Khatry et al., 2019; Graham, 2020) on HEU design. However, recent efforts focus on designing a feasible engine that relies on LEU fuel due to its lower cost and nuclear proliferation risk (Venneri and KIM, 2015a; Venneri and KIM, 2015b; Gates et al., 2018). In LEU design, moderator assembly is employed to cooperate with fuel assembly to support the whole structure and take some heat away. What is more, this design can slow down the fast neutrons produced from the latter and keep the system critical. These new LEU-based generations of NTP require geometrical modifications associated with the assemblies’ thickness and pitch, as well as the core’s length and configuration. These modifications may bring large uncertainties in using empirical models or correlations based on HEU design. Hence, it is essential to understand the influence of moderator assembly on heat transfer, which is required for a core design using a comprehensive analysis method involving neutronics, thermal-hydraulics, and dynamic characteristics of the NTP system.

Several numerical methods have been proposed to model NTP sub-elements over a few decades. Hall et al. (1992) analyzed the thermal-hydraulic performance of the gaseous flow using the KLAXON code. This was among the first computational sequences developed to study an NTP system. Nevertheless, only a simple 1D model was used to predict pressure distribution from the inlet of the reactor core to the exit of the converging-diverging nozzle. Cheng and Yen-Sen (2015) conducted 3D numerical simulations on NERVA-type engine FE using self-developed computational fluid dynamic code (UNIC). The simulation implemented a conjugate heat transfer (CHT) routine to couple the solid and flow regions. Husemeyer (2016) conducted the coupled neutronic-thermal-hydraulic analysis on a Pewee-type engine. The calculation routine relied on a self-developed 3D finite element method to capture temperature distributions across the core. However, most of the work above relied on empirical heat transfer correlations and focused on HEU objects instead of LEU ones. Recently, Wang and Kotlyar (2021) utilized a reduced-order 1.5D semi-analytic solution to implant the legacy heat transfer correlations for a 3D CHT numerical solver on the OpenFOAM platform. Moreover, the in-house 1.5D reduced-order solver, namely, THERMO, uses an inaccurate heat transfer model between fuel assembly and moderator assembly due to the current implementation limitation within THERMO (Krecicki and Kotlyar, 2020).

The purpose of this work is to introduce the heat transfer model of the LEU nuclear thermal propulsion assembly. The flow heat transfer in each coolant channel is simplified into one-dimensional axial heat transfer, characterized by the Newton cooling formula. The radial solid heat conduction is accurately expressed by a three-dimensional heat conduction equation. The heat conduction equation provides heat flux boundary conditions for axial convective heat transfer. The temperature distribution in the fuel assembly and moderating assembly can be obtained by solving the two heat transfer processes. This model can provide thermo-hydraulic compliance for the subsequent conceptual design of the core. Still, the accuracy of the model largely depends on the heat transfer coefficient in the expression. Therefore, referring to the empirical expression of legacy experimental data (Walton, 1992) and utilizing a large number of numerical simulations conducted by STAR-CCM+, the fitted heat transfer coefficient expressions are obtained, which covers as many working conditions as possible to improve the applicability of the correlation.

For NTP reactors, cores with different power levels are usually designed according to the mission to meet different mission requirements, such as a crewed lunar landing or crewed Mars landing. Different tasks correspond to different core power levels, and the fuel assembly as the basic unit of output power may operate at different power. It is essential to ensure that the fuel and moderator assembly operate safely within a reasonable power range. In order to ensure that the temperature does not exceed the limit, different output powers of components correspond to different flow characteristics (i.e., Reynolds number). As a typical case, the inlet Reynolds number of the simulation domain is around 1,00,000. According to this typical reference case (case 2), the other four cases are expanded, as shown in Table 4. In cases 1 and 3, the inlet Reynolds number is decreased and increased by 20%, respectively. In cases 4 and 5, 20% reduction and increment will occur, respectively, when the heating power is modified with other parameters of case 2.

Based on case 2, preliminary simulation can be first carried out by STAR-CCM+ to understand heat transfer between assemblies. Figures 7, 8 give the cross-sectional profiles of velocity (axial direction) and temperature at two assembly heights. They show that the velocity and temperature vary significantly among the channels, ranging from 23 to 537 m/s and from 45 to more than 2,200 K, respectively. The negative values shown in Figure 7 indicates that the fluid flow direction is perpendicular to the paper surface outward. The temperature of the coolant varies significantly in the axial direction. Therefore, the density, pressure, and speed vary greatly, and the speed of the drastic change, in turn, affects the temperature distribution. The influence between them is a two-way intense coupling process.

Significant temperature changes can be confirmed by the heat flow distribution in each component, as shown in Figure 9. It shows axial linear heat power distribution of all channels, and it can be seen that most of the heat is taken away by the coolant in the fuel channel, so its temperature changes most violently, which can also be further confirmed by the axial temperature distribution of the fuel assembly, as shown in Figure 10. Besides, Figure 9 shows that the fuel assembly transfers heat to both the return channel and supply channel of the moderating assembly. The quantitative results show that this part of heat accounts for about 11.6% of the total. Moreover, the heating power of the return channel has some negative values near the outlet section, indicating that the coolant has an outward heat transfer process.

The reverse heat transfer can be seen from the axial distribution of the bulk and wall temperatures in the return channel, as shown in Figure 11, which also shows that the bulk temperature in the return channel increases and then decreases along the flow direction because the coolant flow direction of the return channel is contrary to that of other channels, which adds the complexity of finding the heat exchange law between assemblies. In case of the possible failure of the existing empirical correlation of heat transfer coefficient, it is necessary to establish a new inter-assembly heat transfer model. This model can be regarded as reliable for future core design by accurately predicting the fuel and moderator peak temperature. The above analysis reveals the general phenomenon of heat transfer between assemblies. However, the heat distribution in each channel and peak temperature inside the assembly under different working conditions still need more profound concern. Before further discussion, a new heat transfer model between assemblies must be established. The following part will carefully deduce the heat transfer function based on the simplified axial convective heat transfer model and radial thermal resistance.

Figure 12 displays the equivalent FE and full ME model that preserves the hydraulic parameter and fuel volume. The equivalent FE model draws on the method by Wang and Kotlyar (2021). The method uses 1/19 of the full explicit FE model to get a conservative but quick result, although it is more suited for solving circular and annular pins instead of hexagonal prisms. As for the heat transfer mechanism between assemblies, different from the 1.5D conduction–convection solution (Krecicki and Kotlyar, 2020; Wang and Kotlyar, 2021), the current research focuses on deriving the heat transfer equations between the two assemblies according to the principle of heat flow diversion and superposition (Holman, 2010). Heat flux Q 2 → 1 through the cylinder surface of the fuel channel can be deduced from the heat conduction differential equation and corresponding boundary conditions for the steady-state temperature distribution in the equivalent fuel assembly, shown as Eq. 7. In this equation, L represents the length of the assembly, k is the conductivity of fuel, T 2 means the maximum temperature of the fuel, T 1 is the bulk temperature of the fuel channel, r 1 and r 2 are the equivalent fuel radius (distance from centroid to the position of maximum fuel temperature) and equivalent coolant channel radius, respectively. q v is volumetric heat power. Eq. 7 contains two terms in total. The first term can be regarded as the heat flow generated by the internal and external temperature ( T 1 on L 1 and T 2 on L 2 ) difference acting alone, and the second term can be regarded as the heat flow generated by the interval [ r 1 , r 2 ] uniform heat source acting alone. Combining Eq. 7 with Newton’s cooling theorem (Eq. 8) and eliminating T 1 , Eq. 7 can be expressed as Eq. 9. The calculation Q 2 → 3 , the heat flux from the fuel to the moderator, is more complex than that of Q 2 → 1 because the calculation process of Q 2 → 1 is essentially a steady-state heat conduction problem of concentric rings with an internal heat source, and there are multi-layer structural materials between the return channel of moderator assembly and fuel assembly. Therefore, it is necessary to calculate the heat flux without an internal heat source. The sum of the two parts of heat flux is Q 2 → 3 , as shown in Eq. 10. In Eq. 10, T 3 , is the temperature of the interface between fuel assembly and moderator assembly and r 3 is the equivalent radius of the fuel assembly. Combining Eq. 10 with Eq. 11 and eliminating the intermediate variable T 3 , Eq. 12 can be obtained. In Eq. 12, T b 4 is the bulk temperature of return channel. ∆ r i ,   i = 1,2,3 are the corresponding thickness of graphite, ZrC, and Zr, respectively, and A l i ,   i = 1,2,3,4 are the corresponding heat exchange area of graphite, ZrC, Zr, and ZrH_x, respectively: Q 2 → 1 = 2 π L k ( T 2 − T 1 ) ln r 2 / r 1 + q v π L r 1 2 ln r 2 / r 1 { ( r 2 r 1 ) 2 ⁡ ln r 2 r 1 − 1 2 [ ( r 2 r 1 ) 2 − 1 ] } (7) Q 2 → 1 = h A ( T 1 − T b 1 ) (8) Q 2 → 1 = ( T 2 − T b 1 ) + q v r 1 2 2 k { ( r 2 r 1 ) 2 ⁡ ln r 2 r 1 − 1 2 [ ( r 2 r 1 ) 2 − 1 ] } ln r 2 / r 1 2 π L k + 1 h 1 A 2 (9) Q 2 → 3 = 2 π L k ( T 2 − T 3 ) ln r 3 / r 2 + q v π L r 2 2 ln r 3 / r 2 { ( r 3 r 2 ) 2 ⁡ ln r 3 r 2 − 1 2 [ ( r 3 r 2 ) 2 − 1 ] } (10) Q 2 → 3 = T 3 − T b 4 ∆ r 1 k 1 A l 1 + ∆ r 2 k 2 A l 2 + ∆ r 3 k 3 A l 3 + 1 h 4 A l 4 (11) Q 2 → 3 = ( T 2 − T b 4 ) + q v r 2 2 2 k { ( r 3 r 2 ) 2 ⁡ ln r 3 r 2 − 1 2 [ ( r 3 r 2 ) 2 − 1 ] } ln r 3 / r 2 2 π L k + ∆ r 1 k 1 A l 1 + ∆ r 2 k 2 A l 2 + ∆ r 3 k 3 A l 3 + 1 h 4 A l 4 (12) Q f u e l   e l e m e n t = Q 2 → 3 + Q 2 → 1 (13)

Eqs 9, 12 are only functions of the fuel center and bulk temperatures. If Eq. 13 is simultaneously used when the inlet and outlet temperature of the fuel channel and fuel assembly power (this is usually considered as boundary conditions) are known, the fuel center and bulk temperatures of the return channel can be calculated more conveniently. The maximum temperature limit of the fuel center is one of the criteria for future core design. Moreover, it is also essential to predict the maximum temperature of the moderator assembly, which can release the design margin of the coolant flow rate in the current LEU core design. Therefore, it is necessary to establish the relationship between the return channel’s bulk temperature and the moderator’s maximum temperature through STAR-CCM+ simulation so that the maximum temperature of the moderator can be predicted once the bulk temperature is known. Eqs 9, 12 provide an effective way to predict temperature, but the accuracy of the prediction depends on the selection of heat transfer coefficient. The usual practice is to select the empirical correction based on legacy experimental data. Nevertheless, it may not apply to LEU assembly, especially the flow and heat transfer of the return channel. In this work, the CFD method is used to simulate different working conditions, and the heat transfer coefficients are fitted from the simulation results concerning the form of general legacy empirical correction. The evaluated formula is substituted into Eqs 9, 12 for further calculations. The prediction effect evaluation will also be discussed subsequently.

Firstly, based on the CFD simulation, the Nusselt number distribution in the axial direction of the fuel channel under various working conditions is calculated, as shown in Figure 13. With the increase in height, the Nusselt number in the fuel channel increases gradually. The increase in Nu can be explained by the coolant flowing through the fuel assembly from top to bottom. With the decrease in height, the temperature difference between coolant and fuel gradually decreases, and the heat exchange process slows down. It is worth mentioning that when the coolant approaches the outlet, the Nusselt number first decreases and then increases, which is due to the outlet effect, and the bulk temperature of the fuel channel also presents a similar phenomenon.

Figure 14 shows that the Nusselt number in the return channel decreases slowly along the flow direction. However, when approaching the outlet section of the return channel, there is distortion in the Nusselt number curve due to the propellant flowing in two opposite directions in the fuel and return channels. The coolant in the fuel channel cools the fuel assembly from top to bottom. Thus, its temperature rises from top to bottom, as Figure 10 shows. Because the fuel power presents a cosine distribution in the axial direction, the fuel center temperature increases and decreases from top to bottom. Moreover, the fuel channel’s coolant will make the fuel center’s maximum temperature shift downstream along the axis. The coolant in the return channel cools the moderator assembly from bottom to top, and the temperature rises first and then decreases in Figure 11. Therefore, the so-called reverse heat transfer process will occur when the fuel center temperature drops below the return channel bulk temperature. This process occurs after about z = 0.8. Figure 15 shows that the linear heat flux has a negative value. Once reverse heat transfer occurs, the bulk temperature will continue to decrease, which can protect the moderator on the right from high-temperature damage. Further, for evaluating the heat transfer coefficient, only the heat transfer coefficient in the height range of 0–0.8 m can be considered without affecting the conservative estimation of the maximum temperature of the bulk temperature in the return channel. N u f u e l   c h a n n e l = 0.0095 × R e b 1 0.935 Pr b 1 1.5236   ⁡ ( T w 1 / T b 1 ) − 0.0389 , (14) N u r e t u r n   c h a n n e l = 0.023 × R e b 1 0.6257 Pr b 1 1.3736   ⁡ ( T w 4 / T b 4 ) 0.742 + 1.3085 ∗ D 4 / X (15)

Referring to the Dittus–Boelter correlation and approximating it according to the Nussle number curve in Case 2, the mathematical expression for the Nussle number can be obtained as shown in Eqs 14, 15, where R e 1 and R e 4 are the Reynolds numbers of the coolant in the fuel and return channels, Pr 1 and Pr 4 are the Prandtl numbers of the fuel and return channels, T w 1 and T w 4 are the wall temperature of the fuel and return channels, and T b 1 and n are the bulk temperature of the fuel and return channels, respectively. Because the existing Dittus–Boelter correlation equation or Sieder–Tate equation is inapplicable to the operating conditions discussed in this study, our proposed correlation equation is equivalent to an extension of the Dittus–Boelter equation. These equations are all functions of independent variables, Re and Pr. For the expression of the Nusselt number in the return channel, the correction of the entrance effect as an exponential item is also considered. A Python script can realize this approximate process. With the two equations, the abovementioned closed calculation of Eqs 9, 13 can be carried out, and the assemblies’ bulk temperature, wall temperature, and maximum fuel center temperature can be predicted conveniently and quickly with given boundary conditions. It can provide a convenient and fast calculation tool for reactor conceptual design. Nevertheless, from the Nusselt number curves shown in Figures 16, 17, the model curve should be consistent with the CFD calculation results as much as possible. However, as the exit approaches, the consistency between the two models decreases, whether within the fuel or return channel.

Further analysis should be considered to prove the proposed model’s reliability for temperature prediction. Firstly, the difference between the Nusselt number calculated by the two equations and the Nusselt number obtained by STAR-CCM+ numerical simulation is compared. Figures 18, 19 show that the deviations in the Nusselt number between the CFD model and the proposed model are mostly within 10%. The points with a relatively large coefficient deviation are located in the outlet section in the fuel channel. The height is less than 0.2 m, which can also be seen intuitively in Figure 16. In this area, the bulk and wall temperature difference first decreases and then increases. In other words, the bulk temperature increases slowly from low to high, and the wall temperature decreases slowly from high to low. When the two temperatures are approaching, the flow heat transfer capacity will decrease, leading to the falling of the heat transfer coefficient; the bulk temperature slightly exceeds the wall temperature in the approaching outlet section, where a reverse heat transfer process occurs. The heat transfer coefficient increases slightly, but this range is minimal (0–0.05 m), and the heat transfer coefficient changes little. For the two model coefficient evaluation curves in Figure 18, on the premise of ensuring that most of the deviation range is within 10%, the calculation curve of the heat transfer coefficient of the proposed model is lower than that of CFD, which can ensure a conservative estimation temperature prediction in this channel.

Heat transfer coefficient evaluation in the return channel should be considered as well. Figure 19 shows that the Nusselt number has a large deviation at the exit section. Therefore, the predicted value of the heat transfer coefficient in this narrow area deviates greatly from the CFD result. Nevertheless, it does not affect the maximum temperature prediction of the return channel. Figure 11 shows that the bulk temperature in the return channel begins to drop at a height above 0.7 m. Thus, the prediction of the maximum bulk temperature of the return channel only needs to consider the temperature change at the height below 0.7 m. In addition, just like the evaluation of the heat transfer coefficient in the fuel channel, the return channel also ensures that the model calculation curve is lower than the CFD calculation curve, which can not only balance the deviation of the outlet section but also obtain the conservative prediction of the bulk temperature in the return channel.

The purpose of predicting the maximum bulk temperature of the return channel is to obtain a conservative estimate of the maximum ZrH_x moderator temperature. The limit of the maximum ZrH_x moderator temperature (around 1,000 K) is one of the criteria for the future NTP conceptual design. Figure 20 presents the maximum temperature comparison between the return channel and the ZrH_x moderator and shows that the maximum bulk temperature in the return channel is about 100 K higher than that of the ZrH_x moderator. There is a simple linear relationship between the two curves. Then, it can be concluded that it is appropriate to regard the prediction of the maximum bulk temperature of the return channel by the model as a conservative estimation of the maximum moderator temperature.

A simplified heat transfer model between LEU fuel assembly and moderator assembly is proposed by deriving the thermal resistance assumption between the assemblies. The model’s correction evaluation of the heat transfer coefficient is conducted based on CFD simulations and legacy empirical formula. Furthermore, the calculation deviation of the model is calculated and analyzed under different Reynolds numbers and heating power. The results show that the model’s evaluation curve of the axial distribution of heat transfer coefficient is in good agreement with the CFD calculation results, and the overall deviation of the axial distribution of heat transfer coefficient is less than 10%, whether in the fuel or return channel. The deviation near the outlet is large, but the overall evaluation curve is slightly lower than the CFD calculation curve to ensure that the temperature prediction of the model is conservative. The proposed model can be applied to the LEU core design in the future NTP concept design.