1. Introduction
The study of passive scalars evolving within wall-bounded turbulent flows has great practical importance, being relevant for the behaviour of diluted contaminants and/or as a model for the temperature field under the assumption of low Mach number and small temperature differences (Monin & Yaglom Reference Monin and Yaglom1971; Cebeci & Bradshaw Reference Cebeci and Bradshaw1984). It is well known that measurements of concentration of passive tracers and of small temperature differences are quite difficult and, in fact, available information about even basic passive scalar statistics are rather limited (Gowen & Smith Reference Gowen and Smith1967; Kader Reference Kader1981; Subramanian & Antonia Reference Subramanian and Antonia1981; Nagano & Tagawa Reference Nagano and Tagawa1988), mostly including basic mean properties and overall mass or heat transfer coefficients. Although existing semi-empirical correlations have sufficient accuracy for engineering design (Kays, Crawford & Weigand Reference Kays, Crawford and Weigand1980), their theoretical foundations are not firmly established. Furthermore, assumptions typically made in turbulence models such as constant turbulent Prandtl number are known to be crude approximations in the absence of reliable reference data.
Given this scenario, direct numerical simulation (DNS) is the natural candidate to establish a credible database for the physical analysis of passive scalars in wall turbulence, and for the development and validation of phenomenological prediction formulae and turbulence models. Most DNS studies of passive scalars in wall turbulence have been so far carried out for the prototype case of planar channel flow, starting with the work of Kim & Moin (Reference Kim and Moin1989), at 
$Re_{\tau } = 180$ (here 
$Re_{\tau } = u_{\tau } R / \nu$ is the friction Reynolds number, with 
$u_{\tau } = (\tau _w/\rho )^{1/2}$ the friction velocity, 
$R$ the pipe radius, 
$\nu$ the fluid kinematic viscosity, 
$\rho$ the fluid density and 
$\tau _w$ the wall shear stress), in which the forcing of the scalar field was achieved using a spatially and temporally uniform source term. Additional simulations at increasingly high Reynolds number were carried out by Kawamura, Abe & Matsuo (Reference Kawamura, Abe and Matsuo1999) and Abe, Kawamura & Matsuo (Reference Abe, Kawamura and Matsuo2004), based on enforcement of strictly constant heat flux in time (this approach is hereafter referred to as CHF), which first allowed scale separation effects to be appreciated and a reasonable value of the scalar von Kármán constant 
$k_{\theta } \approx 0.43$ to be deduced, as well as effects of Prandtl number variation (the molecular Prandtl number is here defined as the ratio of the kinematic viscosity to the scalar diffusivity, 
$Pr=\nu / \alpha$). Those studies showed close similarity between the streamwise velocity and passive scalar field in the near-wall region, as after the classical Reynolds analogy. Specifically, the scalar field was found to be organised into streaks whose size scales in wall units, with a correlation coefficient between streamwise velocity fluctuations and scalar fluctuations close to unit. Computationally high Reynolds numbers (
$Re_{\tau } \approx 4000$, with 
$Pr \le 1$) were reached in the study of Pirozzoli, Bernardini & Orlandi (Reference Pirozzoli, Bernardini and Orlandi2016), using spatially uniform forcing in such a way as to maintain the bulk temperature constant in time (this approach is hereafter referred to as CMT). Recent large-scale channel flow DNS with passive scalars using the CHF forcing at 
$Pr=0.71$ (as representative of air) have been carried out by Alcántara-Ávila, Hoyas & Pérez-Quiles (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021).
Flow in a circular pipe is clearly more practically relevant than planar channel flow in view of applications such as heat exchangers, and it has been the subject of a number of experimental studies, mainly aimed at predicting the heat transfer coefficient as a function of the bulk flow Reynolds number (Kays et al. Reference Kays, Crawford and Weigand1980). High-fidelity numerical simulations including passive scalars in pipe flow have been quite scarce so far and mainly include studies at 
$Re_{\tau } \le 1000$ (Piller Reference Piller2005; Redjem-Saad, Ould-Rouiss & Lauriat Reference Redjem-Saad, Ould-Rouiss and Lauriat2007; Saha et al. Reference Saha, Chin, Blackburn and Ooi2011; Antoranz et al. Reference Antoranz, Gonzalo, Flores and Garcia-Villalba2015; Straub et al. Reference Straub, Forooghi, Marocco, Wetzel and Frohnapfel2019). The current knowledge about gross properties and mean scalar profiles across the Reynolds and Prandtl numbers envelope thus currently rests on semi-empirical correlations (Kays et al. Reference Kays, Crawford and Weigand1980; Kader Reference Kader1981), which, although suitable for practical applications, deserve careful scrutiny. Clearly, the state of the art is not as well developed as for planar channel flow.
In this paper, we thus present DNS data of turbulent flow in a smooth circular pipe up to 
$Re_{\tau } \approx 6000$, including a passive scalar field with 
$Pr = 1$, at which some asymptotic high-Reynolds-number effects appear which were not observed in previous studies. Relying on the DNS database, we carry out an analysis of the structure of passive scalars in turbulent pipe flow, revisit current theoretical inferences and discuss implications about possible trends in the extreme-Reynolds-number regime. From a more engineering standpoint, we also revisit formulae for heat transfer prediction, as well as assumptions made in Reynolds-averaged Navier–Stokes (RANS) models for passive scalar turbulence. Although, as previously pointed out, the study of passive scalars is relevant in several contexts, the main field of application is heat transfer and, therefore, from now on we refer to the passive scalar field as the temperature field (denoted as 
$T$), and scalar fluxes will be interpreted as heat fluxes. Details on the velocity statistics from the present DNS database were reported in a separate publication (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021).
2. Numerical dataset
Numerical simulations of fully developed turbulent flow in a circular pipe are carried out assuming periodic boundary conditions in the axial (
$z$) and azimuthal (
$\phi$) directions, as shown in figure 1. The velocity field is controlled by two parameters, namely the bulk Reynolds number (
$Re_b = 2 R u_b / \nu$, with 
$u_b$ the bulk velocity namely averaged over the cross section), and the relative pipe length, 
$L_z/R$. The incompressible Navier–Stokes equations are supplemented with the transport equation for a passive scalar field (hence, buoyancy effects are disregarded), with the same diffusivity as the velocity field (hence, we assume 
$Pr=1$) and with isothermal boundary conditions at the pipe wall (
$r=R$). The passive scalar equation is forced through a time-varying, spatially uniform source term (CMT approach), in the interest of achieving complete similarity with the streamwise momentum equation, with the obvious exclusion of pressure. Although the total heat flux resulting from the CMT approach is not strictly constant in time, it oscillates around its mean value under statistically steady conditions. Differences of the results obtained with the CMT and CHF approaches have been described by Abe & Antonia (Reference Abe and Antonia2017) and Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021), which although generally small deserve some attention.
Figure 1. Definition of the coordinate system for DNS of pipe flow, where 
$z$, 
$r$ and 
$\phi$ are the axial, radial and azimuthal directions, respectively, 
$R$ is the pipe radius, 
$L_z$ the pipe length and 
$u_b$ is the bulk velocity.
The computer code used for the DNS is the spin-off of an existing solver previously used to study Rayleigh–Bénard convection in cylindrical containers at extreme Rayleigh numbers (Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013). That code is, in turn, the evolution of the solver originally developed by Verzicco & Orlandi (Reference Verzicco and Orlandi1996), and used for DNS of pipe flow by Orlandi & Fatica (Reference Orlandi and Fatica1997). A second-order finite-difference discretisation of the incompressible Navier–Stokes equations in cylindrical coordinates is used, based on the classical marker-and-cell method (Harlow & Welch Reference Harlow and Welch1965), whereby pressure and passive scalars are located at the cell centres, whereas the velocity components are located at the cell faces, thus removing odd–even decoupling phenomena and guaranteeing discrete conservation of the total kinetic energy and passive scalar variance in the inviscid limit. The Poisson equation resulting from enforcement of the divergence-free condition is efficiently solved by double trigonometric expansion in the periodic axial and azimuthal directions, and inversion of tridiagonal matrices in the radial direction (Kim & Moin Reference Kim and Moin1985). An extensive series of previous studies about wall-bounded flows from this group proved that second-order finite-difference discretisation yields in practical cases of wall-bounded turbulence results which are by no means inferior in quality to those of pseudo-spectral methods (e.g. Moin & Verzicco Reference Moin and Verzicco2016; Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016). A crucial issue is the proper treatment of the polar singularity at the pipe axis. A detailed description of the subject is reported in Verzicco & Orlandi (Reference Verzicco and Orlandi1996), but, basically, the radial velocity 
$u_r$ in the governing equations is replaced by 
$q_r = r u_r$ (
$r$ is the radial space coordinate), which by construction vanishes at the axis. The governing equations are advanced in time by means of a hybrid third-order low-storage Runge–Kutta algorithm, whereby the diffusive terms are handled implicitly, and convective terms in the axial and radial direction explicitly. An important issue in this respect is the convective time step limitation in the azimuthal direction, due to intrinsic shrinking of the cells size toward the pipe axis. To alleviate this limitation, we use implicit treatment of the convective terms in the azimuthal direction (Akselvoll & Moin Reference Akselvoll and Moin1996; Wu & Moin Reference Wu and Moin2008), which enables marching in time with similar time step as in planar domains flow in practical computations. In order to minimise numerical errors associated with implicit time stepping, explicit and implicit discretisations of the azimuthal convective terms are linearly blended with the radial coordinate, in such a way that near the pipe wall the treatment is fully explicit, and near the pipe axis it is fully implicit. The code was adapted to run on clusters of graphic accelerators (GPUs), using a combination of CUDA Fortran and OpenACC directives, and relying on the CUFFT libraries for efficient execution of fast Fourier transforms (FFTs) (Ruetsch & Fatica Reference Ruetsch and Fatica2014).
From now on, inner normalisation of the flow properties will be denoted with the ‘+’ superscript, whereby velocities are scaled by 
$u_{\tau }$, wall distances by 
$\nu /u_{\tau }$ and temperatures with respect to the friction temperature,
\begin{equation} T_{\tau} = \frac{\alpha}{u_{\tau}} \left\langle \frac{\mathrm{d} {T}}{\mathrm{d} y} \right\rangle_w. \end{equation}
In particular, the inner-scaled temperature is defined as 
$\theta ^+ = (T-T_w)/T_{\tau }$, where 
$T$ is the local temperature and 
$T_w$ is the wall temperature. Capital letters are used to denote flow properties averaged in the homogeneous spatial directions and in time, brackets to denote the averaging operator and lowercase letters to denote fluctuations from the mean. Finally, bulk values of axial velocity and temperature are defined as
\begin{equation} u_b = 2 \left. \int_0^R r \langle u_z \rangle \, \mathrm{d} r \right/ R^2 , \quad T_b = 2 \left. \int_0^R r \langle T \rangle \, \mathrm{d} r \right/ R^2 . \end{equation} A list of the main simulations that we have carried out is given in table 1. The mesh resolution is designed based on the criteria discussed by Pirozzoli & Orlandi (Reference Pirozzoli and Orlandi2021). In particular, the collocation points are distributed in the wall-normal direction so that approximately 30 points are placed within 
$y^+ \le 40$ (
$y=R-r$ is the wall distance), with the first grid point at 
$y^+ < 0.1$, and the mesh is progressively stretched in the outer wall layer in such a way that the mesh spacing is proportional to the local Kolmogorov length scale, which there varies as 
$\eta ^+ \approx 0.8 {y^+}^{1/4}$ (Jiménez Reference Jiménez2018). Regarding the axial and azimuthal directions, finite-difference simulations of wall-bounded flows yield grid-independent results as long as 
${\rm \Delta} x^+ \approx 10$, 
$R^+ {\rm \Delta} \phi \approx 4.5$ (Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016), hence we have selected the number of grid points along the homogeneous flow directions as 
$N_z = L_z / R \times Re_{\tau } / 9.8$, 
$N_{\phi } \sim 2 {\rm \pi}\times Re_{\tau } / 4.1$. According to the established practice (Hoyas & Jiménez Reference Hoyas and Jiménez2006; Ahn et al. Reference Ahn, Lee, Lee, Kang and Sung2015; Lee & Moser Reference Lee and Moser2015), the time intervals used to collect the flow statistics 
$({\rm \Delta} t_{stat})$ are reported as a fraction of the eddy-turnover time (
$R/u_{\tau }$).
Table 1. Flow parameters for DNS of pipe flow. Cases are labelled in increasing order of Reynolds number, from A to F. Suffixes SH and LO indicate DNS in short and long domains, respectively; FF, FR and FZ denote refinement along the 
$\phi$, 
$r$ and 
$z$ directions, respectively.
The sampling errors for some key properties discussed in this paper have been estimated using the method of Russo & Luchini (Reference Russo and Luchini2017), based on extension of the classical batch means approach. The results of the uncertainty estimation analysis are listed in table 2, where we provide expected values and associated standard deviations for the Nusselt number (
$Nu$), mean temperature at the pipe centreline (
$\varTheta ^+_{{CL}}$) and peak temperature variance and its wall distance (
$\langle {\theta ^2}\rangle ^+_{{IP}}$ and 
$y^+_{{IP}}$, respectively). We find that the sampling error is generally quite limited, being larger in the largest DNS, which have been run for shorter time. In particular, in DNS-F, the expected sampling error in Nusselt number, centreline temperature and peak temperature variance is approximately 
$0.5\,\%$. Additional tests aimed at establishing the effect of axial domain length and grid size have been carried out for the DNS-C flow case, whose results are also reported in table 2. We find that even halving the pipe length yields minimal change in the basic flow properties, which is well within the uncertainty bounds. This is in contrast to properties related to the velocity field, which are significantly affected from use of a short domain (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021). The interesting consequence of this observation is the possibility to carry out DNS of scalar fields in more limited domains, as also noted by Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021). In order to quantify uncertainties associated with numerical discretisation, additional simulations have been carried out by doubling the number of grid points in the azimuthal, radial and axial directions, respectively. Based on the data reported in the table, after discarding the short-pipe case, we can thus quantify the uncertainty due to numerical discretisation and limited pipe length to be approximately 
$0.2\,\%$ for the Nusselt number, 
$0.4\,\%$ for the pipe centreline temperature and 
$0.7\,\%$ for the peak temperature variance.
Table 2. Uncertainty estimation study: mean values of representative quantities and standard deviation of their estimates, where 
$Nu$ is the Nusselt number, 
$\varTheta _{{CL}}^+$ is the mean pipe centreline temperature, 
$\langle \theta _z^2\rangle ^+_{{IP}}$ is the peak temperature variance and 
$y^+_{{IP}}$ is its distance from the wall.
3. Results
3.1. General organisation of the temperature field
Qualitative information about the structure of the flow field is provided by instantaneous perspective views of the axial velocity and temperature fields, which we show in figure 2. Although finer-scale details are visible at the higher 
$Re_b$, the flow in the cross-stream planes is always characterised by a limited number of bulges distributed along the azimuthal direction, which correspond to alternating intrusions of high-speed fluid from the pipe core and ejections of low-speed fluid from the wall. Streaks are visible in the near-wall cylindrical shells, whose organisation has clear association with the cross-stream pattern. Specifically, 
$R$-sized low-speed streaks are observed in association with large-scale ejections, whereas 
$R$-sized high-speed streaks occur in the presence of large-scale inrush from the core flow. At the same time, smaller streaks scaling in wall units appear, corresponding to buffer-layer ejections/sweeps. Hence, organisation of the flow on at least two length scales is apparent here, whose separation increases with 
$Re_{\tau }$. As figure 2 shows, the temperature field has the same qualitative organisation as axial velocity, and low-speed streaks correspond to low-temperature thermal streaks. This is not surprising, given the formal similarity of the controlling equations at 
$Pr=1$, and close association of the two quantities pointed out in many previous studies (e.g. Abe & Antonia Reference Abe and Antonia2009; Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016; Alcántara-Ávila, Hoyas & Pérez-Quiles Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2018). It is interesting that this association includes both the large flow scales in the pipe core and the small, near-wall streaks. Zooming in closer (see figure 3), one can nevertheless detect some differences between the two fields, in that temperature tends to form sharper fronts, whereas the axial velocity field tends to be more blurred. As noted by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016), this is due to the fact that the axial velocity is not passively advected, but rather it can react to the formation of fronts through feedback pressure. As a result, whereas the organisation at large scales is similar, smaller features are found in the temperature fields, as clearly highlighted in the corresponding spectral densities.
Figure 2. (a),(c) Instantaneous axial velocity and (b),(d) temperature contours in turbulent pipe flow as obtained from (a),(b) DNS-A and (c),(d) DNS-F. Thirty contours (from zero to the mean centreline value) are shown on a cross-stream plane and on a near-wall cylindrical shell (
$y^+ \approx 15$), in colour scale from blue to red.
Figure 3. (a) Instantaneous axial velocity and (b) temperature contours in a subregion of the pipe cross section for DNS-F.
The spectral maps of 
$u_z$ and 
$\theta$ are depicted in figure 4, for the DNS-F flow case. In order to isolate changes in the typical length scales, in the figure we show the azimuthal spectral densities normalised by the respective variances, defined as
\begin{equation} \hat{E}_{x}(k_{\phi}) = {E}_{x}(k_{\phi}) / \langle x^2 \rangle, \end{equation}
where 
$k_{\phi }=2 {\rm \pi}/ \lambda _{\phi }$ is the relevant wavenumber for the 
$\phi$ direction and 
$x$ is the generic flow property. The axial velocity spectra clearly bring out a two-scale organisation of the flow field, with a near-wall peak associated with the wall regeneration cycle (Jiménez & Pinelli Reference Jiménez and Pinelli1999), and an outer peak associated with outer-layer large-scale motions (Hutchins & Marusic Reference Hutchins and Marusic2007). The latter peak is found to be centred around 
$y/R \approx 0.3$, and to correspond to eddies with typical wavelength 
$\lambda _{\phi } \approx 1.5 R$, consistent with that found by Ahn et al. (Reference Ahn, Lee, Lee, Kang and Sung2015) for pipe flow at 
$Re_{\tau } = 3000$. Secondary peaks corresponding to harmonics of this fundamental wavelength are also observed here, suggesting that the typical outer modes are not purely sinusoidal with respect to the azimuthal direction. Notably, very similar organisation is found in the temperature field, the main difference being a somewhat broader peak at large wavelengths. Both the axial velocity and the temperature field exhibit a prominent spectral ridge corresponding to modes with typical azimuthal length scale 
$\lambda _{\phi } \sim y$, extending over about two decades, which can be interpreted as the footprint of a hierarchy of wall-attached eddies following Tonwsend's hypothesis (Townsend Reference Townsend1976). Hence, we may expect that inferences of the attached-eddy hypothesis regarding the behaviour of the axial velocity field also carry on to the temperature field.
Figure 4. Variation of pre-multiplied, normalised azimuthal spectral densities of 
$u_z$ (
$\hat {E}_{u_z}$, (a)) and 
$\theta$ (
$\hat {E}_{\theta }$, (b)) with wall distance, for flow case DNS-F. Wall distances and wavelengths are reported both in inner units (bottom, left), and in outer units (top, right). The solid diagonal line marks the trend 
$\lambda _{\phi } = 7.16 y$. Contour levels from 0.05 to 0.5 are shown, in intervals of 0.05.
Differences between velocity and scalar spectra are better scrutinised in figure 5, where we show spectral densities at a discrete set of wall distances. Figure 5(a) clearly brings out the bi-modal distribution of energy between the inner and the outer energetic sites. At intermediate wall distances (
$y^+ \approx 100$) there is some mild evidence for a 
$\hat {E}_{x} \sim \lambda _{\phi }$ range which is also predicted from Townsend's theory (Nickels et al. Reference Nickels, Marusic, Hafez and Chong2005). Most importantly, the figure shows extra energy at small wavelengths in the temperature spectra, with exception of the nearest wall location. This difference is emphasised in figure 5(b), which shows spectra at 
$y/R=0.3$ (at which the Taylor microscale Reynolds number is 
$Re_{\lambda } \approx 400$) in the classical log–log, non-pre-multiplied representation. Whereas the 
$u_z$ and 
$\theta$ spectra are very similar at the largest scales of motion, temperature tends to have a more shallow decay in the inertial and dissipative regions. This is well seen in the compensated representations shown in the insets. Whereas the classical 
$k^{-5/3}$ behaviour can be traced in the 
$u_z$ spectra (at least in a tiny range of wavenumbers), the 
$\theta$ spectra seem to feature instead a 
$k^{-4/3}$ range, which is the theoretically expected behaviour for passive scalars in sheared turbulence (Lohse Reference Lohse1994).
Figure 5. (a) Pre-multiplied, normalised spectral densities of 
$u_z$ (solid) and 
$\theta$ (dashed), at 
$y^+=15$ (gold), 
$y^+=50$ (green), 
$y^+=100$ (cyan) and 
$y/R=0.3$ (purple), for the DNS-F flow case. (b) Normalised spectral densities of 
$u_z$ (solid) and 
$\theta$ (dashed) at 
$y/R=0.3$, compensated by 
$k^{5/3}$ (top inset) and by 
$k^{4/3}$ (bottom inset).
Differences between axial velocity and temperature fields are also apparent in the close proximity of the wall. In figure 6 we show the probability density functions (p.d.f.s) for the wall-normal derivatives of 
$u_z$ and 
$\theta$. Both variables seem to tend to limit distributions in the infinite-
$Re$ limit, however whereas 
$\theta$ is mathematically bound to be positive, hence its wall-normal derivative must also be positive, 
$u_z$ can have instantaneously negative values corresponding to local flow reversal. As a result, we find that the p.d.f. of the temperature gradient is well approximated by a log–normal distribution, as resulting from random multiplicative events. On the other hand, the p.d.f. of 
$u_z$ obviously deviates from log–normality near the origin, but also its positive tail seems to be less prominent than for 
$\theta$. The existence of local flow inversion at the wall, although with small probability (about 
$0.1\,\%$ overall) was noted in several previous studies (e.g. Lenaers et al. Reference Lenaers, Li, Brethouwer, Schlatter and Örlü2012), and related to the presence of oblique vortices inducing negative pressure fluctuations. This again corroborates the interpretation of different behaviour of 
$u_z$ and 
$\theta$ as being due to the action of pressure.
Figure 6. Probability density function of wall-normal derivatives of (a) axial velocity and (b) temperature. The colour codes are as in table 1. The dashed grey lines denote a log–normal distribution made to fit the DNS-F data.
3.2. Mean temperature field
The mean temperature profile in turbulent pipes has received extensive attention from theoretical and experimental studies, and the general consensus (Kader Reference Kader1981) is that a logarithmic law well fits the experimental data. Recent studies have instead questioned the validity of the logarithmic law for the mean velocity field at finite Reynolds number (Jiménez & Moser Reference Jiménez and Moser2007; Pirozzoli, Bernardini & Orlandi Reference Pirozzoli, Bernardini and Orlandi2014; Lee & Moser Reference Lee and Moser2015), and corrections to account for the effect of the core flow on the overlap layer have been proposed (e.g. Luchini Reference Luchini2017; Cantwell Reference Cantwell2019; Monkewitz Reference Monkewitz2021). Such corrections mainly amount to addition of a linear term to the logarithmic profile which can be justified as a higher-order term in the asymptotic matching between the inner and the outer layer (Afzal & Yajnik Reference Afzal and Yajnik1973), or as due to the presence of a mean pressure gradient in internal flows (Luchini Reference Luchini2017). Deviations of the profiles of passive scalars from the assumed logarithmic distribution were also observed in DNS of channel flow by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016) and Alcántara-Ávila et al. (Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021), amounting to a linear correction whose inner-scaled slope decreases with 
$Re_{\tau }$. In figure 7, we show a series of temperature profiles computed with the present DNS, fitted with a logarithmic function with inverse slope 
$k_{\theta }=0.459$, determined as described in the next section. The additive constant resulting from best fitting of the DNS data is 
$C_{\theta } \approx 5.78$. As shown in the inset of figure 7(a), the velocity profiles for 
$Re_{\tau } \ge 10^3$ follow this distribution with deviations smaller than 
$0.1$ wall units from 
$y^+ \approx 30$ to 
$y/R \approx 0.1$. Hence, the standard log law is a good approximation of the temperature profile in the overlap layer for most practical purposes. The functional expression proposed by Kader (Reference Kader1981, (9), circles in panel (a)) is also found to provide reasonable approximation of the data throughout the wall layer, even if with somewhat unnatural behaviour in the buffer layer and slight overprediction in the outer wall layer.
Figure 7. (a) Inner-scaled mean temperature profiles and (b) corresponding log-law diagnostic functions. Deviations from the assumed logarithmic wall law, 
$\varTheta _{log}^+ = \log y^+ / 0.459 + 5.78$, are highlighted in the inset of (a). Circles denote the functional approximation proposed by Kader (Reference Kader1981), here evaluated for 
$Re_{\tau }=6019$, 
$Pr=1$. In (b), the dashed horizontal line denotes the inverse of the Kármán constant, 
$1/k_{\theta }$, and the dash-dotted lines in the inset denote the linear fit (3.3), with 
$k_{\theta }=0.459$, 
$\alpha _{\theta } = 1.81$. See table 1 for colour codes.
Very similar considerations apply to the mean axial velocity profiles (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021, figure 3), which are visually well fitted with a logarithmic distribution, with estimated value of the von Kármán constant 
$k \approx 0.387$ (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021), hence much less than 
$k_{\theta }$, as consistently noted in previous studies on the subject. Based on the results presented in § 3.1, this difference can be justified by recalling that in the log layer (if present), 
$k/k_{\theta } \approx \nu _t / \alpha _t$, where 
$\nu _t$ and 
$\alpha _t$ are the turbulent kinematic and thermal diffusivities. As noted previously, the temperature field has a tendency to form sharper fronts, with steeper gradients, hence its effective diffusivity is expected to be larger than for the axial momentum. Accordingly, one may expect 
$k$ to be smaller than 
$k_{\theta }$, and the turbulent Prandtl number to be less than unity in the outer layer. More formal arguments to justify this difference, based on the properties of the similarity solution for the logarithmic mean profile over the inertial (non-diffusive) domain, were offered by Zhou, Klewicki & Pirozzoli (Reference Zhou, Klewicki and Pirozzoli2019).
More detailed scrutiny about the behaviour of the mean temperature profile is carried out in figure 7(b), where we show the logarithmic diagnostic function,
\begin{equation} \varXi_{\theta} = y^+ \, \mathrm{d} {\varTheta}^+{/} \mathrm{d} y^+ , \end{equation}
which is expected to be constant in the presence of a genuine logarithmic layer. As found previously for the axial velocity field (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021, figure 4), no region with flat distribution of this indicator is, in fact, present. Rather, we note the occurrence of a nearly linear distribution from 
$y^+ \approx 100$ to 
$y/R \approx 0.4$, whose slope is approximately constant in outer units, hence the diagnostic function can be expressed as
\begin{equation} \varXi_{\theta} \approx \frac 1{k_{\theta}} + \alpha_{\theta} \frac yR ,\end{equation}
with 
$\alpha _{\theta } \approx 1.81$. In other words, whereas a simple logarithmic profile is a reasonable approximation for engineering estimates, a linear correction yields significant improvement in the representation of the temperature profile, over a wider range of wall distances. Based on (3.3), a genuine log layer in the mean temperature profile would only emerge at infinite Reynolds number.
The structure of the core region of the flow is inspected in figure 8, where the mean temperature profiles are shown in defect form. Disregarding the DNS at the lowest Reynolds number (DNS-A and DNS-B) the scatter across the various temperature profiles for 
$y/R \ge 0.2$ is less than 
$1\,\%$, which suggests that outer-layer similarity is very nearly achieved. As suggested by Pirozzoli (Reference Pirozzoli2014) and Orlandi, Bernardini & Pirozzoli (Reference Orlandi, Bernardini and Pirozzoli2015), the core velocity and temperature profiles can be closely approximated with simple universal quadratic distributions, which one can derive under the assumption of constant eddy diffusivity of momentum and temperature. In particular, we find that the formula
\begin{equation} \varTheta_{{CL}}^+{-} \varTheta^+ = C_O ( 1 - y/R )^2, \end{equation}
with 
$C_O = 5.5$, fits the mean temperature distribution in the pipe core quite well. Closer to the wall, the corrected logarithmic profile sets in at 
$y/R \lesssim 0.44$, here expressed in outer coordinates,
\begin{equation} \varTheta_{{CL}}^+{-} \varTheta^+= - \frac 1{k_{\theta}} \, \log (y/R) - \alpha_{\theta} \frac yR + B_{\theta} , \end{equation}
where data fitting yields 
$B_{\theta }=0.732$. Although more elaborate descriptions of the outer velocity profiles are possible (e.g. Krug, Philip & Marusic Reference Krug, Philip and Marusic2017; Luchini Reference Luchini2018), the composite profile compounding (3.4) and (3.5) yields accurate representation of the whole outer-layer mean temperature profile to within the scatter of the available DNS data.
Figure 8. Mean defect temperature profiles in (a) linear and (b) semi-logarithmic scale. The dashed grey line marks a parabolic fit of the DNS data (
$\varTheta ^+_{{CL}}-\varTheta ^+ = 5.5 (1-y/R)^2$) and the dot-dashed purple line in (b) the corrected outer-layer logarithmic fit 
$\varTheta ^+_{{CL}}-\varTheta ^+ = 0.732 - 1/0.459 \log (y/R) - 1.81 (y/R)$. Only datasets DNS-C to DNS-F are shown here, see table 1 for colour codes.
3.3. Heat transfer coefficients
The primary subject of engineering interest in the study of passive scalar fields is the transfer coefficient at the wall, which can be expressed in terms of the Stanton number,
\begin{equation} St= \frac{\alpha \left\langle \dfrac{\mathrm{d} {T}}{\mathrm{d} y} \right\rangle_w}{u_b ( T_m - T_w )} = \frac 1{u_b^+ \theta_m^+}, \end{equation}
where 
$T_m$ is the mixed mean temperature (Kays et al. Reference Kays, Crawford and Weigand1980),
\begin{equation} T_m = 2 \left. \int_0^R r \langle u_z \rangle \langle T \rangle \mathrm{d} r \right/ ( u_b R^2 ) , \end{equation}
with 
$\theta _m = (T_m-T_w)/T_{\tau }$ or, more frequently, in terms of the Nusselt number,
\begin{equation} Nu = St \, Re_b \, Pr . \end{equation} A predictive formula for the heat transfer coefficient in wall-bounded turbulent flows was derived by Kader & Yaglom (Reference Kader and Yaglom1972), based on the assumed existence of logarithmic layers for the mean velocity and temperature profiles as a function of the wall distance, and on universality of the core layer in defect representation. For the purpose of critically evaluating the assumptions made in the derivation of Kader's formula, we show (in figure 9) the distributions of the bulk and mean centreline values (namely, at 
$r=0$) of velocity (figure 9a) and temperature (figure 9b), as a function of the friction Reynolds number, Consistently with theoretical expectations (e.g. Monkewitz Reference Monkewitz2021), the data suggest logarithmic increase of the bulk and centreline velocity with 
$Re_{\tau }$ according to
\begin{equation} u_b^+= \frac 1k \log Re_{\tau} + B, \quad U_{{CL}}^+= \frac 1{k} \log Re_{\tau} + B_{{CL}}, \end{equation}
with 
$k=0.387$, 
$B= 1.229$ and 
$B_{{CL}}=5.85$ (Pirozzoli et al. Reference Pirozzoli, Romero, Fatica, Verzicco and Orlandi2021). The relative standard deviation of the above formulae with respect to DNS data is approximately 
$0.2\,\%$ for 
$u_b$ and 
$0.5\,\%$ for 
$U_{{CL}}$. Similar trends are here observed for the temperature field, which also exhibits logarithmic growth of the mean and centreline temperature with 
$Re_{\tau }$, according to
\begin{equation} \theta_b^+= \frac 1{k_{\theta}} \log Re_{\tau} + \beta, \quad \varTheta_{{CL}}^+= \frac 1{k_{\theta}} \log Re_{\tau} + \beta_{{CL}}, \end{equation}
with 
$\theta ^+_b = (T_b-T_w)/T_{\tau }$, and fit of the DNS data suggests 
$k_{\theta } = 0.459$, 
$\beta = 2.96$ and 
$\beta _{{CL}}=6.46$. The relative standard deviation of the above formulae with respect to DNS data is approximately 
$0.06\,\%$ for 
$\theta _b$ and 
$0.5\,\%$ for 
$\varTheta _{{CL}}$, hence the quality of the fits is excellent as for the axial velocity field, and the resulting estimates of the flow constants, and especially of the scalar von Kármán constant appear to be quite robust.
Figure 9. Bulk and centreline values of (a) axial velocity and (b) temperature. Bulk values (
$u^+_b$, 
$\theta ^+_b$) are denoted with squares and centreline values (
$U^+_{{CL}}$, 
$\varTheta ^+_{{CL}}$) with circles. Diamonds in (b) denote the mixed mean temperature (
$\theta ^+_m$). The dashed lines denote logarithmic fits of the DNS data. The dash-dotted line in (b) refers to the fit (3.11).
Considering a large number of experimental works, Kader & Yaglom (Reference Kader and Yaglom1972) suggested 
$k_{\theta } \approx 0.47$, and provided empirical formulae for the additive constants as a function of the Prandtl number, 
$\beta (Pr)=12.5 Pr^{2/3} + 1/k_{\theta } \log Pr - 5.3$ and 
$\beta _{{CL}} = \beta + 0.6$. Studies carried out by means of DNS in planar channel flows with CHF show some scatter in the prediction of 
$k_{\theta }$, likely due to low-Reynolds-number effects. For instance, Kawamura et al. (Reference Kawamura, Abe and Matsuo1999) reported 
$0.40 \le k_{\theta } \le 0.42$, Abe et al. (Reference Abe, Kawamura and Matsuo2004) reported 
$k_{\theta } \approx 0.43$, whereas more data at higher Reynolds number suggest 
$k_{\theta } \approx 0.44$ (Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021). Studies carried out with CMT typically tend to yield slightly higher values, namely 
$k_\theta \approx 0.46$ (Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016), thus closer to Kader's prediction. It should be noted that all the above estimates were based on attempting to fit the temperature profiles with a logarithmic law, which as shown previously may not be a good approximation, especially at low Reynolds number. The method herein used to estimate the scalar von Kármán constant from the bulk and centreline temperatures yields greater accuracy and robustness, resulting in a value similar to that suggested by Kader & Yaglom (Reference Kader and Yaglom1972), and to values obtained using a similar approach, based on channel flow DNS data (Abe & Antonia Reference Abe and Antonia2017). It is also worth pointing out that the additive constants 
$\beta, \beta _{{CL}}$ in (3.10a,b) do depend on the molecular Prandtl number (Kader & Yaglom Reference Kader and Yaglom1972), hence the values herein reported are specific to the 
$Pr=1$ case.
Determination of the appropriate Reynolds number trends of the mixed mean temperature is not as straightforward as for the bulk or the centreline temperature, because it involves integrating the product of the mean velocity and temperature distributions. Simple developments (Kader & Yaglom Reference Kader and Yaglom1972) suggest the expected behaviour to be
\begin{equation} \theta_{m}^+= \varTheta_{{CL}}^+{-} \beta_2 + \frac{\beta_3}{u_b^+},\end{equation}
with 
$\beta _2$ and 
$\beta _3$ universal constants, hence, on account of (3.9a,b) and (3.10a,b), deviations from logarithmic dependence on 
$Re_{\tau }$ are to be expected. Figure 9 shows that it is indeed the case, and the mixed mean temperature (diamond symbols) seems to follow a near logarithmic distribution only at the highest Reynolds numbers under consideration. Values of the constants 
$\beta _2=4.92$, 
$\beta _3=39.6$, in (3.11) yield a fair fit of the DNS data. Although the scatter in the analytical fit (see the dash-dotted line) is more significant than for the other flow properties, this is clearly more accurate than a simple logarithmic fit. Corrections to the logarithmic law in the mixed mean temperature were generally disregarded in previous studies (Kader & Yaglom Reference Kader and Yaglom1972; Abe & Antonia Reference Abe and Antonia2017), although they can be included in the analysis without additional difficulty.
Proceeding as proposed by Kader & Yaglom (Reference Kader and Yaglom1972), it is straightforward to derive a predictive formula for the inverse Stanton number,
\begin{equation} \frac 1{St} = \frac k{k_{\theta}} \frac 8{\lambda} + \left( \beta_{{CL}} - \beta_2 - \frac k{k_{\theta}} B \right) \sqrt{\frac 8{\lambda}} + \beta_3, \end{equation}
where the friction factor 
$\lambda = 8 / {u_b^+}^2$ can be obtained from (3.9a,b). Assuming strictly logarithmic variation of the mixed mean temperature with 
$Re_{\tau }$ (hence, setting 
$\beta _3=0$), Kader & Yaglom (Reference Kader and Yaglom1972) arrived at the following expression,
\begin{equation} \frac 1{St} = \frac{2.12 \log ( Re_b \sqrt{\lambda/4} ) + 12.5 Pr^{2/3} + 2.12 \log Pr - 10.1}{\sqrt{\lambda/8}} ,\end{equation}which could also be rearranged to a form more similar to (3.12). Additional correlations which are in wide use in the engineering practice were proposed by Gnielinski (Reference Gnielinski1976),
\begin{equation} Nu = \frac{Pr \lambda /8 ( Re_b - 1000 )}{1 + 12.7 (\lambda/8)^{1/2} ( Pr^{2/3} - 1)}, \end{equation}and by Kays et al. (Reference Kays, Crawford and Weigand1980),
\begin{equation} Nu = 0.022 \, Re_b^{0.8} \, Pr^{0.5}. \end{equation}
Last, direct fitting of the DNS data (at 
$Pr=1$) with a power-law expression yields
\begin{equation} Nu = 0.0219 \, Re_b^{0.804}. \end{equation} All the above predictive formulae are tested in figure 10, showing the predicted inverse Stanton number (figure 10a) and Nusselt number (figure 10b). With little surprise, we find that (3.12) with DNS-informed definition of the coefficients matches the DNS data quite well, with maximum relative error of 
$0.8\,\%$. Despite minor differences in the coefficients with respect to the baseline predictive formula (3.15), the direct power-law fit given in (3.16) is also quite accurate, except at low 
$Re_b$. All other formulae fall short of the DNS data for 
$St$ by up to 
$5\,\%$. This difference may be partly due to inaccuracy of correlations based on old experimental data, to the fact that those are mainly tuned for the 
$Pr=0.71$ case, whereas here 
$Pr=1$, but also to the fact that the CMT setup herein used tends to slightly overpredict the heat flux as compared with the CHF approach (Abe & Antonia Reference Abe and Antonia2017; Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021). It is interesting that differences are levelled off when the popular representation in terms of the Nusselt number is used, as in figure 10(b). This suggests that the 
$1/St$ representation should be used when relatively small differences must be discriminated, or perhaps the compensated Nusselt representation used in the inset of figure 10(b).
Figure 10. Distribution of (a) inverse Stanton number and (b) Nusselt number obtained from DNS (circles), and as predicted from (3.12) (solid line), from Kader's formula ((3.13), dashed), from the power-law data fit (3.16) (dotted), from Gnielinski analogy ((3.14), dot-dot-dashed) and from Kays–Crawford correlation ((3.15), dot-dashed). The inset in (b) shows the Nusselt number in compensated form (
$Nu \times Re_b^{-0.8}$).
3.4. Temperature fluctuations statistics
The distributions of the axial velocity and temperature variances are shown in figure 11, in inner scaling. All the profiles feature a prominent peak in the buffer layer at 
$y^+ \approx 15$, and an outer-layer shoulder which starts to form at sufficiently high Reynolds number. The most notable Reynolds number effect on the temperature variance profiles is sustained increase, as is the case of the axial velocity variance (Marusic & Monty Reference Marusic and Monty2019). According to Townsend's attached eddy model (Townsend Reference Townsend1976), growth of the wall-parallel velocity variances is expected to be logarithmic with 
$Re_{\tau }$, on account of increased influence of ‘distant’ eddies. According to the spectral maps presented in figure 4, the attached-eddy hypothesis is also expected to apply to the temperature field. In fact, logarithmic growth of the peak temperature variance in turbulent planar channels was observed by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016). This is confirmed by the present pipe DNS data, see figure 11(b), which compare the growth of the peak temperature variance with the axial velocity variance. Although the inferred growth rate is the same, the magnitude of the temperature variance peak is larger than for the axial velocity. This is the consequence (Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2014) of near equality of the corresponding production terms in the buffer layer, however with extra energy draining in the axial velocity variance equation from the pressure term which tends to equalise kinetic energy across all the velocity components. No evidence is found based on the present data for saturation of the logarithmic growth, which has been inferred for the axial velocity variance in recent theoretical studies (Chen & Sreenivasan Reference Chen and Sreenivasan2021).
Figure 11. Distribution of (a) temperature variances and (b) corresponding peak value as a function of 
$Re_{\tau }$. The dashed lines in (a) denote the distributions of the axial velocity variance. In (b), circles correspond to the peak temperature variance and squares to the peak axial velocity variance. Dash-dotted and dashed lines correspond to the associated logarithmic fits, namely 
$\langle \theta ^2\rangle ^+_{{IP}} = 0.68 \log Re_{\tau } + 3.9$, 
$\langle u^2_z\rangle ^+_{{IP}} = 0.67 \log Re_{\tau } + 3.3$. Refer to table 1 for colour codes.
The distributions of the turbulent heat flux, 
$\langle u_r \theta \rangle$, shown in figure 12, are visually indistinguishable from the corresponding turbulent shear stresses (reported with dashed lines). This is a further confirmation that the lift-up mechanism which is responsible for the establishment of correlation of 
$u_z$ and 
$\theta$ with vertical velocity fluctuations, 
$u_r$, is nearly the same, and of essentially linear nature (Jiménez Reference Jiménez2013). In both cases, the peak position grows approximately as the square root of 
$Re_{\tau }$, corresponding to the minimal wall distance for a logarithmic layer to develop (Klewicki, Fife & Wei Reference Klewicki, Fife and Wei2009). Slight differences between the two distributions are nevertheless responsible for differences in the mean axial velocity and temperature profiles, on account on mean momentum balance and mean thermal balance (Saha et al. Reference Saha, Klewicki, Ooi and Blackburn2015; Zhou, Pirozzoli & Klewicki Reference Zhou, Pirozzoli and Klewicki2017). These differences are better appreciated in figure 12(b), where we show the peak turbulent heat flux and shear stress as a function of 
$Re_{\tau }$. Based on mean momentum balance, and assuming the presence of a logarithmic layer (which is a bit inaccurate in view of what previously said), Orlandi et al. (Reference Orlandi, Bernardini and Pirozzoli2015) inferred that the peak turbulent shear stress should scale as
\begin{equation} \max \langle u_r u_z \rangle \approx 1 - \frac 2{\sqrt{k Re_{\tau}}} . \end{equation}Similarly, it can be readily shown, based on mean thermal balance and assuming a log layer in the mean temperature distribution, that the peak turbulent heat flux should scale as
\begin{equation} \max \langle u_r \theta \rangle \approx 1 - \frac 2{\sqrt{k_{\theta} Re_{\tau}}} . \end{equation}These two trends are compared in figure 12(b). Whereas the data points coincide at low Reynolds number, some segregation is observed at the highest Reynolds numbers, at which the peaks tend to follow their respective theoretical distributions.
Figure 12. Distribution of (a) turbulent heat flux and (b) corresponding peak value (complement to one and premultiplied by 
$Re_{\tau }^{1/2}$) as a function of 
$Re_{\tau }$. The dashed lines in (a) (barely visible) correspond to the distributions of the turbulent shear stress. In (b), circles correspond to the peak turbulent heat flux and squares to the peak turbulent shear stress. Dashed and dash-dotted lines correspond to the theoretical predictions (3.17) and (3.18), respectively. Refer to table 1 for colour codes.
A quantity of great relevance in turbulence models of scalar transport is the turbulent Prandtl number, defined as (Cebeci & Bradshaw Reference Cebeci and Bradshaw1984)
\begin{equation} Pr_t = \frac{\nu_t}{\alpha_t} = \frac{\langle u_r u_z \rangle}{\langle u_r \theta \rangle} \frac{\mathrm{d} {\varTheta} / \mathrm{d} y}{\mathrm{d} {U} / \mathrm{d} y}, \end{equation}
whose distributions are shown in figure 13. As expected based on its definition, 
$Pr_t \approx k/k_{\theta } \approx 0.843$, through a large part of the outer layer, say from 
$y^+ \approx 100$ to 
$y/h \approx 0.25$. This is quite similar to what was found in channels (Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017; Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Pérez-Quiles2021) and in general agreement with the values 
$Pr_t \approx 0.85$ suggested in reference publications (Kader Reference Kader1981; Cebeci & Bradshaw Reference Cebeci and Bradshaw1984). Closer to the wall, the turbulent Prandtl number tends to exhibit a plateau with nearly unit value within the buffer layer (
$5 \le y^+ \le 40$), as a result of the close similarity of the velocity and temperature fields in that region. At 
$y^+ \lesssim 1$ the eddy viscosity tends to exceed the eddy diffusivity, and as a consequence 
$Pr_t > 1$. Theoretical estimates for the near-wall behaviour of 
$Pr_t$ were proposed by Cebeci (Reference Cebeci1973), based on a mixing length model with van Driest near-wall damping. This results in the estimate
\begin{equation} Pr_t = \frac{k}{k_{\theta}} \frac{1-\exp ({-}y^+{/}A)}{1-\exp ({-}y^+{/}B)},\end{equation}
where 
$A$ and 
$B$ are the damping functions for the velocity and scalar fields, respectively. The original choice of those two parameters, 
$A=26$, 
$B=34.96$ (for 
$Pr=1$), captures the near-wall growth of 
$Pr_t$, although its value is overpredicted. Changing the damping constant for the temperature field to 
$B = 32.2$ (see the dotted grey line in figure 13a) improves the fit significantly, although the buffer-layer plateau is not captured. Regarding the outer layer, figure 13(b) seems to show tendency to universality in outer scaling at sufficiently high Reynolds number. As suggested by Rotta (Reference Rotta1962) and Abe & Antonia (Reference Abe and Antonia2017), the data are well fitted by a quadratic correlation, and we find
\begin{equation} Pr_t = 0.87 - 0.3 (y/R)^2, \end{equation}
for 
$y/R \gtrsim 0.25$.
Figure 13. Distribution of turbulent Prandtl number, in (a) inner and (b) outer coordinates. The dashed line denotes 
$Pr=k/k_{\theta }=0.843$. In (a), the dash-dotted line denotes the prediction of (3.20) with the original set of constants and the dotted line the same formula, with 
$B=32.2$. The dotted line in (b) denotes the fitting function (3.21). Refer to table 1 for colour codes.
4. Concluding comments
We have analysed the dynamics of passive scalars in turbulent pipe flow up to 
$Re_{\tau } \approx 6000$, at unit Prandtl number. Compared with previous studies, the Reynolds number is sufficiently high here that the results are representative of realistic fully developed turbulence. A general expected finding, at least at 
$Pr=1$, is that passive scalars exhibit a behaviour similar to the axial velocity field, being organised into streaks in the buffer layer and featuring low-wavenumber azimuthal modes in the bulk of the flow. The spectral maps also highlight close similarities at the large scales, which at the higher Reynolds number under study (
$Re_{\tau } \approx 6000$), exhibit a prominent ridge with eddy size proportional to the wall distance, as suggested by Townsend's attached-eddy model. Besides overall similarity, differences are found as the passive scalar field has stronger tendency to form steep fronts, which are prevented in the velocity field by the action of pressure. Hence, scalar spectra tend to exhibit shallower 
$k^{-4/3}$ inertial range, for which the present DNS data provide convincing evidence.
Regarding the one-point statistics, we find that the mean scalar profiles in the overlap layer can be conveniently approximated by logarithmic distributions. However, significant improvement is obtained through addition of a linear corrective term scaling in outer units, which yields good fit of the data up to 
$y/R \approx 0.44$. Further away from the wall, the mean scalar profile is approximated with excellent accuracy by a simple quadratic distribution. Notable differences are found in the von Kármán constants, which we estimate to be 
$k_{\theta } \approx 0.459$ for the scalar field, and 
$k \approx 0.387$ for the axial velocity field, which we obtain by fitting the trends of the respective bulk values with the friction Reynolds number, with estimated error much less than 
$1\,\%$. The DNS data help explaining this difference as the result of greater effective diffusivity of the scalar field resulting from the formation of smaller scales. Reynolds number effects are apparent in the distributions of the passive scalar variance, which show sustained logarithmic growth of the inner peak magnitude, as after Townsend's attached-eddy hypothesis. No evidence for saturation of this growth is found, at the Reynolds numbers under scrutiny. Notably, the growth rate is found to be very nearly the same as for the axial velocity variance.
Full insight into the scalar and velocity statistics provided by DNS allows to pinpoint possible limitations of classical analyses of heat transfer in smooth pipes and of classical modelling assumptions for passive scalar transport. In particular, whereas logarithmic dependence of the mixed mean temperature over 
$Re_{\tau }$ was assumed by Kader (Reference Kader1981), we find that finite-
$Re$ deviations should be accounted for, thus obtaining the predictive formula (3.12), which is found to be more accurate than the classical Kader's formula (3.13). It is noteworthy that deviations of empirical formulae employed in the engineering practice are sometimes hidden in the traditional 
$Nu$ versus 
$Re_b$ representation, whereas they show up much more clearly when the inverse Stanton number is reported, as in figure 10(a). Regarding the distributions of the turbulent Prandtl number, which is needed for turbulence modelling, the DNS data show that in the overlap layer the assumption 
$Pr_t \approx 0.843$ is quite appropriate. However, significant deviations are found farther from the wall, with quadratic decrement with the wall distance. Deviations are also found in the near-wall region, where 
$Pr_t$ exceeds unity. This fact is roughly acknowledged in existing engineering correlations (Cebeci Reference Cebeci1973; Kays et al. Reference Kays, Crawford and Weigand1980), but the accuracy is far from acceptable. Whereas the above-noted shortcoming of common modelling assumptions and predictive formulae may be minor to be appreciated in practical contexts, as may be related with the assumed unit Prandtl number, with the particular form of volumetric heating used as forcing and with assumed isothermal wall, the procedure herein outlined nevertheless serves to illustrate a general use of DNS as a way for testing hypotheses in more general situations, for instance at higher Prandtl number, at which the thermal boundary layer is thinner than the kinematic one.
A natural extension of the present study would be, in fact, considering the case of lower and higher Prandtl numbers, which we could not include in the present DNS owing to memory restrictions. That would allow to verify the Prandtl number dependence of the 
$\beta _{{CL}}$ constant in the main predictive equation for the heat transfer (3.12). Equally interesting would be carrying out the present DNS within the CHF approach, whereby strictly constant heat flux is retained in time. That would allow the reasons for the small (but sizeable) deviations observed in figure 10 to be elucidated with respect to engineering correlations in current use.
Acknowledgements
We thank A. Ceci and M. Yu for help in processing the data. We acknowledge that the results reported in this paper have been achieved using the PRACE Research Infrastructure resource MARCONI based at CINECA, Casalecchio di Reno, Italy, under project PRACE no. 2021240112.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at http://newton.dma.uniroma1.it/database/.
Open access
