A Logarithmic Turbulent Heat Transfer Model in Applications with Liquid Metals for Pr = 0.01–0.025

: The study of turbulent heat transfer in liquid metal ﬂows has gained interest because of applications in several industrial ﬁelds. The common assumption of similarity between the dynamical and thermal turbulence, namely, the Reynolds analogy, has been proven to be invalid for these ﬂuids. Many methods have been proposed in order to overcome the difﬁculties encountered in a proper deﬁnition of the turbulent heat ﬂux, such as global or local correlations for the turbulent Prandtl number and four parameter turbulence models. In this work we assess a four parameter logarithmic turbulence model for liquid metals based on the Reynolds Averaged Navier-Stokes (RAN) approach. Several simulation results considering ﬂuids with Pr = 0.01 and Pr = 0.025 are reported in order to show the validity of this approach. The Kays turbulence model is also assessed and compared with integral heat transfer correlations for a wide range of Peclet numbers.


Introduction
The fluid-dynamics of low-Prandtl-number liquid metals is currently being studied due to an increasing interest in the applications of these fluids in engineering problems. Advanced nuclear reactor designs cooled with liquid metals have been proposed in the Generation IV Forum and some liquid metal fast nuclear reactors have been built and operated [1]. Liquid metal fluid dynamics is also studied in other engineering fields as they are considered possible coolants for new solar power systems [2]. These fluids show low dynamical viscosity with high thermal conductivity and therefore low Prandtl numbers. Two main important coolants for Generation IV nuclear reactors are sodium and lead-bismuth eutectic. The coolant Prandtl number depends on the operational temperature of the reactor; in particular, with increasing temperature the Prandtl number decreases. The Prandtl number of lead-bismuth eutectic is in the range of 0.03-0.01, while sodium has a lower Prandtl with values ranging from 0.01 to 0.005 [3]. Despite the increasing interest in reliable computational tools for these fluids for the study of turbulent heat transfer in complex industrial applications with liquid metals, they are still lacking in commercial codes.
The study of reliable and robust computational techniques to deal with this problem is still in the development phase [1, [4][5][6][7][8]. Due to the complexity of the modeling, different approaches have been proposed by different authors. Grötzbach in his work [1] attempted to model the temperature variance equation for turbulent natural convection; meanwhile, Cheng et al. in [4] proposed a modeling of the turbulence Prandtl number based on experimental data in circular tubes. Near-wall modeling and Kays correlation based on Large Eddy Simulations (LES) for forced convection at low Prandtl numbers is described by Duponcheel et al. [5] and Bricteux et al. in [6]. Shams et al. in [7] presented an overview of different algebraic models for prediction of heat transfer turbulent and different operational temperature ranges. Furthermore, we present the model in a logarithmic k-Ω-k θ -Ω θ formulation with the purpose of improving its numerical robustness. In two-equation turbulence models, when state variable limiters are not used, numerical instabilities usually appear due to negative values of the state turbulence variables, especially near boundaries where temperature and velocity field profiles show steep variations. A logarithmic formulation of the problem avoids these numerical undershooting or overshooting errors and improves the smoothness of dynamical and thermal field profiles at the walls for a thinner boundary refinement [23]. Since a four-equation turbulence model k-Ω-k θ -Ω θ could be computationally expensive, with respect to a SED model, we report also the results obtained by using Kays correlation.

Turbulence Model k-Ω-k θ -Ω θ
The Reynolds averaged Navier-Stokes equations for incompressible flows can be written as [8,22] ∇ · u = 0 , (1) ∂T ∂t where u is the average fluid velocity, P the pressure, T is the average temperature and Q represents the volume heat source or sink. In Equations (2) and (3) Reynolds stresses τ R and the turbulent heat flux q R are calculated using the mean velocity u, the mean temperature T and two new variables, the eddy viscosity ν t and the eddy heat diffusivity α t by setting [25,26] where u and T are the values of fluctuating velocity and temperature field, respectively. The hypothesis in Equation (4) is known as the Boussinesq hypothesis, wherein k represents the turbulent kinetic energy. The hypothesis in Equation (5) is known as diffusion gradient hypothesis. When α t is a constant we have simple eddy diffusivity model (SED), but it is possible to consider α t as a function of the velocity field, leading to the generalized diffusion gradient model. In this paper, in order to extend the range of the turbulence model validity with respect to the Prandtl number, we employ a four-parameter logarithmic turbulence model k-Ω-k θ -Ω θ (KLW model). It is obtained through a change of state variables in the k-ε-k θ -ε θ turbulence model [8,[16][17][18]. Let k be the turbulent kinetic energy, ε be its dissipation rate, k θ be the mean-squared temperature fluctuations value and ε θ be its dissipation rate defined as To compute the time scales of turbulence we propose to employ logarithmic values of k and k θ specific dissipation rate ω and ω θ , where C µ is a model constant.
The motivation to use logarithmic variables is well known [8,23]. During numerical computations the state values of turbulence equations, such that k θ , ε θ , k and ε, may take negative values near the boundary. Laplacian operator with negative coefficients looses its stability and convergence of the Navier-Stokes system may not be attained. This is one of the main problems with k-ε models and numerical limiters are very common [23]. Furthermore logarithmic variables take smoother variations at the wall where turbulence state variable profiles take larger changes [23]. By substituting the new definitions in Equation (7) in the k-ε formulation of the turbulence model, we obtain [8,[16][17][18] The terms P k and P kθ are the source terms for turbulent kinetic energy and temperature fluctuations. In our modeling we set A summary of the model's constant values is collected in Table 1, where f ε and C d2 are two functions needed to correct the behavior of the turbulence model near the wall [15][16][17][18]  We define the turbulent Reynolds number R t and the non-dimensional length R δ by using the variables k and ε. Let δ be the wall distance and (ν 3 /ε) 1/4 be the Kolmogorov length scale. Now we have (15) and, when logarithmic variables are used, R t and R δ become

Eddy Viscosity and Heat Diffusivity Model
In order to model the eddy viscosity ν t and heat diffusivity α t , we need to introduce different time scales. For the dynamical case, let τ lu be the dynamical turbulence time scale and the eddy viscosity ν t be defined as [8,[16][17][18]27] where The far wall characteristic dynamical time scale τ u is defined as k/ε = 1/C µ e Ω , f 1µ A 1µ and f 2µ A 2µ are the correction terms for the bulk and the near wall behavior, respectively. As the wall distance increases, f 2µ decreases, and it is significantly different from zero only in the near wall region. This function forces the correct near wall behavior of ν t when the wall distance δ tends to zero; namely, ν t has to be proportional to δ 3 .
For the thermal case, let τ lθ be the dynamical turbulence time scale and the eddy heat diffusivity α t be defined by [8,[16][17][18][19][20]27] where f 2aθ = f 1θ e −(R t /500) 2 , As in the dynamical case the time scale τ lθ is calculated as the sum of a near and a far wall term. B 1θ is the bulk term and depends only on the dynamical time scale τ u . It is the most important term in the region far from the wall, while B 2θ dominates in the near and medium regions. In the near and medium region, another characteristic thermal time scale should be introduced. Let τ θ be the thermal time scale defined by k θ /ε θ = 1/C µ e Ω θ . Let R be the thermal to dynamical time scale ratio τ θ /τ u = k θ ε/(kε θ ) = e Ω−Ω θ . B 2θ , which consists of two contributions, and is the sum of a near wall term proportional to the time scale τ u √ R = √ τ u τ θ and a term proportional to the mixed time scale τ m calculated as the harmonic mean between the time scales τ u and τ θ ; namely, τ m = τ u R/(C γ + R). This setting allows the model to reproduce the near wall behavior of the eddy thermal diffusivity with α t proportional to the cubic power of the distance from the wall δ in the presence of no thermal fluctuations on wall surfaces, and α t proportional to the square power of the distance from the wall δ, when thermal fluctuations are considered [8,16,17,27].
In this study we shall also consider the Kays model which has been used in many RANS simulations [5,13]. This model computes α t as a function of the eddy viscosity ν t . We consider ν t as in Equation (17) with the appropriate time scales and model functions already defined. α t is computed through the definition of a variable turbulent Prandtl number as where ν t /ν is the eddy viscosity ratio. This model has some advantages because it is very simple, not computationally expensive and allows one to compute the Pr t as a local function.
In the Numerical Results section we compare the results obtained with these two models with DNS data and an experimental correlation for liquid-metal flows.

Boundary Conditions for Turbulence Models
Two-equation turbulent models can be solved if appropriate boundary conditions are applied. When a near-wall approach with no wall functions is used, the boundary conditions can be computed by a near wall series expansion, as reported in [22]. In Table 2 the expansion is computed for the mean and fluctuating velocity, and for the mean and fluctuating temperature, by considering a small, plane surface area with x-z plane coordinates and y the wall distance. We remark that the velocity normal to the surface is assumed to vanish like its fluctuating components. Table 2. Expansion for the components of the velocity and temperature for plane coordinates x-z around wall points and y distance from the wall.

Mean Components
Fluctuating Components By using similar expansions we calculate the near wall behaviors of k, ε, k θ and ε θ , which are labeled by the lower-script w ("wall"). Then we combine these expansions and derive those for Ω and Ω θ . Dynamical turbulence variables can be expanded around a wall point as From Equation (24) one can see that the k w tends to zero as y vanishes. The constant ξ is not known but the k w boundary conditions can be written in a closure form by taking the derivative of k. In fact we can write ∂k The form in Equation (27) is a well known result that allows us to determine Dirichlet boundary conditions for fluctuations in terms of their derivatives. We see that the dissipation of turbulent kinetic energy, ε, in the near wall region does not depend on the wall distance. The near wall behavior is where ξ is a constant value which depends on the fluctuating velocity components. This value cannot be determined a priori and it should be calculated as ξ = 2k w /y 2 . On wall surfaces we thus impose a Dirichlet b.c. ε w = νξ with ξ varying iteration by iteration. This could be a reason for the KE model showing poor numerical stability. Ω w depends only on the known kinematic viscosity ν, wall distance y and model constant C µ so that the KLW model is not affected by this instability issue. If we take the Ω derivative in the y direction we obtain so we can also impose a Neumann b.c. on wall surfaces. The boundary conditions for the dynamical turbulent variables of the KLW model are summarized in Table 3.
The issue of boundary conditions on fluctuating variables is still an open question [8,9,16,17,27]. As it concerns the energy equation, in this paper we apply only boundary conditions with uniform heat flux on the wall surfaces. Boundary conditions of this kind have to be satisfied for the total variation which is the sum of an instantaneous and an averaged temperature. For this reason the derivative of the fluctuating temperature tends to vanish as we tend to the wall. In this paper we use the mixed boundary condition (MX) in which the temperature fluctuation and its derivative on the wall-so the terms d 0 and d 1 in the near wall expansion in Table 2-are set to zero.
With these assumptions on the fluctuating temperature field we obtain the following near wall behaviors We observe that, as for the case of dynamical turbulence, the temperature fluctuations k θ are proportional to y 2 while its dissipation ε θw tends to a constant value. The quantity ε θw is thus affected by the same issue of ε w , so we can only apply a Dirichlet b.c. with a value αd 2 2 that varies iteration by iteration. For the variable k θw we impose a Neumann b.c.
For Ω θw we can impose both forms of b.c., as explained for the dynamical turbulence. For thermal turbulence variables one can refer to Table 3. We finally observe that if we impose the exact Dirichlet b.c. on the variables Ω and Ω θ then the time scale ratio R is This implies that the value of R tends to Pr near the wall surface. In some DNS simulations this limit has been observed when MX b.c. are used [9].

Numerical Results
In this part of the paper we investigate fully developed turbulent flows in plane channels and cylindrical ducts for several values of the Reynolds number by using the KLW and Kays models. In particular, the simulations reported are obtained for a fluid with Pr = 0.01 in the plane channel and for fluids with Pr = 0.01 and Pr = 0.025 in the cylindrical pipes.
By reporting these results we aim at extending the validation of our turbulence model to fluids with different Prandtl numbers. In Table 4 we report the values of the physical properties employed for the simulations. The materials used as coolant liquid metals are lead, lead-bismuth eutectic, sodium and some other alloys. We study liquid metals with non-dimensional Pr numbers in the range 0.025-0.01. We use the values in Table 4 since not enough data-related information is reported in many reference papers. However, the model depends only on non-dimensional numbers which are independent of the specific property values. For each case we have specified the non-dimensional results.
We implement and solve the system in Equations (8)- (11) together with the Navier-Stokes Equations (1)-(3) in the code FEMuS [28]. The code uses Taylor-Hood finite elements for the coupled velocity-pressure Navier-Stokes equations, and quadratic finite elements for the turbulence equations that are solved uncoupled at the end of the Navier-Stokes step where the turbulence fields are updated. When dealing with new turbulence models it is very important to have available reference data for the model validations. These data can either be results of direct numerical simulations or results of physical experiments. In the first case the model can be validated after a quantitative comparison of the calculated fields with respect to DNS results, with a great level of detail, since fields can be compared on each point of the computational grid. In the case of experimental results, the model validation is done mainly through integral quantities; namely, pressure drops measured between different sections of the domain or heat transfer coefficients. DNS results are very useful, but are usually restricted to low Reynolds number, because for higher Re number values the required computational cost is still too high. Experimental results are instead available on a wider range of Re numbers. In the following sections DNS results are used to compared the obtained results for the two-dimensional channel case, since DNS data are available for three cases of Re number at Pr = 0.01, while for the study of the cylindrical duct we focus our attention on heat transfer coefficients values and validation against experimental correlations.

Two-Dimensional Channel Tests
Experimental data with forced flows in engineering applications cannot provide velocity and temperature profiles due to the measurement limits of the liquid-metal instrumentation. However, these profiles can be computed by direct numerical simulations which are available mainly in low Prandtl numbers, in plane geometries and in fully developed turbulence regime [9,10]. For the case Pr = 0.025, even if the KWL model is slightly modified, the results are very similar to those reported in [8] and therefore we focus our attention to Pr = 0.01.
Since DNS data in [10] are available for both the dynamical and thermal turbulence for a fluid with Pr = 0.01, we study three KLW model simulation cases for Re τ = 180 (A), 395 (B) and 590 (C). The corresponding bulk Reynolds numbers for these cases are respectively Re = 5760, 14,160 and 22,170. The physical domain consists of a channel of width L = 0.0605 m. On the wall surfaces a uniform heat flux q per unit area is applied with q = 3.6 × 10 5 W/m 2 .
We consider half physical domain because the fully developed turbulent flows present symmetry with respect to the middle channel axis. The boundary conditions we apply are reported in Table 5. We denote with x the direction normal to the wall surfaces and with y the axial direction. The inlet and outlet surfaces are represented with Γ in and Γ out respectively; Γ sym indicates the symmetry plane and Γ w the wall surface. Different mesh refinements have been used to study the turbulent flows for the various Re t cases, in order to obtain a value of non-dimensional wall distance y + smaller than 1 on the first mesh nodes away from the wall. Table 5. Two-dimensional channel tests. Boundary conditions.

Var.
Γ The results are set to non-dimensional values by using the kinematic viscosity ν; the friction velocity v τ = τ w /ρ, with τ w the wall shear stress and ρ the fluid density; and the friction temperature T τ = q/(v τ ρC p ), where C p is the fluid specific heat capacity.
In Figure 1 the non-dimensional axial velocity profiles v + are reported for different Re t together with the DNS results. We see that the linear behavior v + = y + is well reproduced in the viscous region, while the typical logarithmic behavior v + = 1/κ ln(y + ) + 5.2 is better captured by the case of Re t = 590, since for the other two cases the Re t is quite small and the deflection of the velocity profile, close to the channel center, occurs before the complete development of the logarithmic profile. The non-dimensional Reynolds stresses τ + R = u + v + and the non-dimensional total stresses τ + e f f = τ + R + νv ,y /u 2 τ are shown as a function of the non-dimensional wall distance y + in Figure 1b together with the DNS data. The comparison is good in the logarithmic region (y + > 30), while some discrepancies are present in the buffer region (5 < y + < 30) and in a part of the viscous region. A more quantitative comparison between the KLW results and the DNS ones is reported on the bottom plots of Figure 1, where the percentage difference is defined on the base of the effective value. For the non-dimensional velocity, the difference with respect to the DNS results is bounded in the range [−2%; 6%], with slightly higher values obtained in the buffer region. Even if the non-dimensional Reynolds stress component τ + r shows some underestimation with respect to the DNS value in the viscous region, the total stress τ + e f f = τ + + τ + R shows very good agreement with the reference values due to the the turbulent stress. In fact, the viscous stress is predominant in the viscous region, while the turbulent stress becomes significant in the buffer and bulk regions. For the total stress the difference with DNS data is in the range [−0.5%, 2%]. The non-dimensional mean temperature profilesT + are shown in Figure 2a for all the simulated cases. We see that the numerical results match the linear behaviorT + = Pr y + and the relative DNS data.
In Figure 2b we report the profiles of the non-dimensional turbulent heat flux q + R ·n = u + T + , wheren is the unit vector normal to the wall surface, and of the total non-dimensional heat flux q + e f f = (q + R + q + ) ·n. In the viscous region the non-dimensional heat flux profiles agree with the reference behavior, i.e., u + T + ∝ (y + ) 3 , which is obtained in the case of zero temperature fluctuations on the heated wall. In fact, from the near wall behaviors described in Table 2, it follows Having imposed null temperature fluctuations on the heated wall, the value of d 0 is equal to zero, so u + T + = b 2 d 1 y 3 + . . . , with b 2 and d 1 being unknown coefficients. As for the case of non-dimensional velocity and Reynolds stress component, a quantitative analysis is reported with the percentage difference plots shown on the bottom of Figure 2. For the non-dimensional temperature, the agreement is very satisfactory, with a diff value in the range [−0.25%; 0.5%], while for the total non-dimensional heat flux the difference with DNS is in the range [−0.5%, 2.5%]. In particular, in the viscous and buffer regions the computed turbulent heat flux overestimates the reference values. Since the turbulent heat flux is computed as u T = α t T ,x and the temperature profile perfectly matches the reference one, the overestimation of the heat flux is due to the modeled value of the eddy thermal diffusivity. At low Reynolds, in order to capture these components, an appropriate non-isotropic model may be needed. However this difference in the turbulent heat flux does not influence the temperature profile because the molecular heat flux has a great importance for these very low Prandtl and Reynolds numbers, and for very high y + (see [5]), as can be seen from the comparison of q + R and q + e f f .
klw, q + ef f 180 dns, q + ef f 180 10 0 10 1 10 2 10 3 10 −6 10 −4 klw, q + ef f 590 dns, q + ef f 590 10 0 10 1 10 2 10 3 Figure 2. Two-dimensional channel tests (Pr = 0.01). Non-dimensional mean temperatureT + profiles (a) and non-dimensional turbulent heat flux q + R and total heat flux q + e f f (b) as a function of the non-dimensional distance from the wall y + . The results obtained with KLW model are compared with DNS data [10]. On the bottom, quantitative comparisons with DNS data forT + and q + e f f .

Cylindrical Pipe Tests
With the reference to a duct geometry that consists of a cylinder with a diameter D = 0.0605 m, we simulate several cases of fully developed turbulent flows with both Pr = 0.01 and Pr = 0.025. A summary of these cases is reported in Table 6 together with the relative value of the turbulent Reynolds number Re τ and of the Reynolds number Re. On the wall surface, a uniform heat flux q per unit area is applied with value of 3.6 × 10 5 W/m 2 . As for the case of plane channel, we denote with Γ w the wall surface; with Γ in and Γ out the inlet and outlet sections; and with Γ sym the cylinder axis which is the axial symmetry axis of the problem. The boundary conditions are reported in Table 7, where r is the radial coordinate, u its component of the mean velocity field and v the axial one.
The physical properties are shown in Table 4 where we consider two fluids with different Prandtl numbers. We implement the modified KLW model and also the Kays correlation in the same numerical code. The same test cases are carried out for this model and the same comparison with the experimental correlations is given. In this way we intend to compare two different thermal turbulence models using the same numerical code and the same model for the dynamical turbulence, thereby allowing us to highlight the differences obtained in the simulated heat transfer performances, which are then due only to the different thermal turbulence models. In Figure 3a the non-dimensional temperatureT + profiles are reported for the simulated cases. We see that the simulations match the linear behaviorT + = Pr y + , in accordance with the near wall behavior described in Table 2. In Figure 3b the non-dimensional root-mean-squared temperature fluctuation k θ is shown as a function of non-dimensional wall distance y + for all the simulated cases and for both Pr = 0.01 and 0.025. We can see how the different molecular Prandtl numbers affect the temperature fluctuations. In fact, for any given Reynolds number, the simulations with Pr = 0.01 are characterized by a maximum value of T + rms which is smaller than the one obtained with Pr = 0.025 and the peak position is shifted nearer the pipe's middle axis. In Table 8 we report the mean values across the pipe transverse section of non-dimensional mean temperature and non-dimensional RMS temperature fluctuations together with their ratio T + rms,m /T + m . These parameters allow us to compare the mean temperature fluctuations modulus with the average difference between the fluid and the solid temperatures. For low velocities, i.e., with Re τ ≤ 1000, we see that the temperature fluctuations increase for increasing mean velocities and are higher for Pr = 0.025 than for Pr = 0.01. For the cases Re τ > 1000 we observe that the importance of the temperature fluctuations, relative to the mean temperature, decreases for increasing velocities and that it is higher for Pr = 0.01 than for Pr = 0.025.   Table 6.
For forced convection flows and high Reynolds numbers DNS simulations are not available and only experimental data can be compared with the results obtained with the KLW model. As a consequence, temperature profiles actually cannot be captured from liquid-metal experiments, and only integral measurements, i.e., Nusselt number values and heat transfer coefficients, are available. Nusselt numbers usually define the heat transfer coefficients between a wall surface and a fluid flow as where λ is the thermal conductivity, D h the hydraulic diameter and h the convective heat transfer coefficient. If a constant heat flux q is applied on wall surfaces, the convective heat transfer coefficient h can be written as where T w is the wall surface temperature and T b is the bulk temperature of the fluid defined as where the unit vector normal is denoted byn and A is the channel cross section. Experimental data summarize their results in the form of correlations for each type of geometry as a function of the Peclet number Pe = Pr · Re. Many experimental correlations are known for liquid metals, but unfortunately many of them do not agree. One can easily see that there is a disagreement between the different experiments and authors have attempted many explanations for this discrepancy; see for example, [4]. The general expressions of these correlations for the cylindrical pipe geometry take the following form: where a and n are constant assumed to be positive numbers. As said before, there are many different experimental correlations and each correlation comes from a different experiment that brings its own range of validity expressed in terms of Reynolds numbers. Some of these well known correlations are reported as The correlations in Equations (40)-(45) can be seen in Figure 4. Each correlation has a name related to its investigator, as labeled in Figure 4. In this Figure, correlation data from different experiments, such as lead, lead-bismuth, bismuth and sodium, are mixed with the aim to obtain a law valid for all liquid metals. The interested reader can see details for each experiment in [29][30][31][32][33][34][35][36]. For low values of Peclet numbers (Pe < 1000) the Kirillov curve is the correlation that fits better all the experimental points [4,34]. In [4] the authors proposed the following new correlation It is easy to see that this correlation matches the Kirillov correlation in the low Peclet region while it tends to the Stromquist one in the range of high Peclet numbers. In Figure 5 we compare the values of Nu obtained with the KLW and with Kays models for simulations in the range 80 < Pe < 10000 with Pr = 0.025 and 0.01 with Kirillov and Cheng correlations. The same data are shown in Table 9 as functions of different Pr and Re numbers.   The KLW model, for Pr = 0.025, matches close to the Kirillov correlation when Pe > 1000, while some difference can be seen for Pr = 0.01. The KLW model lays between Kirillov and the Cheng correlations, in the high Peclet region, while in the low Peclet region the model overestimates the values. From Figure 5 it can be seen that the Nusselt number values obtained from the KLW model, in both Pr = 0.01 and Pr = 0.025 cases, lay in the range [0, +5%] with respect to the Cheng correlation. As it regards the Kays model, we can see that for low Peclet and Prandtl numbers the results are greater than those obtained with KLW model and predicted by both correlations, while with increasing Peclet the mismatch becomes very large, with Nusselt number values lying in the range [+30%; +50%] with respect to Cheng correlation. When eddy viscosity ratio increases, namely, for high Peclet numbers, the average turbulent Prandtl, given by Equation (23), tends to 0.85, which is known to be inaccurate for predicting this type of flow.

Conclusions
In this work we have extended the investigation of heat turbulence modeling to a wider range of low-Prandtl numbers characteristic of liquid-metal flows. The investigation, which presents a two-equation logarithmic k-Ω-k θ -Ω θ model, includes forced flow studies with both low Peclet numbers, where only DNS simulations are available, and high Peclet numbers where only experimental results are available. Two approaches have been compared, a Kays model and a logarithmic k-Ω-k θ -Ω θ turbulence model, to be used in engineering applications with low-high Peclet numbers. Since only DNS results are available, for low-Peclet-number, in plane and channel geometries, we have compared them for sodium and lead's characteristic Prandtl numbers. The comparison between the numerical results obtained with the four-parameter turbulence model with DNS data and experimental correlations can be considered satisfactory for a wide range of Peclet numbers and for both the Prandtl numbers studied. On the other hand, the Kays model does not seem to be able to adapt to the study of all these different cases, especially for high Peclet numbers, and a two-equation model for the heat turbulent diffusivity seems necessary.