Direct numerical simulation of a zero-pressure-gradient turbulent boundary layer with passive scalars up to Prandtl number Pr = 6

Abstract The objective of the present study is to provide a numerical database of thermal boundary layers and to contribute to the understanding of the dynamics of passive scalars at different Prandtl numbers. In this regard, a direct numerical simulation (DNS) of an incompressible zero-pressure-gradient turbulent boundary layer is performed with the Reynolds number based on momentum thickness ${\textit {Re}}_{\theta }$ ranging up to $1080$. Four passive scalars, characterized by the Prandtl numbers ${\textit {Pr}} = 1,2,4,6$ are simulated using the pseudo-spectral code SIMSON (Chevalier et al., SIMSON : a pseudo-spectral solver for incompressible boundary layer flows. Tech. Rep. TRITA-MEK 2007:07. KTH Mechanics, Stockholm, Sweden, 2007). To the best of our knowledge, the present DNS provides the thermal boundary layer with the highest Prandtl number available in the literature. It corresponds to that of water at ${\sim }24\,^{\circ }{\rm C}$, when the fluid temperature is considered as a passive scalar. Turbulence statistics for the flow and thermal fields are computed and compared with available numerical simulations at similar Reynolds numbers. The mean flow and scalar profiles, root-mean-squared velocity and scalar fluctuations, turbulent heat flux, turbulent Prandtl number and higher-order statistics agree well with the numerical data reported in the literature. Furthermore, the pre-multiplied two-dimensional spectra of the velocity and of the passive scalars are computed, providing a quantitative description of the energy distribution at the different length scales for various wall-normal locations. The energy distribution of the heat-flux fields at the wall is concentrated on longer temporal structures with increasing Prandtl number. This is due to the thinner thermal boundary layer as thermal diffusivity decreases and, thereby, the longer temporal structures exhibit a different footprint at the wall.


Introduction
Wall-bounded turbulence is a phenomenon of huge technological importance in many industrial and environmental applications.Understanding and predicting the behaviour of turbulent flows are critical in various fields, such as energy production, environmental modelling and fluid-dynamics research.However, investigating turbulent flow in complex geometries poses considerable challenges from both numerical and experimental perspectives.For this reason, simpler geometries are chosen when the fundamental physics of the flow is studied.One canonical flow case widely used in the literature is the boundary layer developing on a flat surface (Spalart 1988;Schlatter et al. 2009).The spatially evolving fully turbulent boundary layer has been studied using different experimental techniques (Österlund 1999;Rahgozar, Maciel & Schlatter 2013;Shehzad et al. 2021) with researchers constantly improving the measurement techniques (Örlü & Alfredsson 2010;Bailey et al. 2013;Vinuesa & Nagib 2016) to obtain reliable measurements at high Reynolds numbers (Samie et al. 2018).At the same time, direct numerical investigations of turbulent boundary layers have been performed in several studies (Spalart 1988;Ferrante & Elghobashi 2005;Wu & Moin 2009;Simens et al. 2009;Schlatter et al. 2010) implementing different solution methods for an increasing range of Reynolds numbers.Given the resolution limitations of experimental techniques in capturing the near-wall region of boundary-layer flow, direct numerical simulations (DNSs) have emerged as a valuable tool for investigating the transport phenomena in turbulence (Araya & Castillo 2012).However, the applicability of DNS is currently constrained to low-Reynolds-number flows due to its significant computational costs.Despite this limitation, DNS has proven instrumental in the study of turbulent boundary layers.Further, not limited to canonical flows, DNSs have been extended to more complex geometries such as airfoils (Vinuesa et al. 2017).These simulations have not only led to detailed insights into velocity profiles, turbulence statistics and the near-wall dynamics to name a few, but have also served as a validation to experiments and thereby they have contributed to refining turbulence models, establishing scaling laws and exploring turbulent coherent structures.
In practical engineering applications, the consideration of heat and mass transfer, turbulent mixing, combustion and other related phenomena is crucial for understanding and optimizing various systems (Kozuka, Seki & Kawamura 2009).In these scenarios, scalar quantities such as temperature play a significant role and need to be accurately simulated.Additionally, understanding and predicting the dynamics of passive scalars like air and water pollutants play an important role in local and global environmental problems (Kasagi & Iida 1999;Lazpita et al. 2022), as well as in the design of transport and energy systems (Straub et al. 2019).
Several experimental studies (Kays 1972;Perry & Hoffmann 1976;Simonich & Bradshaw 1978;Subramanian & Antonia 1981;Krishnamoorthy & Antonia 1987) have analysed different aspects of heat transfer of passive scalars in turbulent boundary layers.The investigation by Kays (1972) presented the variation of skin-friction coefficient and Stanton number (which characterizes the ratio of the heat transfer into the fluid to the thermal capacity of the fluid) in the boundary layer over a transpiring wall, with different blowing and suction conditions for a constant free-stream velocity condition.They have also discussed and proposed theoretical models to enable the prediction of heat transfer in a turbulent boundary layer.A zero-pressure-gradient turbulent boundary layer with constant wall-temperature conditions was set up by Perry & Hoffmann (1976), which enabled them to test similarity relations between instantaneous heat and momentum fluxes.Simonich & Bradshaw (1978) investigated the effects of free stream on heat transfer in a turbulent boundary layer and reported an increase in Stanton number with respect to free-stream turbulence.Subramanian & Antonia (1981) studied the effects of Reynolds number in a turbulent boundary layer and reported the Kármán and additive constants in the logarithmic law for velocity and temperature to be independent of Reynolds number.Further, Krishnamoorthy & Antonia (1987) were able to measure the three components of average temperature dissipation very close to the wall in a turbulent boundary layer, in the effort to model turbulence for the computation of temperature fields.Thereby, the heat-transfer behaviour has been a subject of continuous investigation from both engineering and numerical-modelling perspectives.
In numerical investigations, the fluid temperature can be considered as a passive scalar, provided that the buoyancy effects and the temperature dependence of fluid properties are considered as negligible (Chandrasekhar 1961;Monin & Yaglom 1971).Note that the passive scalar as simulated in the DNS is a diffusive contaminant in the fluid flow.Due to its low concentration, it does not have an influence on the fluid but is influenced by the fluid motion.Although the discussion of passive scalars is in the context of a thermal boundary layer, it could very well be considered as a pollutant concentration, in which case the Schmidt number would be the mass transfer analogous to the Prandtl number.Considering the passive scalar to represent temperature, many DNS studies of turbulent scalar transport have been performed to analyse the convective heat transfer between the fluid and solid walls in spatially developing flows.Bell & Ferziger (1993) first performed the DNS for a turbulent thermal boundary layer.Later, Kong, Choi & Lee (2000) performed a DNS at a Prandtl number of Pr = 0.71 with different boundary conditions including isothermal (Dirichlet) and isoflux (Neumann) conditions, for Re θ ranging between 300 and 420 (note that Re θ is the Reynolds number based on momentum thickness where the momentum thickness quantifies the loss in momentum due to the presence of a boundary layer).The Reynolds-number range was extended in the studies by Hattori, Houra & Nagano (2007), who simulated Re θ from 1000 to 1200, at Pr = 0.71.At the same time, the numerical investigations for Prandtl numbers up to Pr = 2 were performed by Tohdoh, Iwamoto & Kawamura (2008) for a relatively lower Reynolds-number range, up to Re θ = 420.The effect of different boundary conditions at Pr = 0.2, 0.71 and 2.0 in the Reynolds-number range of Re θ ∈ [350, 830] was reported by Li et al. (2009).The thermal channel-flow simulations have been conducted at higher Prandtl numbers of 49 and low Re by Schwertfirm & Manhart (2007) and at a Pr of 10 and high Reynolds number by Alcántara-Ávila & Hoyas (2021).However, the thermal turbulent boundary layers have been only partially explored at a medium Prandtl number of 2 by Li et al. (2009) owing to the significant computational cost associated with higher Pr.In this study, we consider higher Prandtl numbers in a turbulent thermal boundary layer, reporting analyses that are currently not available in the literature, according to the authors' knowledge.Thereby, the passive scalars at Pr = 1, 2, 4 and 6 are simulated for Re θ up to 1080 in a zero-pressure-gradient turbulent boundary layer using an isothermal wall boundary condition.
The details of the simulation set-up are provided in § 2. The statistics obtained from the simulation at different Prandtl numbers are compared with the data available in the literature for the fully developed thermal turbulent boundary layer in § 3. Since the thermal channel flow and thermal boundary layer exhibit similar behaviours in the near-wall region, the statistical quantities of the channel flow reported in the literature are also compared with the thermal boundary layer quantities at similar Reynolds number.In § 4 we analyse the distribution of energy in different scales for the wall-heat-flux field and wall-parallel fields at y + = 15, 30 and 50 (where the superscript '+' denotes scaling in terms of friction velocity u τ , see § 2.4).The premultiplied two-dimensional power-spectral density provides additional insight into the scalar transport at different Prandtl numbers.Finally, a short summary of the observations discussed in this work is reported in § 5.

Governing equations
A DNS of the zero-pressure-gradient (ZPG) turbulent boundary layer (TBL) is performed using the pseudo-spectral code SIMSON (Chevalier et al. 2007).The code solves the governing equations in non-dimensional form (here, written in index notation), in particular the flow and scalar variables are non-dimensionalized as where (x 1 , x 2 , x 3 ) = (x, y, z) are the Cartesian coordinates in the streamwise, wall-normal and spanwise directions, respectively, and t denotes the time.The length scale used for the non-dimensionalization is the displacement thickness at x = 0 and t = 0, denoted by δ * 0 .The corresponding instantaneous velocity components are denoted by with the mean quantities identified by ( U , V , W ) and the fluctuations by (u, v, w).Here U ∞ is the undisturbed laminar free-stream velocity at x = 0 and time t = 0.The total pressure is denoted by P and the density and kinematic viscosity of the fluid are represented by ρ and ν, respectively.In this study, four different passive scalars (Θ 1 , Θ 2 , Θ 3 , Θ 4 ) are simulated at different Prandtl numbers (Pr = 1, 2, 4, 6), respectively.Here, Θ i,∞ , Θ i,w correspond to the ith scalar concentration in the free stream and at the wall, respectively, with the mean quantities indicated by Θ i and the corresponding fluctuations by , 2, 3, 4].The superscript • introduced in (2.1a-e) identifies a non-dimensional variable and it shall be dropped in the non-dimensional quantities for the rest of the sections for simplicity.The non-dimensional form of the incompressible Navier-Stokes equation and the transport equation for passive scalars are given by where Re δ * 0 identifies the Reynolds number based on the free-stream velocity (U ∞ ) and the displacement thickness at the inlet (δ * 0 ).The product of the Reynolds number based on free-stream velocity and the Prandtl number results in another non-dimensional number, called the Péclet number (Pe = Re δ * 0 Pr), which measures the ratio between the scalar convective transport and the scalar molecular diffusion.Here, F j and F Θ i correspond to the volume force terms for the velocity and passive scalars, respectively.The velocity-vorticity formulation of the incompressible Navier-Stokes equation is implemented in the solver as the divergence-free condition is implicitly satisfied by the formulation.

Boundary conditions
Having defined the governing equations, the problem definition is completed by providing appropriate boundary conditions.At the wall, the velocity of the fluid is the same as that of the solid surface and is given by the following no-slip and no-penetration boundary conditions: U| y=0 = 0, V| y=0 = 0, W| y=0 = 0. (2.5a-c) From the continuity equation, we also obtain The flow is assumed to extend to an infinite distance perpendicular to the plate, but discretizing an infinite domain is not feasible.Hence, a finite domain has to be considered, for which artificial boundary conditions have to be applied.A simple Dirichlet condition can be considered; however, the desired flow solution generally contains a disturbance that would be forced to zero.This would introduce an error due to the increased damping of the disturbances in the boundary layer (Lundbladh et al. 1999).An improvement to the aforementioned boundary condition can be made by using the Neumann boundary condition given by ∂U j ∂y y=y L = ∂U j ∂y y=y L , (2.7) where y L is the height of the solution domain in the wall-normal direction in physical space and U i is the laminar base flow that is chosen as the Blasius flow for the present study.For the passive scalars a constant (isothermal) wall boundary condition is applied, as given by which corresponds to a vanishing thermal-activity ratio K.The thermal-activity ratio defines the ratio between the fluid density, thermal conductivity and specific heat capacity and the same properties of the boundary surface as defined below (2.9) Here, ρ w , k w and C p,w correspond to the density, thermal conductivity and specific heat capacity of the wall.The isothermal wall boundary condition corresponds to the fluid that exchanges heat with the boundary surface, without modifying the wall temperature.The boundary condition in the free stream is given by (2.10)

Numerical scheme
The DNS is performed with a pseudo-spectral method, where Fourier expansions are used in the streamwise and spanwise directions and Chebyshev polynomials T k (ξ ) (on −1 ≤ ξ ≤ 1) are used in the wall-normal direction employing the Chebyshev-tau method for faster convergence rates.The time advancement is performed using the second-order Crank-Nicholson scheme for linear terms and the third-order Runge-Kutta method for nonlinear terms, with a constant time step t.The maximum Courant number is set to 0.6.The nonlinear terms are calculated in physical space and the aliasing errors in the evaluation of these terms are removed by the 3/2 rule.
Since the TBL is developing in x, the periodic boundary condition cannot be directly used in this particular direction, which requires a specific numerical treatment.In this regard, one approach is to impose an appropriate instantaneous velocity and scalar profile at the inlet for every time step.Assuming self-similarity of the flow in the streamwise direction, Lund, Wu & Squires (1998) proposed a rescaling-recycling method to generate the required inlet profiles based on the solution downstream.An alternative approach is the addition of the fringe region downstream of the physical domain to retain the periodicity in the streamwise direction as described by Bertolotti, Herbert & Spalart (1992); Nordström, Nordin & Henningson (1999).In this method, the disturbances are damped, and the flow is forced from the outflow of the physical domain to the same profile as the inflow.The fringe technique is used in the present study, as the inflow conditions from a laminar profile followed by tripping produce natural instantaneous fluctuations for the velocity and the scalar fields (Araya & Castillo 2012).The fringe region is implemented by adding a volume force (F j ; F Θ i ) to the momentum and scalar transport equations ((2.3) and (2.4)), respectively.The forcing term is given by (2.11) where λ(x) is the strength of the forcing, which is non-zero only in the fringe region.The flow field at the inlet is the laminar Blasius profile U j and for the scalar Θ i it is the linear variation with y ranging from 0 to 1, denoted by Θi .

Computational domain and numerical set-up
The laminar base flow is tripped by a random volume force strip (at x/δ * 0 = 10) to trigger the transition of the flow to a turbulent state.For this simulation, a three-dimensional cuboid is considered with length, height and width equal to x L , y L , z L , respectively.The lower surface of the cuboid is considered as a flat plate with no-slip boundary conditions.The boundary layer grows in the considered computational domain with initial thickness denoted by δ * 0 .In the streamwise direction, the computational domain terminates with the fringe region.The vertical extent of the computational domain includes the whole boundary layer and the domain height is chosen based on the free-stream boundary condition.In particular, the problem set-up in this work is similar to that studied by Li et al. (2009).
The computational domain is discretized by N x , N y and N z grid points in the streamwise, wall-normal and spanwise directions, respectively.The grid spacing is uniform in the streamwise and spanwise directions.For the wall-normal direction, the collocation points follow the Gauss-Lobatto distribution given by (2.13) The computational box has dimensions 1000δ * 0 × 40δ * 0 × 50δ * 0 in the streamwise, wall-normal and spanwise directions, respectively.The number of grid points in each direction is 3200 × 385 × 320, correspondingly.Note that the smallest scale in the scalar fluctuations is inversely proportional to Pr 1/2 (Tennekes & Lumley 1972) and, hence, the Batchelor length scale η Θ i (which is analogous to the smallest scale in turbulent flow, the Kolmogorov scale η) is estimated as ηPr −1/2 (Kozuka et al. 2009).Similarly, the ratio of the largest to the smallest scales is proportional to Re 1.5 Pr 0.5 at high Pr (Batchelor 1959;Tennekes & Lumley 1972).In this study, an adequate grid resolution is adopted to resolve all the physically relevant scales.The Reynolds number based on free-stream velocity and displacement thickness at the inlet is Re δ * 0 = 450 and the friction Reynolds number based on local friction velocity (u τ ) and boundary-layer thickness (δ 99 ) is Re τ = 46.At the outlet, the Reynolds number based on displacement thickness is Re δ * 0 = 1, 580 and Re τ = 396.
Considering the friction velocity u τ at the middle of the computational domain (x/δ * 0 = 500, corresponding to Re θ = 794), the grid resolution in viscous units is x + = 6.6 and z + = 3.3 in the wall-parallel directions.In the wall-normal direction, we have an irregular distribution of collocation points, hence y + varies between 0.01 and 3.5.The results discussed at a particular streamwise position are inner scaled with respect to the local friction velocity, whereas, when the data along the streamwise direction are considered, the reference quantity at x/δ * 0 = 500 is considered for inner scaling.This indicates that if a particular wall-parallel plane (x-z plane) is sampled at a wall-normal distance (for e.g.y + = y s ), considering the friction velocity at x/δ * 0 = 500, the actual inner-scaled location varies along the streamwise direction, which is within ±0.1y s for the range of Reynolds numbers simulated in this work.
In this study, five different realizations of ZPG TBL are performed by introducing different trip forcings through modification of the random seed parameter, to obtain an ensemble average of the statistical quantities.All the different realizations are run for approximately 2400 time units (δ * 0 /U ∞ ) after one flow through of initial transience, which corresponds to 1000 time units (δ * 0 /U ∞ ).The converged statistics are obtained with the data corresponding to 12,000 time units (δ * 0 /U ∞ ), equivalent to 600 eddy-turnover times.

Integral quantities and non-dimensional numbers
The shape factor H 12 , which measures the ratio between the displacement thickness δ * and the momentum thickness θ, is plotted in figure 1.The shape factor is lower in the turbulent region with increasing Re θ and agrees well with the experimental and numerical data for Re θ > 600 where the TBL is fully developed.
Figure 2 depicts the streamwise variation of the skin-friction coefficient C f for the present simulation and indicates that the obtained result is in good agreement with the turbulent skin-friction solution provided by Schoenherr (1932), which is given by The computed skin-friction coefficient is also in good agreement with the correlation proposed by Smits, Matheson & Joubert (1983), which is given as The trip location is at x/δ * 0 = 10, with a strong peak in the skin-friction coefficient followed by the transition to turbulence before x/δ  provided by Erm & Joubert (1991) with wire tripping also closely correspond to the calculated turbulent skin-friction coefficient.It should be noted that Erm & Joubert (1991) also repeated the experiments with different tripping devices and found the influence of tripping to persist until Re θ ≈ 1500.Due to this, Jiménez et al. (2010) found the experimental data to be scattered for Re θ < 1070 and also showed the scatter to decrease with Re θ > 1600.
The Stanton number measures the convective heat transfer into the fluid with respect to the thermal capacity of the fluid.The spatial evolution of the Stanton numbers for different passive scalars scaled with the square root of Pr (as plotted in figure 3) is very similar to the skin-friction profiles plotted in figure 2. There is a difference observed in the laminar Stanton number at x = 0 with the present scaling.However, the individual Stanton-number profiles correspond well to the laminar solution (not shown in the figure) given by Kays & Crawford (1993), which is and to the turbulent solution obtained from the Reynolds-Colburn analogy as given by Lienhard & John (2005).For Θ 1 , due to the Reynolds analogy, the Stanton profile matches with the skin-friction profile scaled by a factor of two.For the passive scalars at higher Prandtl numbers, a generalization of the Reynolds-Colburn analogy can be obtained, as reported in the study by Lienhard ( 2020) with the values of a 1 = 1, a 2 = 12.8 and a 3 = 0.68 provided in Lienhard & John (2005) and C f being expressed as The Stanton-number plot for the passive scalar at Pr = 6 agrees well with the turbulent solution obtained from the Reynolds-Colburn analogy, as shown in figure 3. On the other hand, reference curves for Kays & Crawford and the Reynolds-Colburn analogy are only reported for Θ 4 for clarity.
Using the data for water with Pr = 5.9, Hollingsworth (1989) developed an interpolation for Prandtl numbers from 0.7 to 5.9 assuming the critical thickness of the sub-layer to be a simple function of Prandtl number.The empirical expression is given by which is plotted for the scalar at Pr = 6 in figure 3. We find that the relationship proposed by Hollingsworth (1989) underpredicts the Stanton number at Pr = 6.However, the interpolation relation provides a good match with the calculated Stanton number for Θ 1 , better than the Reynolds-Colburn analogy, which is not indicated in the plot for clarity.
Using the data obtained in the present simulation, a correlation between the Nusselt, Prandtl and Reynolds numbers can be obtained as which yields R 2 = 0.9985.This is similar to the correlation proposed by Kays & Crawford (1993) but for fully developed profiles in circular tubes and computed for gases Nu = 0.021Re 0.8 Pr 0.5 . (3.8)

Mean velocity and scalar profiles
The mean velocity profile obtained at the streamwise location corresponding to Re θ = 670 is shown in figure 4. The streamwise velocity profile is compared with the DNS data from Spalart (1988) at Re θ = 670 and Komminaho & Skote (2002) at Re θ = 666.The comparison of the present data with the existing DNS results shows a good agreement in the inner region.There is a slight deviation of the mean velocity profile reported by Spalart (1988) in the wake region with respect to the present data but it agrees well with the data provided by Komminaho & Skote (2002).
The mean profiles of the various passive scalars Θ i are normalized with the respective Prandtl numbers and represented in inner scaling by dividing with the friction scalar Θ i,τ defined as where C p is the heat capacity of the fluid and q i,w is the rate of heat transfer from the wall to the fluid and is defined by where k is the thermal conductivity of the fluid.The normalized inner-scaled mean scalar profiles are plotted at Re θ = 1070, corresponding to Re τ = 395, as shown in figure 5.
The mean scalar profiles follow the conductive sub-layer relation (Θ + i = Pry + ) and is clearly identified in the plots for y + < 5.The profiles of the passive scalars at Pr = 1, 2 are compared against the channel DNS data provided by Kozuka et al. (2009) at the same Prandtl numbers.The comparison shows a good agreement in the near-wall and overlap regions.It should be noted that the boundary condition used by Kozuka et al. (2009) is the uniform-heat-flux condition as opposed to the uniform scalar boundary condition applied in this study.Note that the study by Kawamura, Abe & Shingai (2000) showed that the Kármán constant k Θ as appearing in (3.11) (3.11) varies between the Dirichlet and Neumann boundary conditions based on low-Reynoldsnumber simulations.Here, A Θ i denotes the additive constant for scalar Θ i .However, based on higher-Reynolds-number simulations, Pirozzoli, Bernardini & Orlandi (2016) reported that the difference in boundary condition affects the mean passive scalar profiles only in small magnitudes in the logarithmic region, although the effect is evident in the scalar fluctuation profiles reported later.
The scalar profile at Pr = 4 is compared against the channel DNS data reported by Alcántara-Ávila et al. (2018).The channel DNS data reported by Kozuka et al. (2009) at Pr = 7 are used for comparison of the passive scalar at Pr = 6, since there are no simulations reported in the literature at exactly the same Prandtl number.This gives us the opportunity to highlight the difference between the profiles at these high Prandtl numbers.Due to the difference in the considered Prandtl numbers, we observe a discrepancy in the mean velocity profile for y + > 40 in the overlap region.Since the channel data are used for the comparison, there is a difference observed near the wake region for all the cases.However, a good agreement of the profiles is observed for the inner region.Based on experimental data, semi-empirical fits were provided by Kader (1981) for a boundary layer with constant heat flux.In this study it was assumed that the overlap layer exhibited logarithmic variation as given in (3.11), and an empirical relation was provided to determine the additive constant B Θ i .The comparison of the mean scalar profiles against the relationships provided by Kader (1981) shows that the value at the wake is slightly overestimated with respect to the DNS data.The deviation for scalar Θ 4 corresponding to Pr = 6 is approximately 3 % for y + ∈ [100, 500].One possible reason for this small deviation could also be the constant scalar boundary condition imposed in our simulations as opposed to the constant-heat-flux boundary conditions considered by Kader (1981).The comparison against Kader (1981) is not presented, instead, the mean scalar profiles are evaluated against the recent correlations proposed by Pirozzoli (2023).These correlations were developed by fitting the functional form of eddy thermal diffusivity to DNS data, resulting in predictive formulae for mean scalar profiles.Our comparison demonstrates good agreement across the different Prandtl numbers.For the passive scalar Θ 4 , the maximum observed deviation is approximately 2 %.
In the overlap region, the mean scalar profiles are nearly logarithmic and are characterized as (3.12) The mean scalar and velocity profiles are plotted in defect form in figure 6.We observe that the slope of the scalar profiles in defect form is nearly parallel indicating a constant Kármán constant for the scalar fields.The Kármán constant for the velocity is observed to be 0.41, consistent with the values reported in Spalart (1988).The diagnostic function for scalar is defined as is plotted in figure 7. The Kármán constant for a scalar field is obtained as k Θ = 0.4, which is in between k Θ = 0.33 (Wikström 1998) and k Θ = 0.47 (Kader 1981).The obtained k Θ from the present simulation is the same as observed by Kawamura, Abe & Matsuo (1999) and close to the value of 0.41 as reported by Li et al. (2009).In the outer region, we observe a collapse of the profiles.

Velocity and scalar fluctuations
As shown in figure 8, the present velocity-fluctuation root-mean-squared (r.m.s.) data show a trend similar to that of the results by Jiménez et al. (2010).The r.m.s.profiles of the three velocity components are in good agreement in the inner region, while a minor difference can be observed in the outer region.Nonetheless, the peaks of the velocity fluctuations match in both position and magnitude.There is a slight offset in the plots of p rms , which can be attributed to the small difference in the considered Re θ for comparison.The r.m.s. of the velocity components are also compared with the channel DNS data provided by Abe et al. (2004).The streamwise r.m.s.agrees well with the present DNS results and the near-wall peak value coincides with the present observations.As expected, there is a difference observed in the outer region of flow, since channel and boundary-layer flows are fundamentally different farther from the wall.Additionally, the r.m.s. of the pressure fluctuations observed in the boundary layer is different compared with the channel flow.From the present DNS data, the peak of the streamwise velocity fluctuation is found at y + = 14, corresponding in outer units to y/δ 99 = 0.035.
The r.m.s. of the scalars at different Prandtl numbers are plotted in figure 9.The scalar r.m.s.profile at Pr = 1 is similar to the streamwise velocity r.m.s. and has a higher (roughly 5 %) near-wall peak comparatively, as expected.The comparison of scalar-fluctuation with R 2 = 0.99, where R 2 corresponds to the coefficient of determination.It is defined as  (2016) have suggested that the attached-eddy arguments support the increase of inner peak of streamwise velocity r.m.s. with respect to Re τ due to the effect of overlying attached eddies and they assumed that the same argument applies to the passive scalars.On the other hand, we find that the inner peak of the passive scalar r.m.s. at high Pr does not follow the previous argumentation.It should also be noted that the range of Re τ in our present simulation is narrow compared with those of the works by Pirozzoli et al. (2016) and Alcántara-Ávila et al. (2021) and a more detailed investigation of this topic is necessary to make any conclusive statements.
The above procedure was also performed for the streamwise heat flux, as shown in figure 12.A similar behaviour was observed with the overall variation in streamwise heat with R 2 = 0.9913.

Turbulent Prandtl number
An important parameter for scalar transport is the turbulent Prandtl number Pr t , which is defined as the ratio between turbulent eddy viscosity and turbulent eddy diffusivity The eddy viscosity and the eddy diffusivity arise from the Boussinesq hypothesis (Boussinesq 1877) for modelling turbulent stresses and the heat-flux vector, respectively.For parallel flows (the ones in which the velocity profile does not vary in the streamwise direction), the turbulent eddy viscosity and the turbulent eddy diffusivity are used to describe the turbulent momentum transfer and heat transfer with respect to the mean-flow conditions, in particular the mean velocity strain and temperature gradients, respectively.They are defined as and It should be noted that the eddy viscosity or diffusivity does not represent a physical property of the fluid, like the molecular viscosity, but rather a property of the flow.The Reynolds analogy introduces the similarity between the turbulent momentum exchange and turbulent heat transfer in a fluid.Reynolds (1874) noted that, for a fully turbulent field, both the momentum and heat are transferred due to the motion of turbulent eddies.This yields a simpler model for the turbulent Prandtl number, where the turbulent eddy viscosity for the momentum exchange and turbulent eddy diffusivity for the scalar transport are equal, such that Pr t = 1.Substitution of the eddy viscosity and diffusivity into the definition of turbulent Prandtl number results in From this definition, evaluating the turbulent Prandtl number at any point in the boundary layer would require the turbulent shear stress, turbulent heat transfer, velocity gradient and temperature gradient.Experimental investigations have limited accuracy in the simultaneous measurement of the Reynolds shear stress and wall-normal turbulent heat flux, in particular close to the wall (Araya & Castillo 2012).For this reason, experimental investigations like Blom (1970), Antonia (1980) and Kays & Crawford (1993) exhibit significant scatter in the data.
The variation of turbulent Prandtl number with the wall-normal distance in inner units at a given Re θ = 1070 is reported in figure 13.We observe that the turbulent Prandtl number varies at the wall and increases with respect to molecular Prandtl number of the scalar.From figure 14, we also see the turbulent Prandtl number decays for y + > 15 and this decay becomes steeper as Re decreases.The turbulent Prandtl number was assumed to approach a constant value close to the wall independent of the molecular Prandtl number (Li 2007).Studies such as the ones by Kestin & Richardson (1963) and Blom (1970) have analysed experimentally the turbulent Prandtl number in order to assess the validity of this assumption.Such experimental investigations have reported a constant value of turbulent Prandtl number as the wall is approached.However, different experimental campaigns provided data that could not show a conclusive interpretation of the behaviour of turbulent Prandtl number close to the wall.Nonetheless, many authors have proposed different correlations based on the experimental data to predict the heat-transfer coefficient through the prediction of turbulent Prandtl number with respect to molecular Prandtl number.
For the different passive scalars considered in the present study, the turbulent Prandtl number approaches a constant value > 1 close to the wall in the viscous sub-layer.The plots of Pr t exhibit a significant difference closer to the wall with respect to various molecular Prandtl numbers.There is a slight decrease in Pr t up to y + ≈ 20 and then the increase is maintained farther from the wall up to y + ≈ 50, the point after which the trend steadily decreases.It is pointed out in Kays (1994) that the peak between y + ≈ 20 and 100 is not observed in the experimental data, an observation which was attributed to the high Reynolds number of the experiments, where DNS data are not available for comparison.
Following the discussions about a constant Pr t in the logarithmic region, Kays & Crawford (1993) proposed a correlation for the turbulent Prandtl number that is applicable to air.In this correlation, the value of Pr t approaches a constant value of 0.85 in the logarithmic region.In the studies by Hollingsworth (1989), a correlation was proposed based on the measurement of the temperature profile of water at Pr ≈ 6.Again, the value of Pr t approaches a value of 0.85 as y + is increased.This observation of constant Pr t in the logarithmic region is not clearly observed with the present DNS data.This can be due to the low-Reynolds-number range considered in the present study.Based on the correlation proposed by Hollingsworth (1989), Kays (1994) suggested a constant value of Pr t = 1.07 for 0 < y + < 5. Indeed, if we consider only the passive scalars at Pr = 1, 2, the turbulent Prandtl number approaches a constant value of 1.07 closer to the wall.It appears that the turbulent Prandtl number close to the wall is independent of the molecular Prandtl number, as shown in the studies by Li et al. (2009).This independence of the turbulent Prandtl number at the wall with respect to the molecular Prandtl number has also been reported in many other studies like Kong et al. (2000) and Jacobs & Durbin (2000) for TBL flow, as well as in Kim & Moin (1989) and Kasagi, Tomita & Kuroda (1992) for turbulent channel flow.However, from the present study, we find that the turbulent Prandtl number indeed depends on the molecular Prandtl number and this observation is based on the increasing value of Pr t at the wall with respect to the scalars with Pr = 4, 6, as shown in figure 13.
In order to verify the plausibility of present observations at Pr = 4, 6, the obtained DNS data were compared with the DNS channel data reported by Alcántara-Ávila & Hoyas (2021).Figure 15 shows the comparison of Pr t at different molecular Prandtl numbers, where the present DNS data were at Re θ = 1070, corresponding to Re τ = 395, and the data from Alcántara-Ávila & Hoyas (2021) were at Re τ = 500.Despite these differences, the turbulent Prandtl numbers close to the wall are in good agreement, confirming that the Pr-scaled wall-normal heat flux decreases with an increase in Pr for Pr 4, as stated by Alcántara-Ávila & Hoyas (2021).Thus, the present observation confirms the constant behaviour of the turbulent Prandtl number very close to the wall and highlights its dependence on the molecular Prandtl number, which has been often ignored in turbulent heat-transfer calculations.
A brief discussion of the Reynolds-stress budget is provided in Appendix A.
3.5.Higher-order statistics, shear stress and heat flux The higher-order statistics (specifically, the values of the third-and fourth-order moments of a quantity) provide information on the non-Gaussian behaviour of turbulence.The thirdand fourth-order moments are also called the skewness and flatness, respectively, and for a statistically stationary variable m, they are defined as The skewness and flatness of the streamwise and wall-normal velocity components are shown in figure 16 where, for the spanwise velocity, we observe a Gaussian behaviour in the overlap region, as expected.The higher-order moments obtained by Vreman & Kuerten (2014) for the turbulent channel flow are compared with the present TBL results at Re τ = 180 (corresponding to Re θ = 420), showing a reasonable agreement.From figure 2, it should be highlighted that Re θ = 420 is roughly at the beginning of the turbulent regime and the deviation of the profiles for y + > 100 as observed in figure 16 is due to the intermittent wake region which decays to the values for a Gaussian distribution as the free stream is approached.There are some differences observed in the overlap region where the skewness of u is roughly 10 % higher than that reported by Vreman & Kuerten (2014).Although the plots indicate an overall good agreement, some drastic differences are observed in the flatness of the wall-normal and spanwise velocity components close to the wall, where the turbulence is highly intermittent.The flatness of the wall-normal velocity component approaches a value of ≈ 29 close to the wall for the data provided by Vreman & Kuerten (2014) whereas, in the TBL, it converges to ≈ 20.This is still lower than the value of 22 obtained by Kim, Moin & Moser (1987).The skewness of the pressure fluctuations is roughly 20 % higher in the overlap region of the TBL compared with the channel.Further, the intermittency (flatness) of the pressure is higher than the velocity components as reported by Kim et al. (1987), whereas the data obtained with TBL show an offset of 15 % throughout and are lower than the flatness behaviour observed in the channel.The flatness factor for pressure fluctuation at the wall approaches a value of 4.5 in the present simulations as compared with the values of 4.7 (Li et al. 2009) and 4.9 (Schewe 1983), whereas a value of ≈5.2 is reported by Vreman & Kuerten (2014).The Reynolds shear stress is plotted in figure 17 and it is compared with the TBL data provided by Jiménez et al. (2010) and the channel DNS data by Kozuka et al. (2009).Even if there is a good agreement between the data in the inner region, there is a clear difference in the overlap region between the channel and the TBL.Such differences were also observed in the plots of the fluctuations in the streamwise and wall-normal directions shown in figure 8.The peak of the Reynolds shear stress is higher for the TBL compared with a channel flow, which indicates a higher momentum transfer by the fluctuating velocity field in the TBL.Looking at the energy budgets for shear stress component in both the channel and TBL corresponding to Re τ = 395 (plots not shown here), we find that there is higher production in the TBL compared with channel flows  (with uniform-heat-flux wall conditions) in the overlap region whereas the dissipation was of similar magnitude.The turbulent streamwise heat flux shows a reasonable agreement with the channel data available at Pr = 1, 2, as shown in figure 18, although the peaks for the channel-flow data are slightly higher and closer to the wall compared with the present observations.For Pr = 4, however, the peak of the streamwise heat flux is higher in the channel data by Alcántara-Ávila et al. (2018), since the heat flux is reported at a higher Reynolds number than our DNS.Furthermore, the comparison at the highest Pr in our simulation was at a slightly lower Pr than that in the channel and hence there is a difference in the slope of the streamwise heat flux close to the wall.Figure 19 shows the wall-normal heat flux for the different simulations.A clear difference can be observed in the outer region when comparing similar Prandtl numbers.On the other hand, the boundary-layer data and the channel data show a better agreement in the inner region.
The correlation coefficients provide more information on the statistical association between the fields, and here the structure of the flow field and scalar fluctuations are analysed in terms of The correlation-coefficient plot in figure 20 corresponding to u − θ i shows a strong correlation of streamwise velocity and scalar fluctuations at Pr = 1 and it decreases with increase in Prandtl number.This is related to the similarity in the momentum and passive scalar transport by the turbulent eddies close to the wall.On the other hand, the v − θ i correlation coefficient also shown in figure 20 exhibits an increasing trend with the Prandtl number.In the conductive sub-layer the correlation coefficients coincide for the various Prandtl numbers under study and then approach different values at the wall.This highlights a similar behaviour of the turbulent wall-normal momentum and passive scalar fields with a caveat, i.e. the differences are present very close to the wall for different Prandtl numbers.Li et al. (2009) for (×, black) U, ( , black) Θ 2 .Inset: spanwise two-point correlation for (--, blue) Θ 1 , (--, orange) Θ 2 , (--, green) Θ 3 , (--, red) Θ 4 at Re θ ≈ 830.

Spanwise two-point correlations
The two-point correlations provide some quantitative information on the turbulent structures near the wall.For example, the streak spacing near the wall is of interest and can also be observed in an experimental setting.In order to identify the mean streak spacing, spanwise two-point correlations of the velocity components and passive scalars were obtained at five different positions along the streamwise direction at different wall-normal positions.Overall, the obtained results at Re θ = 830 were compared with the data reported by Kim et al. (1987) and show good agreement (not depicted here).The obtained two-point correlations for different scalars at Re θ = 830 are shown in figure 21 to assess the differences for varying Pr.The two-point correlation becomes negative and reaches a minimum at an inner-scaled correlation length of δz + ≈ 50.The length at which the minimum occurs provides an estimate of the half-mean separation between the streaks in the spanwise direction, i.e. (λ + s /2).The streak spacing is plotted along the streamwise direction at a wall-normal position of y + ≈ 7, as shown in figure 21.Overall, the velocity and scalar streak spacings increase with Re θ , as reported in the works of Li et al. (2009) and Schlatter et al. (2009).Note that the correlations are available only at five streamwise locations (solid line for U in figure 21 is for clarity), but the comparison with Li et al. (2009) is in reasonably good agreement for the velocity and scalar streaks at Pr = 2.The velocity streak spacing increases from 102 to 115 and appears to saturate for Re θ > 830.Note that such a saturation of streak spacing was also pointed out by Li et al. (2009) for Re θ > 1500.From figure 21, we also observe that the streak spacing decreases with increasing Pr and that the streak spacings for the scalars at Pr = 4, 6 are indistinguishable, although the rate of decay of the two-point correlations for the scalars is different.A higher grid resolution might be necessary to quantify the possible differences in the scalar streak spacing at higher Pr.

Scaling of wall-heat-flux fields
The wall-shear and heat-flux fields at different Prandtl numbers are normalized by subtracting the mean and dividing by the corresponding r.m.s.quantities.Although the distribution of data is different in the various fields, after normalization the distributions become identical.This indicates the uniformity in the distribution of the fluctuations of shear and heat flux when scaled with the corresponding r.m.s.quantities.This observation is useful for certain applications, for instance, in the prediction of fluctuating flow quantities from the wall, as discussed in the study by Guastoni et al. (2022).

Spectral analysis
The analysis of thermal boundary-layer statistics reveals that the scalar fluctuations and heat flux are strongly affected by the Prandtl number of the scalar in the flow and the corresponding scalings were reported in § 3. Additional insight can be obtained by analysing the energy distribution at the different length scales for the scalars at the simulated Prandtl numbers.In this regard, time series of the wall-shear, wall heat flux, streamwise velocity and different scalars were sampled at different wall-normal locations with a sampling time (∼ t + s = 1 with the reference friction velocity at x/δ * 0 = 500) 974 A49-26 corresponding to 12 000 time units (δ * 0 /U ∞ ).The two-dimensional (2-D) premultiplied power-spectral density (PSD) k z k t φ(λ + z , λ + t ) is obtained in the spanwise direction and in time based on the time series, for a total sampled time of about 11.5 flow-through times.Note that, here, k z , k t denote the spanwise, temporal wavenumber, respectively, and φ is the PSD defined for the particular quantity under study.
For the calculation of PSD, the procedure outlined in Pozuelo et al. (2022) is utilized.After the mean subtraction of the sampled time series to obtain the turbulent quantities, first a 1-D spectrum E[t, x](λ + z ) is obtained by performing a fast Fourier transform (FFT) in the spanwise direction, due to the periodicity condition imposed along z.As a result, a spectral decomposition of the energy content into different wavenumbers k z is obtained with the corresponding wavelengths λ z given by 2π/k z .Note that the local friction velocity is used to obtain the inner-scaled quantities.It should be pointed out that the flow is developing along the streamwise direction and hence the FFT along the streamwise direction is not applicable.In the present spectral analysis, we consider the Re θ range between 470 and 1070 and use Welch's overlapping-window method to address the non-periodicity in the temporal As a next step, the spectrum in time is obtained using Welch's method with 15 bins in total, where 8 of them are independent.A Hamming window is used for imposing the periodicity in the bins, and the 2-D spectrum is obtained as E[x](λ + z , λ + t ) by using FFT along z and Welch's method in t, with λ t denoting the temporal wavelength (i.e. the period), which is defined as λ t = 2π/k t .The obtained spectrum is divided by k t k z and premultiplied with k t k z .Finally, an averaging along x is performed to yield the 2-D premultiplied PSD k t k z φ(λ + z , λ + t ).The 2-D PSDs of the streamwise wall shear and wall heat flux at different Prandtl numbers are provided in figure 25, along with the PSD of the streamwise velocity and scalar fluctuations at y + = 15.The obtained 2-D PSD at y + = 15 agrees well with the results reported by Pozuelo et al. (2022), although the latter are at higher Reynolds number.At the wall, the spectral peak is observed at λ + z ≈ 100 and λ + t ≈ 100.Furthermore, at y + = 15, we observe that the maximum of spectral-energy distribution is at λ + z ≈ 120, which corresponds to the characteristic streak spacing in wall turbulence (Smith & Metzler 1983).It is observed that the PSD for the streamwise wall-shear stress is very similar to that for the wall heat flux at Pr = 1, which is an expected result for the reasons outlined in § 3.1.However, with increasing Pr, the PSD shifts to the right, indicating that the energy is not spread over a wider range of scales, and instead is concentrated in longer temporal structures.One could argue that, with a shorter boundary layer, the structures at Pr = 6 can become larger than the ones developing at Pr = 1.Because of this, the larger structures have a different footprint at the wall.Additionally, the plots in figure 25 also exhibit a slight trend downwards for higher Pr, a fact that indicates the presence of smaller spanwise scales, in agreement with the discussion in § 3.6.Overall, the temporal wavelength range at which we have the most energetic structures in the wall-heat-flux fields increases for larger Prandtl numbers.From the above observations, considering the dominant energetic structures to be composed of streaks at the wall, the scalar at Pr = 6 (in general for higher Pr) might exhibit longer and thinner scalar streak structures at the wall compared with the case at Pr = 1.
The PSDs calculated at y + = 30 and 50 are shown in figure 26.From this figure, we observe that the similarity in the distribution of energy for the scalar at Pr = 1 and the streamwise velocity is lost as we move farther from the wall.In contrast to the observation of large streamwise structures at the wall, we find an increasing concentration of energy in smaller scales as we increase the Prandtl number.Further, the range of scales in which  the energy is distributed also increases with Pr.Focusing on the most energetic structures, we find that these are concentrated in a region of smaller temporal and spanwise scales with increasing Pr at both y + = 30 and 50.For the scalar at Pr = 6, the spectral peak is observed at λ + z ≈ 100 and λ + t ≈ 20 for y + = 30 and 50.974 A49-28

Summary and conclusions
In the present study, a DNS of the thermal TBL is performed up to Re θ = 1080 with passive scalars at Pr = 1, 2, 4 and 6, which to authors' knowledge is the highest Prandtl number simulated for the thermal boundary layer.Various statistical quantities for the flow and scalars were computed and compared against the reported channel and TBL data in the literature.Overall, the statistical quantities are in good agreement with the existing data whenever a comparison is possible.For higher Pr, we also observed that the peak of the scalar fluctuations decreases when Re τ increases, which is different from the trend reported in the literature.Further, we showed that the variation of the peak in scalar fluctuation has a weak dependence with Re τ compared with the Prandtl number.Similarly, the peak in the heat flux also exhibits a weak dependence with Re τ compared with Pr of the scalar and the heat flux scales as ∼Pr 0.4 .
In the present study, we also highlighted the behaviour of the turbulent Prandtl number Pr t , which does not approach a constant value of 1.07 as the wall is approached for higher Prandtl numbers.In addition, we also found the corresponding Pr t to increase with Pr, confirming the findings of Alcántara-Ávila & Hoyas (2021).Finally, a brief description of the energy distribution in the scales for different Pr at different wall-normal locations is presented by analysing the 2-D pre-multiplied PSD.
The analysis and data provided in this work are expected to serve as a database for the research community to assess the validity of new turbulence models and validate other numerical and experimental results.The scalar-flux equation is written (in index notation) as where P Θ i j is the production term due to both the mean gradients of the velocity and scalar Θ i , Θ i j is the dissipation term, C Θ i j is the scalar-pressure-gradient correlation term (which can be split into a pressure scalar-gradient correlation term and a divergence of the pressure-scalar correlation term), T Θ i j is the turbulent diffusion and D Θ i j is the molecular diffusion term.The terms are defined as T Θ i j := − ∂ ∂x l u j u l θ i , (A11) A detailed description of the scalar-flux budget terms can be found in Kozuka et al. (2009) and Li et al. (2009).Due to the lack of TBL data at higher Prandtl numbers, the budgets for the scalar at Pr = 2 are compared against the channel DNS data from Kozuka et al. (2009), as shown in figure 28.Overall, there is a good comparison obtained for the different terms in the scalar-flux budgets.Also, as discussed in § 3.5, we observe a higher production compared with the channel-flow case for the vertical passive scalar-flux budget in the overlap region.In addition, the scalar-pressure diffusion term also exhibits the same behaviour as discussed above and it should be noted that the discrepancy not only stems from the different problem set-up but also the wall boundary condition for the scalar, which is Dirichlet in the present study and Neumann in the works of Kozuka et al. (2009).Further, a comparison of the data obtained at Pr = 6 with the data obtained by Kozuka et al. (2009) at Pr = 7 for a channel flow is also provided in figure 29.Overall, the comparison of the stress budgets at different parameter points shows a good agreement with the data available in the literature.Note that the small discrepancy observed in figure 29 is due to the different Prandtl numbers.

Figure 22 .
Figure 22.Instantaneous normalized streamwise wall-shear and heat-flux fields at different Prandtl numbers.
Pirozzoli et al. (2016)orresponding to the fitted value of x t .The observed correlation shows a weak dependence on Re compared with Pr.Further, from (3.14), a decaying trend of θ + i rms,max with Re is obtained, although such a trend has not been observed in the literature except for Θ 3 , Θ 4 in the present study.The studies byPirozzoli et al. (2016)and Alcántara-Ávila et al. (2021) have reported the increasing trend of θ + i rms,max with respect to Re but for a lower Pr than in the present work.Pirozzoli et al.