The accurate prediction of supercritical water behavior especially heat transfer requires not just high mesh fidelity but also sophisticated turbulence modeling that captures the evolving and often anisotropic nature of turbulent structures under pseudo-critical conditions. Traditional turbulence and heat-transfer models, calibrated on uniform or simple geometries, tend to underpredict heat transfer in rod bundles. According to Rowinski et al. (2017) non-uniform heating or buoyancy effects are strong conditions that become especially severe near the pseudo-critical point. As observed by Palko and Anglart (2008), empirical models fail to account for complex axial and radial property gradients, which dominate supercritical water behavior and influence the onset of heat transfer deterioration. They emphasized the need for geometry-sensitive modeling approaches and localized correction factors, especially in tight lattice configurations. SST k-ω models, which combine the robustness of k-ε in the bulk flow with the wall-resolving capacity of k-ω, have shown good performance in simulating near-wall heat transfer and capturing laminarization effects near deteriorated zones. Podila and Rao (2016) found the SST model reliable for low Reynolds number flows in heated bundles, especially near pseudo-critical conditions where viscosity increases sharply. Reynolds Stress Models (RSMs) which are associated with significantly higher computational costs and require fine-tuned wall functions have been used to resolve three-dimensional turbulence anisotropy and capture core flow structures more effectively than eddy-viscosity models. Hybrid turbulence models, such as those developed by Zhang et al. (2024) combine the SST and RSM approaches to retain near-wall accuracy while resolving large-scale turbulent structures in bundle cores. Their application to wire-wrapped 2×2 bundles resulted in improved predictions of both wall temperature distributions and onset of heat transfer deterioration thresholds. Despite advances in turbulence modeling, heat transfer remains challenging to predict due to its sensitivity to temperature-dependent properties, flow acceleration, and thermal stratification. In et al. (2015) stressed the importance of using real-fluid property libraries with fine temperature steps to capture the sharp changes in specific heat, density, and viscosity near the pseudo-critical point.

The choice of an appropriate turbulence model is crucial when modeling flow of supercritical water in rod bundles due to the presence of sharp property gradients and flow instabilities near the pseudo-critical point. Zou et al. (2025) concluded that Reynolds Stress Transport (RST) provided the most accurate prediction of heat transfer when the wall temperature crossed the pseudo-critical threshold. They conducted a comparative study using five turbulence models including SST k‑ω and Reynolds Stress Transport models. The SST k‑ω model, though more stable, underpredicted buoyancy and turbulent fluctuations in near-wall regions under low mass flux conditions. Xiong et al. (2015) also analyzed the performance of advanced turbulence models such as BSL-RSM and ω-RSM in a 2×2 rod bundle. Their study confirmed that models incorporating full Reynolds stress transport better captured the wall heat transfer behavior in regions experiencing flow acceleration and pseudo-critical gradients. Similarly, Han and Zhang (2023) emphasized the need for incorporating buoyancy models to accurately resolve wall temperature and turbulence fields in Canadian SCWR geometries. In a related study, Tripathi and Basu (2020) observed that enhanced near-wall turbulence helped delay the onset of heat transfer deterioration, reinforcing the role of turbulence modeling in mitigating flow-induced instability. Silas (2018), using STAR‑CCM+ in a square rod bundle simulation, found that low-Reynolds SST k‑ω model with detailed near-wall treatment yielded the best compromise between convergence and physical accuracy, particularly for simulations near the pseudo-critical region. In the present study, numerical simulation of heat transfer for flow of supercritical water in a 2×2 rod bundle was conducted to determine the extent of reproducibility of the experimental data obtained at the SWAMUP-II test facility at Shanghai Jiao Tong University and to also characterize thermal hydraulics parameters.

The basis for the present research activity was the experimental study of flow of supercritical water through a 2 × 2 rod bundle test section with two channels, carried out at the Shanghai Jiao Tong University by Gu et al. (2015). The bundle consisted of four heated rods with an Outer Diameter (OD) of 10 mm and a pitch-to-diameter ratio of 1.18 and a heated length of 833 mm. The 2×2 rod bundle was inserted into a square assembly box (25.4 mm inner width, 2 mm wall thickness) with rounded corners. The test section was separated into two channels. A downward flow of water was enabled in the first channel between the pressure tube and the assembly box which then turns upward in the second channel inside the assembly box to cool the four heated rods, which are directly heated by the supply of a DC Power. The SWAMUP-II test facility at Shanghai Jiao Tong University is shown in Fig. 2.

The SWAMUP-II test facility consists of a main loop, a cooling loop, a purification loop, and an I&C system. The main loop has a plunger pump, a pre-heater, a re-heater, a mixing chamber, a heat exchanger and the test section. This loop is designed for pressure of up to 35 MPa, temperature of up to 550 °C, and mass flow rate of up to 2.8 kg/s. Every rod has four thermocouples. The inlet and outlet fluid temperatures are measured by the sheathed thermocouples positioned at the inlet and outlet of the test section, respectively. Fluid temperature at the inlet of the second channel is measured by the sheathed thermocouple in the mixing chamber at the bottom of the test section. Table 2 shows the parameters of the experimental test section and Table 3 presents the experimental parametric range with inlet temperature varying (200 °C to 420 °C) such that the outlet temperature range (350 °C to 500 °C) covers the pseudocritical values where unusual thermal hydraulic behavior is expected.

Pressure, mass flow rate and power to the heater rods were kept constant for each experimental case. The fluid temperature at the inlet of the test section was increased gradually by controlling the heating power of the pre-heater. After the steady state condition was reached, pressure at the inlet, mass flow rate, power to the heated rods, fluid temperature and inner wall temperature data were recorded at each step. The test ended when the wall temperature at any location exceeds 580 °C or when the inlet temperature is higher than 420 °C. The pressure in the test loop was controlled by the pressure regulator valve at the outlet of the loop once the piston pump is operating. The mass flow rate was adjusted by the control valves in the bypass loop. The power to the heating rods was increased step by step to the expected value.

STARCCM+ was employed to model and discretize the portion of the experimental test session relevant to the present study. Implicit unsteady model for the computation of unsteady flows was employed to account for time dependent flow and thermal fluctuations that are expected under the conditions specified. Reynolds-Averaged Navier-Stokes (RANS) equations which describe the conservation of mass, momentum, and energy for Newtonian fluids adopted in the solution scheme of STARCCM+ were used. For the flow solution, the 3-dimensional momentum conservation equation were solved together with the continuity and energy equations expressed by Han and Zhang (2023) as follows:

Continuity equation

∂ρu¯i∂xi=0 (Eq. 1)

Momentum conservation equation

∂∂xj(ρuiuj―)=−∂p¯∂xi+∂∂xj(μ∂u¯i∂xj−ρui′uj′―)+ρgi (Eq. 2)

Energy conservation equation

∂∂xi(ui′ρcpT)=∂∂xi[(λ+cpμtPrt)∂T∂xi]+ϕ (Eq. 3)

where, u is velocity, ρ is density, μ is viscosity, C_p fluid heat capacity at constant pressure,T is temperature, λ is molecular thermal conductivity and Pr_t is turbulent Prandtl number.

The heat transfer coefficient behavior was experimentally determined by Gu et al (2015) in the heated rods based on the cross-sectional average parameters in the channel. The local heat transfer coefficient α was defined as the ratio of the heat flux (q) in kW/m² to the difference between the wall outer temperature T_w and the bulk fluid T_b and expressed as,

α=qTw−Tb (Eq. 7)

The electrical heat flux supplied was obtained by dividing the power (P) by the heated area and given as,

q=P4πDL=UI4πDL (Eq.8)

U and I are the electrical voltage and current respectively.

In the present numerical analysis, the heat transfer from the fluid to the assembly box and also from the fuel rod to the fluid were determined using numerical simulation data obtained. The rod to fluid heat transfer coefficient was determined as,

α=qTRW−TBF (Eq. 9)

where, q is the heat flux in kW/m², T_BF is the bulk fluid temperature in °C obtained over constrained planes positioned along the length of the flow domain and T_RW is the fuel rod wall surface temperature in °C taken along the length of the domain. The fluid to assembly box average heat transfer coefficients were also determined in the present work using the expression,

α=qTBF−TTW (Eq. 10)

where, q is the heat flux in kW/m², T_BF is the bulk fluid temperature in °C and T_RW is the tube (assembly box) wall surface averaged temperature in °C.

To select a working mesh that supports the available computational resource and timely convergence, plots of wall and flow temperature profiles at a mass flux of 228 kg/m²s were compared for the base, coarse and fine meshes developed and are presented in Figs 12, 13.

The wall temperature profiles in Fig. 12 were obtained over 16-line probes (4 per rod), each of which runs from the inlet to the outlet with 17 nodes and positioned on the wall surfaces of the 4 rods as shown in Fig. 11. The average temperature value at each position for the 16 rods is plotted for the three (3) mesh sizes considered. While the data in Fig. 12 shows some deviation in the wall temperature prediction for the mesh sizes, such deviations are totally not present for the bulk fluid temperature profile. The bulk flow parameters are surface averaged values obtained over constrained planes positioned in the bulk flow domain from inlet to outlet whiles the wall temperature data are not surface averaged value but taken at specific points on the wall boundary. The observed deviations in the wall temperature graphs could be attributed to the averaging of point-specific values along the line probes instead of a surface average value obtained over the defined wall surface. Data obtained over the base mesh are much closer to those obtained over the fine mesh than the coarse mesh and the deviations of the wall surface average parameters are observed to be negligible and therefore insensitive to the mesh refinement than seen in the point-specific averaging on the line probes. Jin et al. (2024) and Kukanov et al. (2024) confirmed that base meshes capturing local wall temperature distributions within ±0.2% of fine-mesh results are sufficient for SCW simulations.

Hence the base mesh was chosen for use in the present work, more especially to meet the available computational resource and enhance timely computation and convergence.

Researchers have performed many numerical studies using CFD simulations of flow and thermal field in the supercritical water channels. As observed by Han and Zhang (2022) most of the CFD studies applied the Reynolds-averaged Navier – Stokes approach. A comparison of the k-epsilon and k-Omega models prediction of the near wall and bulk flow temperatures to the experimentally determined data were made for the continuum flow at 400 kg/m²s and heat flux of 300 kW/m². Fig. 14 shows the extent of prediction of the experimental data by the k-ω and k-ε models and Fig. 15 presents a comparison of the experimental and simulation fuel rod surface temperature prediction.

As observed in Fig. 14 the deviation in the near wall temperature prediction by the k-epsilon turbulence model is much appreciably wider than that of the k-Omega prediction. Studies by Launder and Spalding (1974); Yakhot and Orszag (1986) established that k-ε performs reliably in core flows but tends to over-diffuse turbulence near walls. Menter’s k-ω SST correction addressed these weaknesses, improving accuracy in near-wall thermal layers consistent with the superior near-wall performance seen in Fig. 14. The closeness of the k-Omega prediction at near wall therefore better represents the experimental observation. In the mean flow, it is observed that the k-epsilon models underpredict the experimental bulk temperature data much more than the curve of the k-Omega model. While at the near wall, the k-Omega turbulence model slightly under-predicts the experimental data, in the mean flow, the model almost fits the experimental data. The extent of k-Omega turbulence model prediction as observed in Fig. 14 was considered appreciable and hence, the k-Omega turbulence model was selected and employed in the present work. As observed in Fig. 15, the simulation captures the near-inlet behavior reasonably well but fails to reproduce both the continued downstream temperature rise and the pressure sensitivity observed in the experimental data. This points to a likely gap in the data acquisition methods adopted, particularly near supercritical pseudo-critical conditions. The observed rod surface temperature deviation from the experimental data will be consequentially propagated in the heat transfer coefficient determination by equation 6.

The physical or numerical origin of the peak observed in Fig. 14 especially in the k-ω model prediction could be attributed to certain possible limitations in the CFD analysis and could include non-uniform local mesh quality and or limited availability of high-resolution experimental data for local wall temperature under supercritical conditions. This may be as a result of the practical difficulty of instrumenting high-pressure, high-temperature test loops. The observed peak could also be a classic situation for heat transfer deterioration where strong property variations suppress near wall turbulent mixing and standard eddy viscosity RANS model with fixed turbulent Prandtl number as documented do struggle to predict. The RANS-based k-ω and k-ε closures rely on the Boussinesq isotropic eddy-viscosity hypothesis, which is less accurate under strong buoyancy forces and steep near-wall property gradients near the pseudo-critical point; thus, k-ω’s superior performance should be viewed as relative, not absolute, with LES/DNS offering more rigorous benchmarks beyond the scope of the present study. Sharp thermophysical property variations demand fine property-table resolution, and coarseness may introduce localized inaccuracies. In the prediction of supercritical flows, the RANS model may either overshoot (ie. k-ω) the transition of heat transfer deterioration, smear or even miss it completely without capturing it (ie. k-ε). This is mainly due to strong property variations which are sometimes compounded by buoyancy and flow acceleration effects.

The prediction discrepancies between SST k-ω and Realizable k-ε trace directly to the fact that ω possesses an exact analytical wall boundary condition while ε requires empirical modeling, creating near-wall inaccuracy in k-ε even when the realizability constraint is satisfied. SST resolves the viscous sublayer directly using a mathematically exact wall boundary condition, while Realizable k-ε relies on empirical wall functions that are only valid away from the wall. This makes SST inherently more accurate for separation, reattachment, and skin friction predictions. In strongly non-equilibrium flows, SST k-ω dynamically reduces turbulent eddy viscosity through its Bradshaw-based limiter, physically bounding turbulent shear stress and preventing over-estimation of turbulent mixing, while Realizable k-ε keeps eddy viscosity elevated producing excessive mixing. This systematically over-predicts the turbulent eddy viscosity by relying on mathematical consistency alone rather than the actual stress-carrying capacity of turbulence. However, this advantage reverses in buoyancy-dominated configurations. The Realizable k-ε directly couples buoyancy to its length-scale equation through the C₃ε term, capturing stratification and natural convection more accurately, while SST’s absent ω-buoyancy term causes it to underestimate stratification effects and under-predict natural convection heat transfer by 5–15%. For supercritical fluid applications in the context of the present research, wall-bounded and buoyancy-driven turbulence mechanisms are simultaneously active, indicating that neither model is unconditionally superior. However, given the geometric complexity and critical importance of near-wall accuracy in supercritical flow simulations, SST k-ω with explicit buoyancy correction modifications represents the more physically justified modeling choice. This choice preserves SST’s superior near-wall resolution and stress-bounding capability while addressing its deficiency in buoyancy-scale coupling

Pioro and Duffey (2005) reported that the working fluid employed in supercritical water-cooled reactors, without going through the distinct phase of boiling, transitions from liquid water to saturated steam and becomes supercritical above the critical temperature of 374 °C as presented in Fig. 1. In the present work, STARCCM+ simulations were carried out at flow pressures of 23 MPa and 25 MPa, keeping continuum inlet temperature constant at 310 °C. Imposition of heat flux on the fuel rod surface and varying coolant flow led to corresponding exit temperatures of the continuum. The effect of pressure on the continuum exit temperatures and fuel rod surface average temperatures are presented in Figs 16, 17.

Supercritical fluid at pressures of 23 MPa and 25 MPa are associated with temperatures of 377.5 °C and 384.9 °C respectively as shown in Table 1. For the imposed conditions of heat and mass flux of 300 kW/m² and 400 kg/m²s respectively in the present numerical study, the fluid exit temperatures realized (ie. 353.98 °C at 23 MPa and 352.16 °C at 25 MPa). Bulk flow exit and rod surface average temperatures realized at 23 MPa as observed were higher than those obtained at 25 MPa, indicating that, under the conditions specified (ie. 300 kW/m² and 400 kg/m²s), a much higher and desirable cooling of the fuel rod surface is attained at 25 MPa. This is because the coolant at 23 MPa and specified conditions as observed in Fig. 16 is of higher heat content than at 25 MPa and hence unable to remove much more heat relative to the extent of heat removal observed in at 25 MPa in Fig. 17.

The effect of pressure variation on the fuel rod surface average temperature and bulk flow temperature profiles along the flow direction are as presented in Figs 18, 19.

Fuel rod surface temperature and bulk flow temperature along the mean flow direction are higher at 23 MPa than at 25 MPa. A rise and fall in fuel rod surface average temperature were noticed, some few distance away from the entry to the domain. Whiles a peak temperature of 58 °C at 23 MPa was attained at a position of 0.21 m away from the entry to the domain, at 25 MPa, a relatively lower fuel rod peak temperature of 540 °C was observed 0.26 m away from the entry to the domain. The fluctuation in fuel rod surface average temperature was observed at 23 MPa between point 0.11 m at 400 °C and 0.26 m at 518 °C whereas at 25 MPa, the fluctuation occurred between points 0.11 m and 0.37 m at 381 °C and 407 °C respectively. A higher magnitude of temperature fluctuation was observed at 23 MPa whiles the extent of temperature decay from the peak was much more pronounced at 25 MPa.

Higher turbulence kinetic energy was observed for the operation at 23 MPa, 300 kW/m² and 400 kg/m²s than at 25 MPa as shown in the plot of turbulent kinetic energy (TKE) along the flow direction in Fig. 20.

An abrupt rise in turbulence was observed between 0.16 m and 0.37 m away from the entry to the flow domain followed by a steady rise in turbulence. Turbulent kinetic energy at mass fluxes of 800 kg/m²s and 1200 kg/m²s for operating pressures of 23 MPa and 25 MPa showed a marked rise in intensity at virtually equal magnitude and trend as presented in Fig. 21. Variations in trends and magnitude were observed at mass fluxes below 800 kg/m²s.

It is worth to note that supercritical condition was not established under the conditions plotted in Fig. 21 as justified in the bulk exit temperature values reported in Tables 8, 9 which are all lower than the required corresponding flow exit temperature of 384.9 °C for supercritical fluid at a pressure of 25 MPa and 377.5 °C at a pressure of 23 MPa reported in Table 1.

Turbulent kinetic energy for cases when supercritical conditions was established at varying heat flux and constant mass flux of 228 kg/m²s is presented in Fig. 22 with the setting in of supercritical condition justified by the bulk exit temperature values presented in Tables 10, 11 under the operating pressures of 25 MPa and 23 MPa respectively.

As observed in Fig. 22, higher turbulent kinetic energy is realized at higher pressure and Table 11 shows that at a lower pressure, supercritical condition only set in at a much higher heat flux than observed at 25 MPa. Fig. 23 presents contour plots of turbulent kinetic energy showing the presence of higher turbulence under supercritical conditions than compressed conditions of the fluid considered. The increase at 25 MPa was observed at a ratio of 312 % for a bulk exit temperature range of 377 °C to 427 °C and 354% at 23 MPa for a bulk exit temperature range of 354 °C to 391 °C.

Peak TKE values increase from 0.0021 J/kg (a) and 0.0023 J/kg (c) in compressed fluid to 0.0094 J/kg (b) and 0.0103 J/kg (d) in supercritical fluid representing 4.48 times increase in TKE when transitioning from compressed to supercritical fluid. From the figure, minimum TKE values remain comparable across all four cases (~5–9×10^–9 J/kg), meaning the increase in TKE is entirely concentrated at the high-energy end, expanding the dynamic range from ~10⁵ in the compressed fluid to ~10⁶ in the supercritical fluid. At 23 MPa, peak TKE is approximately 9.5% higher than that at 25 MPa within both thermodynamic states, confirming pressure as a weak secondary influence. The higher TKE at lower pressure (23 MPa) directly reflects proximity to the pseudo-critical temperature. At this pressure, the Pseudocritical temperature is lower, meaning the fluid’s thermophysical property gradients are steeper at the simulated temperatures, generating stronger buoyancy production and baroclinic vorticity. High near-wall TKE in supercritical fluids as observed does not translate into high near-wall turbulent mixing due to the suppression of heat transport by turbulence resulting from thermal laminarization in the near wall layer, thickening of the viscous sub-layer which elevate molecular conduction a the wall creating a barrier that outer layer turbulence cannot overcome and elevated turbulent Prandtl number which decouples turbulent momentum transport from thermal transport. The observed increase in turbulence (TKE) signals altered flow structure and thermophysical coupling that can produce heat transfer deterioration, sharp near-wall enthalpy absorption analogous to boiling, and spatial redistribution of thermal energy. The near wall thermophysical state, not bulk turbulence energy, is the true governing parameter of convective heat transfer performance in the supercritical regime.

Sharabi et al. (2008) in their prediction of three dimensional turbulent heat transfer to CO₂ at supercritical pressure of 8 MPa, mass flux of 314 kg/m²s and heat flux of 30 kW/m², observed a reduced production of turbulence along the symmetry line BC shown in Fig. 24 and along axial nodes of a triangular channel, as a consequence of the effect of buoyancy in redistributing the flow by increasing the velocity close to the wall and decreasing it in the core. This led to low turbulent kinetic energy at low flow velocities with the maximum reduction in turbulence occurring where the wall temperature reaches their maximum value.

At the prevailing supercritical condition of CO₂, in the investigation conducted by M. Sharabi et al., the magnitude of turbulent kinetic energy attained was in the order of 10^-3 as presented in Fig. 25 which agrees with the present findings. This is indicative of the occurrence of low turbulence or laminarization in the flow of supercritical fluids through bare rod bundles.

The intensity of turbulence observed in the present work at exit increased with increase in mass flux with relatively negligible effect of pressure at 300 kW/m² and much higher effect when the heat flux is doubled. These observations are as presented in Figs 26, 27.

Heat transfer from fuel rod to the coolant was analyzed using the numerical data obtained in the present work. The heat transfer coefficient (HTC) from rod to fluid were determined using equation 6 (Eq. 6). Fig. 28 shows the profiles of heat transfer coefficient for two distinct conditions of the coolant, specifically, existing as compressed fluid (defined in Fig. 1) and occurring at 300kW/m² and also as supercritical fluid (defined in Fig. 1) occurring at 600 kW/m² for two different pressures of 23 MPa and 25 MPa respectively.

It could clearly be observed per the results of the present analysis as presented in Fig. 28 that, the coolant existing as compressed fluid at pressures of 23 MPa and 25 MPa with corresponding temperatures of 354 °C and 377 °C respectively is more efficient is fuel rod heat removal than that existing as supercritical fluid at pressures of 23 MPa and 25 MPa with corresponding temperatures of 391 °C and 427 °C respectively. Interestingly for each condition of the fluid, the heat transfer at the higher pressure is higher (25 MPa) than that at the lower pressure (23 MPa). This result was obtained at a constant mass flux of 228 kg/m²s. It suffices to state that, the relatively high levels of turbulence recorded in supercritical fluid as shown in Fig. 22 could not translate into high rod to fluid heat removal by the supercritical fluid as presented in Fig. 28.

As presented in Fig. 29, when mass flux was increased to 800 kg/m²s from 228 kg/m²s, at a fixed heat flux of 600 kW/m², the rod to fluid heat transfer coefficient for the compressed fluid (800 kg/m²s, 600 kW/m² at 23 MPa and 25 MPa) still exceeded that of the supercritical fluid (228 kg/m²s, 600 kW/m² at 23 MPa and 25 MPa) with higher heat transfer occurring at 25 MPa for the compressed fluid and almost the same levels of heat transfer at both 23 MPa and 25 MPa for the supercritical fluid. The turbulent kinetic energy profile for the cases presented in Fig. 29 is as shown in Fig. 30.

Clearly, relatively higher turbulent kinetic energy is present in the flow of the supercritical fluid (228 kg/m²s, 600 kW/m² at 23 MPa and 25 MPa) than that of the compressed fluid (800 kg/m²s, 600 kW/m² at 23 MPa and 25 MPa) as presented in Fig. 30. However, the higher extent of turbulence observed in the supercritical fluid does not lead to more effective heat transfer as higher heat transfer coefficients were rather associated with the compressed fluid which exhibited relatively lower extent of turbulence.

Turbulent Kinetic Energy (TKE) and the corresponding Heat Transfer Coefficient (HTC) at exit of the domain for both pressures of 23 MPa and 25 MPa for the cases presented in Figs 22, 30 are be summarized in Tables 12, 13.

In Table 12, the highest turbulent kinetic energy (TKE) of 5.4×10^-3 J/kg was obtained at 23 MPa in a supercritical fluid with a corresponding lowest heat transfer coefficient (HTC) of 1.0, whiles the lowest TKE of 2.3×10^-3 J/kg occurred in a compressed fluid at 25 MPa with a corresponding HTC of 3.9 which is relatively higher than that obtained in the supercritical fluid. It could also be noted that, whiles the highest TKE (5.4×10^-3 J/kg) occurred under supercritical fluid condition at 23 MPa, the highest HTC (11.4) occurred under compressed fluid condition at 25 MPa. These observations were made under constant heat flux of 600 kW/m² and varying heat flux and pressure. Table 13 present the scenario for constant mass flux of 228 kg/m²s, varying heat flux and pressure.

In Table 13, the highest TKE of 2.9×10^-3 J/kg occurred in supercritical fluid at 23 MPa with a corresponding low HTC of 1.1 whiles the lowest TKE of 1.3×10^-3 J/kg occurred under compressed fluid condition at 25 MPa with a corresponding HTC of 1.6, relatively higher than that observed under the supercritical condition. Whiles the highest TKE (4.7×10^-3 J/kg) occurred under supercritical fluid condition at 25 MPa, the highest HTC (1.61) occurred under compressed fluid condition at 23 MPa. In both cases considered at constant mass flux of 228 kg/m²s and constant heat flux of 600 kW/m², the highest TKE predominantly occurred under supercritical fluid condition whiles the highest HTC are prevalent in compressed fluid condition indicating that the high extent of turbulence realized in supercritical fluids does not translate into effective heat removal from the rod to the coolant under supercritical condition. Largely, heat transfer deterioration from compressed fluid to supercritical fluid is observed in the present study. Several researchers, Yamagata et al. (1972); Koshizuka et al. (1995); Cheng et al. (2007); Zhang et al. (2010); Wen and Gu (2010) have conducted experimental and numerical studies on fluid flow and heat transfer of supercritical water flow in circular tube. Results of these investigations showed that the contributions of buoyancy and the thermal acceleration due to the sharp variations of thermophysical properties of the supercritical water near the pseudo-critical point might be the main reasons of the abnormal heat transfer phenomenon.

Whiles the analysis based on bulk-averaged or wall-averaged quantities as presented is entirely relevant, it is required to also discuss the local thermal-hydraulic behavior, which is particularly critical in supercritical water flows, where near-wall behaviors appreciably influence heat transfer. In rod bundle configurations, strong spatial non-uniformities are expected at near rod bundle walls. To significantly enhance the physical insight of the presentation, local near wall temperature and turbulent kinetic energy distributions obtained over probes positioned, 4 per each rod are presented in Figs 32, 33, 34, 35. Fig. 31 shows the location of the probes, circumferentially on the rod wall surface.

To understand how well heat is being removed from each rod, temperature and turbulence intensity were measured at four points around the surface of each rod. The fluid enters at a lower temperature (~630 K) and heats up rapidly near the inlet, with wall temperatures peaking between 1149 K and 1170 K. Rod 3 recorded the highest value of 1170.6 K across the bundle. These sharp temperature rises near the inlet are characteristic of supercritical flows, driven by the steep decline in heat capacity of the fluid near the pseudo-critical temperature. This behavior is widely reported in recent studies ie. Zhang et al. (2021); Gu et al. (2022).

The key finding is that not all surfaces of a rod cool equally. Probe positions facing the open gaps between rods benefit from stronger fluid mixing, which carries heat away more effectively; these locations show turbulence levels of 1.5×10^–4 to 2.4×10^–4 m²/s² and moderate outlet temperatures of 1035–1082 K. In contrast, probe positions facing the confining wall where fluid movement is restricted show very low turbulence (1.3×10^–6 to 6.0×10^–6 m²/s²) and persistently high surface temperatures of 1134–1154 K at the outlet. The direct relationship between reduced turbulence and diminished heat removal in wall-confined subchannels of rod bundles has been consistently demonstrated in recent numerical and experimental investigations such as in Chen et al. 2023) and Wang et al. (2024).

The temperature difference around the circumference of each rod ranging from 18.7 K (Rod 1, most uniform) to 27.1 K (Rod 3, least uniform) confirms that heat transfer is uneven around each rod. Rod 1 performs most consistently, while Rods 3 and 4 show the greatest imbalance between their gap-facing and wall-facing sides. According to Shen et al. (2023), in supercritical water reactor designs, such circumferential non-uniformity is recognized as a key indicator of localized thermal stress and potential heat transfer deterioration and the wall-confined surfaces identified warrant careful attention in safety and thermal design assessments.