Large Eddy Simulation Study on Forced Convection Heat Transfer to Water at Supercritical Pressure in a Trapezoid Annulus

The supercritical water-cooled reactor (SCWR) is one of the six prospective generation IV nuclear reactors that are currently under development. One of the key aspects of the SCWRs is the unique features of the thermal hydraulics of such designs. The thermal–physical properties of the fluid at supercritical pressure vary dramatically at the vicinity of pseudocritical temperature. This often leads to heat transfer deterioration or enhancement which has been reported by many researchers. Predicting the detailed flow and the temperature distributions for the design development is a challenge in this area. The fuel assemblies of the SCWR designs are in a form of rod bundles similar to many other nuclear reactor designs. The coolant flows through the rod bundles via the subchannels formed by the fuel rods array. Various gaps are formed between the fuel rods, which connect these subchannels to each other and are characterized by the pitch-to-diameter ratios. Instability and interchannel mixing have been found to exist in such narrow gap regions. Therefore, the knowledge of the flow and heat transfer at supercritical pressure in a rod bundle is an essential input to the SCWR core design.

Most experimental studies of supercritical water-cooled reactors are in pipe or channel configurations. There is limited literature on the experimental study at supercritical pressure in rod bundles. One of the early studies of the flow of fluid at supercritical pressure in a rod bundle was performed by Dyadyakin and Popov. The bundle consisted of seven-rods with helical fins on them and a local heat transfer correlation was produced based on the experimental data [1]. Significant pressure oscillations were observed at high heat fluxes and large mass fluxes.

There are some recent experimental studies on 2 × 2 rod bundles with supercritical water as working fluid [2–4]. Several correlations were assessed against the test data from the experimental studies, and it has been concluded that the best correlations compared are those of the Jackson and Fewster [5] and Bishop et al. [6]. Reference [2] also showed that the buoyancy parameter developed by Cheng et al. [7] characterizes well the dependence of heat transfer coefficient (HTC).

Computational fluid dynamics (CFD) simulations are also a tool that has been widely used to study the flow behaviors in rod bundles. The Reynolds-averaged Navier–Stokes/unsteady Reynolds-averaged Navier–Stokes approaches are found to be most popular to study the flow and heat transfer of fluids at supercritical pressure in rod bundles [8–11]. Hooper and Rehme [12] first demonstrated the existence of some periodic azimuthal variations of the turbulent-velocity component through the gaps of the closely spaced rod arrays. The azimuthal turbulent-velocity component is not associated with the secondary-flow velocities driven by the Reynolds-stress gradients. In addition, flow instability has also been observed in other nonuniform channel flows [13–19]. An experimental validation study of 2 × 2 rod bundle conducted by Xiong et al. [20] used the ω-Reynolds stress model and baseline-Reynolds stress model. The models were found to be capable of predicting the cross flow pattern in between subchannels. They found the discrepancies in turbulent mixing prediction due to the strong cross flow between subchannels. However, these did not lead to a large discrepancy in bulk parameters of the subchannels. Although numerical studies using large eddy simulation (LES) are still limited due to the excessive requirements of computational time and mesh density, it has been demonstrated by some authors to be superior to the Reynolds-averaged Navier–Stokes in predicting the characteristics of the flow such as the vortex street and flow pulsations in the narrow gaps in fuel rod bundles [18,21–23]. In particular, it has been demonstrated in Refs. [22] and [23] that the buoyancy force has non-negligible effects on the behavior of the large flow structures in the narrow gaps of the tight-packed rod bundle.

This paper focuses on the forced convection of supercritical water flow in a trapezoid annulus, which is a simplified model for the tightly packed triangular-array rod bundle in the nuclear reactor. The LES with wall adapting local eddy viscosity (WALE) subgrid-scale model embedded in an open-source CFD code, “code_saturne,” is used to simulate the heated supercritical water flow. It is well known that the properties of the supercritical water change rapidly around the pseudocritical temperature, including, for example, the density and heat capacity as illustrated in Fig. 1. This behavior may affect the effectiveness of the heat transfer through affecting the turbulence production, which causes heat transfer enhancement or deterioration. The thermal–physical properties of water at supercritical pressure of 25 MPa are used, which are integrated in the code_saturne software for the simulations using a subroutine. It should be mentioned here that the pseudocritical temperature (T_pc) of the considered fluid is ∼358 °C. The purpose of this study is to develop a better understanding of the heat transfer phenomenon and flow structure in the tightly packed rod bundles.

The considered configuration of the flow passage is a trapezoid annulus which is similar to that studied experimentally by Wu and Trupp [13] and also that adopted by Duan and He study [22,23]. The cross-sectional view of the channel is illustrated in Fig. 2, which has a narrow gap and a wide gap between the outer wall and the rod located at the bottom and top of the channel, respectively. The “main channels” are located on either sides of the rod connected by these two gaps. The enclosed rod has a diameter D of 0.0508 m and is located S = 0.004 m above bottom wall. Consequently, the ratio of S/D for this gap is 0.079. The height of the trapezoid channel is 0.066 m, and the widths of the top and bottom wall are 0.0508 m and 0.127 m, respectively. The hydraulic diameter, D_h, of the channel is 0.0314 m. Referring to Fig. 3, the total length of the domain is divided into two parts, namely, the inflow generator and the test section. The inflow generator is an isothermal domain with a hydraulic diameter length of 10D_h at the upstream. At the 5D_h length of the isothermal domain, the velocity field at every time-step is copied and fed back to the inlet as a boundary condition for the next time-step. This aims to produce a fully developed inlet flow. Meanwhile, at the test section, a constant heat flux is applied on the inner rod wall, and the length of this heated section is 30D_h. Therefore, the total length of the simulated domain is 40D_h. The values of various parameters used for the simulations are summarized in Table 1.

In this study, a constant heat flux is imposed on the rod wall in the heated section. Three different heating cases have been considered to study the effects of the variations of thermal properties on the behavior of the large flow structures. The heat fluxes considered are 10 kW/m², 50 kW/m², and 75 kW/m². The inlet temperature is chosen to be 360 °C which is below the pseudocritical temperature of the working pressure of 25 MPa. The thermal properties were obtained from the National Institute of Standards and Technology database. The Reynolds number is 10,540 based on the bulk velocity of 0.0388 m/s.

A fully structured mesh containing 33 × 10⁶ elements with finer mesh grids near the walls was created. An overview of the mesh cross section can be seen in Fig. 4. The first mesh nodes next to the wall are in the range of 5 ≤ Δx⁺ ≤ 17, 0.13 ≤ y⁺ ≤ 0.2, and Δz⁺ ≈ 23. There are about 13 cells located within y⁺ ≤ 20 with ∼9 cells located in the range 0 ≤ y⁺ ≤ 10. The time-step is 0.001 s for all cases, with the Courant–Friedrich–Lewy number below 0.5. The flow domain was solved by using LES with the WALE subgrid-scale model embedded in code_saturne. The velocity and pressure coupling is obtained using semi-implicit method for pressure linked equations-consistent algorithm. A second-order central difference scheme in space and time is adopted.

To aid the discussion, it is necessary to introduce the locations used to present the results. The points at the center of the narrow gap and wide gap at the symmetric line of the channel are referred to as NG and BG, respectively, see Fig. 2. The histories of the instantaneous velocities at these points are recorded for several axial locations for the spectrum and correlation analyses. An equidistant plane (inset of Fig. 2) is also used to present the instantaneous flow field at the narrow gap.

There are no experimental data or direct numerical simulation data which can be used to directly validate the current LES. Similar to the work [22], we quantify of LES simulations using the contours s and $LES_IQv$ that have been suggested by Geurts and Fröhlich [24] and Celik et al. [25], respectively. Results are shown in Figs. 5 and 6 for several distances down the channel for half cross sections of the trapezoid annulus thanks to the symmetry of the geometry. To represent the quality of the results, the case of the highest heat flux of 75 kW/m² is chosen here. $μsgs$ is a function of the element's size and appears in both criteria. Evidently, it can be seen that the value of the parameter s in the channel largely reflects the mesh distribution in both figures. From Fig. 5, the values of s are very low in the near wall regions, while in this region the sizes of mesh elements are really small. The value of s is maintained at around $∼$0.15 in this region. The overall values of s increase slightly down the channel. The region of the values of s > 0.45 spreads more over the cross section down the channel. However, in the region near the heated rod, the values of s can be seen to decrease down the channel, which can be clearly seen across the half top of the rod pin. Inversely, the values of $LES_IQv$ increased in the region near the heating rod. Overall, this parameter in the channel is maintained at a value around 0.97–0.98. The value of s can be seen as an estimation of the fraction of the turbulent kinetic energy modeled by LES. According to Pope [26] states that to obtain a good LES simulation, the turbulence kinetic energy has to be resolved over 80%. For the current simulations, the values of s are reasonably small in all wall regions, but are rather high in some central regions. On the other hand, Celik et al. [25] commented that an LES is considered to be good when $LES_IQv$> 0.8. In this case, the value of $LES_IQv$ is found to be around 1 near the wall and >0.97 across the channel.

Figure 7 shows the general behaviors of the mean values of the flow and temperature in the channel. As expected, the bulk fluid temperature and wall temperature increase monotonically down the channel in each of the cases. The temperature increases more rapidly in the case with a higher heat flux. From these results, it can be noted that for the first heating case of 10 kW/m², the average wall temperature, the bulk temperature of the flow, and the pseudocritical temperature are related in the T_b < T_w < T_pc manner. Meanwhile, for the other two higher heating cases, the average wall temperature exceeds the pseudocritical temperature so that T_b < T_pc < T_w. Consequently, that the fluid temperature in the latter two cases goes across the T_pc at some point, which results in a rapid decrease in the density and strong variations in other thermal properties. Since the density is negatively correlated to the fluids temperature, the bulk velocity is also expected to accelerate down the channel, which is clearly shown in Fig. 7(c). It also shows that the acceleration rate is positively related to the heat flux. In comparison, the acceleration is relatively linear in the case of 10 kW/m², especially in the initial region of the heated section. The slower density change in this case is the reason. The averaged HTCs however show an opposite trend, that is, it decreases in value with increase of heat flux. The value of HTC also decreases along the heated channel, as can be observed in Fig. 7(d). The reduction of HTC is the strongest in the case of the lowest heat flux. The higher the heat flux, the lower the HTC is. Hereafter, the case with a heat flux 10 kW/m² is referred to as the “low heating case,” and the cases of 50 kW/m² and 75 kW/m² are referred to as the “high heating cases” when the phenomena of the large flow structures are discussed.

In the study of Duan and He [22], the general features of the large flow structure were visualized in the narrow gap of the channel. The existence of swinging large flow structures in the cases studied here is demonstrated in Fig. 8, using the contours of the instantaneous velocity and temperature field at the equidistant plane between the rod wall and the bottom wall of the trapezoid channel (see the inset of Fig. 2). The general structures for the swinging flow are very similar for all heating cases. Similar to the bulk velocity, the streamwise velocity at the center of the narrow gap increases with the increase of heat flux. Similarly, the temperature reaches a maximum at the center of the narrow gap, and this value is higher for a case with a higher heat flux. The large flow structures, bringing the high-temperature fluids in and out of the narrow gap region, are expected to increase the mixing and help the heat transfer.

Figure 9 shows the instantaneous temperature on the rod surface facing the narrow gap and the wide gap. A difference in the distribution of the temperature can be seen on the rod wall facing the two gaps. Clearly, the swinging flow structures in the narrow gap region cause an oscillation in the temperature distribution on rod surface facing the narrow gap in all cases. Such oscillations of the rod temperature may increase the thermal fatigue, a risk that needs to be considered in design even though generally speaking the large flow structures are good for heat transfer reducing the mean temperature in small gapping regions. Meanwhile, there is no obvious temperature oscillations observed on the rod wall facing the wide gap.

The flow structures in the narrow gap can be further inspected using the time history of the spanwise velocity. Figure 10 shows the time histories of normalized fluctuating spanwise velocity at NG down the heated channel for each of the heating cases. Strong quasi-periodic oscillation of the spanwise flow at the narrow gap can be seen in all cases. The dominant periods and the trends of the oscillations remain to be similar down the heated channel. To quantify the dominant frequency of the quasi-periodic flow structures, the power spectrum density (PSD) of the u′ is evaluated down the heated channel and is illustrated in Fig. 11. To assist the presentation of data, the PSD values for the 50 kW/m² and 75 kW/m² are multiplied by a factor 10² and 10⁴, respectively. There are visible peaks shown in the distribution of the PSD of u′ at the narrow gap for each of the heating cases at different heights. The dominant frequencies of the structures where the peaks are located are referred to as the peak frequencies (f_p). Throughout the channel, f_p remains the same for each heating case with a similar value of ∼0.21 Hz. The dimensionless frequency Strouhal number which is evaluated using the bulk velocity, U_b, at the unheated section is around St⁻¹ ≅ 5.62 for all three cases. This value is close to the value obtained by Wu and Trupp [13] and Duan and He [22], although the current Reynolds number is only 20% of the former and doubled that in the latter. This result shows that the Strouhal number is only dependent on the geometric configuration when the buoyancy force is not considered. Furthermore, stronger dissipations of energy can be observed at the frequencies above 4 Hz for all cases. The noise level (i.e., the oscillations) at the high-frequency end of the spectra varies significantly from case to case. For example, the oscillations of PSD at higher frequency above 10 Hz at z/D_h = 25 and 35 in the case of 50 kW/m² are significantly lower than in other cases. It appears that such noise level is dependent on the samples used. Some samples seem to bear more noises than others. This might be linked to the history of the fluid path of that sample. The noise/oscillations can be smoothed out with ensemble average using a large number of samples. Only three samples are used for each spectrum presented here, since each such spectrum requires a time history of 20 s and this is resources demanding.

It is useful to inspect the flow structures at the wide gap, noting that Refs. [13] and [22] have shown the existence of the large flow structures not just in the narrow gap but also in the wide gap. Figure 12 shows the representation of the time history of the normalized fluctuating spanwise velocity down the heated channel for each case. For all heating cases, the signal shows strong turbulent fluctuations with no obvious regular oscillation observed. It is probable that the swing flow structures are weak, and the turbulent noises overwhelm their existence in the time wave form. The PSDs of the u′ at BG at several locations are shown in Fig. 13 for each heating case. There are some dominant peaks seen for the PSD at some axial location for each heating case. These peaks are clearly seen at the length of z/D_h = 15 and reduced down the channel and hardly identifiable at z/D_h = 35. There are two sets of peaks at the beginning of the heating section of z/D_h = 15. The first is in the range of 0.2–0.25 Hz, and the second set is at around 0.5 Hz. Note that the first set of peaks is very similar to f_p of the flow structures in the narrow gap, but the magnitude of the peaks is lower than those in NG. These results show the existence of weaker large flow structures in the wide gap, and these flow structures are potentially linked to the swing flow structures in the narrow gap. The second set of peaks may be related to the existence of a second large flow structure which disappears quickly down the channel.

To investigate the relationship between the flow structures in the narrow and wide gaps, the cross-correlation functions of u′ between flows at NG and BG are calculated, which are illustrated in Fig. 14. The correlations show that the flows at the two gaps are inversely correlated at 0 s lag for the heat flux of 10 kW/m². This correlation is very strong at upstream but reduces down the heated channel. For both of the higher heating cases, the flow structures between the gaps are weakly correlated near the heating inlet; then is the negative correlation increases somewhat before reducing again. This may be due to the changes in local temperature in the gaps where strong changes of properties occur, consequently leading to changes in the flow structures in these gaps. These inverse correlations of the flow structures at 0 s lag between the gaps are similar to the forced convection case in the study [22]. It is interesting to mention that the study carried out in Ref. [22] was based on constant properties. However, at high heat fluxes where the wall temperature could exceed the pseudocritical temperature, the correlation of the flow structure between gaps appears to be reduced.

Thanks to the swing behavior of the large flow structures, the fluctuating velocity u′ is a good quantity to describe the flow behavior as shown earlier. The scale of the flow structures can be statistically measured using two-point correlations [26]. The axial scale of the large flow structures is studied using cross-correlation between u′ at the point near the entrance of the heating section (z/D_h = 11.25) and several positions down the heated channel in the narrow gap and the wide gap. As shown in Fig. 15(a), the large flow structure passing through the narrow gap is similar to each other for all cases. Some elongation of the structure is seen for the higher heating cases where the size starts to increase at z > 0.6 m. This may be the result of the accelerated flow in that region due to thermal expansion. Unlike the flow structures in the narrow gap, the flow structures in the wide gap are weakly correlated. Nevertheless, a pseudo-oscillation behavior of the peaks could be seen for z < 0.75 m for the low heating case before the structures are elongated. The peaks are weak in the high heating cases, but there are some additional peaks around 0.45 m and 0.6 m. These additional peaks are in between the main peaks, which are similar to the peaks of the narrow gap. These two sets of peaks may be related to the two sets of peaks exhibit in the PSD in Fig. 13.

Figures 16 and 17 present correlation coefficients of Ruu as suggested by Guellouz and Tavoularis [15] between two points at several axial locations at the symmetry lines of the wide and narrow gaps of the test section. To perform the calculation of the correlation coefficient of the flow in gaps, instantaneous velocities are recorded at a point downstream from the heated inlet (z = 25D_h) as well as additional points further downstream (Δz/D_h = 0.2, 0.4, 0.6). Reference [15] estimated the average convection speed, U_c, of the structures using these space–time correlations. The convection speed, U_c, is the ratio of the streamwise probe separation Δz over the time delay for maximum correlation Δt_max. Table 2 shows the calculated convection speed for each case in both gaps. The convection speed changes when the heat flux is changed. At the narrow gap, the convection speed is inversely proportional to the increase of the heat flux, while the opposite is true at the wide gap. At the wide gap, the convection speed increases gently as the heat flux is increased.

In Fig. 18, u_eff normalized with U_b down the heated channel are presented for both the narrow gap and wide gap. The trends of u_eff/U_b are different for the different cases at the two gaps. The variation of u_eff/U_b for the low heating case of 10 kW/m² is largely unchanged along the heated domain, while for the higher heating cases u_eff/U_b reduces along the heated channel. This result shows that the mixing between the subchannels reduces with distance from the start of heating when the temperature follows T_b < T_pc < T_w. It is also observed that the variations of u_eff/U_b are not monotonic in most cases, and there is a minimum point before the end of the inspected length. The location of the minimum moves downstream with the increase of heat flux. The turning point of the trend can be observed at z/D_h = 20, 25, and 30 for heat flux 10 kW/m², 50 kW/m², and 75 kW/m², respectively. This can be explained as the development of the mixing taking longer distance for higher heat flux. Longer domain may be required to further study the mixing down the heated channel after the recovery from the reduction, especially at higher heating.

In the wide gap, an opposite behavior can be seen, that is, u_eff/U_b increases down the heated channel for the higher heating cases, while u_eff/U_b gently decreases in the low heat flux case. The value of u_eff/U_b is much lower than in the narrow gap, since the unsteady flow structures at the wide gap are weaker in comparison. It is interesting to compare u_eff/U_b obtained here with that of Ref. [22], which is 0.135. This value is quite close to the value obtained at z/D_h = 35. In Table 3, u_eff obtained here using LES results at $x/δij$ = 0.18 is compared with the u_eff calculated using Eq. (5). The u_eff from the correlation by Rehme does not correlate well with the u_eff obtained using LES. The ratio between the results of the correlation and those of LES is around 0.46, which is far lower than the results obtained in Ref. [22]. This suggests that the variable properties considered in this study have a strong influence on the mixing velocity.

The mixing coefficients, Y, calculated at NG are shown in Fig. 19. These results show a similar trend to that of u_eff down the heated channel. The results obtained here can be compared with the result of forced convection case in Ref. [22]. The values at z/D_h = 35 are about 15–20% higher than the value of 113.50 in Ref. [22] for the heat fluxes 10 kW/m² and 75 kW/m², while the value is almost the same with that of Ref. [22] for the heat flux of 50 kW/m². Moreover, the differences between the two high heating cases, 50 kW/m² and 75 kW/m², for both u_eff/U_b and Y and that of Ref. [22] could be due to the differences in the local temperature at the gap region, where the temperature in this region is exceeding the pseudocritical temperature for some distances in the heated domain. However, as mentioned earlier, it is difficult to determine trend of the mixing after the recovery due to the limited heated domain length in this study especially for the 75 kW/m². This mixing has shown a complex behavior in the high heating cases where the temperature is following T_b < T_pc < T_w order.

Large eddy simulations of turbulent forced convection of water at supercritical pressure in a trapezoid annulus have been performed using code_saturne. The effects of the thermal expansion and the variation of other properties on the large flow structures in the channel have been discussed. The results obtained clearly demonstrate the existence of the large flow structures in the narrow gap and the wide gap. This agrees with the finding by Wu and Trupp [13] and Duan and He [22]. There are differences between the behaviors of the large flow structures in the higher heating cases and the low heating case. The swinging flow pulsation observed in the narrow gap led to the temperature variations on the rod wall. This may increase the thermal fatigue of the rod, especially on the side of the narrow gap. The flow structures in the wide gap are inversely linked to those in the narrow gap, but they are much weaker. The St⁻¹ of the channel is found to be similar to that of Ref. [22] and close to Ref. [14]. The convection speed in the narrow gap reduces with the increase of heat flux, but the trend is reversed in the wide gap.

The effective mixing velocity, u_eff, is also investigated at the gaps and a point away from the gap. u_eff obtained using LES away from the gap is compared with the Rehme correlation in Ref. [27]. The Rehme correlation is unable to predict the u_eff of the results obtained using LES, which implies that the effect of variable property is strong. Indeed, the effect of the variations of the properties on the behavior of large flow structure in the gaps has been observed as discussed earlier. The mixing factor, Y, and u_eff/U_b are also compared to the forced convection result in Ref. [22]. The variation of the properties due to the heating influences the differences in mixing velocity at the gaps between this study with Ref. [22].

An earlier version of this paper was presented at the Eighth International Symposium on Supercritical Water-Cooled (ISSCWR-8) held on Mar. 13–15, 2017 in Chengdu, China.

• Engineering and Physical Sciences Research Council through UK Turbulence Consortium, which provides access to the facilities of the UK national supercomputer ARCHER (EP/L000261/1).

• Majlis Amanah Rakyat (MARA).

D =

diameter of the rod, m

D_h =

hydraulic diameter, m

E_uu(f) =

power spectrum density function

f* =

friction factor

f_p =

frequency for the peak, Hz

$LES_IQv$ =

LES quality parameter [25]

q =

wall heat flux, kW/m²

Ruu =

two-point correlation coefficient

s =

LES quality of $s≈μsgs/(μsgs+μ)$

S =

the size of the narrow gap, m

t_max =

time of maximum correlation

T_b =

bulk temperature, °C

T_pc =

pseudocritical temperature, °C

T_w =

wall temperature, °C

$u′,v′,w′$ =

fluctuating velocity component, m/s

u_eff =

effective mixing velocity, m/s

U_b =

bulk velocity, m/s

U_c =

mean convection velocity of coherent structure, m/s

x, y, z =

spanwise direction, wall normal direction, and streamwise direction

Y =

mixing coefficient

Greek Symbols

Δ =

difference or delay

δ_ij =

distance between two subchannels

$ν$ =

kinematic viscosity

Subscripts or Superscripts

b =

bulk

c =

convective

eff =

effective

h =

hydraulic

ij =

two adjacent subchannels

max =

maximum correlation

p =

peak

pc =

pseudocritical

sgs =

subgrid scale

w =

wall

0 =

initial or isothermal

Nondimensional Numbers

Re =

Reynolds number, $UbDh/ν$

St =

Strouhal number, $fpDh/Ub$

Acronyms or Abbreviations

HTC =

heat transfer coefficient

LES =

large eddy simulation

PSD =

power spectrum density

SCWR =

supercritical water-cooled reactor

WALE =

wall adapting local eddy viscosity