International Journal of Heat and Mass Transfer 127 (2018) 345–358 Contents lists available at ScienceDirect International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt The influence of low Prandtl numbers on the turbulent mixed convection in an horizontal channel flow: DNS and assessment of RANS turbulence models Dante De Santis a, Andrea De Santis a, Afaque Shams a,⇑, Tomasz Kwiatkowski b a b Nuclear Research and Consultancy Group NRG, Westerduinweg 3, 1755 LE Petten, The Netherlands National Center for Nuclear Research, Andrzeja Sołtana 7, 05-400 Otwock-Świerk, Poland a r t i c l e i n f o Article history: Received 15 March 2018 Received in revised form 29 July 2018 Accepted 30 July 2018 Keywords: Turbulent heat transfer DNS Mixed convection Algebraic heat flux model Low Prandtl number fluids a b s t r a c t In this paper, the turbulent mixed convection in a channel flow with differentially heated walls is considered for a fixed Richardson number (Rib ¼ 0:5) and three different Prandtl numbers (Pr ¼ 1; 0:1 and 0:01). Numerical simulations are conducted assuming constant fluid properties and the effect of the buoyancy is taken into account by means of the Boussinesq approximation. Direct Numerical Simulations (DNS) are performed first and the effect of the buoyancy on the first and second order statistics of the fluid and thermal fields is highlighted. Furthermore, it is found that in mixed convection the Prandtl number has a much larger effect on the results than in the case of forced convection. The obtained DNS results are then used as a validation database for two different RANS turbulent heat flux models, i.e. the classical Reynolds analogy and a recently proposed three-equation algebraic heat flux models called AHFM-NRG+. It is observed that, as expected, the Reynolds analogy fails to predict the thermal field even for unitary Prandtl numbers fluids. On the other hand, it is shown that the AHFM-NRG+ is in a reasonable agreement with the reference DNS over the entire range of Prandtl numbers considered in the study. Ó 2018 Elsevier Ltd. All rights reserved. 1. Introduction The mixed convection is a process in which momentum and heat are transferred within the fluid flow due to the combined effects of shear and buoyancy. This phenomenon represents an important mechanism of heat transfer which can be found in several applications such as heat exchangers, nuclear reactors and cooling systems for electronic components. Practically speaking, the mixed convection regime is bound by the two extreme cases: the forced convection and the natural convection regimes. Nevertheless, the results in mixed convection regime are far from being just the combination of those in the two extreme cases. With respect to the forced convection regime, the pressuredriven turbulent flow in a plane channel (i.e., the Poiseuille flow) has been used extensively as a representative configuration. Starting from the seminal work of Kim et al. , several DNS of planar channel flows in forced convection have been performed. When thermal effects are considered, the temperature is simply assumed to be a passive scalar transported within the fluid flow and the parameters that define the problem are the Reynolds and the ⇑ Corresponding author. E-mail address: email@example.com (A. Shams). https://doi.org/10.1016/j.ijheatmasstransfer.2018.07.150 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved. Prandtl numbers. In addition, simulations with imposed thermal boundary conditions [2–5] as well as simulation with conjugate heat transfers [6–8] were considered. On the other hand, the Rayleigh-Bénard convection (RBC) has been widely considered as the prototypical case for the natural convection regime. It corresponds to a free convection between two differentially heated walls and it has been extensively studied both numerically and experimentally [9–11]. In this case, the controlling parameters are represented by the Grashof number (Gr) and the Prandtl number. Both experimental and numerical investigations of the mixed convection regime have been less frequent so far; this is due to: (i) the challenging nature of the mixed convection regime, which includes features typical of both forced and natural convection, (ii) the lack of an single representative prototypical case for this regime. With respect to the latter point, the Poiseuille-RayleighBénard (PRB) flow is an attractive option since it represents a combination of the main features of the prototypical cases of forced and natural convection. In fact, the PRB is a pressure-driven Poiseuille flow where buoyancy effects are induced through an imposed temperature difference between the bottom and the top walls of the channel. The most interesting PRB configuration is the one in which the bottom wall of the channel is kept at a higher 346 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 temperature than the top wall since, due to the effect of the buoyancy forces, it results in an unstable thermal stratification. This PRB configuration has been the subject of early experimental studies [12–14] carried out at low and moderate Richardson numbers (Rib ) and at unitary Prandtl number. It was observed that, compared to the forced convection regime, the buoyancy significantly alters the near-wall turbulence mechanism, leading to the onset of large scale convective thermal plumes. With respect to numerical investigations, an early simulation of the RBC with superimposed shear was performed in  in order to assess the combined effects of the natural convection and the mean shear on the flow. It was observed that the unstable stratification condition enhances the diffusion of the turbulent kinetic energy between the near-wall and the bulk regions. Subsequently, DNS of the PRB flow at friction Reynolds number Res ¼ 150 and at low and moderate Richardson numbers have been carried out by Iida and Kasagi : the statistics for both the fluid and the thermal fields were studied and it was observed that the buoyancy effects have a substantial impact on the near-wall turbulence, resulting in an increase in both the skin friction and the wall heat transfer. In the DNS study of Zonta and Soldati , the PRB flow was investigated at different friction Reynolds numbers (Res from 110 to 180) and different friction Richardson numbers (Ris from 926 to 346) and for the Prandtl number Pr ¼ 3. It was found that, with respect to the mixed convection case, in the PRB flow the wall-normal transport of momentum and heat is enhanced due to the presence of large scale convective structures similar to those observed in RBC. In the same work, the effects of temperaturedependent fluid properties on the PRB flow were also investigated. In the work of Sid et al. , a DNS of the PRB flow was performed at two different friction Reynolds numbers (Res ¼ 180 and 395) and two different Richardson numbers (Rib ¼ 0:1 and 1) and at unitary Prandtl number. The study confirmed that both the skin friction and the wall heat transfer are enhanced under unstable stratification for sufficiently large Richardson numbers. The authors also showed the breakdown of the constant turbulent Prandtl number (Prt ) assumption under the considered conditions. In a very recent work, Pirozzoli et al.  performed DNS of several mixed convection cases in an horizontal channel flow for a wide range of Reynolds and Rayleigh numbers and they showed that the mixed convection regime is characterized by the presence of longitudinal rollers within the channel height; these rollers are characterized by a large spanwise aspect-ratio. In addition to the previous observations, it should be pointed out that the RANS approach still represents the workhorse for the modelling of the turbulent heat transfer in complex flow configurations of industrial relevance, mainly due to its reduced computational costs. In this respect, a large number of models have been proposed and assessed for the closure of the turbulent momentum flux, whilst relatively little attention has been paid to the modelling of the turbulent heat flux (THF) term . In particular, in most commercial CFD codes the so-called Reynolds analogy is often the only option available for the latter purpose. This approach assumes similarity in the turbulent transport of momentum and energy, and is extremely popular due to its robustness and simplicity. Nevertheless, it presents well-known limitations especially in applications involving low-Prandtl fluids and/or non-negligible buoyancy effects. This has led to a growing interest towards the development of more advanced THF closures which would allow to overcome, at least to some extent, the shortcomings of the Reynolds analogy [21,22]. On the other hand, it is widely recognized that one of the main hampering factors for the development of such models is represented by the lack of reliable reference data [21,23,24]. In the present work, several DNS of the PRB case at moderate Reynolds number and at unitary and lower Prandtl numbers are considered. The focus of the work is on assessing the impact of the Prandtl number on the fluid and thermal fields under the effects of the buoyancy forces. To this purpose, DNS calculations at fixed bulk Reynolds number Reb ¼ 5639, corresponding to Res ¼ 180 in forced convection, and Richardson number Rib ¼ 0:5 are performed for different values of the Prandtl number: Pr ¼ 1; 0:1, and 0:01. The choice of low values of the Prandtl number is particularly interesting for nuclear applications and has not been considered yet in the existing literature. Furthermore, the DNS database generated in this work can be regarded as a valuable tool to develop, calibrate and assess turbulence models in the RANS framework. Consequently, the present DNS results are employed to assess the performance of two different RANS closures for the THF term. In particular, the classical Reynolds analogy and a recently proposed algebraic closure, the so-called AHFM-NRG+ [23,25], are considered. The former closure has been chosen, despite its well-known limitations, since it still represents the most common approach employed in the RANS framework. On the other hand, the latter model has been considered since it has been developed with a focus on low-Prandtl fluids and buoyant flows. The remaining of the paper is organised as follows: Section 2 presents a description of the considered cases. The DNS methodology and the related results are discussed in Section 3. Then, the set-up of the RANS simulation and the results obtained for the different cases are described in Section 4. Finally, conclusive remarks are given in Section 5. 2. Test case description In this work, a turbulent planar channel flow with differentially heated walls and constant fluid properties is considered. The streamwise direction is parallel to the x-axis, the wall-normal direction is parallel to the y-axis and the spanwise direction is parallel to the z-axis, as shown in Fig. 1. The fluid flow is driven by an imposed mass flow rate and the bulk Reynolds number of the channel flow is Reb ¼ U b Ly =m ¼ 5639, where U b is the bulk velocity and m is the kinematic viscosity of the fluid. For a forced convection case, this corresponds to a friction Reynolds pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ number Res ¼ ut d=m ¼ 180, where us ¼ sw =q, with sw the wall shear stress. The gravity (g) acts downward along the y-axis. The bottom wall is kept at the uniform high temperature (T h ¼ 1) and the top wall is kept at the uniform low temperature (T c ¼ 0). As a consequence, the wall-normal temperature difference (DT ¼ T h T c ) between the bottom and the top walls is responsible for the buoyancy effects in the flow field. For the thermal field, a fixed Richardson number Rib ¼ gbDTd=U 2b ¼ 0:5 is imposed, with b the thermal expansion coefficient. Three different Prandtl numbers (Pr ¼ mqcp =k) are considered, namely Pr ¼ 1; 0:1 and 0:01 with, q the fluid density, cp the specific heat and k the thermal conductivity. Fig. 1. Sketch of the three-dimensional computational domain for the threedimensional mixed convection simulations. 347 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 On the walls, the no-slip boundary condition is imposed for the velocity and the iso-thermal boundary condition is imposed for the temperature. The flow is considered periodic in the streamwise and spanwise directions, hence periodic boundary conditions for the velocity and the temperature are applied at the boundary faces normal to the periodic directions. walls along the wall-normal direction. The number of elements for the three DNS is reported in Table 1. For each element, a grid refinement with a polynomial degree 9 in each direction is used to achieve the desired accuracy. In terms of friction units, the average mesh resolution is Dxþ ¼ 9:8 and Dzþ ¼ 4:2 in the streamwise and spanwise directions, respectively. The wall-normal mesh resolution near the wall and in the bulk region is Dyþ wall ¼ 0:34 and Dyþbulk ¼ 3:6, respectively. 3. DNS approach 3.2. Results and discussions 3.1. Governing equations and numerical method The obtained DNS results for the turbulent flow field are presented first, followed by the thermal field analysis. For the sake of completeness, the results of the mixed convection cases are compared with those for a forced convection case. The DNS results of the forced convection are taken from the DNS database of Kasagi and Iida , since this is one of the few public databases for the forced convection case with differentially heated walls. In this database, Res ¼ 150 and Pr ¼ 0:71 which are different from the parameters adopted in the present DNS, hence the results of the DNS database of Kasagi and Iida are considered only for the purpose of a qualitative comparison. In all the simulations, in order to achieve statistically steady solutions, after the flow field has sufficiently developed the statistics are collected every time step for approximately 50 flow through time. All the profiles reported are obtained by averaging also in space along the two homogeneous directions (i.e., streamwise and spanwise); the time and space average operator is indicated as h i. For all the results presented in this section, the quantities are made non-dimensional with the wall friction quantities and are indicated with the superscript þ . Furthermore, since results are symmetric around the middle plane of the channel, only the profiles along half the channel height from the heated wall are presented. 2.1. Boundary conditions The governing equations are the incompressible Navier-Stokes and energy equations coupled together through the Boussinesq hypothesis, namely $ u ¼ 0; @u 1 2 r u þ Rib Tey ; þ ðu $Þu ¼ $p þ @t Reb @T 1 r2 T; þ u $T ¼ @t Pr Reb ð1Þ where u is the velocity vector with components u; v and w, in the streamwise, wall-normal and spanwise directions, respectively, p is the pressure, T is temperature and ey is the unit vector parallel to the y-axis. The simulations are performed using the spectral element code NEK5000  which solves the governing equations on hexahedral elements by means of local approximations based on high-order orthogonal polynomial basis and using the Gauss-LobattoLegendre node distribution. The velocity and the pressure spaces are represented by the same polynomial degree. Furthermore, over-integration and filtering stabilization are adopted. The equations are advanced in time using a third order mixed Backward Difference/Extrapolation (BDF3/EXT3) scheme and using a constant CFL number equal to 0:2. The dimensions of the computational domain have been carefully chosen in order to minimize their effects on the numerical results. Initial simulations were performed using a standard domain size commonly used for forced convection DNS at unitary Prandtl number, namely Lx Ly Lz ¼ 4pd 2d 4=3pd. The a posteriori analysis of the obtained results showed that a larger domain was needed as a consequence of the large thermal structures which are present in the case of mixed convection, especially at lower Prandtl number. The two-point spatial correlations computed for a coarse DNS (not reported here) on a large domain Lx Ly Lz ¼ 15pd 2d 15pd have been used to estimate the dimensions of the computational domain to be employed for the final simulations. On the basis of this study, different optimal domain sizes for the three simulations have been identified and are reported in Table 1. These dimensions are then verified with the a posteriori analysis using the two-point spatial correlations, see Appendix A. The same spatial resolution is enforced in the three DNS, hence a different number of elements is used in each case. A uniform spatial distribution of the elements is used in the streamwise and spanwise directions, while the elements are clustered near the Table 1 Domain sizes and number of elements for the three DNS cases. Case Pr ¼ 1 Pr ¼ 0:1 Pr ¼ 0:01 Dimensions (Lx Ly Lz ) Number of elements (N x N y N z ) 5pd 2d 5pd 7:5pd 2d 7:5pd 10pd 2d 10pd 58 20 25 90 24 40 120 24 56 3.2.1. Flow field analysis As mentioned previously, numerical simulations are performed with a fixed mass flow rate corresponding to Res ¼ 180 in the forced convection case. As a consequence, in the mixed convection regime the wall shear stress will be different from that observed in the forced convection and, in general, it will be a function of the Prandtl number. In Table 2, the computed Res for the different values of Pr is reported. It is worth noticing that, for the selected configuration, in forced convection Res ¼ 180, while here it can be seen that the effect of the buoyancy is to increase the Res with respect to the forced convection case, and this effect is more evident for lower values of Pr. The phenomenon is a consequence of the fact that, under unstable thermal stratification, the large scales rolls induced by the thermal plumes tend to confine low-speed streaks near the wall region which are intensified by buoyancydriven instabilities. Hence, as observed also in other works [17,18], at a sufficiently high Richardson number the buoyancy effects in the PRB case enhance the wall skin friction with respect to the forced convection case. In Fig. 2, the mean velocity profiles along the wall-normal direction for the three PRB cases are reported together with the profile in forced convection. It can be noted that, with respect to the forced convection, in the mixed convection cases the velocity has Table 2 Friction Reynolds number (Res ) computed in the DNS of the PRB channel flow for the different Prandtl numbers (Pr) considered. The nominal Reynolds number in the case of forced convection is Res ¼ 180. Pr ¼ 1 Pr ¼ 0:1 Pr ¼ 0:01 Res ¼ 184:5 Res ¼ 192:4 Res ¼ 203:2 348 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Fig. 2. Mean velocity profiles along the height of the channel calculated with the DNS for the mixed convection case with different Pr and in the case of forced convection. a rather flatter profile in the bulk region as a consequence of the higher turbulence levels. Furthermore, this effect is accentuated at the low Prandtl numbers. Indeed, as indicated by the increment of the Res with the Pr in mixed convection, the main effect of the buoyancy is that of increasing the level of turbulence together with the associated wall-normal transport of momentum. The profiles along the wall-normal direction of the channel of the root mean square (RMS) of the velocity components fluctuation (hu0 irms ; hv 0 irms and hw0 irms ) and of the velocity covariance (huv i) are reported in Fig. 3 for the different Prandtl numbers together with the profiles obtained in the case of the forced convection. It can be observed that in the mixed convection regime the values of the velocity fluctuation are significantly higher than those in the forced convection regime, except for the case at unitary Prandtl number, and the difference increases as the Prandtl number decreases. Furthermore, compared to the forced convection case, the streamwise velocity fluctuation (hu0 irms ) has a flat profile in the bulk region while the profile of hu0 irms in the forced convection regime decreases in the bulk region. The profile of hv 0 irms in the mixed convection regime is significantly different from that in the forced convection and tends towards the profile observed in natural convection; more specifically, it can be noticed how the Fig. 3. Profiles of the RMS of the velocity components and of the covariance huv i along the height of the channel for the mixed convection case with different Pr and for the case of forced convection. 349 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 fluctuation of the wall-normal component of the velocity increases monotonously along the height of the channel, hence the maximum is located in the bulk region due to the enhanced turbulent level induced by the buoyancy effects. The spanwise component of the RMS in mixed convection is higher than that observed in forced convection as a consequence of the ejection of plumes in the near-wall region that enhances the turbulent fluctuation in the spanwise direction. Also the velocity covariance huv i is significantly higher with respect to the forced convection case, and higher velocity covariance values are observed for lower Prandtl numbers. To further highlight the impact on buoyancy on the structure of the turbulence, the anisotropy invariant maps for the mixed convection case at the three Prandtl numbers are reported in Fig. 4, together with the forced convection case (Res ¼ 150). The anisotropy invariant map is analyzed along the wall-normal direction of the channel by plotting the second and third invariants of the deviatoric part of the Reynolds stress tensor bij , defined respectively as II ¼ bij bji and III ¼ bij bjk bki , where bij ¼ hui uj i 1 dij ; 2k 3 i; j ¼ 1; 2; 3 with k the turbulent kinetic energy and dij Kronecker delta. In the case of forced convection, the turbulence has a two-component behavior in the near-wall region while, in the logarithmic region the Reynolds stress tensor has an axi-symmetric configuration with one large eigenvalue, and finally almost an isotropic condition is reached in the bulk region of the channel. In the cases in mixed convection regime, a two-component axi-symmetric state is still observed near the wall as a consequence of the fact that here the wall-normal component of the velocity fluctuation is much lower than the other two components, but the level of anisotropy is lower than that observed in forced convection. In the case of Pr ¼ 1, a fast return to the isotropic state is observed, while in the cases at lower Prandtl numbers, the behavior of the anisotropy invariant map is significantly different from that observed at unitary Prandtl. More specifically, outside the viscous sub-layer the Reynolds stress tensor initially tends toward a disc-like configuration and in the logarithmic regions the Reynolds stress tensor has a rod-like configuration similar to that observed in the case of forced convection. 3.2.2. Thermal field analysis In Fig. 5a and b, the temperature profile along the height of the channel is reported for the different Prandtl numbers respectively in bulk units, y ¼ y=Ly ; T ¼ T=DT, and in friction units hþ ¼ ðT T w Þ=T s , with T w the wall temperature, T s the friction temperature defined as T s ¼ qw =ðqcp us Þ, with qw the wall heat flux. From Fig. 5a, it is clear that the main effect of increasing the Prandtl number is that of increasing the wall-normal temperature gradient at the wall and to decrease it in the bulk region. Furthermore, it can be observed that the temperature profile exhibits an inversion for Pr ¼ 1 while at the lower Pr the sign of the temperature gradient does not change through the height of the channel. This type of inversion of the temperature has been already observed in the case of Rayleigh-Bénard convection, in particular in  it was suggested that this behavior might be caused by an excessive amount of cold plumes brought by the flow pattern. In Fig. 5b, the temperature profiles are reported in friction units for the different cases of mixed convection and, as a reference, the temperature profile computed in the case of forced convection of Kasagi and Iida  is reported as well. It can be observed that the effect of the buoyancy is significant both on the values and the profiles of hþ . In particular, the temperature in the bulk region has a flat profile in the mixed convection regime while in the 0.6 0.6 0.5 0.5 y 0.3 0.05 0.1 0.15 0.2 0.5 0.1 0 0 1 0.2 0.5 0.1 * y 0.3 1 0.2 0 −0.05 0.4 * II II 0.4 0 −0.05 0.25 0 0 0.05 III 0.1 0.15 0.2 0.25 III 0.6 0.5 0.4 II y* 0.3 1 0.2 0.5 0.1 0 −0.05 0 0 0.05 0.1 0.15 0.2 0.25 III Fig. 4. Anisotropy invariant map for the mixed convection case at the three different Prandtl numbers (symbols) as function of the wall-normal distance: Pr ¼ 1 (a), Pr ¼ 0:1 (b) and Pr ¼ 0:01 (c). In the figure the profile for the forced convection case at Res ¼ 150 is reported for reference in dashed line. 350 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Fig. 5. Mean temperature profiles along the height of the channel for the mixed convection case with different Pr. The profiles computed in the case of forced convection are reported as well for sake of comparison. Profiles in bulk units (a), profiles in friction units (b). forced convection it keeps increasing along the height of the channel. Furthermore, due the reduction of the wall heat flux with the decrease of the Prandtl number, the non-dimensional temperature decreases significantly with the Pr. The flattening of the temperature profile in the bulk region for the mixed convection case is due to an enhanced turbulent mixing with respect to the forced convection case. As a consequence, the temperature fluctuation in the core region of the channel is lower in the mixed convection case (see Fig. 6) than that observed in forced convection. The profiles of the streamwise and wall-normal turbulent heat flux are reported in Fig. 7 for the different Prandtl numbers. Compared to the case of the forced convection, it can be observed that a lower value of the streamwise turbulent heat flux is observed in the case of the mixed convection. In particular, as the Prandtl decreases the streamwise component of the turbulent heat flux is significantly reduced in the near wall region. The wall-normal component of the turbulent heat flux shows a similar trend with the Prandtl number, and its values are lower near the walls than that in forced convection except for the case at Pr ¼ 1. One interesting quantity that can be evaluated from the DNS simulation is the turbulent Prandtl number which is defined as the ratio of the momentum eddy diffusivity (mt ) to the thermal eddy diffusivity (jt ), namely Prt ¼ mt hu0 v 0 i @hTi=@y ¼ : jt hv 0 T 0 i @hui=@y ð2Þ In the classical Reynolds analogy, the Prt is assumed to be constant with values in the range 0:85—1:0. It is well known that this represents a substantial simplification which breaks down in several cases, such as in buoyant flows and when dealing with lowPrandtl fluids [21,28]. Due to the temperature inversion observed for Pr ¼ 1 (see Fig. 5), the wall-normal derivative of the temperature changes in sign along the height of the channel resulting in a negative value of the Pr t in the bulk region, as observed also in other works . Nevertheless, since in the context of the Reynolds analogy Prt is assumed to be a positive number, the absolute value of this quantity is reported in Fig. 8. This figure shows the profiles of the calculated turbulent Prandtl number along the height of the channel for the different values of Pr. It can be observed that in the mixed convection regime, Prt deviates significantly from the constant value that is commonly used in the Reynolds analogy approach (i.e., Prt 0:9) and it shows a noticeable variation along the height of the channel. In the cases with Pr ¼ 0:1 and Pr ¼ 0:01 the profiles of Pr t are qualitatively comparable within the bulk region, although a higher near wall value is observed in the case Pr ¼ 0:01. 4. RANS approach The DNS database described in the previous section has been employed as a reference for the assessment of two different RANS turbulence models. Details about the RANS numerical setup and the results obtained in the simulations are discussed below. 4.1. Numerical setup Fig. 6. RMS of the temperature fluctuation in mixed convection regimes for different Prandtl numbers. For reference, the profile in forced convection is reported as well. All the RANS simulations have been performed using the commercial code STAR-CCM . The computational domain is twodimensional, with Ly ¼ 2d and Lx ¼ 4pd the wall-normal and streamwise lengths of the domain, respectively. It consists of 15,000 structured quadrilateral cells (N x N y ¼ 100 150). The D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 351 Fig. 7. Streamwise (a) and wall-normal (b) turbulent heat flux profiles along the height of the channel for the mixed convection case with different Pr and for the forced convection case as well. and the mean rate of strain. For the closure of the THF term, two different approaches have been considered. The first model relies on a Simple Gradient Diffusion Hypothesis (SGDH), in which the THF is expressed as u0 T 0 ¼ jt $T ð3Þ where jt is the unknown thermal eddy diffusivity. This quantity is evaluated by invoking the Reynolds analogy, which results in u0 T 0 ¼ mt Prt $T ð4Þ where mt is the turbulent viscosity and the turbulent Prandtl number Pr t is a model constant taken equal to 0:9. The second THF closure is represented by the so-called AHFMNRG+ model. This model expresses the THF through an implicit algebraic expression having the form  Fig. 8. Turbulent Prandtl number profiles. size of the mesh in the wall-normal direction is such as to ensure that the condition yþ < 1 is satisfied in all the calculations. Details on the mesh sensitivity study are given in Appendix B. The governing equations have been solved in their steady-state formulation, using the SIMPLE algorithm  for pressure-velocity coupling and a second-order spatial discretisation. In terms of boundary conditions, following the DNS settings, all the calculations have been performed at a fixed bulk Reynolds number Reb ¼ 5639, corresponding to Res ¼ 180 in forced convection and Richardson number Rib ¼ 0:5. The walls are treated as no-slip boundaries at a fixed temperature, whilst periodicity has been enforced in the streamwise direction. For the transport equation of the temperature variance which is solved by the AHFM-NRG+ model, a Dirichlet condition enforces T 0 ¼ 0 at the walls. 2 4.2. Turbulence models Two different turbulence models have been assessed against the reference DNS data. The two models employ the same closure for the turbulent momentum flux, i.e. the low-Reynolds linear k e proposed in  and implemented in STAR-CCM+ , with a linear constitutive relationship between the Reynolds stress tensor R u0 T 0 ¼ C t0 k 2 C t1 R $T þ C t2 $u u0 T 0 þ C t3 bT 0 g þ C t4 A u0 T 0 e ð5Þ where A is the Reynolds stress anisotropy tensor and C t0 ; . . . ; C t4 are model coefficients. The previous expression has been originally proposed for its application to the natural and mixed convection regimes at unitary Prandtl number in the work of , which is henceforth referred to as AHFM-2005. This closure has been successively extended to the forced convection regime and to a broad range of Pr values in the AHFM-NRG+ model [23,25,34]. This has been achieved by expressing the model coefficients C t1 and C t3 as a function of the Reynolds and the Prandtl numbers and of the Rayleigh and the Prandtl numbers, respectively, as: C t1 ¼ 0:053lnðReb PrÞ 0:27 with Reb Pr P 180 ð6Þ C t3 ¼ 4:5 109 ðlog10 ðRaPrÞÞ þ 2:5 with 1 6 RaPr 6 1017 7 ð7Þ It is worth to point out that, given the imposed value of Reb for the present calculations, the Pr ¼ 0:01 case results in a bulk Peclet number Peb ¼ Reb Pr of about 56. This value is well below the limit of 180 reported in Eq. (6) and, more in general, can be regarded as rather low in the RANS framework. Although it would be desirable to have validation data available at higher Peb , in the case of low Prandtl fluids this is hampered by the high value of Reb , and therefore by the elevated computational cost that would be required to 352 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Table 3 Coefficients for the AHFM model. AHFM-2005 AHFM-NRG+ C t0 C t1 C t2 C t3 C t4 R 0.15 0.2 0.6 Eq. (6) 0.6 0.6 0.6 Eq. (7) 1.5 0.0 0.5 0.5 Fig. 9. RANS results at Pr ¼ 1: mean velocity (a), mean temperature (b) and temperature RMS (c) profiles compared with the reference DNS. generate such a database. Therefore, at the present stage, most of the high-fidelity numerical data for low Prandtl fluids are available at low Peb . For the current calculations, the value of C t1 corresponding to Peb ¼ 180 has been employed for the Pr ¼ 0:01 case. In the AHFM-NRG+ the unknown temperature variance T 0 present in Eq. (5) is evaluated by solving a suitable transport equation for this quantity. The temperature variance dissipation appears in this equation in an unclosed form, and this latter parameter is evaluated in the AHFM-NRG+ model by assuming a constant thermal to mechanical time-scale ratio R , resulting in a three-equation 2 k e T 0 model. The assumption of a constant value for the time-scale ratio allows to simplify the problem significantly; more in detail, such an assumption has been found to be reasonable in a broad range of flow configurations, especially at unitary Prandtl number and away from solid walls [35,36]. Nevertheless, it has 2 to be pointed out that it has also been observed that R is not constant and shows a marked dependency on the wall-distance and on controlling parameters such as Reb and Pr in several flow configurations, in particular when buoyancy effects are non-negligible [33,37]. A summary of the model coefficients for the AHFM-2005 and the AHFM-NRG+ models is reported in Table 3. 4.3. Results and discussions This section discusses the results obtained in the steady-state RANS simulations performed with the two different models for the THF described above, i.e. the classical Reynolds analogy and the AHFM-NRG+. The DNS database presented in Section 3 has been used to assess the RANS results and to compare the performance of the two turbulence models. Since the DNS and the RANS D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 calculations have been performed at the same fixed bulk Reynolds number, all the variables in the plots presented in this section are normalised by the corresponding bulk quantities, and this is indicated by the superscript ‘⁄’. Fig. 9 shows a comparison between the DNS and the two RANS models at Pr ¼ 1. For the mean velocity, at this value of the molecular Prandtl number, it can be seen that the two RANS closures give similar results. In addition, both RANS calculations are in a good agreement with the reference DNS in the near-wall region as well as in the central part of the channel. With respect to the mean temperature, the Reynolds analogy seems to overpredict the temperature gradient in the thermal boundary layer region. This model overestimates the thickness of the thermal boundary layer significantly and this, in turn, results in the inability o reproduce the mean temperature plateau which is observed in the reference DNS in the central part of the channel. On the other hand, the AHFM-NRG+ is in a reasonable agreement with the reference data over the entire height of the channel, although this model results in a slightly higher mean temperature value with respect to the reference data between y 0:05 and y 0:35, which corresponds to an underprediction of the temperature in the specular region of the channel on the cold wall side. None of the two RANS closures is able to correctly reproduce the presence of the temperature 353 inversion observed in the reference data at this Pr value. Since the algebraic model solves an additional transport equation for the temperature variance, the predicted temperature fluctuation can be compared against those obtained in the reference DNS simulations. Fig. 9 shows that the AHFM-NRG+ is able to reproduce the profile of the temperature fluctuations observed in the DNS results from a qualitative point of view, although the value of the nearwall peaks is underestimated. In the bulk region of the channel the deviation with respect to the DNS calculations is less evident, even if the RANS model is not able to reproduce the flat hT 0 irms profile observed between y 0:35 and y 0:65. The results obtained for Pr ¼ 0:1 are reported in Fig. 10. The two RANS models are in a good agreement with each other with respect to the mean velocity field. The RANS results for this variable are also in a reasonable agreement with the reference DNS, although it appears that both models overpredict the value of the friction Reynolds number corresponding to the imposed Reb value. Also, a flatter velocity profile is observed in the central part of the channel with respect to the reference DNS. For the mean temperature profiles, the discrepancy between the Reynolds analogy approach and the DNS, already observed at Pr ¼ 1, is even more evident here. The Reynolds analogy overpredicts the thickness of the boundary layer substantially, resulting in a significant discrepancy with respect to Fig. 10. RANS results at Pr ¼ 0:1: mean velocity (a), mean temperature (b) and temperature RMS (c) profiles compared with the reference DNS. 354 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Fig. 11. RANS results at Pr ¼ 0:01: mean velocity (a), mean temperature (b) and temperature RMS (c) profiles compared with the reference DNS. the reference DNS also in the bulk region of the channel. Conversely, the AHFM-NRG+ is in a good agreement with the reference data in the central part of the channel, whilst it seems to slightly underpredict the temperature gradient in the thermal boundary layer. Lastly, it can be noted that the latter model is able to provide a reasonable prediction of the temperature fluctuations at this Pr value. More in detail, the AHFM-NRG+ results in a lower value of hT 0 irms over the entire height of the channel, but the deviation from the DNS results appear smaller with respect to the Pr ¼ 1 case. A comparison for the Pr ¼ 0:01 case is shown in Fig. 11. For this value of the molecular Prandtl number the overprediction of the friction Reynolds number by both RANS models, already noted at Pr ¼ 0:1, can be clearly observed. This results in a higher velocity gradient at the wall and flatter velocity profiles in the central part of the channel in comparison with the reference DNS. These deviations with respect to the reference DNS appear slightly more evident in the AHFM-NRG+ results. For the mean temperature, the plateau in the centre of the channel can still be observed in the DNS profile, although its extent is reduced with respect to the higher Pr cases due to the increased thickness of the thermal boundary layer. The Reynolds analogy fails to predict the presence of this plateau region and also largely overpredicts temperature gradient in the proximity of the wall. Similarly to what has been observed in the Pr ¼ 0:1 case, the AHFM-NRG+ predicts a slightly lower temperature gradient with respect to the DNS results in the thermal boundary layer and it is in a good agreement with the reference data in the central region of the channel. Moreover, the algebraic THF closure is able to qualitatively predict the change in the temperature fluctuations profile which is observed in the DNS calculations at this low Pr value; in particular, it is in good agreement with the DNS in the bulk region of the channel, whilst it is overestimating the value of hT 0 irms moving closer to the walls. 5. Conclusions Mixed convection is an extremely important phenomenon which is crucial in many practical applications, yet its mechanisms are not fully understood and only a limited number of studies are available. In this work, DNS of the turbulent channel flow with unstable thermal stratification effects have been performed. The value of the Prandtl number has been varied systematically, in order to obtain new insights on the effects of this parameter on the turbulent statistics in this flow configuration. A particular focus has been put on the case of low Prandtl numbers, since the low values of Pr are of particular relevance in nuclear reactor applications based on the use of liquid metals as coolant. It has been observed that: D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 The main effect of the buoyancy is to increase the turbulence levels in the near wall region as indicated by the increase in the value of the Res with respect to forced convection for the same value of the bulk Reynolds number. This effect is increased as the Prandtl number decreases. This phenomenon is a consequence of the fact that under unstable thermal stratification the large scales rolls induced by the thermal plumes tend to confine low-speed streaks near the wall region which are intensified by buoyancy-driven instabilities. In the bulk region, the velocity profile tends to level off in the presence of buoyancy effects. Compared to the forced convection regime, the velocity fluctuations are consistently larger in the bulk region and this effect is accentuated at reduced values of the Prandtl number. The thermal field is affected by the combination of the enhanced turbulence in the flow field and the Prandtl number. It has been observed that the temperature levels off in the bulk region compared to the forced convection case. By increasing the Prandtl number, the gradient of the temperature increases at the wall, decreases in the remaining of the domain and temperature inversion is observed for Pr ¼ 1. The fluctuation of the temperature is less intense in the mixed convection regime and it is reduced as the Prandtl number decreases. This effect is due (1) to a better thermal mixing attained with a more turbulent flow field and hence with reduced temperature fluctuations, (2) to the more effective conduction filtering resulting from the larger thermal conductivity of the fluid. It has been observed that large thermal structures with high spanwise aspect-ratio appear in the mixed convection regime. Therefore, the DNS requires significantly larger domains compared to the forced convection case. Since the size of the thermal structures becomes larger by reducing the Prandtl number, the size of the computational domain must increase accordingly. The DNS results have been used as a reference for the assessment of two different RANS turbulence models, i.e. the classical Reynolds analogy based on the use of a constant turbulent Prandtl number and a recently proposed algebraic closure named AHFMNRG+. It is has been found that: The breakdown of the Reynolds analogy in mixed convection leads to a significant deviation from the DNS data in the results obtained with this closure for the mean temperature, even at 355 unitary Prandtl number. The limitations of this approach for the considered case appear even more evident at low values of the molecular Prandtl number. The algebraic model improves the prediction of the thermal field significantly and a fairly good agreement with the DNS results is observed over the entire range of the Prandtl numbers considered in this study. Furthermore, within the limitation of the steady-state RANS approach, the AFHM-NRG+ is also able to provide reasonable qualitative results for the temperature fluctuations. In the near future this work could be extended to investigate the effect of the other non-dimensional parameters, i.e. Richardson and Reynolds numbers, on the flow and the thermal fields. Such an extension in the parameter space would also be beneficial for the development and validation of the RANS models. In particular, for the AHFM-NRG+ model, a more accurate prediction of the temperature fluctuations could be attained by further refining the formulation of the model to account for a variable mechanical-tothermal time-scale ratio. Conflict of interest The authors declared that there is no conflict of interest. Acknowledgments The work described in this paper is funded by the Dutch Ministry of Economic Affairs and by the Euratom research and training program 2014–2018 under grant agreements Nos. 662186 (MYRTE) and 654935 (SESAME). Part of the simulations presented in this paper were performed on the Świerk Computing Centre in the Department of Complex System in the National Centre for Nuclear Research, Poland. The authors acknowledge the insightful comments of the anonymous referees that helped to improve the quality of the work. Appendix A. Two point correlations In this appendix, the spatial two point correlation functions in the spanwise and streamwise directions are reported for the mixed convection simulations. The correlation functions are computed for the spanwise and streamwise component of the velocity as well as for the temperature, and are evaluated at the wall-normal location y=d ¼ 0:5, with d the half height of the channel. In Figs. 12–14 the Fig. 12. Spatial two point correlation function in the streamwise (left) and spanwise (right) directions for the case at Pr ¼ 1. 356 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Fig. 13. Spatial two point correlation function in the streamwise (left) and spanwise (right) directions for the case at Pr ¼ 0:1. Fig. 14. Spatial two point correlation function in the streamwise (left) and spanwise (right) directions for the case at Pr ¼ 0:01. two point correlation functions are reported for the cases at Pr ¼ 1; Pr ¼ 0:1 and Pr ¼ 0:01, respectively. It can be observed that the two-point correlation functions show a periodic oscillation behavior instead of approaching asymptotically the zero value, as typical in many accurate forced convection DNS. The oscillations are more significant at lower Prandtl numbers and they have been observed also on the preliminary coarse DNS simulation performed on a much large domain. This suggests that these oscillations of the two-point correlation functions are the results of the physics of the problem instead of the numerical approximation. The periodic oscillations of the two-point spatial correlation function might be due to the presence of the large roll-up vortices in the domain. As the Prandtl number decreases these vortices become larger, consequently large oscillations of the two-point correlation functions are observed. given in Table 4. A refinement factor of two in both the streamwise and the wall-normal directions has been applied between two consecutive meshes. The same factor has been employed to refine the size of the first cell at the wall in the normal direction (Dywall ). The resulting yþ ; Res and Nusselt number Nu evaluated at the bottom wall are summarised in Table 5 for Pr ¼ 1; 0:1 and 0:01. In addition, the relative deviations () of the calculated values of Res and Nu with respect to the most refined mesh are reported in the same tables. It can be noticed that all the meshes satisfy the condition yþ < 1 over the entire range of Prandtl numbers considered in the study, although the yþ value tends to increase with decreasing Pr due to the higher values of Res associated with low Prandtl numbers. In addition, it can be seen that Mesh 2 results in small deviations (less than 0:3% for the considered quantities) with respect to Mesh 3 for all the values of the molecular Prandtl Appendix B. Mesh sensitivity for the RANS approach Table 4 Parameters defining the three grids employed for the mesh sensitivity study. This appendix reports the details of the mesh sensitivity study performed for the RANS calculations. The mesh sensitivity has been assessed for all the Pr values considered in this study and employing the Reynolds analogy for the THF closure. The parameters defining the three grids considered for the sensitivity study are Nx Ny Dywall (m) N tot Mesh 1 Mesh 2 Mesh 3 50 75 0.006 3750 100 150 0.003 15000 200 300 0.0015 60000 357 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358 Table 5 Calculated yþ ; Res and Nu at the bottom wall on the three meshes for the three Prandtl numbers. Pr ¼ 1 yþ Res Res ð%Þ Nu Nu ð%Þ Pr ¼ 0:1 Pr ¼ 0:01 Mesh 1 Mesh 2 Mesh 3 Mesh 1 Mesh 2 Mesh 3 Mesh 1 Mesh 2 Mesh 3 0.608 203.21 0.80 47.39 1.17 0.302 201.95 0.18 46.96 0.26 0.151 201.59 – 46.84 – 0.664 222.15 0.55 18.63 0.16 0.331 221.21 0.13 18.61 0.05 0.166 220.93 – 18.60 – 0.711 238.24 0.36 7.39 0.40 0.356 237.62 0.10 7.41 0.13 0.178 237.39 – 7.42 – Fig. 15. Mesh sensitivity at Pr ¼ 1: mean velocity (a), turbulent kinetic energy (b). number. The profiles of the mean velocity and turbulent kinetic energy obtained with the three different meshes for Pr ¼ 1 are reported in Fig. 15. It can be observed that the three meshes give practically identical results for the mean velocity, whilst Mesh 2 is closer to the most refined grid when considering the near-wall peak of the turbulent kinetic energy; on the other hand, the value of this peak is visibly overestimated with Mesh 1. Therefore Mesh 2 has been selected to perform all the remaining simulations in the present work. Appendix C. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer. 2018.07.150. References  J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developed channel flow at low Reynolds number, J. Fluid Mech. 177 (1987) 133–166.  M. Kozuka, Y. Seki, H. Kawamura, DNS of turbulent heat transfer in a channel flow with a high spatial resolution, Int. J. Heat Fluid Flow 30 (3) (2009) 514– 524.  N. Kasagi, O. Iida, Progress in direct numerical simulation of turbulent heat transfer, in: Proceedings of the 5th ASME/JSME Joint Thermal Engineering Conference, San Diego, California, 1999.  S. Pirozzoli, M. Bernardini, P. Orlandi, Passive scalars in turbulent channel flow at high Reynolds number, J. Fluid Mech. 788 (2016) 614–639.  I. Tiselj, E. Pogrebnyak, C. Li, A. Mosyak, G. Hetsroni, Effect of wall boundary condition on scalar transfer in a fully developed turbulent flume, Phys. Fluids 13 (4) (2001) 1028–1039.  C. Flageul, S. Benhamadouche, E. Lamballais, D. Laurence, DNS of turbulent channel flow with conjugate heat transfer: effect of thermal boundary conditions on the second moments and budgets, Int. J. Heat Fluid Flow 55 (2015) 34–44.  I. Tiselj, C. Cizelj, DNS of turbulent channel flow with conjugate heat transfer at Prandtl number 0.01, Nucl. Eng. Des. 253 (2012) 153–160.  I. Tiselj, R. Bergant, B. Mavko, I. Bajsicc, G. Hetsroni, DNS of turbulent heat transfer in channel flow with heat conduction in the solid wall, ASME J. Heat Transfer 123 (5) (2001) 849–857.  I. Otić, G. Grötzbach, Statistical Analysis of Turbulent Natural Convection in Low Prandtl Number Fluids, Springer Berlin Heidelberg, 2005.  G. Grötzbach, Direct numerical simulation of laminar and turbulent Bénard convection, J. Fluid Mech. 119 (1982) 27–53.  G. Silano, K.R. Sreenivasan, R. Verzicco, Numerical simulations of RayleighBénard convection for Prandtl numbers between 101 and 104 and rayleigh numbers between 105 and 109 , J. Fluid Mech. 662 (2010) 409–446.  T. Mizushina, F. Ogino, N. Katada, Ordered motion of turbulence in a thermally stratified flow under unstable conditions, Int. J. Heat Mass Transfer 25 (9) (1982) 1419–1425.  K. Fukui, M. Nakajima, Unstable stratification effects on turbulent shear flow in the wall region, Int. J. Heat Mass Transfer 28 (12) (1985) 2343–2352.  K. Fukui, M. Nakajima, H. Ueda, Coherent structure of turbulent longitudinal vortices in unstably-stratified turbulent flow, Int. J. Heat Mass Transfer 34 (9) (1991) 2373–2385.  J.A. Domaradzki, R.W. Metcalfe, Direct numerical simulations of the effects of shear on turbulent Rayleigh-Bénard convection, J. Fluid Mech. 193 (1988) 499–531.  O. Iida, Kasagi, Direct Numerical Simulation of unstably stratified turbulent channel flow, J. Heat Transfer 119 (1) (1997) 53–61.  F. Zonta, A. Soldati, Effect of temperature dependent fluid properties on heat transfer in turbulent mixed convection, J. Heat Transfer 136 (2) (2014) 1–11.  S. Sid, Y. Dubief, V.E. Terrapon, Direct numerical simulation of mixed convection in turbulent channel flow: on the Reynolds number dependency of momentum and heat transfer under unstable stratification, in: Proceedings of the 8th International Conference on Computational Heat and Mass Transfer, ICCHMT, 2015.  S. Pirozzoli, M. Bernardini, R. Verzicco, P. Orlandi, Mixed convection in turbulent channels with unstable stratification, J. Fluid Mech. 821 (2017) 482– 516.  G. Grötzbach, Challenges in low-Prandtl number heat transfer simulation and modelling, Nucl. Eng. Des. 264 (2013) 41–55.  A. Shams, (U)RANS: turbulence modelling, SESAME Project: Lecture Series on Thermohydraulics and Chemistry of Liquid Metal Cooled Reactors, von Karman Institute for Fluid Dynamics, 10th-14th April 2017, 2017. 358 D. De Santis et al. / International Journal of Heat and Mass Transfer 127 (2018) 345–358  F. Roelofs, A. Shams, I. Otic, M. Böttcher, M. Duponcheel, Y. Bartosiewicz, D. Lakehal, E. Baglietto, S. Lardeau, X. Cheng, Status and perspective of turbulence heat transfer modelling for the industrial application of liquid metal flows, Nucl. Eng. Des. 290 (2015) 99–106.  A. Shams, F. Roelofs, E. Baglietto, S. Lardeau, S. Kenjeres, Assessment and calibration of an algebraic turbulent heat flux model for low-Prandtl fluids, Int. J. Heat Mass Transfer 79 (2014) 589–601.  A. De Santis, A. Shams, Application of an algebraic turbulent heat flux model to a backward facing step flow at low Prandtl number, Ann. Nucl. Energy 117 (2018) 32–44.  A. Shams, Assessment and calibration of an algebraic turbulent heat flux model for high Rayleigh number flow regimes, in: 8th Conference on Severe Accident Research, May 16th-18th, Warsaw, Poland, 2017.  J.W.L. Paul, F. Fischer, S.G. Kerkemeier, nek5000 Web page, 2008. <http:// nek5000.mcs.anl.gov>.  J. Wang, K. Xia, Spatial variations of the mean and statistical quantities in the thermal boundary layers of turbulent convection, Eur. Phys. J. B – Condens. Matter Complex Syst. 32 (1) (2003) 127–136.  G. Grötzbach, Anisotropy and buoyancy in nuclear turbulent heat transfer critical assessment and needs for modelling, Forschungszentrum Karlsruhe FZKA 7363 (2007).  SIEMENS PLM Software, STAR-CCM+ User Guide - Version 11.06, 2016.  S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (10) (1972) 1787–1806.  F.S. Lien, W.L. Chen, M.A. Leschziner, Low-Reynolds number eddy-viscosity modelling based on non-linear stress-strain/vorticity relations, Eng. Turbul. Model. Exp. (1996) 91–100.  S. Kenjereš, K. Hanjalić, Convective rolls and heat transfer in finite-length Rayleigh-Bénard convection: a two-dimensional numerical study, Phys. Rev. E 62 (6) (2000) 7987–7998.  S. Kenjereš, S. Gunarjo, K. Hanjalić, Contribution to elliptic relaxation modelling of turbulent natural and mixed convection, Int. J. Heat Fluid Flow 26 (4) (2005) 569–586.  A. Shams, Towards the accurate numerical prediction of thermal hydraulic phenomena in corium pools, Ann. Nucl. Energy 117 (2018) 234–246.  H. Kawamura, H. Abe, Y. Matsuo, DNS of turbulent heat transfer in channel flow with respect to Reynolds and Prandtl number effects, Int. J. Heat Fluid Flow 20 (3) (1999) 196–207.  K. Hanjalić, One-point closure models for buoyancy-driven turbulent flows, Annu. Rev. Fluid Mech. 34 (1) (2002) 321–347.  I. Otić, G. Grötzbach, M. Wörner, Analysis and modelling of the temperature variance equation in turbulent natural convection for low-Prandtl number fluids, J. Fluid Mech. 525 (2005) 237–261.