Geometry model

The SIMONS energy conversion system is based on a closed Brayton cycle, where the core serves as the heat source. The core was designed with a solid structure to achieve system miniaturization. Figure 1 shows a schematic illustration of the core. The solid core comprised staggered coolant channels and fuel rods. The nuclear fuel used was uranium carbide, the moderator was graphite, and the cladding was composed of a titanium zirconium molybdenum (Mo-TZM) alloy.

Fig. 1Full-screen View

(Color online)  Longitudinal cross-sectional view of the core  

In Fig. 1, the dashed hexagonal line on the left represents a radially repeatable element located in the core. Because this element exhibits radial symmetry, a radial 1/6 structure was utilized for modeling. The details of the constructed geometric model are shown on the right side of Fig. 1. This model comprised solid regions for the fuel, moderator, and cladding, whereas the coolant channel was represented by a fluid region.

Under normal operating conditions, the low-temperature helium–xenon gas mixture enters the solid core from the channel inlet at the bottom of the core, is driven by the compressor, absorbs the heat released by the nuclear fuel, and flows out from the outlet at the top of the core. An adiabatic section is incorporated before the inlet to achieve a fully developed flow velocity distribution at the inlet of the heating section. The geometric parameters of the model are listed in Table 1.

Thermophysical properties

For helium–xenon cooled microreactors that utilize a closed Brayton cycle, the operating pressure of the helium–xenon gas mixture is typically approximately 2 MPa, and the operating temperature is between 400 and 1300 K[21]. To balance the thermal power, efficiency, miniaturization, and lightweight requirements of the system, a helium–xenon gas mixture with a xenon volume fraction of 12% was selected as the working medium for the cycle.

The temperature-dependent properties of the mixture at a pressure of 1.9 MPa are shown in Fig. 2. Physical properties, such as specific heat, thermal conductivity, and viscosity, are introduced into the CFD software using polynomials that are dependent only on temperature. Density was calculated using the ideal gas equation of state. The physical properties of the solid regions were simplified to constants, as listed in Table 2.

Fig. 2Full-screen View

(Color online) Properties of the helium–xenon gas mixture over the operation temperature

Numerical models

This study employs the widely utilized CFD software STAR-CCM+ to conduct numerical simulations of flow and heat transfer in nuclear power systems[40, 41]. To balance computational time and resources, the Reynolds-averaged Navier–Stokes method was employed to accommodate turbulence fluctuation terms. A turbulence model was introduced to solve the Navier–Stokes equations. Currently, k-epsilon and shear stress transfer (SST) k-omega models are widely used in high Reynolds number (Re) calculations. The SST k-omega model has been utilized in related studies[21, 22, 42] and yielded satisfactory results; thus, it was employed in this study. To enhance the accuracy and convergence speed, an implicit coupling solver was employed for both the fluid and solid regions.

Pr_t is a dimensionless parameter which affects the turbulent thermal conductivity. As mentioned previously, the helium–xenon mixture used in this study is a low Pr fluid that possesses properties distinct from those of conventional fluids. A study[7] showed that the Pr_t model developed by Kays[42] can be used to accurately calculate the Pr_t of a helium–xenon gas mixture, as expressed in Eq. 1. $P{r}_{\text{t}}={\left\{\frac{1}{2P{r}_{\text{t},\text{b}}}+CP{e}_{\text{t}}\sqrt{\frac{1}{P{r}_{\text{t},\text{b}}}}-{\left(CP{e}_{\text{t}}\right)}^2\left[1-\text{exp}\left(-\frac{1}{CP{e}_{\text{t}}\sqrt{P{r}_{\text{t},\text{b}}}}\right)\right]\right\}}^{-1},$ (1) $P{e}_{\text{t}}=\frac{Pr\rho {\epsilon }_{\text{m}}}{\mu },$ (2) where Pr_t,b is the Pr_t of the bulk region (0.85), C is a constant (0.3), and Pe_t is the turbulent Peclet number calculated using Eq. 2. Under the conditions investigated, the temperature did not significantly affect Pr.

Using Nu as a metric to quantify the degree of convective heat transfer is a widely recognized practice. Nu is expressed in Eq. 3 as follows: $Nu=\frac{Lh}{k},$ (3) where h is the convection heat transfer coefficient, which is expressed as shown in Eq. 4: $h=\frac{q}{t_{\text{w}}-t_{\text{b}}}.$ (4)

Here, t_b is the bulk temperature, which can be calculated using Eq. 5 as follows: $t_{\text{b}}=\frac{\int Cp\rho ut\text{d}A}{mCp}.$ (5)

Model validation

The model was validated by comparing the simulation and experimental results of Taylor et al. ^[11]. A schematic diagram of the experimental setup is shown in Fig. 3; the adiabatic and heating sections of the tube were 328.72 and 352.20 mm long, respectively, and the inner diameter of the tube was 5.87 mm.

Fig. 3Full-screen View

(Color online) Schematic diagram of Taylor’s experimental section

To verify the simulation results, data from two runs were selected: Run-689H, which featured gas Pr and Re ranges closest to those in the current study; and Run-715H, which featured gas mole mass and temperature ranges closest to those in the current study. The turbulence model, solver, and Pr_t model used in the simulation were consistent with those described in Sect. 2. The boundary conditions for the two experiments are listed in Table 4, and a comparison of the simulation and experimental results is presented in Fig. 4.

Fig. 4Full-screen View

(Color online) Comparison between simulation and experimental results. (a) Run-689H; (b) Run-715H. 

As shown in Fig. 4, both experiments showed a significant decrease in the wall temperatures near the outlet, which is likely due to the axial heat conduction of the tube wall[22]. Under two separate operating conditions, the relative error between the experimental and simulated data was less than 5%. The simulated Run-689H and Run-715H exhibited maximum relative errors of 4.4% and 4.2%, respectively, whereas their average relative errors were 2.8% and 1.9%, respectively, except for two points situated near the outlet.

In addition, both the aforementioned model and the core channel model satisfied the criteria for grid independence. Figure 5 illustrates the grid partitioning of the core channel model along the axial direction. The grid was refined in the vicinity of the wall within the fluid region to ensure that the center of the first layer of the grid cells was positioned within the viscous sublayer.

Fig. 5Full-screen View

(Color online) Detailed grid view of the core flow channel model.

Flow and temperature field of standard condition

To analyze the flow and heat transfer characteristics of the helium–xenon gas mixture in the core environment, one must examine the velocity and temperature distributions in both the flow channel and the surrounding solid region. Figure 6 shows the radial velocity distribution at various axial positions within the flow channel.

Fig. 6Full-screen View

(Color online) Flow velocity distribution at different axial positions.

Based on Fig. 6, the fluid flowed gradually from the inlet to the outlet and the gas expanded owing to the heating of the fuel elements, thus resulting in a velocity increase. Under normal operating conditions, the maximum local velocity in the flow channel reached approximately 199 m/s, which corresponded to approximately 20% of the local speed of sound. Additionally, as shown in Fig. 6, a significant velocity gradient appeared in the near-wall region. This is attributed to the fluid experiencing resistance from the adjacent wall, which resulted in a lower velocity compared with that at the mainstream flow region.

Figure 7 presents the temperature distribution in both the flow channel and the surrounding solid region. Figures 7(a)–(c) illustrate the radial temperature distributions at different axial positions, whereas Fig.7(d) depicts the axial temperature distributions of all four regions.

Fig. 7Full-screen View

Temperature distribution in solid and fluid regions

Based on Fig.7, the fuel region demonstrates unsatisfactory thermal uniformity, primarily because of its relatively low thermal conductivity. By contrast, the moderator and cladding regions exhibited better thermal uniformity. Additionally, a significant temperature gradient was observed in the near-wall region of the fluid, where the cladding temperature significantly exceeded that of the helium–xenon gas mixture.

The axial distributions of the radial average temperature for each solid and fluid region is illustrated in Fig.8. As shown, the temperature peak points of the three solid regions were located in the latter half of the flow channel. Furthermore, the proximity of a region to the fluid region corresponded to the proximity of its peak temperature point to the outlet. This phenomenon arises because of the gradual increase in fluid temperature along the flow channel and the cosine distribution of power in the axial direction. Additionally, the helium–xenon gas mixture experienced a rapid temperature increase in the middle of the flow channel and a slower temperature increase at both ends, which is attributable to the cosine distribution of the axial power.

Fig. 8Full-screen View

(Color online) Axial distribution of radial mean temperature in different regions

Effects of operating parameters on Nu

3.3.1

Power distribution

To render the current study more applicable to practical operating conditions, the axial power distribution for four different operating conditions, namely Run-D1 to Run-D4, was designed based on the actual configuration of the axial reflector in the reactor core and previous heat transfer experiments. Table 3 presents the boundary conditions for each operating condition and Fig. 9 shows the normalized power distributions.

Fig. 9Full-screen View

(Color online) Power distribution under different conditions

Figure 9 shows that the addition of the reflector significantly increased the normalized wall heat flux density near the inlet of Run-D1 and outlet of Run-D2. Additionally, although the normalized wall heat flux distributions for Run-D3 and Run-S were almost identical, the former demonstrated a significantly higher wall heat flux density near the inlet and outlet of the channel than the latter.

Figure 10 shows the axial distribution of Nu under different operating conditions. Table 5 lists the heat transfer characteristic parameters for each operating condition. To analyze the effects of different factors on the axial change rate of Nu, ΔNu was defined as 10% of the average Nu obtained in each run. Additionally, two dimensionless distances were introduced: γ, which represents the distance required downstream of the first measuring point to observe a decrease of ΔNu, and δ, which represents the distance required upstream of the last measuring point to observe an increase of ΔNu. These values are also marked in Fig. 10.

Table 5Full-screen View

Characteristic parameters under different axial power distributions

Run No.	Nu			γ (z/L)	δ (z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	16.76	4.22
D1	93.74	48.13	71.59	9.57	5.16
D2	88.92	56.34	72.96	23.78	17.20
D3	92.04	56.60	72.62	11.97	17.34
D4	89.34	66.88	72.59	5.81	82.96

Show more

Fig. 10Full-screen View

(Color online) Axial Nu distribution under different axial power distributions

As shown in Fig. 10, Nu decreased along the flow direction, and the decreasing trend was consistent across all operating conditions, except for the case of uniform power distribution. In the case of uniform power distribution, Nu decreased significantly near the inlet, whereas the rate of decrease at other positions decelerated considerably. Similar phenomena were reported by Sparrow[34] and Vitovsky[16].

Table 5 indicates that the addition of a reflector at the inlet increased the axial decline rate of Nu near the inlet. By contrast, incorporating a reflector at the outlet can reduce the axial decline rate of Nu, thus improving Nu near the outlet and contributing positively to heat transfer. Hence, incorporating reflectors near the core outlet can potentially provide supplementary thermal safety margins for the core. However, in the case of microreactors, one must scrutinize the tradeoff between the benefits of enhanced thermal safety and the associated increase in core volume and weight resulting from the addition of reflectors.

3.3.2

Power value

To investigate the effect of power on the channel heat transfer, four sets of operating conditions, namely Run-Q1 to Run-Q4, with different power values were designed, and their Nu distributions were compared with those of Run-S. The boundary conditions for each operating condition are listed in Table 3. The axial distribution of Nu and the heat transfer characteristic parameters for each operating condition are presented in Fig. 11 and Table 6, respectively.

Table 6Full-screen View

Characteristic parameters under different power values

Run No.	Nu			δ (z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	4.22
Q1	78.52	39.97	69.96	3.18
Q2	85.10	44.34	72.87	3.71
Q3	89.03	44.96	72.70	4.04
Q4	94.22	43.60	70.51	4.44

Show more

Fig. 11Full-screen View

(Color online) Axial Nu distribution under different power values

As shown in Fig. 7, power exerted a positive effect on Nu near the inlet. However, as the power increased, the effect weakened. Furthermore, the Nu in the axial middle section of the channel decreased gradually as the total power increased. This can be attributed to the lower fluid temperature at lower power levels, which results in a lower viscosity and higher Re under similar flow velocities. Based on the preceding discussion, one can infer that under certain circumstances, such as low-power operation or transient conditions, the phenomenon in which the cladding temperature at the reactor inlet exceeds the expected value must be considered.

In terms of the rate of change of Nu, the distribution of Nu under certain working conditions exhibited a significant nonlinear trend. Hence, only the rate of change of Nu near the outlet is presented in Table 6. As shown in the table, the power level did not significant affect the rate of decrease in Nu near the outlet.

3.3.3

Inlet temperature

Operating conditions Run-T1 to Run-T4 were imposed to investigate the effect of the inlet temperature on the axial distribution of Nu. The boundary conditions, axial distribution of Nu, and heat transfer characteristic parameters for each case are presented in Table 3, Fig. 12, and Table 7, respectively.

Table 7Full-screen View

Characteristic parameters under different inlet temperatures

Run No.	Nu			γ (z/L)	δ (z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	16.76	4.22
T1	104.67	50.75	82.61	20.79	4.20
T2	97.48	47.52	76.85	20.13	4.24
T3	85.90	41.90	67.25	18.68	4.31
T4	81.01	39.55	63.27	17.99	4.33

Show more

Fig. 12Full-screen View

(Color online) Axial Nu distribution under different  inlet temperatures 

Based on the results presented in Fig. 12 and Table 7, as the inlet temperature decreased, the overall Nu in the channel increased, and the rate of increase in the average Nu was approximately –0.094/K. This phenomenon can be attributed to a decrease in the inlet temperature, which resulted to an increase in the fluid density and a decrease in the viscosity, thus causing an increase in Re.

Although reducing the inlet temperature of the core can improve the overall heat transfer efficiency, it may decrease the outlet temperature, which consequently affects the thermal-electric conversion efficiency of the system.

As shown in Table 7, the effect of the inlet temperature on the rate of change of Nu near the inlet and outlet was insignificant.

3.3.4

Inlet velocity

The results of the simulations for Runs-U1 to Run-U4 can be used to investigate the effect of the flow velocity on the axial distribution of Nu based on a comparison with the results of Run-S. The boundary conditions, axial distribution of Nu, and heat transfer characteristic parameters for each case are presented in Table 3, Fig. 13, and Table 8, respectively.

Table 8Full-screen View

Characteristic parameters under different inlet velocities

Run No.	Nu			γ (z/L)	δ (z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	16.76	4.22
U1	78.65	36.76	56.80	9.47	4.94
U2	85.89	40.95	64.46	12.85	4.62
U3	97.63	47.25	78.62	30.22	3.95
U4	103.29	49.14	84.89	47.34	3.67

Show more

Fig. 13Full-screen View

(Color online) Axial Nu distribution under different inlet velocities

Based on Fig. 13 and Table 8, Nu increased with the flow velocity in the axial middle section of the channel at an average increase rate of 0.468 s/m. This is because an increase in flow velocity improves the mass flow rate and reduces the average fluid temperature, thus resulting in an overall increase in Re.

Furthermore, as the velocity increased, the positive effect of high velocity on Nu near the inlet weakened gradually, thus resulting in a lower axial decrease rate of Nu near the inlet. Additionally, the effect of the inlet velocity on the axial change rate of Nu near the outlet was insignificant.

Based on the preceding discussion, one can infer that increasing the inlet flow velocity can enhance the overall heat transfer efficiency of the core. However, the rate of improvement in the heat transfer efficiency at the inlet may not be as significant as that in other regions. Consequently, this increase in the flow velocity can result in higher cladding temperatures at the inlet than anticipated.

3.3.5

Outlet pressure

The results of Run-P1 to Run-P4 were compared with those of Run S to investigate the effect of the outlet pressure on the axial distribution of Nu. During the modification of the outlet pressure, the physical properties of the helium–xenon gas mixture were adjusted accordingly. The boundary conditions, axial distribution of Nu, and heat transfer characteristic parameters for each case are presented in Table 3, Fig. 14, and Table 9, respectively.

Table 9Full-screen View

Characteristic parameters under different outlet pressures

Run No.	Nu			γ (z/L)	δ (z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	16.76	4.22
P1	75.71	35.39	57.50	16.10	4.38
P2	84.55	39.93	64.81	16.48	4.32
P3	98.32	48.89	78.51	21.99	4.22
P4	105.38	53.03	84.96	24.58	4.18

Show more

Fig. 14Full-screen View

(Color online) Axial Nu distribution under different outlet pressures

As shown in Fig. 14 and Table 9, Nu increased with the outlet pressure, and the average Nu increased at a rate of approximately 0.031/kPa. The increase in Nu can be attributed to the increase in fluid density caused by the increase in pressure, which results in an overall increase in Re and an improved heat transfer performance.

Caution must be exercised when increasing the inlet pressure to enhance the heat transfer efficiency in the core because it can affect the volume and mass of the compressor, thus potentially affecting the flexibility of the micronuclear reactor power source.

As shown in Table 9, the effect of the outlet pressure on the axial change rate of Nu near the inlet and outlet was insignificant.

3.3.6

Inlet velocity under fixed inlet and outlet temperatures

One approach to determine the operating range of a closed Brayton cycle nuclear power system is to vary the core power while maintaining a constant temperature increase in the core and then observe changes in parameters such as the system thermoelectric conversion efficiency. Using this method, operating conditions Run-UQ1 to Run-UQ4 were established by adjusting the total power to maintain a stable outlet temperature based on operating conditions Run-U1 to Run-U4. The simulation results of Run-UQ1 to Run-UQ4 were compared with those of Run-S to evaluate the effects of varying the total power. The boundary conditions, axial distribution of Nu, and heat transfer characteristic parameters for each case are presented in Table 3, Fig. 15, and Table 10, respectively.

Table 10Full-screen View

Characteristic parameters under different outlet pressures with fixed inlet/outlet temperatures

Run No.	Nu			γ (z/L)	δ(z/L)
	Maximum	Minimum	Average
S	92.68	44.20	71.73	16.76	4.22
UQ1	76.70	38.19	58.60	13.20	4.72
UQ2	84.63	41.49	65.14	15.41	4.52
UQ3	98.67	47.10	78.25	24.55	4.04
UQ4	106.78	49.19	84.43	27.86	3.87

Show more

Fig. 15Full-screen View

(Color online) Axial Nu distribution under different inlet velocities with fixed inlet/outlet temperatures

As shown in Fig. 15, Nu increased with the inlet velocity, and the average Nu increased at a rate of approximately 0.431 s/m, which was lower than the result presented in Sect. 3.3.4.

As shown in Table 10, similar to the conclusion in Sect. 3.3.4, as the inlet velocity increased, the positive effect of the high velocity on Nu near the inlet gradually weakened. However, under the operating conditions investigated in this section, this weakening trend became less pronounced owing to the change in the power. Furthermore, the inlet velocity did not significantly affect the decreasing rate of Nu near the outlet.

Theoretical derivation

To depict the axial variation in Nu in the channel in the core environment precisely, a theoretical framework that can capture the changing behavior of the thermal boundary layer is necessary. Sparrow et al. presented a theoretical derivation based on the theory of turbulent boundary layers and eddy viscosity. According to their theory, the wall-to-bulk temperature difference along the axial direction of a circuit channel under turbulent heat transfer conditions of forced convection with a uniform axial power distribution can be expressed as shown in Eq. 6[34]. $t_{\text{w}}-t_{\text{b}}=\frac{qL}{2k}\left\{G\left(r_0^{+}\right)+{\displaystyle \sum _{n=0}^{\infty }{F}_n}\exp \left[-\frac{4{\beta }_n^2}{Re}\frac{z}{L}\right]\right\},$ (6) where $G\left(r_0^{+}\right)$ denotes the dimensionless wall-to-bulk temperature difference in the fully developed temperature field; F_n is the coefficient of the series expansion of the wall temperature, which is a negative value that depends only on Re, Pr, and G; ${\beta }_n^2$ is the eigenvalue and is related only to Re and Pr. The derivation of Eq. 6 is based on the following six assumptions:

I. The fluid properties are constant.

II. Compared with the radial heat diffusion, the axial diffusion of both molecular and turbulent heat is negligible.

III. The mean value of the radial velocity is 0.

IV. Viscous dissipation is negligible.

V. The turbulence velocity distribution is fully developed throughout the channel.

VI. The turbulent Pr can be approximated as 1.

This equation can be used to precisely calculate the distribution of the axial wall-to-bulk temperature difference, including that at the thermal entrance region. Additionally, $G\left(r_0^{+}\right)$ can be expressed as $G\left(r_0^{+}\right)=-{\displaystyle \sum _{n=0}^{\infty }{F}_n.}$ (7)

Based on Eq. 6, when z = 0 (i.e., at the entrance), the wall-to-bulk temperature difference is 0; when z = ∞ (i.e., where the temperature is fully developed), the temperature difference is a constant.

For operational conditions with an arbitrary axial power distribution, the power distribution can be regarded as a combination of multiple power steps, with each power step affecting the downstream temperature distribution[35]. Based on Eq. 6, if no heat flux appears on the wall before position z and a uniform heat flux of dq appears after z, then the temperature difference downstream $\tilde{z}$ can be calculated as $t_{\text{w}}-{t_{\text{b}}|}_{\tilde{z}}=dq\cdot \frac{L}{2k}\left\{G\left(r_0^{+}\right)+{\displaystyle \sum _{n=0}^{\infty }{F}_n}\exp \left[-\frac{4{\beta }_n^2}{Re}\frac{\tilde{z}-z}{L}\right]\right\} \left(z\le \tilde{z}\right)$ (8)

Integrating the equation above yields the following expression for the wall-to-bulk temperature difference at $\tilde{z}$ under an arbitrary axial power distribution[35]: $t_{\text{w}}-{t_{\text{b}}|}_{\tilde{z}}=-\frac{2}{kRe}{\displaystyle \sum _{n=0}^{\infty }{F}_n}{\beta }_n^2\underset{0}{\overset{\tilde{z}}{{\displaystyle \int }}}q\left(z\right)\text{exp}\left[-\frac{4{\beta }_n^2}{\text{Re}}\frac{\tilde{z}-z}{L}\right]\text{d}z.$ (9)

Without axial reflectors, the axial power distribution in the core can be approximated using a cosine function as follows: $q\left(z\right)=\alpha \text{cos}\left[\pi \left(z-\frac{1}{2}\right)\right].$ (10)

After simplification, it can be expressed as $q\left(z\right)=\alpha \text{sin}\left(\pi z\right).$ (11)

By substituting Eq. 11 into Eq. 9, the equation for the wall-to-bulk temperature difference at $\tilde{z}$ under a cosine power distribution can be obtained as follows: $t_{\text{w}}-{t_{\text{b}}|}_{\tilde{z}}=-\frac{2\alpha }{k\text{Re}}{\displaystyle \sum _{n=0}^{\infty }{F}_n}{\beta }_n^2\frac{{{\displaystyle \int }}_0^{\tilde{z}}\text{sin}\left(\pi z\right)\text{exp}\left[\frac{4{\beta }_n^2}{Re}\frac{z}{L}\right]\text{d}z}{\exp \left[\frac{4{\beta }_n^2}{Re}\frac{\tilde{z}}{L}\right]}.$ (12)

Let $a_n=\frac{4{\beta }_n^2}{ReL}$. Therefore, the temperature difference at $\tilde{z}$ can be expressed as $t_{\text{w}}-{t_{\text{b}}|}_{\tilde{z}}=\frac{2\alpha }{k\text{Re}}{\displaystyle \sum _{n=0}^{\infty }{F}_n}{\beta }_n^2\frac{\pi \text{cos}\left(\pi \tilde{z}\right)-\frac{\pi }{\exp \left(a_n\tilde{z}\right)}-a_n\text{sin}\left(\pi \tilde{z}\right)}{a_n^2+{\pi }^2}.$ (13)

The axial distribution of Nu^–1 under the cosine power distribution is represented as shown in Eq. 14. $\begin{array}{c}N{u}^{-1}=\frac{t_{\text{w}}-t_{\text{b}}}{q\left(\tilde{z}\right)}\frac{k}{L}=\frac{1}{2}{\displaystyle \sum _{n=0}^{\infty }\frac{F_n{a}_n}{a_n^2+{\pi }^2}}\\ \left\{\pi \left[\text{cot}\left(\pi \tilde{z}\right)-\frac{1}{\text{sin}\left(\pi \tilde{z}\right)\text{exp}\left(a_n\tilde{z}\right)}\right]-a_n\right\}\end{array}$ (14)

In Eq. 14, when $\tilde{z}$ approaches 1 m from the left side, the following equation is satisfied. $\underset{\tilde{z}\to 1^{-}}{\text{lim}}\left\{\pi \left[\text{cot}\left(\pi \tilde{z}\right)-\frac{1}{\text{sin}\left(\pi \tilde{z}\right)\text{exp}\left(a_n\tilde{z}\right)}\right]-a_n\right\}=-\infty$ (15)

Thus, in the region near the channel outlet, Nu^–1 approaches infinity, whereas Nu approaches zero (a_n > 0, F_n < 0). As $\tilde{z}$ approaches 0 m from the right side, it satisfies the following equation: $\underset{\tilde{z}\to 0^{+}}{\text{lim}}\left\{\pi \left[\text{cot}\left(\pi \tilde{z}\right)-\frac{1}{\text{sin}\left(\pi \tilde{z}\right)\text{exp}\left(a_n\tilde{z}\right)}\right]-a_n\right\}=0.$ (16)

Hence, near the inlet, Nu^–1 approaches 0, whereas Nu approaches infinity.

As shown in Eq. 14, F_n and a_n are critical parameters that affect the Nu. In Sparrow and Siegel’s series of studies[34, 35], the values of F_n and a_n were obtained by solving the eigenvalue system of the Sturm–Liouville type, which is based on the six assumptions mentioned above. In this study, the high flow velocity of the helium–xenon gas mixture renders viscous dissipation a significant factor that affects the Nu distribution under certain operating conditions. Moreover, the gas is heated vigorously in the channel, thus causing changes in the physical properties and a pronounced acceleration effect. This renders it challenging to satisfy the conditions of constant properties and a fully developed turbulent velocity distribution. Therefore, the six assumptions mentioned earlier may not be entirely valid, whereas solving the eigenvalue equation to determine F_n and a_n may introduce a significant deviation.

In addition, coefficients F_n and a_n in Eq. (14) are related to the characteristic values, which result in an infinite number of unknown coefficients that must be determined, thereby rendering the equation impractical for engineering applications. To reduce the number of unknown coefficients, F_n and a_n in Eq. 14 were simplified as φ and ω, respectively, and existing simulation data were used to determine the values of φ and ω for various operating conditions. The subsequent analysis indicates that this simplification affected the accuracy of the relationship at a reasonable level. Hence, Eq. 14 can be simplified as follows: $Nu=2\frac{{\omega }^2+{\pi }^2}{\phi \omega }{\left\{\pi \left[\text{cot}\left(\pi \tilde{z}\right)-\frac{1}{\text{sin}\left(\pi \tilde{z}\right)\text{exp}\left(\omega \tilde{z}\right)}\right]-\omega \right\}}^{-1}.$ (17)

Using the MATLAB Genetic Algorithm Toolbox, the axial distribution of Nu for all operating conditions in Sect. 3, which presents axial power distributions in the form of a cosine function (excluding Run-Q1 and Run-Q2, which exhibited clear nonlinear trends), was fitted using Eq. 17. The variations in -φ and ω with the average Re of each operating condition are presented in Fig. 16. The average Re of each condition was calculated as follows: $R{e}_{\text{avg}}=\frac{R{e}_{\text{in}}+R{e}_{\text{out}}}{2}=\frac{2m\left({\mu }_{\text{in}}+{\mu }_{\text{out}}\right)}{\pi L{\mu }_{\text{in}}{\mu }_{\text{out}}}$ (18)

Fig. 16Full-screen View

Relationship between the coefficients and average Re (dashed line represents the function graph of the fitted equation, and the average Re is within the range of 5.3 × 10⁴ to 1 × 10⁵)

When Pr = 0.264 and the average Re was between 5.3 × 10⁴ to 1 × 10⁵, the relationship between φ and the average Re, as well as that between ω and the average Re, were fitted using power functions. The results are expressed as shown in Eqs. 19 and 20. $\phi =-90.72R{e}_{\text{avg}}^{-0.72}$ (19) $\omega =1075.65R{e}_{\text{avg}}^{-0.31}$ (20)

Here, the Nu at any axial position $\tilde{z}$ in the flow channel can be calculated conveniently using Eq. 17. In the following section, $\tilde{z}$ is uniformly substituted with z in the relevant correlations to facilitate understanding.

Verification and analysis

As shown in Fig. 17, the axial Nu distribution calculated using Eq. 17 was compared with the simulation results for all the operating conditions listed in Table 3, where the axial power distribution was cosine (excluding Run-Q1 and Run-Q2).

Fig. 17Full-screen View

Comparison of calculations using Eq. 17 with simulation results (relative error within red and black dashed lines is less than 5% and 10%, respectively)

As shown in Fig. 17, the most significant difference was indicated in the region near the inlet, with a maximum relative error of 75.3%, whereas the average relative error was approximately 5.2%. Approximately 94% of the data points showed a relative error within 10%, of which and 78% indicated a relative error within 5%. Moreover, in the latter half of the channel (z/L ≥ 62.5), almost all the data points showed a relative error within 5%.

Table 11 presents the boundary conditions of the four newly established cases, which are distinct from those in Table 3, used to assess the validity of Eq. 17.

Based on the calculations performed using Eq. 18, the average Re of conditions Run-V1 to V4 were approximately 5.9 × 10⁴, 6.9 × 10⁴, 8.6 × 10⁴, and 9.3 × 10⁴, respectively.

To evaluate the accuracy of Eq. 17, the calculation error and the errors of several existing and widely used correlations were compared. Previous studies[18, 21, 11] showed that the Dittus–Boelter[46], Churchill[47], Kays[48], and Pickett[49] methods are more accurate for calculating the Nu of a helium–xenon gas mixture compared with other correlations. The expressions and applicabilities of the four relationships are listed in Table 12. The accuracies of the four methods above and Eq. (17) were compared, as shown in Fig. 18. Notably, the Nu values obtained from the correlations in Table 12 were calculated using the local Re, whereas used Re_avg was used in Eq. 17 for the calculations.

Fig. 18Full-screen View

(Color online) Comparison of calculations using different equations based on simulation results of Run-V1 to Run-V4 (relative errors within red and black dashed lines are less than 5% and 10%, respectively).  (a) Run-V1, (b) Run-V2, (c) Run-V3, (d) Run-V4

As shown in Fig. 18, Eq. 17 demonstrated high accuracy in most axial positions, with a lower relative error compared with the Dittus–Boelter, Churchill, Kays, and Pickett correlations. Most of the observations indicated a relative error of less than 10%, except for two points near the inlet. In the latter half of the duct, almost all data points indicated a relative deviation of less than 5%. Table 13 lists the maximum and average relative errors of each correlation for all the data points under the four operating conditions listed in Table 11.

As presented in Table 13, the Dittus–Boelter correlation showed the highest maximum and average relative errors among the five correlations, thus indicating that its predictive performance was the least ideal. However, the accuracy of Eq. 17 was the lowest, with only approximately one-half of the error indicated in the Pickett correlation, thus rendering it the most accurate among the five correlations. Furthermore, based on the calculation results for all four operating conditions, the Kays correlation yielded the most accurate prediction for the entrance area (z/L ≤ 18.75), with an average relative error of 4.5%. In the outlet area (z/L ≥ 106.25), Eq. 17 yielded the highest prediction accuracy, with an average relative error of 2.2%. The results above validated the theoretical derivation and confirmed the feasibility of Eq. 17 for turbulent heat transfer calculations in a helium–xenon gas mixture in a core environment.

Furthermore, Eq. 17 showed high accuracies for calculations near the outlet but exhibited low accuracy near the inlet, in contrast to other correlations. The improved accuracy of Eq. 17, particularly near the outlet, can be attributed to the inclusion of the thermal boundary layer changes caused by power fluctuations, which were not considered in the other correlations. As the power decreases rapidly near the outlet, using a correlation that disregards changes in the thermal boundary layer may result in significant distortions in the calculations. Furthermore, consistent with the observed low accuracy of the calculations near the inlet, previous studies by Sparrow et al.[34] showed similar challenges in achieving accurate calculations in this region. They attributed the cause partially to changes in the fluid properties and the uncertainty of eddy diffusivities. However, in the present study, in addition to the reasons above, the simplification of formulas in the derivation process and the viscous dissipation caused by the high-speed fluid may result in inaccurate predictions near the inlet.

The preceding discussion highlights that although Eq. 17 provides satisfactory accuracy for most axial regions in the channel, further improvements can be realized. In this regard, the Kays correlation can be combined with Eq. 17 to derive a new segmented relationship, as shown in Eq. 21. To evaluate the effectiveness of the approach, Eq. 21 was applied to calculate the four operating conditions specified in Table 11, and the resulting calculation accuracy is presented in Fig. 19. $\begin{array}{l}Nu=\{\begin{array}{l}0.022{\mathrm{Re}}^{0.8}P{r}^{0.6} \left(\frac{z}{L}\le 18.75\right)\hfill \\ 2\frac{{\omega }^2+{\pi }^2}{\phi \omega }{\left\{\pi \left[\cot \left(\pi z\right)-\frac{1}{\sin \left(\pi z\right)\exp \left(\omega z\right)}\right]-\omega \right\}}^{-1} \left(\frac{z}{L}>18.75\right)\hfill \end{array}\\ \phi =-90.72R{e}_{\text{avg}}^{-0.72} \omega =1075.65R{e}_{\text{avg}}^{-0.31}\end{array}$ (21)

Fig. 19Full-screen View

(Color online) Comparison of Nu between simulation data and Eq. 21 (relative errors within red and black dashed lines are less than 5% and 10%, respectively) 

Equation 21 employs the Kays correlation to calculate Nu in the vicinity of the inlet (z/L ≤ 18.75). In this correlation, Re can be approximated using the principle of energy conservation. The Nu values in the remaining regions of the channel were calculated using Eq. 17. The data presented in Fig. 19 show a significant improvement in the prediction accuracy achieved via the optimized correlation (Eq. 21). Specifically, the maximum and average relative errors reduced to 13.3% and 2.9%, respectively, while a satisfactory level of accuracy was maintained throughout the channel axial range.