DHFM is a full second-moment differential model for the transport of Reynolds heat fluxes. Compared with the constant P r 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 P r 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 elliptic blending model. The new model was utilized to study the fully developed square duct flow, rotating, and nonrotating channel flow. The numerical results show good agreement with the LES and DNS results.

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- P r fluids were implemented by Shams et al. (2014), Shams (2018), Shams et al. (2018), and De Santis et al. (2018) 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.

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- P r 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 (2014b); Manservisi and Menghini (2015), Cerroni et al. (2015), Cervone et al. (2020), Chierici et al. (2019), and Da Via et al. (2016); Da Vià and Manservisi (2019); and Da Vià et al. (2020) proposed a two-equation k θ − ε θ model suitable for the LBE turbulent heat transfer simulation, and improved new numerical near-wall boundary conditions for turbulence variables. The numerical results obtained by the literature (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- P r 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 two-equation 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 two-equation 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 open-source 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 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 two-equation k θ − ε θ model was presented.

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. S u b 1 , S u b 2 , and S u b 3 are defined as inner sub-channels, while S u b 4 and S u b 5 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 , b u n and D h , s u b 1 are the hydraulic diameter of the bundle and hydraulic diameter of S u b 1 , respectively. The length of the whole computational domain is set to 15 D h , b u n to eliminate the effect of employing periodic inlet boundary conditions. The flow parameters of LBE are reported in Table 2, where P e b u n = P r U b u n D h , b u n ρ / μ is the Peclet number of the bundle.

For forced convection, the incompressible RANS equations with constant physical properties and no gravity are considered ∂ u i ∂ x i = 0 (1) ∂ u i ∂ t + u j ∂ u i ∂ x j = − 1 ρ ∂ P ∂ x i + ∂ ∂ x j ( v ∂ u i ∂ x j ) − ∂ ∂ x j u i ′ u j ′ ¯ (2) 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 eddy-viscosity model can be adopted as follows: u ′ i u ′ j ¯ = − v t ( ∂ u i ∂ x j + ∂ u j ∂ x i ) + 2 k 3 δ i j (3) 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): ∂ h ∂ t + u j ∂ h ∂ x j + ∂ K ∂ t + u j ∂ K ∂ x j = 1 ρ ∂ P ∂ t + ∂ ∂ x j ( ( α + α t ) ∂ h ∂ x j ) (4) 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).

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 t = C μ f μ k 2 ε (5) where C μ is a constant. f μ is the model function, defined as follows: f μ = ( 1 − exp ( − y ∗ 14 ) ) 2 ( 1 + 5 R t 3 / 4 exp ( − ( R t 200 ) 2 ) ) (6) y ∗ = u ε δ ν (7) 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: ∂ k ∂ t + u j ∂ k ∂ x j = ∂ ∂ x j ( ( v + v t σ k ) ∂ k ∂ x j ) + P k − ε (8) P k = − u i ′ u j ′ ¯ ∂ u i ∂ x j = v t ( ∂ u i ∂ x j + ∂ u j ∂ x i ) ∂ u i ∂ x j − 2 k δ i j 3 ∂ u i ∂ x j (9) ∂ ε ∂ t + u j ∂ ε ∂ x j = ∂ ∂ x j ( ( v + v t σ ε ) ∂ ε ∂ x j ) + C 1 ε ε k P k − C 2 ε f ε ε 2 k (10) f ε = ( 1 − exp ( − y ∗ 3.1 ) ) 2 ( 1 − 0.3 ⁡ exp ( − R t 2 6.5 ) ) (11)

The model constants utilized in the Abe k − ε turbulence model are reported in Table 3.

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 P r 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: α t = C θ k τ l θ (12) where C θ is the constant empirical coefficient and τ l θ is the local thermal characteristic time, modeled as follows: τ l θ = f 1 θ B 1 θ + f 2 θ B 2 θ (13) with the appropriate functions set as follows: f 1 θ = ( 1 − exp ( − R δ 19 P r ) ) ( 1 − exp ( − R δ 14 ) ) (14) B 1 θ = 0.9 τ u = 0.9 k ε (15) f 2 θ B 2 θ = τ u ( f 2 a θ 2 R R + C γ + f 2 b θ 2 R P r 1.3 P r R t 3 / 4 ) (16) f 2 a θ = f 1 θ ⁡ exp ( − ( R t 500 ) 2 ) (17) f 2 b θ = f 1 θ ⁡ exp ( − ( R δ 200 ) 2 ) (18) 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: ∂ k θ ∂ t + u j ∂ k θ ∂ x j = ∂ ∂ x j ( ( α + α t σ k θ ) ∂ k θ ∂ x j ) + P θ − ε θ (19) P θ = α t ( ∂ T ∂ x j ) ( ∂ T ∂ x j ) (20) ∂ ε θ ∂ t + u j ∂ ε θ ∂ x j = ∂ ∂ x j ( ( α + α t σ ε θ ) ∂ ε θ ∂ x j ) + ε θ k θ ( C P 1 P θ − C d 1 ε θ ) + ε θ k ( C P 2 P k − C d 2 ε ) (21) C d 2 = ( 1.9 ( 1 − 0.3 ⁡ exp ( − ( R t 6.5 ) 2 ) ) − 1 ) ( 1 − exp ( − R δ 5.7 ) ) 2 (22)

The constant empirical coefficients used in the Manservisi k θ − ε θ turbulent heat transfer are reported in Table 4.

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 19-rod bundle in a wide range of P e b u n . 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: Θ = ( T − T b , b u n ) λ D h , b u n ⋅ q w (24) where T b , b u n 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 P e b u n = 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.

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 S u b 1 is picked for comparison. The N u s u b 1 is calculated as follows: N u s u b 1 = q w D h , s u b 1 λ ( T w , s u b 1 − T b , s u b 1 ) (25) where T w , s u b 1 and T b , s u b 1 are the mean wall temperature and mean coolant temperature of S u b 1 , respectively. Correspondingly, P e s u b 1 = P r U s u b 1 D h , s u b 1 ρ / μ is the Peclet number of the inner sub-channel S u b 1 . 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 N u s u b 1 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.

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 N u s u b 1 calculated by different turbulence models with heat transfer experimental correlations is displayed in Figure 11. As demonstrated in this figure, the N u s u b 1 calculated by four turbulence models is pretty close when P e s u b 1 < 1000 , mainly because the molecular heat conduction is dominant in this Peclet number range. It should be noted, however, that the deviations of N u s u b 1 obtained by each turbulence model gradually increase as the Peclet number grows. As indicated in Figure 11, all N u s u b 1 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 N u s u b 1 lies between the Graber and Subbotin correlations when P e s u b 1 > 1000 . Although the N u s u b 1 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 N u s u b 1 calculated by the Manservisi k θ − ε θ model and the P r t = 0.85 model is reported in Figure 12. It is evident that compared with the experimental correlations plotted in Figure 12, the N u s u b 1 obtained by the P r 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 N u s u b 1 .