A Comprehensive Evaluation of Turbulence Models for Predicting Heat Transfer in Turbulent Channel Flow across Various Prandtl Number Regimes

: Turbulent heat transfer in channel flows is an important area of research due to its simple geometry and diverse industrial applications. Reynolds-Averaged Navier–Stokes (RANS) models are the most-affordable simulation methodology and are often the only viable choice for investigating industrial flows. However, accurate modelling of wall-bounded flows is challenging in RANS, and the assessment of the performance of RANS models for heated turbulent channel flow has not been sufficiently investigated for a wide range of Reynolds and Prandtl numbers. In this study, five RANS models are assessed for their ability to predict heat transfer in channel flows across a wide range of Reynolds and Prandtl numbers ( Pr ) by comparing the RANS results with respect to the corresponding Direct Numerical Simulation data. The models include three Eddy Viscosity Models (EVMs): standard k − ϵ , low Reynolds number k − ϵ LS , and k − ω SST , as well as two Reynolds Stress Models (RSMs): Launder–Reece–Rodi and Speziale–Sarkar–Gatski models. The study analyses the Reynolds number effects on turbulent heat transfer in a channel flow at a Pr of 0.71 for friction Reynolds number values of 180, 395, 640, and 1020. The results show that all models accurately predict velocity across all Reynolds numbers, but the accuracy of mean temperature prediction drops with increasing Reynolds number for all models, except for the k − ω SST model. The study also analyses the Pr effects on turbulent heat transfer in a channel flow with Pr values between 0.025 and 10.0. An error analysis is performed on the results obtained from different turbulence models, and it is shown that the k − ω SST model has the smallest error for the predictions of the mean temperature and Nusselt number for high-Prandtl-number flows, while the low Reynolds number k − ϵ LS model shows the smallest errors for low-Prandtl-number flows at different Reynolds numbers. An analytical solution is utilised to identify Pr effects on forced convection in a channel flow into three different regimes: analytical region, transitional region, and turbulent diffusion-dominated region. These regimes are helpful to discuss the validity of the models in relation to the Pr . The findings of this paper provide insights into the performance of different RANS models for heat transfer predictions in a channel flow.


Introduction
Heat transfer in turbulent channel flow is widely used for the fundamental understanding and modelling of forced convection within ducts due to its simple geometry [1].Under the condition where the fluid properties remain constant, the energy equation in this configuration becomes decoupled from the momentum equations, and the temperature can be considered as a passive scalar carried by the background fluid motion.The study of scalar transport and its modelling have gained considerable attention in recent decades due to the importance of turbulent transport of passive scalars (not only the temperature, but also humidity, pollutants, or other chemical species) in many industrial applications [2].
However, accurate experimental measurements of near-wall statistics are relatively scarce due to the difficulty of making precise measurements in the presence of large temperature gradients [3].Numerical simulation provides an alternative means of obtaining useful information on these flows.There are three types of simulations depending on the level of the resolved scales of turbulence: Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), and Reynolds-Averaged Navier-Stokes (RANS) simulations.
DNS resolves the full range of spatial and temporal scales of turbulence by numerically solving the Navier-Stokes equations without any turbulence model.However, the grid resolution of the discretised domain must be of the order of the Kolmogorov length scale η = (ν 3 /ϵ) 1/4 , where ϵ is the dissipation rate and ν is the kinematic viscosity.The LES approach is a compromise between accuracy and computational cost compared to the DNS and RANS methodologies.In LES, the grid resolution only captures the energycontaining large eddies, while the small-scale physics is modelled using sub-grid scale (SGS) models.However, even LES is often computationally expensive for many industrial and environmental flows with very high Reynolds numbers.By contrast, the RANS methodology models all scales of motion and is the most-computationally affordable technique.However, RANS models often exhibit limited accuracy for some flows, especially those involving heat transfer, where turbulence is anisotropic close to solid boundaries.Additionally, RANS models may not capture large-scale unsteady and coherent turbulent structures that significantly affect heat transfer in some cases.The following is a summary of three different simulation methods used to study turbulent heat transfer in channel flow.
The first DNS of a thermal flow was conducted by Kim and Moin [4] at a friction Reynolds number, Re τ = u τ h/ν (where h is the channel half-height), of 180 and for Prandtl numbers Pr of 0.1, 0.71, and 2.0.The friction Reynolds number Re τ is a measure of the turbulent behaviour of the flow, and Pr is the ratio of the momentum diffusivity to the thermal diffusivity.In the simulations performed by Kim and Moin [4], Dirichlet boundary condition was specified at the channel walls.Kasagi and Ohtsubo [5,6] carried out DNS for Pr = 0.71 and 0.025 at a Reynolds number of Re τ = 150.They used constant heat flux boundary conditions for both the top and bottom walls instead of a Dirichlet boundary condition.Kawamura et al. [7] extended DNS to a wide range of Prandtl numbers from Pr = 0.025 to 5.0 at Re τ = 180.They examined the Prandtl number dependence of a fully developed turbulent channel flow and reported that the turbulent Prandtl number, Pr t , is independent of Pr except for Pr < 0.1.Kawamura et al. [1] also studied the Reynolds number dependence by carrying out fully developed turbulent channel flow simulations at Re τ = 180 and 395 for Pr = 0.025, 0.2, and 0.71.Statistical quantities such as the temperature variance, turbulent heat flux, and turbulent Prandtl number were obtained, and the effects of the Reynolds and Prandtl numbers were analysed.Abe et al. [8] carried out the turbulent channel flow simulations at Reynolds number values of Re τ = 1020 with the Pr = 0.025 and 0.71 by Kawamura et al. [1].
LES has recently contributed significantly to turbulent heat transfer in channel flows due to the advancement in sub-grid scale (SGS) modelling, such as the dynamic eddy viscosity model proposed by Germano et al. [9].Moin et al. [10] conducted the first LES with a dynamic SGS model for compressible turbulent flow with scalar transport in a channel flow configuration at Re τ = 165, considering three different Prandtl numbers (Pr = 0.1, 0.71, and 2.0).Later, Wang and Pfletcher [11] used the dynamic SGS model to simulate turbulent channel flows with heat transfer for Re τ = 162 and 187.However, they used isothermal walls with a large temperature difference between the bottom and top walls, instead of the constant heat flux boundary conditions.Wang et al. [12] conducted heat transfer analysis in turbulent channel flows using the dynamic SGS model with both isothermal walls and constant heat flux boundary conditions for Prandtl numbers ranging from 1 to 100 at Re τ = 180.Subsequently, an LES-based lattice Boltzmann framework was utilized to analyse heat transfer characteristics in turbulent channel flows at Re τ = 180 for Prandtl numbers of 0.025 and 0.71 [13].Currently, heat transfer in turbulent channel flow is a benchmark case for the assessment of model performances and the development of new SGS models.
While LES is a promising method for analysing turbulent heat transfer in ducts, it may not always be affordablein industrial calculations due to its high computational cost.In contrast, RANS is a more-computationally affordable option for such applications, but the accuracy of the modelling of wall-bounded flows poses a significant challenge for RANS models.Therefore, it is important to thoroughly assess the capabilities of RANS models and identify their limitations in terms of the prediction of heat transfer characteristics in turbulent channel flows.
To the best of the authors' knowledge, there has been little research into the assessment of RANS models for heat transfer predictions in turbulent channel flows across various Prandtl number regimes.An algebraic closure for Reynolds-Averaged turbulent scalar-flux was proposed by Abe and Suga [14] for heat transfer in turbulent channel flows, but the model was only assessed at Re τ = 180 while using different wall boundary conditions.Later, Pozorski et al. [15,16] used heated turbulent channel flow as a wall-bounded test case for velocity scalar probability density function at Re τ = 180.A detailed investigation of the variation of the Nusselt number with the Reynolds and Prandtl numbers was performed by Wei [17], and the predictions of the correlations proposed in the literature have been compared with the DNS results.Recently, Mathur et al. [18] evaluated the performance of different scalar flux models in the case of low-Prandtl-number flows along with the influence of varying the turbulent Prandtl number on the accuracy of the model predictions.Mollik et al. [19] assessed several RANS modelling methodologies at various Reynolds numbers (Re τ = 180, 395, 640, and 1020).However, these studies did not systematically analyse the model performance with respect to variations in the Reynolds and Prandtl numbers.
This paper aims to assess the performance of turbulence models in terms of the predictions of the heat transfer characteristics in turbulent channel flow configurations for incompressible fluids with walls subjected to constant heat fluxes for a wide range of Re τ and Pr.The information obtained from this exercise is utilised to identify the relative merits and demerits of different RANS turbulence models.In Section 2, the governing equations and flow setup for a channel flow with constant heat flux on the bottom and top walls are presented.This is accompanied by the summary of the parameters (Re τ , Pr) of the DNS database, which is utilised to assess the performance of the models.In Section 3, the formulations of the Reynolds-Averaged closures for two classes of RANS models, Eddy Viscosity Models (EVMs) and Reynolds Stress Models (RSMs), are presented, as are the wall functions for the velocity and temperature fields.Section 4 presents the predictions of different turbulence models in terms of the first and second moments of the relevant quantities for different values of Re τ and Pr, and these predictions are compared with respect to the corresponding quantities extracted from explicitly Reynolds-Averaged DNS data.Finally, conclusions are drawn from the findings of the above analysis in Section 5.

Governing Equations and Flow Setup
This study considers fully developed channel flow for incompressible fluids, where the bottom and top walls are uniformly heated by a constant heat flux q w , as illustrated in Figure 1.The sides along the z-axis are taken to be periodic, and the flow direction is aligned with the x-axis; a periodic boundary condition is also imposed in the x direction.The governing equations are non-dimensionalised by the channel half-width h (i.e., x * i = x i /h), the friction velocity u τ = |τ w |/ρ, the kinematic viscosity ν, and the friction temperature T τ = q w /ρc p u τ .The non-dimensionalised governing equations are shown as follows (where t * = h/u τ ).Continuity Equation: Momentum Equation: (2) In this case, periodic boundary conditions are applied in the x-direction, which can result in energy accumulation in the enclosed heated channel flow system.The bulk mean temperature increases linearly with respect to x for the constant heat flux boundary condition, which contrasts with the periodic boundary condition in the x-direction.To address this issue, a temperature transformation is used, where the energy equation is established based on the local temperature with the removal of temperature at the wall, defined as Θ = T − T w .After this transformation, the sink term −u + 1 /⟨u + ⟩ arises, which balances the energy input from the bottom and top walls.Here, ⟨u + ⟩ represents the average velocity over the channel section.More details of the energy equation transformation can be found in [20,21].With this transformation, the energy equation becomes: The boundary conditions in this case are u + i = 0 and Θ + = 0 at y = 0 and y = 2h.In the momentum equation, the normalized Reynolds stress components −u +′ i u +′ j , while in the energy equation, the Reynolds heat fluxes −u +′ i Θ +′ are the unclosed terms and, thus, need modelling.A detailed introduction to the models used for these unclosed terms is provided in the next section.The DNS database of heated turbulent channel flow was established by Kawamura et al. [1,7] over a wide range of Pr (from 0.025 to 10.0) and Re τ (from 180 to 1020).The parameters of the DNS database are listed in Table 1, with the available setups marked by a check mark.The RANS simulations are performed using an open-source code called Code Saturne, which is an unstructured finite-volume code designed to compute turbulent flows in industrial applications [22].This code uses a fractional step method based on a prediction-correction algorithm for pressure/velocity coupling (SIMPLEC) to solve the Navier-Stokes equations for Newtonian incompressible flows, with a second-order central differencing scheme for spatial gradients and an Euler-implicit scheme for time integration [23].

Model Formulations for Reynolds-Averaged Closures
As discussed in the Introduction, RANS has the advantage of being more-computationally affordable compared to DNS and LES, and this makes the RANS framework the mostwidely used approach in industrial Computational Fluid Dynamics (CFD).It is often the only viable choice for parametric investigations of industrial applications over a wide range of parameters, such as the Reynolds number and Prandtl number.Currently, there are numerous RANS models available that vary in the way Reynolds stresses are modelled.The two main categories are Eddy Viscosity Models (EVMs) and Reynolds Stress Models (RSMs) [24], which are introduced individually in the following subsections.

Eddy Viscosity Models
Eddy-viscosity-based models, both linear and nonlinear, have become the most-widely used models in turbulence modelling due to their simplicity and efficiency.Among these models, the linear eddy viscosity methodologies are the simplest.These models algebraically evaluate the Reynolds stress using Boussinesq's hypothesis, which relates the Reynolds stress tensor to the mean velocity gradients in a turbulent flow.Boussinesq's hypothesis can be expressed as follows: Here, k is the turbulent kinetic energy, ν t is the eddy viscosity, which unlike the kinematic viscosity, is not a fluid property and can only be scaled using appropriate velocity and length scales u turb and l turb [25].In general, u turb is expected to scale as k 1/2 , while l turb can be approximated as k 3/2 /ϵ, where ϵ represents the turbulent dissipation rate.As a result, twoequation models are considered for turbulence modelling.Among these, the k − ϵ and k − ω models are the most-widely used due to their physical rationale, simplicity, and numerical robustness.In these models, the eddy kinematic viscosity scales as ν t = C µ k 2 /ϵ and ν t = C ω k/ω, respectively, where ω is the specific turbulence dissipation rate, which is closely related to ϵ and provides a measure of the local dissipation rate per unit turbulent kinetic energy.However, both the k − ϵ and k − ω models have limitations in RANS simulations.The k − ϵ model is not accurate for complex flows such as swirling flows, recirculating flows, flows with secondary flow features, and flows in non-circular channels.On the other hand, the k − ω model predictions are highly sensitive to freestream turbulence, and this model without any special treatment can produce unrealistic high values of the turbulent kinetic energy production rate at the stagnation point.To address these limitations, the k − ω Shear Stress Transport (SST) model was proposed by Menter [26], which combines the advantages of both the k − ϵ and k − ω models.The use of the k − ω formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sublayer; hence, the k − ω SST model can be used as a low-Reynolds-number turbulence model.The k − ω SST formulation replicates the k − ϵ behaviour in the freestream and, therefore, circumvents the common k − ω problem of being too sensitive to the freestream turbulence.
In the case of scalars, the turbulent scalar flux can be modelled through the Standard Gradient Diffusion Hypothesis (SGDH) [14], which is given as: The SGDH provides an algebraic expression for the turbulent scalar flux as the product of the turbulent diffusivity and the scalar gradient.In this expression, the turbulent diffusivity is represented by ν t /Pr t with Pr t being the turbulent Prandtl number.For heated turbulent channel flow, Pr t typically approaches a value of around 1.0 near the wall and varies between 0.5 and 1.0 away from the wall, except for very low Prandtl numbers (Pr ≤ 0.1) [7].
In our RANS simulations, we adopted a constant value of Pr t = 0.71, which is commonly used as an empirical estimation.The variation of turbulent Prandtl numbers extracted from DNS data [7] shows that Pr t varies with the wall-normal distance, but remains of the order of unity and depends on the Prandtl number.Therefore, a single choice of Pr t is not straightforward; hence, Pr t = 0.7 is chosen for the RANS simulations conducted here because this is the default value used in most commercial codes.This compromise is made because the effects of Pr t on the simulation predictions are not the main objective of this paper.The results can be made to match better with DNS by modifying Pr t for a specific case, but this is hardly helpful for practical applications because the optimum value of Pr t will not be known a priori.To assess the different turbulence models, we combined the SGDH with various closure techniques for Reynolds stresses.This combination facilitates the comparison of the different models and allows for a more-accurate evaluation of their effectiveness in predicting turbulent heat transfer in a channel flow configuration.

Reynolds Stress Transport Models
This methodology involves solving modelled transport equations for As all components of Reynolds stress are obtained as a part of the solution, RSM models demand a greater computational effort than the previously discussed EVMs.This methodology, nevertheless, can replicate the return to isotropy under decaying turbulence and the limiting behaviour of turbulence under the rapid distortion limit, where turbulence acts as an elastic medium.This category of models is expected to be valid for complex flows.
The modelled transport equation for R ij takes the following form [25]: The first term on the right-hand side of the above equation represents the production term, indicating a production rate of R ij .The second term D ij refers to the molecular diffusion of R ij .These two terms D ij and P ij are exact and given by The remaining three terms are the transport of R ij due to pressure strain rate interactions (i.e., Π ij ), the transport of R ij due to vorticity (i.e., Ω ij ), and finally, the molecular dissipation of R ij (i.e., −ϵ ij ), respectively.The terms P ij and D ij are closed and do not require models in the context of second-moment closure.However, other terms need modelling.Among these, the pressure strain rate interaction term is usually the biggest challenge in the modelling process, where most models differ.The pressure strain correlation term Π ij is typically divided into a rapid part, Π R ij , and a slow part, Π S ij .The slow part is responsible for the return to isotropy and reduces the anisotropy of the Reynolds stresses.This part represents self-interaction between fluctuating fields: The rapid part captures the interaction of the fluctuating velocity field with the mean velocity field and is dependent on the mean strain rate and rotation rate.The simplest form of the rapid part of most pressure rates of strain models can be expressed as [25,27]: Here, b ij = R ij /2k − δ ij /3 in Equation ( 7) represents the anisotropy tensor for Reynolds stresses, where δ ij is the Kronecker delta.The quantities S ij and Ω ij in Equation ( 8) are the normalized mean rate of stain and normalised mean rate of rotation, respectively, defined as Several pressure-rate-of-strain models have been proposed for modelling turbulence in fluid flows.One of the earliest and most-widely used pressure-rate-of-strain models is the Launder-Reece-Rodi (LRR) model, which was introduced by Launder et al. [28].Another example of a Reynolds stress model is the Speziale-Sarkar-Gatski (SSG) model [27].Table 2 lists the coefficients for the LRR and SSG models.The slow part is closed by using the model proposed by Rotta [29].The LRR model takes the Rotta coefficient to be constant, whereas SSG model incorporates a dependence of the coefficient on P/ϵ, where P = P ii /2.Additionally, the SSG model has a non-zero value of f (2) , which leads to the non-linear return to isotropy.The SSG model also considers one of the anisotropy invariants η = (b ij b ji /6) of the Reynolds stress tensor.It is worth noting that both the LRR and SSG models are derived models of the original model proposed by Hanjalic and Launder [30].
Table 2. Specifications of the coefficients f (n) for pressure rate of strain models: LRR and SSG.

Near-Wall Treatment
The standard k − ϵ model is a high-Reynolds-number model that is not solved up to the wall.Instead, wall functions are adopted, and the first grid point in the wall-normal direction is typically located in the log-law region.Thus, it can avoid the computational expense of resolving the viscous sublayer.The demarcation point between the buffer layer and the log-law layer is around y + ≈ 30 [25].Moreover, the turning point for the thermal boundary layer between the intermediate layer and the log-law layer is around 24.32.Therefore, the first grid point for the current analysis for the standard k − ϵ model is chosen at y + = 25.The mean velocity profile in Figure 2a shows the existence of the logarithmic region, and the slope of the logarithmic curve is related to the von Karman constant (κ ≈ 0.4) [31].Similarly, the mean temperature profile in Figure 2b exhibits a logarithmic region with the slope of the logarithmic curve of Pr t /κ.
Instead of relying on wall functions, a low-Reynolds-number model, such as k − ϵ LS [32] and k − ω SST [26] models, can be used to more-accurately predict nearwall turbulence.These models are solved all the way to the wall, and the near-wall grid resolution needs to be fine enough to resolve the viscous sublayer (y + < 5.0).In low-Reynolds-number k − ϵ modelling, damping functions are incorporated into the expression for the eddy viscosity and the terms of the transport equations for k and ϵ to account for the near-wall damping.The impact of the choice of the first grid point adjacent to the wall is also illustrated in Figure 2. The standard k − ϵ model employs wall functions and uses a first grid point at y + = 25, while k − ϵ LS model does not use wall functions and is set up with the first grid point at y + = 0.3.
For the thermal field, the three-layer thermal wall function proposed by Arpaci and Larsen [33] is used instead of resolving the near-wall region.The exact expressions for the three different layers are provided in Table 3.The thermal boundary layer in this heat transfer problem is divided into three layers, which are separated by two turning points: y + 1 = (1000/Pr) 1/3 and y + 2 = (1000κ)/Pr t .For the current setup with κ = 0.42 and Pr t = 0.71, the value of y + 2 is 24.32 and is independent of Pr.The values of the log-law layer for velocity and temperature suggest that the log-law layers for these quantities are nearly aligned.This alignment can also be observed in Figure 2.  To conclude, for the heat transfer analysis, a three-layer wall function is used, as the region of the diffusion sublayer is influenced by the Prandtl number.Since a wide range of Pr was considered in this analysis, a mesh that can resolve the diffusion sublayer may not resolve the viscous sublayer and vice versa.The standard k − ϵ, LRR, and SSG Reynolds stress closures do not integrate up to the wall by design and rely on wall functions by ensuring y + = 25 for the wall-adjacent grid points, whereas k − ϵ LS and k − ω SST are low-Reynolds-number models, which integrate up to the wall by ensuring y + < 1 adjacent to the wall and, thus, do not need wall functions.

A Summary of Selected Turbulence Models
Five different turbulence models, categorised by their type (EVM or RSM) and their requirements in terms of near-wall treatment, are considered here, and these models are summarised in Table 4.The standard k − ϵ, k − ϵ LS, and k − ω SST models are all EVMs, while LRR and SSG are RSM models.The k − ϵ LS and k − ω SST models are low-Reynoldsnumber models that do not need wall functions, while the standard k − ϵ, LRR, and SSG models are high-Reynolds-number models that do require the use of a wall function.All of the turbulence models used in this work employ a thermal wall function.Since it covers the entire thermal boundary layer, it allows for flexibility in the location of the first grid point.
The values of y * = y/h for the grid points in the wall-normal direction N y up to the channel half-height for models with a velocity wall function (e.g., standard k − ϵ, LRR, and SSG models) and without a velocity wall function (e.g., k − ϵ LS and k − ω SST models) are shown in Figure 3 for different values of Re τ , which shows that the computational requirement (i.e., number of grid points in the wall-normal direction) is greater for models without wall functions than those with wall functions.Moreover, the computational cost increases with increasing Re τ .Table 5 provides a summary of the computational cost associated with various turbulence models.To demonstrate the differences in the time costs between these models using one core of the CPU (Intel(R) Xeon (R) Silver 4208 CPU @ 2.10 GHz), an example is provided in Table 5 at Re τ = 1020 and Pr = 0.71.Compared to the k − ϵ LS and k − ω SST models, the standard k − ϵ, LRR, and SSG models require less computational cost in the near-wall region due to their adoption of the wall function.The standard k − ϵ model is the least expensive, requiring 0.76 h for the simulation to converge, while k − ω SST is the most-computationally expensive, requiring 21 h.The mesh requirement is widely different between the standard k − ϵ and k − ω SST models.The wall-adjacent grid spacing for the standard k − ϵ model is taken to satisfy y + = 25, whereas the wall-adjacent grid spacing needs to satisfy y + < 1.0 for the k − ω SST model, and therefore, a mesh that provides y + = 0.3 for the wall-adjacent grid points is used.This means the grid points in the wall-normal direction for the k − ω SST model are 12.5-times those used in the standard k − ϵ model for the grids used in this analysis.Moreover, the k − ω SST model represents a stiffer system of differential equations than that in the case of the standard k − ϵ model.Therefore, the k − ω SST model takes longer to reach the steady state than the standard k − ϵ model.The combination of the longer simulation time and the finer grid is responsible for the higher computational cost of the k − ω SST model simulation in comparison to that using the standard k − ϵ model.Although the k − ω SST model is morecomputationally demanding than the other models, it still offers a significant reduction in computational cost compared to DNS.For example, the number of grid points required for DNS at Re τ = 1020 and Pr = 0.71 is 150,000-times larger than that needed for RANS without wall functions and 17.6 million-times larger than that required for RANS models with wall functions.While the computational time of DNS is not available, comparing the number of grid points needed makes it clear that DNS demands significantly more computational resources than RANS.
The boundary conditions for the selected turbulence models are applied as follows.The turbulent kinetic energy is identically zero at the wall due to the no-slip boundary condition.In the case of high-Reynolds-number models that require wall functions and the first grid point to be in the log-layer (e.g., standard k − ϵ, LRR, and SSG Reynolds stress closures), the cell-averaged value is specified for ϵ as u 3 k /κy w , where u k = C 1/4 µ k 1/2 and y w is the wall-normal distance of the first grid point adjacent to the wall.For the k − ωSST model, the wall value of ω is taken to be 10 × 6ν/[β 1 (y w ) 2 ], where β 1 = 0.075 is a model parameter [26].For the low-Reynolds-number k − ϵ model, the wall value of ϵ is specified as 2ν(∇ √ k) 2 .

Heat Transfer Rate in Turbulent Channel Flow for Different Values of Reynolds Number
To investigate the influence of the Reynolds number on the performance of different turbulence models, four different friction Reynolds numbers of 180, 395, 640, and 1020 are considered, each having a Prandtl number of Pr = 0.71.The obtained mean velocity and mean temperature profiles are compared with the corresponding profiles obtained from the DNS [8,14].Furthermore, the Reynolds stress distributions are analysed and compared with the corresponding quantities extracted from DNS data to gain a deeper understanding of the performance of the models considered here.

Mean Streamwise Velocity and Temperature Profiles
The velocity statistics in this configuration are expected to exhibit self-similarity, which is independent of the Reynolds number.It has indeed been found that the mean behaviour of the streamwise velocity and temperature remains unchanged despite variations in Re τ .This can be substantiated by Figure 4, where self-similarity is demonstrated by the variations of the mean profiles of the streamwise velocity and temperature obtained from both DNS and turbulence model predictions.
Turbulence models encounter difficulties in precisely predicting intricate fluid dynamics spanning a broad spectrum of Reynolds numbers, especially in the context of low-Reynolds-number flows undergoing transitional states.Nevertheless, these models perform well across a wide range of Re τ values.To evaluate the qualitative performance of different turbulence models and their trends with increasing Re τ , the Root-Mean-Squared (RMS) deviation is utilised to measure the differences between the values obtained from DNS and the values predicted by the five models.Thus, the RMS of the relative deviation of the mean velocity between the DNS and RANS predictions based on the turbulence model Error(u + ) and the RMS of the relative deviation of the mean temperature between the DNS and RANS results based on the turbulence model Error(Θ + ) are given as follows: Error(Θ In the above equations, u + DNS represents the mean non-dimensional streamwise velocity obtained from DNS and u + RANS represents the mean non-dimensional streamwise velocity obtained from RANS.Similarly, Θ + DNS and Θ + RANS represent the mean non-dimensional temperatures obtained from DNS and RANS, respectively.Since the data points are available only across the log-law layer in the high-Reynolds-number turbulence models and the major discrepancy also occurs in this region, all five models are evaluated in this region.To achieve this, the DNS and low-Reynolds-number models are sampled in the region where the high-Reynolds-number turbulence models remain valid.Figure 5 provides quantitative insight into how the performance of the five turbulence models varies with Re τ .All of the models performed well in predicting the mean streamwise velocity at different Re τ , with all five models having Error(u + ) below 10% at all Reynolds numbers.Among these models, the standard k − ϵ model performs well in predicting the mean streamwise velocity distribution, while the k − ω SST model performs best in predicting the mean temperature distribution.However, except for k − ω SST, significant errors are obtained in the prediction of the mean temperature, with Error(Θ + ) becoming larger than 10% for the k − ϵ LS, standard k − ϵ, and LRR models, with the error of the SSG model being larger than 20%.The discrepancy between the temperature profiles predicted by the turbulence models and those obtained from DNS can also be seen in Figure 4, where a visible gap exists between DNS and the different turbulence models, except for the k − ω SST model prediction.

Reynolds Stress Prediction
The normalised Reynolds stresses also exhibited self-similar behaviour for the range of Re τ considered here.As shown in Figure 6, the Reynolds shear stress away from the wall collapses with different Re τ .As Re τ increases, the peak value shifts toward the wall, while the magnitudes of the peak Reynolds shear stress remains almost the same.A similar observation can be made for the turbulent kinetic energy in Figure 6.All five models show excellent performance in predicting the Reynolds shear stress and turbulent kinetic energy, as they all align with the DNS data.However, differences can be observed in the prediction of the Reynolds normal stresses in Figure 6.The RSM closures outperform the EVM closures in this regard.The EVM models assume an isotropic distribution of eddy viscosity.In contrast, the RSM closures explicitly solve for the anisotropic Reynolds stress tensor, yielding more-accurate results in certain scenarios, such as flows with complex geometries and strong pressure gradients.Both the LRR and SSG modelling framework include pressure-rate-of-strain effects, but the modelling approach for the pressure-strain correlation term Π ij is different.The SSG model takes into account Reynolds stress anisotropy and non-equilibrium effects, while LRR does not account for these effects.Consequently, the SSG model yields a Reynolds normal stress prediction that is slightly closer to DNS compared to the LRR model.However, in the current study, no major improvement in the prediction of the mean velocity and mean temperature is observed with the use of the RSM models, which is likely due to the canonical nature of the flow being studied.The DNS data for different values of Prandtl number are available only for Re τ = 180 and Re τ = 395, and thus, the RANS predictions of the mean temperature with different turbulence models are compared to the corresponding results obtained from DNS for these Re τ values at all the Prandtl numbers available in the DNS database [1].In Figure 7, the RANS model predictions of mean non-dimensional temperature Θ + = T/T τ are compared to the corresponding results extracted from the DNS data for Re τ = 180 over different values of Pr.Both k − ϵ LS and k − ω SST models exhibited a similar trend as the DNS data, but k − ω SST model shows better performance, especially when the Prandtl numbers are high.However, the turbulence models with wall functions (i.e., standard k − ϵ, LRR, and SSG) exhibits some unexpected results in the very-low-Prandtl-number region (Pr = 0.025, 0.05), where the models significantly over-predicts Θ + .For Pr ≥ 0.1, all the models relying on wall functions provide comparable satisfactory predictions.The temperature field exhibits distinctive characteristics across three regimes with respect to the values of the Prandtl number, which is discussed later in detail in Section 4.2.2.
In the analytical regime, denoted by Pr = 0.025 and Pr = 0.05 for the thermal field, molecular diffusivity solely influences the thermal behaviour, which is distinctly different from the higher values of Pr, where both the fluid flow and temperature fields are affected by turbulent diffusion.The application of turbulence models with standard thermal wall functions is rendered ineffective in this context, as they were derived based on the assumption of eddy thermal diffusivity dominating over the molecular thermal diffusivity.The RMS of the relative error of the mean non-dimensional temperature between DNS and turbulence model predictions for different values of Pr at Re τ = 180 and 395 are shown in Figure 8.The temperature profiles predicted by the standard k − ϵ, LRR, and SSG models have significant errors when the Prandtl number is below 0.1, which can also be seen in Figure 7 for Re τ = 180 and Re τ = 395.For Prandtl numbers greater than or equal to 0.1, the Error(Θ + ) profiles of the standard k − ϵ, LRR, and SSG models do not exhibit monotonicity with increasing Prandtl number.The Error(Θ + ) remains around 10% for both the standard k − ϵ and LRR models.The error of the SSG model is slightly higher than that of these two models, assuming around 20%, and its error variation with the Prandtl number is very similar to that of the LRR model.9a,b that all the models that rely on velocity wall functions (e.g., standard k − ϵ, LRR, and SSG models) significantly underpredict the Nusselt number for Pr < 0.1, but all the turbulence models provide comparable predictions of Nu, which are of the correct order of magnitude in comparison to the DNS results.The percentage deviations of the Nu predictions of the RANS models with respect to the DNS results (i.e., (Nu RANS − Nu DNS )/Nu DNS × 100%) are shown in Figure 9c,d for Re τ = 180 and 395, respectively.It can be seen from Figure 9c,d that the k − ϵ LS and k − ω SST models overpredict Nu for Pr = 0.025, 0.05, and 0.1, whereas all other models considered here significantly underpredict Nu for these Prandtl numbers.For Pr > 0.1, the extent of the underpredictions of Nu by the standard k − ϵ, LRR, and SSG models do not exhibit a monotonic trend with increasing Pr.The k − ϵ LS model shows significant underpredictions of Nu for Pr > 0.1, and this tendency strengthens with the increasing Prandtl number.The Nusselt number prediction by the k − ω SST model remains within ±5% of the DNS results for Pr > 0.1, which is consistent with the fact that the k − ω SST model provides satisfactory predictions of the mean streamwise velocity and temperature distributions (see Figures 5 and 8).

Different Regimes Based on the Ratio of Eddy Thermal Diffusivity to Molecular Thermal Diffusivity
The thermal field in a fully developed channel flow experiences a transition from molecular diffusivity dominance to turbulent diffusivity dominance as the Prandtl number increases from 0.025 to 10.0.Different types of thermal boundary layers are obtained during this transition, which can impact the temperature distribution and heat transfer rate within the channel.Consequently, the performance of models used to describe these flows is also influenced by these characteristics.
In a fully developed thermal flow field in a channel flow configuration, there are three layers: the diffusion sublayer, intermediate layer, and log-law layer.The existence of the intermediate layer depends on the values of the Prandtl number, with y + 1 being influenced by Pr, and y + 2 depends only on Pr t and κ, as defined in Section 3.3.With the decrease in the Prandtl number, the thermal flow field is gradually dominated by molecular diffusion effects.For extremely small values of the Prandtl numbers, the thermal field may only display molecular diffusion effects, and turbulence models with wall functions depending on the log-law may become invalid, which leads to significant errors in the prediction of the temperature fields.To validate this idea, an analytical solution is proposed.Assuming the thermal field is only influenced by molecular diffusivity, the passive scalar transport equation in a one-dimensional channel flow can be reduced as follows: By integrating the Ordinary Differential Equation given above (i.e., Equation ( 11)), an analytical solution can be obtained.The bulk velocity, denoted by ⟨u + ⟩, is equal to 15.736 at Re τ = 180.While there is no analytical solution for u + (y * ) in fully turbulent channel flow, approximate solutions can be obtained using polynomial fitting.This method suits different Re τ , but the velocity profiles have different polynomial expressions at different Re τ .An eighth-order spline is used to fit the mean velocity profile, resulting in the expression u + (y * ) = −94.7yThe analytical solution can be used to divide the thermal field into distinct regions based on the value of Pr.These regions include the analytical solution region (henceforth also called the low Pr region), transitional region, and turbulent diffusivity-dominant region.Within the turbulent diffusivity-dominant region, there are two sub-regions: the Pr ≈ 1.0 region and the high-Pr (i.e., Pr >> 1) region.The reason why these two sub-regions are defined is that the high-Pr region can account for variations in the thickness of the velocity and temperature boundary layers.Table 6 provides an overview of these three regions and their corresponding Pr values in the current configuration.Details along with the characteristics are discussed in Figure 10 for Re τ = 180, but the same practice can be followed for other values of Re τ ; however, the polynomials will be different.In Figure 10a, the DNS result is in perfect agreement with the analytical solution, which means the thermal field is only influenced by molecular diffusivity at Pr = 0.025.At Pr = 0.05 in Figure 10b, the DNS result shows a slight deviation from the analytical solution, indicating that turbulent diffusivity starts to play an influential role.In this regime, the thermal field is dominated by molecular diffusivity, rendering turbulence models with wall functions depending on the log-law invalid.This low-Pr regime can also be verified by using a three-layer thermal wall function and setting y + 1 = y + 2 .In this case, the turbulent Prandtl number is approximately equal to 0.71, and κ = 0.42.As a result, the critical Prandtl number is approximately 0.07, and there is no existence of an intermediate layer or log-law layer when the Prandtl number is less than 0.07.
The transitional regime can be observed at Pr = 0.1 and 0.2 in Figure 10c,d, respectively, where a diffusion sublayer and an intermediate layer exist.The DNS results aligned well with the analytical solution only in the diffusion sublayer, while deviations are seen in the intermediate layer, indicating that turbulent diffusivity effects became important.There is no log-law layer observed because the slope of the DNS data and the standard k − ϵ model deviate from the log-law layer auxiliary line (Θ + = 2.5ln(y + ) + C Θ ).However, as Pr increases, the mean temperature profiles of DNS and the standard k − ϵ in the log-law layer gradually approached the log-law layer's auxiliary line.At Pr = 0.4 and 0.6 in Figure 10e,f, respectively, three layers can be observed in the temperature profiles: a diffusion sublayer, an intermediate layer, and a log-law layer.This behaviour is characteristic of the regime where turbulent diffusivity plays a dominant role, and the effects of turbulent diffusivity strengthen with increasing Pr.Additionally, the intermediate layer becomes larger as Pr increases because the endpoint of the diffusion sublayer, y + 1 , decreases as Pr increases, while the starting point of the log-law layer is independent of Pr.

Conclusions
In this study, five different RANS modelling strategies are assessed for their ability to predict the heat transfer rate in a fully developed channel flow across a wide range of Re τ and Pr.Specifically, these modelling strategies include the standard k − ϵ, k − ϵ LS, and k − ω SST models, which are all categorised as Eddy Viscosity Models (EVMs), as well as Launder-Reece-Rodi (LRR) and Speziale-Sarkar-Gatski (SSG) models, which are classified as Reynolds Stress Models (RSMs).
Firstly, the predictions of the mean velocity and thermal fields obtained from these turbulence models are compared with the corresponding quantities extracted from the DNS data.Four setups with varying friction Reynolds numbers Re τ of 180, 395, 640, and 1020 at Pr = 0.71 are considered to analyse the Reynolds number dependence of the model performance.The mean velocity and temperature profiles obtained from both DNS and RANS for different turbulence models are found to exhibit self-similar behaviour.The Root-Mean-Squared (RMS) deviation from the DNS data is used to quantify the difference between the turbulence model predictions and the corresponding DNS data.All models perform well in the prediction of the mean streamwise velocity, with the RMS of the relative error of the velocity Error(u + ) less than 10% at all values of Re τ considered in this work.However, the performance of these models is less satisfactory for the prediction of the mean temperature, with the RMS of relative error of mean temperature Error(Θ + ) being larger than 10% for all models, except for the k − ω SST model.Additionally, the model predictions of the second moments are also compared with the corresponding quantities obtained from the DNS data.All models perform well in the prediction of the Reynolds shear stress and turbulent kinetic energy, but the RSM closures perform better than the EVM closures in terms of the predictions of normal Reynolds stresses.Among the RSM models, the SSG model outperforms the LRR model due to an improved representation of the pressure strain correlation term.Secondly, the influence of the Prandtl number on the prediction of the thermal field by different turbulence models in fully developed channel flow is considered.Nineteen setups with Prandtl numbers varying from 0.025 to 10.0 at Re τ = 180 and 395 are considered.The turbulence models with the wall function are found not to predict the mean temperature and Nusselt number adequately for small values of Pr.However, except for the very-low-Pr regime, these models performed comparably well for different choices of Pr.The performances of the k − ϵ LS and k − ω SST models show opposite trends in the prediction of the mean temperature and Nusselt number, with k − ϵ LS model being more accurate in the molecular diffusivity-dominated regime (i.e., small values of Pr) and k − ω SST model performing better in the turbulent diffusivity dominated regime (i.e., high values of Pr).Finally, an analytical solution is introduced to distinguish the Prandtl number dependence on thermal field predictions by the RANS models into three different regimes: the analytical regime, also called the low-Prandtl-number regime, the transitional regime, and the turbulent diffusivity-dominant regime.There are two sub-regions in the turbulent diffusivity regime: the Pr ≈ 1.0 regime and the large-Pr regime.These regimes can help understand the performance variation of a model with the variation in Pr.
Author Contributions: L.L.: conceptualisation, methodology, software, data generation, writingoriginal draft preparation and editing.U.A.: conceptualisation, methodology, supervision, software, testing, writing-original draft preparation, reviewing, and editing.N.C.: corresponding author, conceptualisation, methodology, supervision, writing-original draft preparation, reviewing, and editing.All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the Engineering and Physical Sciences Research Council (Grants: EP/V003534/1 and EP/R029369/1).

Figure 1 .
Figure 1.The configuration of turbulent heat transfer in heated channel flow.

Figure 2 .
Figure 2.An illustration of the law of the wall for the variation of mean streamwise velocity (a) and mean temperature (b) with different y + values chosen for turbulence models for low-Reynoldsnumber models and high-Reynolds-number models at Re τ = 180 and Pr = 0.71.

Figure 3 .
Figure 3. Variations of the normalised wall-normal coordinate of grid points (i.e., y * = y/h) with the number of the grid points in the wall-normal direction N y up to the channel half-height for turbulence models with and without the velocity wall function for different values of Re τ .Note that each symbol on the plots represents the cell centre location of the grid point.

Figure 4 .
Figure 4.The velocity and temperature profiles at Pr = 0.71 with different Re τ .

Figure 5 .
Figure 5. Error estimation of mean velocity (a) and temperature (b) of different turbulence models at Pr = 0.71 with different Re τ .

Figure 6 .
Figure 6.The profiles of the Reynolds shear stress, turbulent kinetic energy, Reynolds normal stress in the direction parallel to the wall, and normal stress in the wall-normal direction at Pr = 0.71 with different Re τ .4.2.Heat Transfer Rate in Turbulent Channel Flow for Different Values of Prandtl Number 4.2.1.Mean Temperature Profiles

Figure 7 .
Figure 7.A comparison of mean non-dimensional temperature profiles obtained from RANS models with DNS data at Re τ = 180 for different Prandtl numbers.

Figure 8 .
Figure 8. Error estimation of the mean temperature of different turbulence models with different Pr at Re τ = (a) 180 and (b) 395.For models without a wall function, both the k − ϵ LS and k − ω SST models provides satisfactory predictions for small values of Pr.The Error(Θ + ) of the k − ϵ LS model for Pr ≤ 0.1 is smaller than 10% at both Re τ = 180 and Re τ = 395.As Pr increases, the Error(Θ + ) of the k − ϵ LS model gradually increases and becomes larger than 20% for Pr ≥ 5. Therefore, the k − ϵ LS model is more accurate when the thermal field is dominated by molecular diffusivity and becomes less accurate when the thermal field is dominated by turbulent diffusivity at both Re τ = 180 and Re τ = 395.On the other hand, the k − ω SST model shows the opposite trend.It is more accurate when the thermal field is dominated by turbulent diffusivity and becomes less accurate when the thermal field is dominated by molecular diffusivity.The Error(Θ + ) of the k − ω SST model is the smallest among all the models considered, with an error lower than 10%, except for very small values of Pr.Further analysis regarding the model performance variation with Pr is provided in the following part.The variations of Nu = H • h/k = (Re τ Pr)/(Θ + w − Θ + b ) with Pr for Re τ = 180 and 395 are shown in Figure 9a,b, respectively, where H is the heat transfer coefficient, Θ + w is the nondimensional temperature at the wall, and Θ + b = 2h 0 u(y)T(y)dy/[T τ 2h 0 u(y)dy] is the nondimensional bulk temperature.The Nusselt numbers obtained from different turbulence models are also compared to the Nu values obtained from DNS data in Figure 9a,b.It can be seen from Figure 9a,b that all the models that rely on velocity wall functions (e.g., standard k − ϵ, LRR, and SSG models) significantly underpredict the Nusselt number for Pr < 0.1, but all the turbulence models provide comparable predictions of Nu, which are of the correct order of magnitude in comparison to the DNS results.The percentage deviations of the Nu predictions of the RANS models with respect to the DNS results (i.e., (Nu RANS − Nu DNS )/Nu DNS × 100%) are shown in Figure 9c,d for Re τ = 180 and

Figure 9 .
Figure 9. Variation of Nu with Pr for Re τ = (a) 180 and (b) 395.The percentage deviations of Nu predictions of RANS models with respect to DNS results (i.e., (Nu RANS − Nu DNS )/Nu DNS × 100%) for different values of Pr for Re τ = (a) 180 and (b) 395.

2 Figure 10 .
Figure 10.Three different regions in the thermal field concerning Pr. (a,b) are located in the analytical region; (c,d) are located in the transitional region; (e,f) are located in the turbulent diffusivitydominant region.

Table 3 .
Three-layer thermal wall function.

Table 4 .
An overall summary of the Reynolds-Averaged closure models that were selected for assessment.

Table 5 .
A comparison of the different turbulence models in terms of their computational time at Re τ = 1020 and Pr = 0.71.

Table 6 .
Three different regimes and the corresponding Pr.