Numerical Study of Low Pr Flow in a Bare 19-Rod Bundle Based on an Advanced Turbulent Heat Transfer Model

Compared to assuming a constant turbulent Prandtl number model, an advanced four-equation model has the potential to improve the numerical heat transfer calculation accuracy of low–Prandtl number ( P r ) fluids. Generally, a four-equation model consists of a two-equation k − ε turbulence model and a two-equation k θ − ε θ heat transfer model. It is essential to analyze the influence of dissimilar turbulence models on the overall calculation accuracy of the four-equation model. The present study aims to study the effect of using different turbulence models on the same k θ − ε θ heat transfer model. First, based on the open-source computational fluid dynamics software OpenFOAM, an advanced two-equation k θ − ε θ heat transfer model was introduced into the solver buoyant2eqnFoam, which was developed based on the self-solver buoyantSimpleFoam of OpenFOAM. In the solver buoyant2eqnFoam, various turbulence models built into OpenFOAM can be conveniently called to close the Reynolds stress and an advanced two-equation heat transfer model can be utilized to calculate the Reynolds heat flux of low- P r fluids. Subsequently, the solver buoyant2eqnFoam was employed to study the fully developed flow heat transfer of low- P r fluids in a bare 19-rod bundle. The numerical results were compared and analyzed with the experimental correlations and other simulation results to validate the effectiveness and feasibility of the solver buoyant2eqnFoam. Furthermore, the influence of combining different turbulence models with the same two-equation k θ − ε θ heat transfer model was also presented in this study. The results show that the turbulence model has a considerable influence on the prediction of turbulent heat transfer in the high Peclet number range, suggesting that it should be prudent when picking a turbulence model in the simulations of low- P r fluids.


INTRODUCTION
Nuclear energy is playing an increasingly irreplaceable role in the future energy structure as the demand for energy increases rapidly (Gu and Su, 2021). Lead-cooled fast reactor (LFR) is one of the six types of innovative nuclear power systems proposed by the Generation IV International Forum (Abram and Ion, 2008;Pacio et al., 2015). Benefiting from the excellent performance in chemical inertness, neutron economy, and thermohydraulic properties, lead-bismuth eutectic (LBE) is considered as one of most promising coolants for LFR. It is indispensable to research the thermohydraulic behaviors of the LBE inner fuel assembly, which influences the security and economic performance of LFRs but is poorly understood (Martelli et al., 2017;Pacio et al., 2017).
Since it is expensive, parlous, and complicated to conduct an experiment with LBE under a high-temperature state, computational fluid dynamics (CFD) methods are widely employed to study the thermohydraulic characteristics of LBE. The CFD methods can be subdivided into three categories: direct numerical simulation (DNS), large eddy simulation (LES), and Reynolds-averaged Navier-Stokes simulation (RANS). Despite the high calculation accuracy of DNS and LES, they have a high demand for computational resources, and as a result, they are only suitable for some specific and straightforward geometric models (Kawamura et al., 1999). Since the computational cost of the RANS approach is much lower than that of the DNS and LES, the RANS approach is the most widely adopted CFD method in engineering calculation. In the RANS method, the linear eddyviscosity k − ε or k − ω turbulence model is ordinarily sufficient to accurately predict the momentum transport of various fluids (Nagano, 2002). On the other hand, for reproducing the heat transfer, the Reynolds-analogy hypothesis assuming a constant turbulent Prandtl number Pr t 0.85~0.9 is adopted in almost all commercial codes (Manservisi and Menghini, 2014a). For the simulation of ordinary fluids like water and air, having a relatively high Prandtl number, the rational results can be obtained with a constant Pr t (He et al., 2021). However, LBE is characterized by high thermal diffusivity and low viscosity values, resulting in low Prandtl numbers (Pr ≈ 0.01~0.03). Consequently, the Reynolds-analogy hypothesis is no longer appropriate to be employed to study the thermohydraulic characteristics of LBE by the CFD methods (Cheng and Tak, 2006). For this reason, some advanced turbulent heat transfer models which can reproduce the heat transfer behaviors of LBE with high precision are highly desirable.
In the past four decades, to improve the calculation accuracy of heat transfer for low-Pr fluids, various heat flux models to close the energy conservation equation in the framework of RANS have been developed.

Differential Heat Flux Model
DHFM is a full second-moment differential model for the transport of Reynolds heat fluxes. Compared with the constant Pr t model, DHFM fully considers the convection, diffusion, generation, and dissipation terms of Reynolds heat flux in the differential equations. Carteciano (1995), Carteciano et al. (1997), Carteciano et al. (2001), and Carteciano and Grötzbach (2003) developed a kind of DHFM named turbulence model for buoyant flows (TMBF). The simulations of two-dimensional forced convection and mixed convection with different fluids were carried out to evaluate the accuracy of TMBF. The numerical results obtained by TMBF demonstrate that stratified flows and buoyant effects were well reproduced compared with the constant Pr t model, especially in mixed convection conditions. Based on a summary of the various DHFM models developed in recent years, Shin et al. (2008) proposed a new set of DHFM models with an

Algebraic Heat Flux Model
AHFM is a simplified second-moment form of DHFM, which transports Reynolds heat flux by establishing algebraic equations. Hanjalić et al. (1996), Kenjereš et al. (2005), Otić et al. (2005), and Otić and Grotzbach. (2007) developed and analyzed an implicit algebraic transport equation for the Reynolds heat flux term to close the energy equation. Evaluations and calibrations of AHFM for low-Pr fluids were implemented by Shams et al. (2014), Shams (2018), Shams et al. (2018), and and De Santis and Shams (2018). AHFM has been validated in their works by comparison with the DNS data for turbulent flows in forced, mixed, and natural convection of different fluids. The numerical results obtained by AHFM illustrate that temperature, heat flux field, buoyant effects in-plane, backward-facing step, corium pool, rod bundle, etc. are well predicted.

A two-equation k θ − ε θ model for Reynolds heat flux
The k θ − ε θ model is a first-order 2-equation model for the calculation of Reynolds heat flux, which can be developed in a way similar to that of a first-order 2-equation k − ε model for the turbulent transport of momentum formulated. Compared with DHFM and AHFM, the two-equation model has been widely applied in recent years because of its lower calculation cost. To precisely reproduce the heat transfer of low-Pr fluids, extensive contributions of model coefficient and function, wall boundary conditions, and near-wall thermal turbulence effect of a two-equation k θ − ε θ model had been made by Nagano and Kim (1988); Nagano and Shimada (1996); Nagano et al. (1997);Nagano (2002), Sommer et al. (1992), and Youssef et al. (1992); and Youssef (2006), Abe et al. (1994), Hattori et al. (1993), Hwang and Lin (1999), Deng et al. (2001), and Karcz and Badur (2005). In recent years, based on previous study, Manservisi and Menghini (2014a) (Manservisi and Menghini, 2014a;Manservisi and Menghini, 2014b;Manservisi and Menghini, 2015;Su et al., 2022) indicate that turbulent heat transfer statistics such as in-plane, tube, backward-facing step, triangular rod bundle, square lattice bare rod bundle, and hexagonal rod bundle of forced convection of LBE are well reproduced based on their k θ − ε θ turbulent heat transfer model.
In recent years, the interest in reliable CFD methods used to investigate the turbulent heat transfer of low-Pr fluids in complicated industrial configurations has increased dramatically. Nevertheless, commercial codes are still lacking, except for an AHFM model available on software STAR-CCM+ (Simcenter, 2016). A two-equation k − ε model utilized to calculate Reynolds stress with a twoequation k θ − ε θ model used to calculate Reynolds heat flux is usually called a four-equation k − ε − k θ − ε θ model, which is expected to improve the numerical CFD accuracy of turbulent heat transfer for LBE. However, different turbulence models will have a certain impact on the time-scale transport of a twoequation k θ − ε θ model. It is necessary to evaluate the sensitivity of various turbulence models to a two-equation k θ − ε θ model.
Thus, in the present study, an improved CFD solver buoyant2eqnFoam, which introduces a two-equation k θ − ε θ model to calculate the Reynolds heat flux and can directly call different turbulence models to calculate the Reynolds stress, was first developed based on the solver buoyantSimpleFoam of opensource CFD program OpenFOAM. The fully developed turbulent heat transfer results of LBE inner flow and a bare 19-rod bundle geometry with different Peclet numbers were investigated and Frontiers in Energy Research | www.frontiersin.org June 2022 | Volume 10 | Article 922169 compared with experimental relations to verify the effectiveness of the solver buoyant2eqnFoam and the numerical algorithm. Finally, the heat transfer sensitivity of different turbulence models to the twoequation k θ − ε θ model was presented.

Physical Model
Thermal-hydraulic phenomena in a 19-rod bundle geometry are an essential research topic. In the past decades, numerous experimental and simulation researches have been conducted to precisely obtain flow characteristics and heat transfer correlations of coolant (Pacio et al., 2014;Martelli et al., 2017). In the present study, a bare 19-rod bundle with a fully developed turbulent LBE flow is considered. Figure 1A displays the cross section of the bare 19-rod bundle. Since the cross-flow in the bare 19-rod bundle is negligible and its construction is symmetrical, one-twelfth of the whole bundle is selected to carry out simulation for the sake of economizing computational cost. The computational domain is sketched in Figure 1B, together with the definitions of sub-channels and boundary regions. Sub1, Sub2, and Sub3 are defined as inner subchannels, while Sub4 and Sub5 are the edge sub-channel and corner sub-channel, respectively. The more detailed geometric parameters are summarized in Table 1, which are consistent with Pacio's experiment (Pacio et al., 2015), except that there are no grid spacers in the current study. D h,bun and D h,sub1 are the hydraulic diameter of the bundle and hydraulic diameter of Sub1, respectively. The length of the whole computational domain is set to 15 D h,bun to eliminate the effect of employing periodic inlet boundary conditions. The flow parameters of LBE are reported in Table 2, where Pe bun PrU bun D h,bun ρ/μ is the Peclet number of the bundle.

Conservation Equations
For forced convection, the incompressible RANS equations with constant physical properties and no gravity are considered where ], u i , and P are the molecular viscosity, Reynolds-averaged velocity, and the so-called average pressure, respectively. To obtain the unknown Reynolds stress u′ i u′ j , the linear eddyviscosity model can be adopted as follows: where k and ] t , both derived from the turbulence model, represent the turbulent kinetic energy and the turbulent viscosity, respectively. It should be noted that, in OpenFOAM, the energy conservation equation can be expressed in terms of enthalpy (Darwish and Moukalled, 2021): where h, K |U| 2 /2, α, and α t are the enthalpy per unit of mass, the kinetic energy per unit of mass, molecular thermal diffusivity, and the turbulent thermal diffusivity, respectively. The unknown α t needs to be derived from a two-equation k θ − ε θ turbulent heat transfer model. After solving Eq. 4, the distribution of h can be obtained. Subsequently, the Reynolds-averaged temperature T can be calculated by using the function Thermo.T () coming with OpenFOAM. The derivation of periodic momentum and energy equations in OpenFOAM can be found in the reference (Ge et al., 2017).

Turbulence Model for Momentum Field
Benefiting from replacing the friction velocity u τ with Kolmogorov velocity u ε , the turbulence model proposed by Abe et al. (1994), which can well reproduce the low Reynolds number and near-wall effects of both separated and attached flows, was widely adopted in the calculation of LBE (Manservisi and Menghini, 2014a;Manservisi and Menghini, 2014b;Cerroni et al., 2015;Manservisi and Menghini, 2015;Da Via et al., 2016;Chierici et al., 2019;Da Vià and Manservisi, 2019;Cervone et al., 2020;Da Vià et al., 2020). However, the Abe k − ε turbulence model does not exist in the current turbulence model library of OpenFOAM. Therefore, in the current study, the Abe k − ε turbulence model is compiled into the turbulence model library that comes with OpenFOAM so as to utilize the wall functions of OpenFOAM in this self-compiled turbulence model. In the Abe k − ε turbulence model, the turbulent viscosity ] t is computed as follows: v where C μ is a constant. f μ is the model function, defined as follows: where R t k 2 /]ε and u ε (]ε) 1/4 . Moreover, δ is the distance from the wall. The equations for k and its dissipation ε can be written as follows: The model constants utilized in the Abe k − ε turbulence model are reported in Table 3.

Two-Equation Model for Thermal Field
In the current work, the k θ − ε θ turbulent heat transfer model developed and improved by Manservisi and Menghini (2014a), Manservisi and Menghini (2014b), and Manservisi and Menghini (2015), which introduces the average square temperature fluctuation k θ and its dissipation ε θ in order to well reproduce the near-wall turbulent heat transfer behaviors of LBE having Pr in the range of 0.01-0.03, is adopted to calculate the turbulent thermal diffusivity α t . In the Manservisi k θ − ε θ model, α t is computed as follows: where C θ is the constant empirical coefficient and τ lθ is the local thermal characteristic time, modeled as follows: (13) with the appropriate functions set as follows: where R δ δε 1/4 /] 3/4 , τ u k/ε, and R τ θ /τ u with the thermal turbulent characteristic time τ θ k θ /ε θ . In addition, τ u k/ε represents the dynamical turbulent characteristic time. The equations for k θ and its dissipation ε θ can be written as follows: The constant empirical coefficients used in the Manservisi k θ − ε θ turbulent heat transfer are reported in Table 4.

SOLVER AND BOUNDARY CONDITIONS
To calculate the thermal-hydraulic characteristics of LBE, a CFD solver named buoyant2eqnFoam was developed on the OpenFOAM platform having user-friendly programming language features based on the turbulence model and the aforementioned turbulent heat transfer model. The SIMPLE algorithm is adopted to handle pressure-velocity coupling equations and the coupled multigrid iterations technique is utilized for matrix solutions. All calculations were performed using double precision on OpenFOAM and the convergence conditions of residual error are set as follows: Max Q n+1 Q n − 1 < 10 −6 (23)  where Q stands for u i , T, P, k θ , ε θ , k, and ε. The index i represents the number of iterations. The framework of the buoyant2eqnFoam solver is presented in Figure 2. The boundary condition data, mesh data, physical property data, calculation control, discrete format of each differential operator, algebraic equation solver, and relaxation factor required by buoyant2eqnFoam to perform calculation are included in the 0 folder, constant/polyMesh, constant/ thermophysicalProperties, system/controlDict, system/ fvSchemes, and system/fvSolutions.
In the computational domain, periodic boundary conditions are set on the region of inlet and outlet, considering the fully developed turbulent inner flow in the bundle. It is worth noting that the energy source term needs to be added to the energy Eq. 4 in order to apply periodic boundary conditions to temperature variables. The calculation method of energy source term refers to this literature (Ge et al., 2017). For k θ and ε θ , the boundary condition zeroGradient is employed on the wall under the uniform heat flux condition, according to the research of Deng et al. (2001). The boundary conditions imposed on each boundary are summarized in Table 5. Since the wall functions kLowReWallFunction for k and epsilonWallFunction for ε are both suitable for the low-Reynolds number turbulence model and can well reproduce the near-wall turbulence behaviors when y + is very low (Darwish and Moukalled, 2021), they are employed in this study. Given that there are no wall functions accessible for k θ and ε θ in OpenFOAM, y + must be less than or equal to 1 in order to accurately reproduce the thermal turbulent behaviors near the wall (Manservisi and Menghini, 2014a).

Mesh Independence Analysis
In this section, the buoyant2eqnFoam, which utilizes the Abe k − ε turbulence model for turbulence fields and uses the Manservisi k θ − ε θ model for thermal fields, is employed to investigate the thermohydraulic characteristics of LBE inner flow in the bare 19rod bundle in a wide range of Pe bun . As shown in Figure 3, the computational domain was discretized by GAMBIT unstructured meshes (tetrahedral and hexahedral mesh blending). The first layer grid was set with a height of 0.001 mm in order to satisfy the criterion of the low-Reynolds number turbulence model for y + ≤ 1. A total of 15 layers of boundary grids with a height ratio of 1.3 were designed. Three sets of mesh with different mesh numbers of 2.05 million, 2.64 million, and 3.11 million were adopted to analyze the mesh sensitivity. The dimensionless coolant temperature Θ is defined as follows: where T b,bun is the bulk temperature of the computational domain. The dimensionless temperature profiles along line ab (shown in Figure 1B) of three sets of mesh are displayed in Figure 4 under Pe bun 1500. It is evident that the difference in the dimensionless temperature profile between three sets of mesh is negligible. Consequently, the mesh with a mesh number of 2.05 million is selected, taking the calculation cost into consideration.

Solver Verification
The fully developed turbulent heat transfer characteristics of LBE inner flow in the bare 19-rod bundle were studied by Chierici et al. (2019), using a four-equation model in logarithmic specific dissipation form k − Ω − k θ − Ω θ , which was developed based on the Abe k − ε turbulence model and the Manservisi k θ − ε θ model. The numerical results of Chierici et al. (2019) can provide some reference for developing a CFD solver of LBE turbulent heat transfer. Therefore, the simulation results of Chierici et al. (2019) and some experimental data are picked for comparison to verify the validity of the solver buoyant2eqnFoam. The Nusselt number is selected for comparison since it is a critical parameter in engineering. Table 6 presents some Nusselt number experimental correlations of the triangular  rod bundle channel cooled by liquid metal, obtained by Subbotin et al. (1965), Mikityuk (2009), Gräber and Rieger (1972), and Mareska and Dwyer (1964), respectively.
Because these reported correlations were developed for triangular lattices, the Nusselt number of inner sub-channel Sub1 is picked for comparison. The Nu sub1 is calculated as follows: where T w,sub1 and T b,sub1 are the mean wall temperature and mean coolant temperature of Sub1, respectively. Correspondingly, Pe sub1 PrU sub1 D h,sub1 ρ/μ is the Peclet number of the inner sub-channel Sub1. Figure 5 displays the comparison of the Nusselt number with Chierici simulation and experimental correlations. From this figure, it can be clearly observed that the tendency of Nu sub1 is consistent with the simulation results of Chierici and shows good agreements with experimental data in a specific Peclet number range, illustrating that the rational prediction of LBE turbulent heat transfer can be obtained by the self-compiled solver buoyant2eqnFoam which can use the Abe k − ε turbulence model with wall functions for turbulence fields and the Manservisi k θ − ε θ model with zero-gradient boundary for thermal fields.

Velocity Field
The profiles of velocity magnitude on the computational domain of the bare 19-rod bundle are reported in Figure 6, with Pe bun 500 and Pe bun 1500, respectively. It is obvious that the mean velocity in the inner sub-channels is higher than that in the edge sub-channel Sub4 and corner sub-channel Sub5 because the hydraulic diameter of inner sub-channels is relatively higher, resulting in much lower flow resistance. For the same reason, the average velocity of Sub5 is much lower compared with that of other sub-channels, as summarized in Table 7. Figure 7 shows the distribution of dimensionless temperature from where it can be seen that the maximum temperature is located in the corner sub-channel Sub5, which is mainly due to the lower mean coolant velocity of the corner sub-channel  Sub5. Owing to the larger hydraulic diameter of the edge subchannel Sub4 and the adiabatic boundary condition applied on the outer casing wall, the coolant with the lowest temperature can be found in the area of the edge sub-channel Sub4 near the outer casing wall. Comparing the four cases reported in Figure 7, it can be found that the convective heat transfer of the coolant in each sub-channel is enhanced as the Reynolds number increases, leading to the decrease of maximum temperature and the increase of bulk coolant temperature.

Dimensionless Temperature and Hot Spot Factor Distributions
The dimensionless hot spot factor characterizing the inhomogeneity of wall temperature is defined as follows: where θ wmax,sub , θ b,sub , and θ wb,sub are the maximum wall temperature of the sub-channel, the bulk temperature of the sub-channel, and the mean wall temperature of the sub-channel, respectively. As ϕ 1, it means that the maximum wall temperature of the sub-channel is equal to the average wall temperature. The calculated results of the hot spot factor of each sub-channel are summarized in Table 8, from where it can be deduced that the wall temperature distribution of the inner sub-channel Sub1 is the most homogeneous. On the other hand, due to the coexistence of the heated rod wall and adiabatic wall in the edge sub-channel Sub4 and the corner sub-channel Sub5, the phenomenon of nonhomogeneous wall temperature distribution in these channels is more dramatic.

Dimensionless Thermal Diffusivity Distribution
To analyze the dependence of heat transfer on Pe bun , defining the dimensionless thermal diffusivity α + as follows: α + is the ratio between turbulent thermal diffusivity and molecular thermal diffusivity. Figure 8 reports the calculated dimensionless thermal diffusivity distribution on the computational domain for Pe bun 500 and Pe bun 1500. From this figure, it can be clearly seen that the α + in the center of each sub-channel is higher than that near the wall,  where the heat is mainly derived by the molecular heat conduction. Moreover, in the center of the edge sub-channel Sub4, the maximum dimensionless thermal diffusivity can be found, indicating that the thermal diffusion caused by turbulent flow reaches its peak in this region. It should be mentioned that when Pe bun 500, the α + in the whole computational domain is less than 1, suggesting that the molecular heat conduction affects the entire computational domain dominantly. In addition, with FIGURE 11 | Comparison of Nu sub1 calculated by different turbulence models.
FIGURE 12 | Comparison of Nu sub1 calculated by different turbulent heat transfer models.
Frontiers in Energy Research | www.frontiersin.org the increase of Pe bun , a region where turbulent thermal diffusion is stronger than molecular heat conduction begins to appear.

Turbulent Prandtl Number Distribution
The distribution of Pr t in the computational domain is plotted against different Pe bun numbers and displayed in Figure 9. As revealed in this figure, the Pr t is higher in the district close to the wall of the heated rod. Moreover, the overall turbulent Prandtl number in the computational domain decreases with the increase of Pe bun . In particular, the turbulent Prandtl number in the turbulent core region decreases significantly with the increase of Pe bun . The mean turbulent Prandtl number Pr tm of the computational domain is defined as follows: In order to investigate the influence of Pe bun on the Pr tm , Figure 10 plots the Pr tm against different Pe bun . From this figure, it can be concluded that the Pr tm tends to decrease as the Pe bun increases. However, the rate of decline also decreases with the increase of Pe bun . Furthermore, the Pr tm of the computational domain is higher than 1, suggesting that the analogy about Pr t 0.85 is not appropriate in such low-Pr fluids.

Assessment of Different Turbulence Models and Turbulent Heat Transfer Models
To analyze the effect of the turbulence model on the simulation of heat transfer, in this sub-section, various turbulence models of OpenFOAM are also employed in buoyant2eqnFoam, including standard k − ε (Launder and Spalding, 1972;1983), SST k − ω (Menter and Esch, 2001;Menter et al., 2003;Hellsten, 2012), and LaunderSharma k − ε (Launder and Sharma, 1974). The specific definitions of these turbulence models can be found in the literature (Launder and Spalding, 1972;Launder and Sharma, 1974;Launder and Spalding, 1983;Menter and Esch, 2001;Menter et al., 2003;Hellsten, 2012). A comparison of Nu sub1 calculated by different turbulence models with heat transfer experimental correlations is displayed in Figure 11. As demonstrated in this figure, the Nu sub1 calculated by four turbulence models is pretty close when Pe sub1 < 1000, mainly because the molecular heat conduction is dominant in this Peclet number range. It should be noted, however, that the deviations of Nu sub1 obtained by each turbulence model gradually increase as the Peclet number grows. As indicated in Figure 11, all Nu sub1 predicted by the standard k − ε turbulence model is located between Mareska and Mikityuk correlations. However, for Abe k − ε, LaunderSharma k − ε, and SST k − ω, the calculated Nu sub1 lies between the Graber and Subbotin correlations when Pe sub1 > 1000. Although the Nu sub1 results obtained by Abe k − ε, LaunderSharma k − ε, and SST k − ω are very similar, in general, Abe k − ε is the most conservative. In addition, the maximum deviation of the Nusselt number obtained by the Abe k − ε and the standard k − ε is close to 19%. It is worth mentioning that due to the significant deviation between the various Nusselt number experimental correlations, the quality of these turbulence models cannot be evaluated. Therefore, great care and caution should be exercised when selecting a turbulence model in simulation. More precise experimental and analytical studies are required in the future to identify the thermohydraulic characteristics of heavy liquid metals like LBE.
Moreover, the Nu sub1 calculated by the Manservisi k θ − ε θ model and the Pr t 0.85 model is reported in Figure 12. It is evident that compared with the experimental correlations plotted in Figure 12, the Nu sub1 obtained by the Pr t 0.85 heat transfer model is higher under almost all Peclet numbers. Oppositely, the Manservisi k θ − ε θ heat transfer model provides the more conservative results of Nu sub1 .

CONCLUSION
In the current study, the Abe k − ε turbulence model was compiled into the turbulence model library coming with OpenFOAM. A CFD solver buoyant2eqnFoam, which introduces the Manservisi k θ − ε θ turbulent heat transfer model, was developed. Subsequently, the Abe k − ε turbulence model with wall functions and Manservisi k θ − ε θ turbulent heat transfer model with zero-gradient boundary were employed to analyze the thermohydraulic characteristics of LBE inner flow in the bare 19-rod bundle. In addition, the influence of the turbulence model on the prediction of turbulent heat transfer was investigated by employing various turbulence models in the self-compiled solver buoyant2eqnFoam, including Abe k − ε, standard k − ε, SST k − ω, and LaunderSharma k − ε. Based on the aforementioned discussions, conclusions obtained from the present work can be summarized as follows: 1) The Nusselt numbers obtained by the self-compiled solver buoyant2eqnFoam are in good agreement with experimental correlations and Chierici simulation research, indicating the validity and reliability of the self-compiled solver. 2) In the bare 19-rod bundle with P/D 1.4, the flow resistance of the corner sub-channel is higher than that of other subchannels due to the smaller hydraulic diameter, leading to the appearance of higher temperature distribution and larger hot spot factor in this region. 3) Although the turbulent Prandtl number of LBE inner flow in the bare 19-rod bundle will decrease as the Peclet number increases, the overall turbulent Prandtl number is higher than 0.85, revealing that the Reynolds-analogy hypothesis about Pr t 0.85 is not appropriate for low-Pr number fluids like LBE. 4) The turbulence model has a considerable influence on the calculation of turbulent heat transfer of low-Pr number fluids in the high Peclet number range, suggesting that it should be prudent and rigorous when picking a turbulence model in the simulations. Moreover, compared with the k θ − ε θ turbulent heat transfer model, the Reynolds-analogy hypothesis about Pr t 0.85 may give the much higher Nusselt numbers in the simulation of low-Pr number fluids.
The applicability of the solver developed in the present study for the more complicated geometry like fuel assembly with grid spacer or wire-wrapped configurations requires further verification.