Abstract
Because of their high molecular heat conductivity, low-Prandtl number liquid metal is a promising candidate coolant for various designs of advanced nuclear systems such as liquid metal–cooled fast reactors and accelerator-driven sub-critical system (ADS). With the fast-growing computational capacity, more and more attention has been paid to applying computational fluid dynamics (CFD) methods in thermal design and safety assessment of such systems for a detailed analysis of three-dimensional thermal–hydraulic behaviors. However, numerical modeling of turbulent heat transfer for low-Prandtl number liquid metal remains a challenging task. Numerical approaches such as wall-resolved large eddy simulation (LES) or direct numerical simulation (DNS), which can provide detailed insight into the physics of the liquid metal flow and the associated heat transfer, were widely applied to investigate the turbulent heat transfer phenomenon. However, these approaches suffer from the enormous computational consumption and are hence limited only to simple geometrical configurations with low to moderate Reynolds numbers. The Reynolds-averaged Navier–Stokes (RANS) approach associated with a turbulent Prandtl number accounting for the turbulent heat flux based on Reynolds analogy is still, at least in the current state in most of the circumstances, the only feasible approach for practical engineering applications. However, the conventional choice of in the order of 0.9∼unity in many commercial computational fluid dynamics codes is not valid for the low-Prandtl number liquid metal. In this study, LES/DNS simulation results of a simple forced turbulent channel flow up to a friction Reynolds number of 2000 at of 0.01 and 0.025 were used as references, to which the Reynolds-averaged Navier–Stokes approach with varying was compared. It was found that the appropriate for the RANS approach decreases with bulk Peclet number and approaches a constant value of 1.5 when becomes larger than 2000. Based on this calibrated relation with , a new model for used in the RANS approach was proposed. Validation of the proposed model was carried out with available LES/DNS results on the local temperature profile in the concentric annulus and bare rod bundle, as well as with experimental correlations on the Nusselt number in a circular tube and bare rod bundle.
Highlights
1) Proposal of a new model for the turbulent Prandtl number in the RANS approach for engineering applications.
2) Validation of the proposed model with available LES/DNS results of the local temperature profile in the concentric annulus and bare rod bundle.
3) Validation of the proposed model with experimental correlations on the Nusselt number in a circular tube and bare rod bundle.
1 Introduction
Low-Prandtl number liquid metals are considered as promising candidate coolants in various innovative nuclear systems, such as the sodium-cooled fast reactor (SFR), the liquid lead–cooled fast reactor (LFR), and the accelerator-driven sub-critical system (ADS), due to their high molecular heat conductivity in favor of the reactors’ reliability and safety (Roelofs, 2019). With the fast development of computational power, more and more attention has been paid to applying the computational fluid dynamics (CFD) methods for detailed analysis of three-dimensional thermal–hydraulic behavior, especially in nuclear fuel assembly design (Cheng and Tak, 2006b; Marinari et al., 2019; Chai et al., 2019). However, modeling of turbulent heat flux in low-Prandtl number liquid metal flow remains a challenging task when using CFD methods to solve the turbulent flow heat transfer behavior (Shams, 2019). Advanced high-fidelity numerical approaches represented by the wall-resolved large eddy simulation (LES) method and the direct numerical simulation (DNS) method can provide detailed insight into the physics of the flow and the associated heat transfer (Tiselj et al., 2019). However, these approaches are limited to simple flow configurations and low to moderate Reynolds numbers due to their rather high requirement of computational capacities. Hence, Reynolds-averaged Navier–Stokes (RANS) method is still at least in the current state, in most cases, the only practically feasible approach to deal with high Reynolds industrial flows, especially those in complex geometries such as those encountered in typical nuclear fuel assemblies (Shams et al., 2019). It is thus important to assess and improve the accuracy of the RANS approach and the associated models for turbulent heat transfer. With this objective in mind, simulation results of the velocity and temperature field obtained by LES or DNS are, hence, very valuable references to which RANS models can be compared and calibrated.
Compared to conventional fluids with the Prandtl number in the order of unity, such as water or air, heat transfer in the low-Prandtl number liquid metal is characterized by the high contribution of molecular heat conduction even in high turbulent flow. Consequently, temperature change in the boundary layer is much smoother, even in the very thin laminar sub-layer. In the classical RANS approach with conventional turbulence models, turbulent heat transfer is often predicted solely from the knowledge of turbulent momentum transfer based on the concept of the so-called Reynolds analogy (Cheng and Tak, 2006a), in which the turbulent heat conductivity (or eddy conductivity) is given by the ratio between turbulent momentum conductivity (or eddy diffusivity) and a turbulent Prandtl number () according to:
It is well acknowledged that an accurate prediction of is of crucial importance in modeling turbulent heat transfer in the low-Prandtl number liquid metal flows. For instance, as pointed out by Cheng and Tak (2006b), a decrease in the Nusselt number of about 30% can be obtained by increasing from 0.9 to 1.5. However, the conventional choice of as a constant value of 0.85–0.9, with which satisfying results of the turbulent heat flux can be obtained for most of the engineering purposes with fluid of Prandtl number in the order of unity, is not valid for a low-Prandtl number liquid metal.
In the past, extensive studies were carried out to derive appropriate expressions for for the low-Prandtl number liquid metal. In general, two types of models can be found in the open literature, the first type specifies as a global value depending on the bulk flow parameters such as the bulk Reynolds number or/and the bulk Peclet number represented by the early model of Aoki (1963).
The models of Reynolds (1975).
Also, the model of Jischa and Rieke (1979).
The model by Cheng and Tak (2006b).
in which the constant A is given as:
The second type gives as a local varying value depending on the local flow parameters such as the local eddy diffusivity, represented by the model of Kays (1994).
where the turbulent Peclet number Pet is defined as:
Figure 1 compares the calculated with the model of Aoki (1963), Reynolds (1975), Jischa and Rieke (1979), and Cheng and Tak (2006b) for bulk Peclet number up to 2,800 at of 0.01 (corresponding to a up to 280,000). Also included in the figure is the proposed model for in the current study given by Eq. 21. It is observed that first, generally decreases with increasing . Second, tends to approach a constant value between 1 and 2 as becomes larger than 2000 except for the model of Cheng and Tak (2006b). predicted by the model of Cheng and Tak (2006b) is generally larger than all the other models. Finally, a rather scattering distribution exists in the predicted values of obtained by different models proposed in the literature, especially in the range of lower than 1,000. A possible reason for the scattering is due to the lack of reliable and consistent experimental data on turbulent heat transfer in liquid metal flows.
FIGURE 1

Comparison of turbulent Prandtl numbers obtained with different models (molecular Prandtl Pr = 0.01).
With the fast-growing computational capacities, advanced high-fidelity numerical approaches such as wall-resolved LES method and DNS method become more and more attractive to provide detailed insight into the physics of the flow and the associated heat transfer. Due to their high demands of computational capacities, LES and DNS simulations are often, at least in the current state, limited to rather simple geometrical configurations. Nevertheless, simulation results of LES and DNS are very valuable references to which RANS models can be compared and calibrated. In the current study, LES and DNS results of a simple fully developed forced turbulent channel flow of low-Prandtl number fluid and ) up to of were used as references to assess RANS modeling of the turbulent heat transfer by means of the turbulent Prandtl number concept. RANS simulations of the forced turbulent channel flow with boundary conditions following those in LES/DNS settings were performed by systematically varying from 0.9 up to 8.0. Based on the comparison of the dimensionless temperature field and/or bulk Nusselt number calculated by the RANS approach with those obtained by LES/DNS simulations, an appropriate choice of the associated with the RANS approach can be obtained and a new model for was then proposed and validated.
2 Large eddy simulation and direct numerical simulation results on fully developed forced turbulent channel flow of low-Prandtl number fluid
As depicted in Figure 2, a fully developed forced turbulent channel flow of low-Prandtl number fluid heated by a uniform heat flux on both walls is a simple scenario which has been extensively investigated in the open literature by means of LES and DNS. Table 1 summarizes the representative LES and DNS simulations in the ascending order of the bulk Peclet number . The channel flow is characterized by the friction Reynolds number which is defined with the friction velocity and the channel half height according to:
FIGURE 2

Sketch of the fully developed forced turbulent channel flow heated by a uniform heat flux on both sides of the walls.
TABLE 1
| Source | Pr [-] | Reτ [-] | Reb [-] | Peb [-] | Nub [-] | CFD method |
|---|---|---|---|---|---|---|
| Bricteux et al., (2012) | 0.01 | 180 | 5,600 | 56 | Not available | DNS |
| Kawamura et al., (1999) | 0.025 | 180 | 5,600 | 140 | Not available | DNS |
| Bricteux et al., (2012) | 0.01 | 590 | 22,000 | 220 | 6.02 | LES |
| Kawamura et al., (1999) | 0.025 | 395 | 13,500 | 337.5 | Not available | DNS |
| Duponcheel et al., (2014) | 0.01 | 2000 | 87,000 | 870 | 8.44 | LES |
| Abe et al., (2004) | 0.025 | 1,020 | 41,000 | 1,025 | Not available | DNS |
| Duponcheel et al., (2014) | 0.025 | 2000 | 87,000 | 2,175 | 14.39 | LES |
Summary of LES/DNS simulation on fully developed turbulent channel flow for low-Prandtl number fluid.
in which the friction velocity is defined with relation to the wall shear stress as:
where and stand for the fluid density and dynamic viscosity, respectively. The bulk velocity is an average channel flow velocity defined as:
which is then used to define the bulk Reynolds number as:
The bulk Peclet number is defined as:
Also included in Table 1 is the global heat transfer performance characterized by the bulk Nusselt number defined as:
where the bulk temperature is the average channel flow temperature which is defined as:
For the purpose of assessing RANS approach, DNS results of Kawamura et al. (1999) with and , DNS results of Abe et al. (2004) with and , as well as DNS results of Bricteux et al. (2012) with and , LES results of Duponcheel et al. (2014) with and were chosen as references. The corresponding bulk Reynolds number covers a wide range from 5,600 up to 87,000, while the bulk Peclet number is varied in the range from 56 up to 2,175.
3 Reynolds-averaged Navier–Stokes approach associated with turbulent Prandtl number concept
3.1 Geometrical setups and boundary conditions
As depicted in Figure 3, a quasi-2D RANS simulation of the fully developed channel flow was performed in the current study with the commercial CFD software package Ansys CFX. An arbitrary channel half height of 0.01 m (y-direction) was chosen, while the channel streamwise (x-direction) length was set as 80 times of . For a quasi-2D simulation, only 1 cell element was extruded in the spanwise direction (z-direction) and symmetrical conditions were specified on the two boundaries of the spanwise direction. A translational periodic boundary condition was specified on the two boundaries of the streamwise direction, in order to obtain the fully developed flow conditions. A constant uniform heat flux was given for the top and bottom boundaries following the LES/DNS settings.
FIGURE 3

Sketch of the mesh structure and boundary conditions used in RANS simulation of the fully developed forced turbulent channel flow.
Table 2 summarizes the most important boundary parameters of the RANS simulations performed in the current study. The flow is assumed to be incompressible with constant thermal–physical properties as in consistency with the LES/DNS settings. The governing equations are the RANS equations for incompressible flow with constant thermal–physical properties and the energy conservation equation for temperature. Due to the periodic boundary condition in the streamwise direction and the constant wall heat flux on the top and bottom walls, the temperature will continuously increase in the streamwise direction and a fully developed flow condition will never be achieved. Therefore, a negative volumetric heat source was specified in the RANS simulation, which takes exactly the same amount of energy away as given to the domain by the heat flux on the top and bottom walls. High-resolution advection scheme and turbulence numerics were chosen to ensure a second order accuracy. The buoyancy effect due to the temperature difference in the channel height direction (y-direction) is negligibly small, hence not considered in the RANS simulations. This is justified by the fact that the bulk Richardson number for the RANS setups according to the following equation is quite small (<<1):
TABLE 2
| Turbulent Prandtl number | 0.9, 1.5, 2.0, 4.0, 6.0, 8.0 |
|---|---|
| Turbulence model | , , SST |
| Mesh size: | 0.1, 0.2, 0.6 |
| Boundary type: streamwise (X) | Translational periodic |
| Boundary type: channel height (Y) | No-slip wall with uniform heat flux |
| Boundary type: spanwise (Z) | Symmetric |
Summary of boundary parameters for RANS simulations of fully developed forced turbulent channel flow for low-Prandtl number fluid.
where and are the gravitational acceleration and the thermal expansivity, respectively.
3.2 Presentation of the Reynolds-averaged Navier–Stokes results: normalization in friction units
The mean streamwise velocity profile and mean temperature profile in the channel height direction (y-direction) are normalized with the friction velocity and friction temperature , respectively. The friction velocity is already defined in Eq. 10 with relation to the wall shear stress . The friction temperature is defined with relation to the wall heat flux as follows:
where stands for the specific heat capacity at constant pressure.
The mean streamwise velocity profile is hence characterized by the normalized dimensionless mean velocity given as:
The mean temperature profile is characterized by the normalized dimensionless mean temperature according to:
where is the wall temperature of the heated top or bottom wall. Both the dimensionless velocity and dimensionless temperature profile are specified in terms of a dimensionless universal wall distance defined as:
where y is the actual distance to the wall.
3.3 Effect of mesh refinement and turbulence model
As also included in Figure 3, a block-structured mesh consisting of only hexahedra was used in RANS simulations. The mesh nodes are uniformly distributed among the streamwise direction (x-direction), while local refinement was specified in the y-direction when approaching the top and bottom wall, in order to ensure at least a dimensionless universal wall distance of the first interior node as required for the wall-resolved low-Reynolds number turbulence model in RANS simulations, that is, the model and the shear stress transport (SST) model. The mesh sensitivity study was performed for the LES case of and (Duponcheel et al., 2014) (corresponding and are 87,000 and 870, respectively). The wall-resolved turbulence model and a constant turbulent Prandtl number of 2.0 were used for all the tested meshes. For the mesh sensitivity study, three different levels of mesh refinement were studied with the same block structure as displayed in Figure 3. The total cell numbers of the three meshes were 3,713, 5,841, and 15,721, respectively, meaning that the highest density mesh is 4.23 times finer than the smallest-density mesh. The corresponding values of were 0.6, 0.2, and 0.1, respectively.
Figure 4 compares the normalized mean streamwise velocity profile () and the normalized mean temperature profile () obtained with the three tested meshes. It is observed that the mean streamwise velocity calculated with mesh 1 is slightly higher than that obtained with mesh 2 and mesh 3 in the near wall region with lower than 5, while the difference among the velocity profile obtained with mesh 2 and mesh 3 is negligibly small. Furthermore, the mean temperature profile obtained with the three tested meshes agrees well with each other, despite a large difference in the total cell number and mesh refinement. The bulk Nusselt number obtained by RANS simulations with the three different meshes are 8.81, 8.81, and 8.84, respectively. In conclusion, mesh 2 with of 0.2 was hence chosen for the RANS simulations henceforth in the current study.
FIGURE 4

Mesh sensitivity study: comparison of the (A) mean velocity profile and (B) mean temperature profile normalized in a friction unit obtained with different meshes.
Two types of turbulence models are available in CFX, that is, the wall-resolved low-Reynolds number turbulence model with which the boundary layer is fully resolved, and the wall-unresolved high-Reynolds number turbulence model with which the boundary layer is approximated relying on a logarithmic wall function. For the present analysis, the wall-resolved and turbulence model as well as the wall-unresolved turbulence model were tested. Figure 5 shows the mean streamwise velocity profile and the mean temperature profile normalized in fiction unit for the case of and . A constant turbulent Prandtl number of 2.0 was specified for all the RANS simulations. It is observed that the temperature profile is hardly influenced by the turbulence model. In the logarithmic region where becomes larger than 60, the velocity profile calculated by the different turbulence model agrees well with each other. The difference in the velocity profile in the near wall region () is to be expected, since the model applies the logarithmic wall function for the complete boundary layer. The bulk Nusselt numbers obtained with , , and SST model are 8.72, 8.81, and 8.84, respectively. The turbulence model was chosen for the RANS simulations henceforth in the current study.
FIGURE 5

Effect of different turbulence models on (A) mean velocity and, (B) mean temperature profile in friction unit.
3.4 Effect of turbulent Prandtl number
To study the effect of the turbulent Prandtl number
, the same case of
at
as in the mesh sensitivity study was further investigated. In the RANS simulations,
was systematically varied from 0.9, 1.2, 1.5, 2.0, 2.5, 3.0, 3.5 up to 4.0. The same mesh with
of 0.2 (mesh 2) and the same
turbulence model were employed. Since the velocity field is hardly affected by
for incompressible flow with constant thermal–physical properties,
Figure 6Ashows only the normalized mean temperature profile obtained by the RANS simulations. For a better demonstration of the results, only RANS results of
= 0.9, 2.0, and 4.0 are shown in the figure (the RANS simulation results with the optimum value of
= 2.3 are also included in
Figure 6Aand will be discussed later.). The corresponding LES results obtained by
Duponcheel et al. (2014)are shown as a straight line in the subfigure for comparison with RANS results. Also included in the figure is the theoretical linear law of the temperature profile
(dashed line) (
Kader, 1981). It is shown as follows:
1) Compared to conventional fluid with the Prandtl number in the order of unity, for which the linear law is valid only in the very thin laminar sub-layer of generally much smaller than 30, the theoretical linear law is still valid at a wall distance of up to 60∼70 for low-Prandtl number fluid. This indicates that the heat transfer in this region is dominated by the molecular effect of heat conduction. Consequently, no difference can be observed in the temperature profiles obtained with RANS of different and they all agree well with the temperature profile obtained by LES.
2) As becomes larger than 100, the theoretical linear law is no longer valid. The temperature profile given by the linear law is much steeper than that calculated by LES, which indicates that the turbulent diffusivity also becomes important in the heat transfer in addition to the molecular effect of heat conduction. It is clearly seen that the temperature profile now depends strongly on . The temperature profile obtained with the conventional choice of = 0.9 is much smoother than that calculated by LES, which indicates that the global heat transfer, that is, the heat transfer contribution due to turbulent diffusivity is obviously overestimated. With increasing from 0.9 to 1.5 and 4.0, the temperature profile becomes steeper and an optimum value of for the RANS approach exists apparently in the range of 1.5–4.0.
FIGURE 6

RANS approach with varying turbulent Prandtl numbers: (A) effect on normalized mean temperature; (B) effect on bulk Nusselt number.
In order to obtain the optimum value of for the investigated case of Pr = 0.01 at Reτ = 2000, the global heat transfer behavior was also studied. Figure 6B shows the bulk Nusselt number obtained in the RANS simulations with different (shown as line with circular symbols), compared with that obtained by LES (red straight line) of Duponcheel et al. (2014). It is observed that with = 0.9, the bulk Reynolds number was overestimated by about 40%. The conventional choice of is clearly not valid for low-Prandtl number fluid. The bulk Nusselt number decreases with increasing turbulent Prandtl number and an optimum value of is apparently located in the interval between 2.0 and 2.5. Therefore, further RANS simulations with = 2.1, 2.2, 2.3, and 2.4 were carried out and the optimum value was finally found at (shown as the diamond symbol in Figure 6B), for this value; the obtained by RANS simulation is equal to that obtained by LES. The temperature profile obtained with was also included in Figure 6A with triangular symbols and it agrees well with that obtained by LES (straight line).
4 Development of a new turbulent Prandtl number model for Reynolds-averaged Navier–Stokes approach
4.1 Proposal of a new turbulent Prandtl number model for Reynolds-averaged Navier–Stokes approach
The aforementioned approach to determine the optimum value of for the RANS approach was also conducted for further selected LES/DNS cases listed in Table 1, for which the bulk Nusselt numbers were reported in the respective sources, that is, for the case with at (Bricteux et al., 2012) and for the case with at (Duponcheel et al., 2014). For the purpose of developing a suitable model for , RANS simulations were also performed for the case with at (Bricteux et al., 2012) and for the case with at (Abe et al., 2004). Although no values of were provided in the DNS simulations of the aforementioned two cases, the optimum value of for the RANS simulation was found solely by comparing the temperature profile. Overall, RANS simulations were performed for five cases of different and (corresponding and covering the range from 5,600 to 87,000 and from 140 to 2,175, respectively). By comparing the temperature profile or/and bulk Nusselt number obtained in RANS simulations with those obtained in their respective LES/DNS simulations, the optimum values of for the RANS approach were determined for the five cases, as shown with the five circular symbols in Figure 7.
FIGURE 7

Development of a new model of the turbulent Prandtl number for RANS with LES/DNS results in turbulent channel flow as a reference.
It is observed that
1) Apparently, the turbulent Prandtl number can be correlated solely with the bulk Peclet number independent of Prandtl number .
2) Furthermore, decreases with increasing and tends to approach a constant value of 1.5 when becomes larger than 2000. This behavior is in general agreement with other models proposed in the literature as shown in Figure 1.
3) More importantly, the fact that approaches a constant value for large bulk Peclet number (corresponding also to a large bulk Reynolds number ) means that no further LES/DNS simulations are required for those large Reynolds numbers when developing a suitable model for .
4) Therefore, the following simple model for was proposed for the usage in the RANS approach, in which was assumed as a function of the bulk Peclet number (shown with the red straight line in Figure 7).
5) Since the proposed new model in Eq. 21 was derived based on LES/DNS results, the validity range of the new model should be limited to the range as specified in LES/DNS simulations as summarized in Table 1. Therefore, the validity range of should be from 56 to 2,175 with a validity range of from 0.025 to 0.01.
For assessment of the proposed model for the turbulent Prandtl number, RANS simulation was performed with calculated from the new model for the case of at (corresponding bulk Reynolds number is 5,600 and bulk Peclet number is 56) for which DNS results of the temperature profile are given in Bricteux et al. (2012), and for the case of at (corresponding is 13,500 and is 337.5), for which DNS results of the temperature profile are given in Kawamura et al. (1999). In the RANS simulations, was then calculated with the proposed model given by Eq. 21. For the case of at , the calculated with the proposed model is 4.15, while for the case of at the calculated is 7.98. Figure 8 compares the two cases of the normalized temperature profiles calculated by RANS simulations with those obtained by their respective DNS simulations. It is observed that, despite the fact that only five points were used to develop the new model for , the agreement between the RANS and DNS results is acceptable, which justifies the adequacy of the proposed model.
FIGURE 8

Assessment of the proposed model for the turbulent Prandtl number with DNS results for Pr = 0.025 at Reτ = 395, and for Pr = 0.01 at Reτ = 180: comparison of the mean temperature profile normalized in a friction unit.
4.2 Validation of the proposed turbulent Prandtl number model for Reynolds-averaged Navier–Stokes approach
The general strategy for validation of the proposed new model for the turbulent Prandtl number for the RANS approach is divided into the following two phases:
1) Phase I: taking LES/DNS results on the local temperature profile in the concentric annulus and bare rod bundle heated by constant uniform heat flux as a reference, with which RANS simulations results with the new model for will be compared. In the current study, we chose the LES/DNS simulation performed in a concentric annular channel by Lyu et al. (2015) and Marocco (2018) with a bulk Reynolds number of 8,900 at and the LES simulation performed in a bare rod bundle by Shams et al. (2018) with a bulk Reynolds number of at as a reference.
2) Phase II: LES/DNS simulations in the concentric annulus and bare rod bundle can provide detailed information about the local temperature profile, but are only available to a rather low-Reynolds number and Peclet number due to the enormous computational requirement. From the engineering point of view, the most important parameter to be considered is the heat transfer behavior in terms of the bulk Nusselt number . Therefore, in the second phase of the validation process, experimental correlations on the bulk Nusselt number developed for a circular tube and bare triangular rod bundle were taken as a reference, with which RANS simulations results with the new model for will be compared. The reason for selecting a circular tube and bare triangular rod bundle lies mainly in the fact, that most of the experimental investigations on turbulent heat transfer with the low-Prandtl number liquid metal were actually conducted with these two kinds of geometries and various experimental correlations were available in the open literature as summarized in OECD (2015).
4.2.1 Phase I: validation based on large eddy simulation/direct numerical simulation results for a local temperature profile in the concentric annulus and bare rod bundle
As the first step of the validation process phase I, a concentric annular channel heated on both walls as depicted in Figure 9A was considered due to its closeness to the sub-channel in a bare rod bundle (Ma et al., 2012). LES/DNS simulations on the fully developed state of turbulent heat transfer in the annular channel were performed at with a bulk Reynolds number by Lyu et al. (2015) and Marocco (2018). The ratio of the outer diameter to the inner diameter is 2. In the LES simulation, the annular channel was heated with a constant uniform heat flux on both walls. Since the flow and heat transfer situation on the inner wall is closer to that in a rod bundle, only the LES results normalized based on the parameters of the inner wall were taken as a reference.
FIGURE 9

(A) Sketch of a concentric annular channel heated uniformly on both walls and (B) mesh used in the RANS simulations.
Following the LES setting, RANS simulation was performed with the given by the new model proposed in the current study. With at , the corresponding bulk Peclet number is 231 and the turbulent Prandtl number calculated with the new model is hence 5.2. Lead–bismuth eutectic (LBE) with constant thermal–physical properties taken from the OECD Handbook on a heavy liquid metal (OECD, 2015) was chosen as a working fluid. An arbitrary inner diameter of 10 mm was chosen for the RANS simulation, while the outer diameter is 20 mm following the setting in the LES simulation. The hydraulic diameter of the annular channel is then 10 mm. The channel length in the flow direction was set as 250 times of the channel hydraulic meter in the RANS simulations, in order to assure the establishment of a fully developed state in the channel. Constant and uniform velocity and temperature were given to the inlet boundary, while a constant pressure was specified at the outlet boundary. Both inner and outer walls of the annular channel were modeled as no-slip walls for the velocity field and a constant uniform heat flux was imposed for the temperature field. Also included in Figure 9B is the detail of the mesh structure in RANS simulation. As recommended in the mesh sensitivity study in Section 3.3, the block-structured hexahedral mesh was used with local refinement close to both heated walls.
Figure 10 shows the fully developed normalized mean temperature profile in the radial direction of the annular channel. The wall distance was determined relative to the inner wall and the mean temperature was also normalized with friction temperature and wall temperature on the inner wall. It is observed that the temperature profile obtained by the RANS simulation with the new model for agrees well with that given by LES simulation in the literature (Lyu et al., 2015; Marocco, 2018). RANS simulation with the Kays model, which gives local varying depending on local flow parameters (Eq. 7), was also performed. As depicted in Figure 10, temperature profiles predicted by the Kays model and by the new model proposed in this study agree well with each other.
FIGURE 10

Validation of the proposed model on turbulent Prandtl number: a heated concentric annular channel.
A further validation of the proposed new model for the turbulent Prandtl number was performed in a loosely spaced bare rod bundle, for which wall-resolved large eddy simulation (LES) on the fully developed state of turbulent heat transfer with liquid lead at an operating temperature of 440°C (corresponding to a ) was performed at a bulk Reynolds number of 54,650 (Shams et al., 2018). As depicted in Figure 11A, the diameter (D) of the rod is 10.5 mm and the pitch (P) between the rods is 13.86 mm, corresponding to a pitch-to-diameter ratio of 1.32. The computational domain consists of two adjacent sub-channels with a hydraulic diameter of 9.67 mm. Liquid lead with constant thermal–physical properties at 440°C is considered as a working fluid. In the RANS simulation, the channel flow length was set at 300 times of the hydraulic diameter for the establishment of a fully developed state. All the rods were modeled as no-slip walls imposed with a constant uniform wall heat flux. Following the LES settings, the arrows of the same color (Figure 11A) indicate that the side planes have been connected with the periodic boundary conditions. Also included in Figure 11B is the block-structured mesh used in the RANS simulation, following the recommendation derived from the mesh sensitivity study in Section 3.3.
FIGURE 11

(A) Sketch of a hexagonal rod bundle heated uniformly on walls and (B) mesh used in the RANS simulations.
Figure 12 compares, then, the fully developed temperature profile in the gap region of the two sub-channels (as indicated with a red arrow in Figure 11A). It is observed that the temperature profile obtained by RANS simulation with the new model for agrees well with that provided by LES (Shams et al., 2018). Similar to the case in the annular channel, temperature profiles predicted by the Kays model and by the new model proposed in this study agree well with each other.
FIGURE 12

Validation of the proposed model on the turbulent Prandtl number: a hexagonal rod bundle.
4.2.2 Phase II: validation based on experimental correlations for the Nusselt number in a circular tube and bare rod bundle
The first step of validation phase II was performed in a circular tube, for this geometry has been studied intensively from an experimental point of view. Many correlations on the turbulent heat transfer behavior (in terms of the bulk Nusselt number
) are available in the open literature. However, a common agreement is not yet available, as often contradictory conclusions were reported by different authors (
OECD, 2015). We take, as recommended in the OECD Handbook on the heavy liquid metal (
OECD, 2015) based on a critical review of the available literature for the experimental data and proposed correlations, the following three correlations as references for comparison with RANS results. All the correlations predict the heat transfer behavior of liquid metal in a similar way as represented by the bulk Nusselt number
and bulk Peclet number
.
1) Correlation of Lyon (1949), Lyon (1951) as an upper limit for the bulk Nusselt number.
2) Correlation of Kutateladze et al. (1959) as a lower limit for the bulk Nusselt number.
3) Correlation of Notter and Sleicher (1972) as the best general applicable correlation for all the investigated liquid metals including sodium (Na) and sodium–potassium alloys (NaK), pure lead (Pb), and lead–bismuth eutectic (LBE) as well as pure mercury (Hg).
It should be noted that all the aforementioned correlations were proposed for the thermal boundary condition of the uniform wall heat flux and have been developed for fully developed turbulent flow conditions with bulk Reynolds numbers in the range between and . As depicted in Figure 13, RANS simulations were, hence, performed in a circular tube heated with a uniform wall heat flux. Constant uniform velocity and inlet temperature were given at the inlet boundary, while a constant pressure boundary condition was set at the outlet boundary. In order to achieve the fully developed state, the streamwise length of the tube was set as 250 times of the tube diameter. RANS simulations were performed for a lead–bismuth eutectic (LBE) with constant thermal–physical properties at . The bulk Reynolds number determined with the bulk flow velocity of the tube, covers a wide range between and , while the corresponding bulk Peclet number varies between 250 and 1,000. The block-structured mesh with local refinement toward the heated tube wall was used in all the RANS simulations, in accordance with the recommendation given by the mesh sensitivity study in Section 3.3.
FIGURE 13

(A) Sketch of a tube heated uniformly on the wall; (B) and (C) mesh used in the RANS simulations.
Figure 14 compares then the bulk Nusselt number obtained in RANS simulations with that calculated with the experimental correlations of Eqs. 22–24. With the conventional choice of , the bulk Nusselt number is largely overpredicted. RANS simulations with the Kays model for still overpredict the bulk Nusselt number, since a good agreement with the upper limit defined by the correlation of Lyon (1949), Lyon (1951) is observed. With the proposed new model for , on the other hand, a much better agreement with the experimental correlations is achieved. The bulk Nusselt number obtained with the proposed new model for lies within the range of the upper and lower limit defined by the correlation of Lyon (1949), Lyon (1951), and Kutateladze et al. (1959), respectively, and it agrees well with that given by the correlation of Notter and Sleicher (1972), which is recommended as the best general applicable correlation for all the investigated liquid metals including sodium (Na) and sodium–potassium alloys (NaK), pure lead (Pb) and lead–bismuth eutectic (LBE) as well as pure mercury (Hg).
FIGURE 14

Validation of the proposed model on the turbulent Prandtl number: heated tube.
In the second step of validation phase II, a bare rod bundle in the hexagonal arrangement was investigated, for which most of the experiments were performed in hexagonal arrays or bundles (
OECD, 2015). Based on a review of available experimental data, the following three correlations were chosen as references. All the correlations express the bulk Nusselt number
in a similar way in terms of the bulk Peclet number
and the pitch-to-diameter ratio (
).
1) Correlation of Gräber and Rieger (1972).
2) Correlation of Ushakov et al. (1977).
3) Correlation of Mikityuk (2009).
It should be noted, that the aforementioned Eqs. 25–27 were developed based on the experimental data of alkali metals such as liquid sodium (Na) or sodium–potassium alloy (NaK). However, as reviewed by Mikityuk (2009), the three correlations were recommended in the OECD Handbook on the heavy liquid metal (OECD, 2015) as the most relevant engineering heat transfer correlations for heavy liquid metal flows, that is, LBE or liquid lead. Similar to the correlations proposed for a circular tube, the aforementioned three correlations for the bare rod bundle are defined for the fully developed state. Table 3 summarizes the validity range of the respective correlations in terms of the pitch-to-diameter ratio and bulk Peclet number. However, it should be noted that in the review conducted by Mikityuk (2009), all the three correlations were recommended for use in the range of and bulk Peclet number of .
TABLE 3
| Source | P/D [-] | Peb [-] |
|---|---|---|
| (Gräber and Rieger, 1972) | 1.2–2.0 | 150–4,000 |
| Ushakov et al., (1977) | 1.3–2.0 | 1–4,000 |
| Mikityuk, (2009) | 1.1–1.95 | 30–5,000 |
Summary of experimental correlations on the bulk Nusselt number in a bare rod bundle of the hexagonal arrangement for the low-Prandtl number liquid metal.
Figure 15A defines the geometry investigated in the current study, where a triangular bundle configuration is shown. One-sixth of a sub-channel (a colored region in the figure) was defined as a computational domain, in which a block-structured hexahedral mesh was defined (Figure 15B). The definition of the boundary conditions was specified in Figure 15C. The rod was defined as a no-slip wall with a constant uniform heat flux, while the symmetry boundary condition was imposed on the other three side planes. In order to achieve a fully developed state, a domain length of 300 times of the hydraulic diameter was given in the flow direction. Constant uniform velocity and temperature were specified in the inlet boundary and a constant pressure condition was defined at the outlet boundary. LBE with constant thermal–physical properties (corresponding Prandtl number is 0.025) was chosen as a working fluid. In the current study, four pitch-to-diameter ratios , and were investigated. For each pitch-to-diameter ratio, RANS simulations were performed for bulk Reynolds numbers , and , respectively. The corresponding bulk Peclet numbers are , and , respectively. The turbulent Prandtl numbers for RANS simulations calculated with the new model are then , and , respectively.
FIGURE 15

(A) Sketch of a bare rod bundle uniformly on the wall; (B) mesh used in the RANS simulations. (C) boundary conditions is the RNAS simulations; (D) sketch of the streamwise domain length.
compare, then, the bulk Nusselt number obtained by RANS simulations with that given by the experimental correlations
Eqs. 25–27for
, and
, respectively. It is shown as follows:
1) The same trend of the bulk Nusselt number increasing with the bulk Peclet number but decreasing with the ratio is predicted by all the three experimental correlations, as well as by the RANS simulations.
2) For the tight-spaced bundle of , a large scattering exists among the experimental correlations. Although the values of calculated with the correlation of Gräber and Rieger (1972) and Ushakov et al. (1977) are still comparable, they are both much higher than those predicted by the correlation of Mikityuk (2009). obtained by the RANS approach with is higher than that obtained by the RANS approach with the new model for . But they all lie within the range defined by the three experimental correlations.
3) However, the situation is much different for the loosely spaced rod bundle of ratio larger than 1.2. predicted by the three correlations agrees well with each other. In general, for the loosely spaced rod bundle, the bulk Nusselt number obtained by the RANS approach with the conventional choice of is much higher than that predicted by the experimental correlations. With the new proposed model for , however, a much better agreement with the experimental correlations was able to be achieved. Compared to the RANS simulations with the Kays model, RANS simulations with the new model proposed in this study show a comparable yet slightly better agreement with the three experimental correlations.
FIGURE 16

Validation of the proposed model on turbulent Prandtl number: bare rod bundle of P/D ratio equals (A) 1.1; (B) 1.2; (C) 1.3 and (D) 1.5.
5 Conclusion and outlook
In this study, LES/DNS simulation results of a fully developed forced turbulent channel flow up to
at
and
were used as a reference, to which the RANS approach with the simple turbulent Prandtl number concept was compared and calibrated. Based on the calibrated relation with bulk Peclet number
, a new model for turbulent Prandtl number
used in the RANS approach was proposed and validated. The main conclusions derived are summarized as follows:
1) Heat transfer characteristics of the low-Prandtl number liquid metal depend strongly on the turbulent Prandtl number . An accurate estimation of is of essential importance when applying the RANS approach to simulate the turbulent heat transfer of liquid metal. The conventional choice of is proven to be not suitable for the liquid metal, with which the bulk Nusselt number could be overestimated by up to 40%.
2) Taking the LES/DNS simulation results of the fully developed forced turbulent channel flow as a reference, the RANS approach with the turbulent Prandtl number concept was assessed. It was found that the optimum value of used in the RANS approach for the low-Prandtl number liquid metal decreases with bulk Peclet Number and tends to approach a constant value of 1.5 as becomes larger than 2000.
3) Based on the aforementioned relation between and , a new model expressing solely in dependence on was proposed in this study. Validation of the proposed model was carried out with available LES/DNS results of a local temperature profile in a concentric annulus and a loosely spaced bare rod bundle, as well as with experimental correlations on bulk Nusselt number in a circular tube and bare rod bundle in a triangular arrangement. An overall good agreement of the RANS simulation results with the LES/DNS results, as well as with the experimental correlations was achieved, which demonstrates the adequacy of the proposed new model for in this study.
4) To summarize, it could be concluded that the RANS approach with the simple concept of a constant global turbulent Prandtl number is a suitable tool for simulating the forced convective turbulent heat transfer with the low-Prandtl number liquid metal, provided an appropriate is specified. The RANS approach with the turbulent Prandtl number concept is, and will still be the only feasible approach when dealing with industrial application of turbulent heat transfer with the low-Prandtl number liquid metal. Therefore, it is believed that the model developed in the current study will have a promising perspective for engineering applications.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
XH—writing–original draft preparation; methodology; investigation; data curation; visualization. BP—conceptualization; methodology; investigation; writing–original draft preparation; writing–review and editing; project administration; funding acquisition. XC—software; investigation; writing–review and editing YY—methodology; project administration; writing–review and editing.
Funding
This work was financially supported by the Natural Science Foundation of Guangdong Province (2021A1515010391).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
AbeH.KawamuraH.MatsuoY. (2004). Surface Heat-Flux Fluctuations in a Turbulent Channel Flow up to Reτ=1020 with Pr=0.025 and 0.71. Int. J. Heat Fluid Flow25 (3), 404–419. 10.1016/j.ijheatfluidflow.2004.02.010
2
AokiS. (1963). A Consideration on the Heat Transfer in Liquid Metal. Tokyo, Japan: Bulletin of the Tokyo Institute of Technology, 63–73.
3
BricteuxL.DuponcheelM.WinckelmansG.TiseljI.BartosiewiczY. (2012). Direct and Large Eddy Simulation of Turbulent Heat Transfer at Very Low Prandtl Number: Application to Lead-Bismuth Flows. Nucl. Eng. Des.246, 91–97. 10.1016/j.nucengdes.2011.07.010
4
ChaiX.LiuX.XiongJ.ChengX. (2019). Numerical Investigation of Thermal-Hydraulic Behaviors in a LBE-Cooled 19-pin Wire-Wrapped Rod Bundle. Prog. Nucl. Energy119, 103044. In press, corrected proof, Available online 1 May 2019, Article 103044. 10.1016/j.pnucene.2019.103044
5
ChengX.TakN.-I. (2006). Investigation on Turbulent Heat Transfer to Lead-Bismuth Eutectic Flows in Circular Tubes for Nuclear Applications. Nucl. Eng. Des.236, 385–393. 10.1016/j.nucengdes.2005.09.006
6
DuponcheelM.BricteuxL.ManconiM.WinckelmansG.BartosiewiczY. (2014). Assessment of RANS and Improved Near-Wall Modeling for Forced Convection at Low Prandtl Numbers Based on LES up to Reτ=2000. Int. J. Heat Mass Transf.75, 470–482. 10.1016/j.ijheatmasstransfer.2014.03.080
7
GraberH.RiegerM. (1972). Experimentelle Untersuchung des Wärmebergangs an Flüssigmetalle (NaK) in parallel durchströmten Rohrbündeln bei konstanter und exponentieller Wärmeflussdichteverteilung. Atomkernenergie19 (1), 23–40.
8
JischaM.RiekeH. B. (1979). About the Prediction of Turbulent Prandtl and Schmidt Numbers from Modeled Transport Equations. Int. J. Heat Mass Transf.22, 1547–1555. 10.1016/0017-9310(79)90134-0
9
KaderB. A. (1981). Temperature and Concentration Profiles in Fully Turbulent Boundary Layers. Int. J. Heat Mass Transf.24, 1541–1544. 10.1016/0017-9310(81)90220-9
10
KawamuraH.AbeH.MatsuoY. (1999). DNS of Turbulent Heat Transfer in Channel Flow with Respect to Reynolds and Prandtl Number Effects. Int. J. Heat Fluid Flow20, 196–207. 10.1016/s0142-727x(99)00014-4
11
KaysW. M. (1994). Turbulent Prandtl Number-Where Are We?J. Heat Transf.116, 284–295. 10.1115/1.2911398
12
KirillovP.UshakovP. A. (2001). Heat Transfer to Liquid Metals: Specific Features, Methods of Investigation, and Main Relationships. Therm. Eng.48 (1), 50–59.
13
KutateladzeS.BorishanskiiV.NovikovI. (1959). Heat Transfer in Liquid Metals. J. Nucl. Energy, Part B, React. Technol.9 (1-4), 214–229. 10.1016/0368-3265(59)90177-x
14
LyonR. (1949). Forced Convection Heat Transfer Theory and Experiments with Liquid Metals. Oak Ridge, Tennessee: Oak Ridge National Laboratory.
15
LyonR. (1951). Liquid Metal Heat-Transfer Coefficients. Chem. Eng. Prog.47 (2), 75–79. 10.1080/00223131.2014.980349
16
LyuY.GeZ.ZhaoP. (2015). Large Eddy Simulation for the Turbulent Heat Transfer of Liquid Metal in an Annulus. J. Univ. Sci. Technol. China45 (11), 917–922. 10.3969/j.issn.0253-2778.2015.11.006
17
MaZ.WuY.QiuZ.TianW.SuG.QiuS. (2012). An Innovative Method for Prediction of Liquid Metal Heat Transfer Rate for Rod Bundles Based on Annuli. Ann. Nucl. Energy47, 91–97. 10.1016/j.anucene.2012.04.023
18
MarinariR.Di PiazzaI.TarantinoM.AngelucciM.MartelliD. (2019). Experimental Tests and Post-test Analysis of Non-uniformly Heated 19-pins Fuel Bundle Cooled by Heavy Liquid Metal. Nucl. Eng. Des.343, 166–177. 10.1016/j.nucengdes.2018.12.024
19
MaroccoL. (2018). Hybrid LES/DNS of Turbulent Forced and Aided Mixed Convection to a Liquid Metal Flowing in a Vertical Concentric Annulus. Int. J. Heat Mass Transf.121, 488–502. 10.1016/j.ijheatmasstransfer.2018.01.006
20
MikityukK. (2009). Heat Transfer to Liquid Metal: Review of Data and Correlations for Tube Bundles. Nucl. Eng. Des.239 (4), 680–687. 10.1016/j.nucengdes.2008.12.014
21
NotterR. H.SleicherC. A. (1972). A Solution to the Turbulent Graetz Problem-III Fully Developed and Entry Region Heat Transfer Rates. Chem. Eng. Sci.27 (11), 2073–2093. 10.1016/0009-2509(72)87065-9
22
OECD (2015). Handbook on Lead-Bismuth Eutectic Alloy and Lead Properties, Materials Compatibility, Thermal-hydraulics and Technologies. Paris, France No: OECD/NEA, 7268.
23
ReynoldsA. J. (1975). The Prediction of Turbulent Prandtl and Schmidt Numbers. Int. J. Heat Mass Transf.18, 1055–1069. 10.1016/0017-9310(75)90223-9
24
RoelofsF. (2019). “Chapter 1: Introduction to Liquid Metal Cooled Reactors,” in Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors. Editor RoelofsF. (Cambridge: Woodhead Publishing, Elsevier Ltd).
25
ShamsA. (2019). “Chapter 6.2.1 Turbulent Heat Transport,” in Thermal Hydraulics Aspects of Liquid Metal Cooled Nuclear Reactors. Editor RoelofsF. (Cambridge: Woodhead Publishing, Elsevier Ltd).
26
ShamsA.De SantisA.KoloszarL. K.Villa OritzA.NarayananC. (2019). Status and Perspectives of Turbulent Heat Transfer Modelling in Low-Prandtl Number Fluids. Nucl. Eng. Des.353, 110220. 10.1016/j.nucengdes.2019.110220
27
ShamsA.MikužB.RoelofsF. (2018). Numerical Prediction of Flow and Heat Transfer in a Loosely Spaced Bare Rod Bundle. Int. J. Heat Fluid Flow73, 42–62. 10.1016/j.ijheatfluidflow.2018.07.006
28
TiseljI.FlageulC.OderJ. (2019). Direct Numerical Simulation and Wall-Resolved Large Eddy Simulation in Nuclear Thermal Hydraulics. Nucl. Technol.206, 164–178. 10.1080/00295450.2019.1614381
29
UshakovP.ZhukovA.MatyukhinN. (1977). Heat Transfer to Liquid Metals in Regular Arrays of Fuel Elements. High. Temp. (USSR)15 (10), 1027–1033.
Summary
Keywords
low-Prandtl number liquid metal, turbulent heat transfer, turbulent Prandtl number, RANS, CFD
Citation
Huang X, Pang B, Chai X and Yin Y (2022) Proposal of a turbulent Prandtl number model for Reynolds-averaged Navier–Stokes approach on the modeling of turbulent heat transfer of low-Prandtl number liquid metal. Front. Energy Res. 10:928693. doi: 10.3389/fenrg.2022.928693
Received
26 April 2022
Accepted
27 June 2022
Published
18 July 2022
Volume
10 - 2022
Edited by
Wenzhong Zhou, Sun Yat-sen University, China
Reviewed by
Yaou Shen, Laboratory of Reactor System Design Technology (LRSDT), China
Hui Cheng, Sun Yat-sen University, China
Updates
Copyright
© 2022 Huang, Pang, Chai and Yin.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Bo Pang, bo.pang@szu.edu.cn
This article was submitted to Nuclear Energy, a section of the journal Frontiers in Energy Research
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.