Direct numerical simulation of forced thermal convection in square ducts up to ${{Re}}_{\tau } \approx 2000$

Abstract

We carry out direct numerical simulations (DNS) of flow in a turbulent square duct by focusing on heat transfer effects, considering the case of unit Prandtl number. Reynolds numbers up to ${{Re}}_{\tau } \approx 2000$ are considered that are much higher than in previous studies, and that yield clear scale separation between inner- and outer-layer dynamics. Close similarity between the behaviour of the temperature and the streamwise velocity fields is confirmed as in previous studies related to plane channels and pipes. Just like the mean velocity, the mean temperature is found to exhibit logarithmic layers as a function of the nearest wall, but with a different slope. The most important practical implication is the validity of the traditional hydraulic diameter as the correct reference length for reporting heat transfer data, as we show rigorously here. Temperature and velocity fluctuations also have similar behaviour, but apparently logarithmic growth of their inner-scaled peak variances is not observed here, unlike in canonical wall-bounded flows. Analysis of the split contributions to the heat transfer coefficient shows that mean cross-stream convection associated with secondary motions is responsible for about $5\,\%$ of the total. Finally, we use the DNS database to highlight shortcomings of traditional linear closures for the turbulent heat flux, and show that substantial modelling improvement in principle may be obtained by retaining at least the three terms in the vector polynomial integrity basis expansion.

1. Introduction

Heat transfer in internal flows is a subject of utmost relevance in mechanical and aerospace engineering applications. A large amount of experimental and numerical studies have been carried out in the past, and a variety of analytical and semi-empirical prediction tools have been developed, which are reported extensively in classical textbooks (Kays & Crawford 1993; Rohsenow, Hartnett & Cho 1998). Most studies have been carried out for the canonical case of ducts with circular cross-section or planar channels, whereas much less is known about the case of ducts with more complex geometry, which also have great practical relevance, for instance in water draining or ventilation systems, nuclear reactors, heat exchangers, space rockets and turbomachinery. In that case, the typical engineering approach is to use the same correlations established for the case of circular pipes, by replacing the pipe diameter with the hydraulic diameter of the duct (Kays & Crawford 1993; White & Majdalani 2006). Although this approach is found to be rather successful in practice, it lacks solid theoretical foundations. Furthermore, the large scatter in experimental data makes it difficult to quantify the actual accuracy of semi-empirical prediction formulas, which are reported to have $\pm$9 % uncertainty for smooth ducts with uniform heating (Shah & Sekulib 1998).

Heat transfer in square ducts was first studied experimentally by Brundrett & Burroughs (1967), who considered air flow at bulk Reynolds number ${{Re}}_b = H u_b/\nu$ (where $H=2h$ is duct side length and also duct hydraulic diameter, $u_b$ is the bulk velocity, and $\nu$ is the fluid kinematic viscosity) values 33 000 and 67 000. Close similarity between the wall shear stress and the heat flux distributions was shown (with quoted discrepancy $\pm$2 %), which the authors connected with similar mixing action of the secondary currents on momentum and temperature. As a result, they observed that the ratio of the average friction and heat transfer coefficients for a square duct is approximately the same as for a circular pipe. Measurements of the wall-normal temperature profiles highlighted close universality when the local wall heat flux is used for normalization, and the presence of a sizeable logarithmic layer, with extrapolated value of the scalar Kármán constant $\kappa _{\theta } \approx 0.51$. Those findings were supported qualitatively from later experiments by Hirota et al. (1997), who also analysed temperature fluctuations and velocity/temperature fluctuations correlations, and found significant distortions over the cross-section associated with the secondary motions. In that study, the inferred scalar von Kármán constant was $\kappa _{\theta } \approx 0.46$, hence more similar to the values generally quoted for circular pipe flow, namely $\kappa _{\theta } \approx 0.47$ (Kader & Yaglom 1972).

Early computational studies of heat transfer in square ducts relied mostly on the use of Reynolds-averaged Navier–Stokes models. For instance, Launder & Ying (1973) developed a full Reynolds stress closure to predict flow and heat transfer in a square duct. Although the global heat transfer coefficient showed general underprediction by about $10\,\%$, the distribution of the wall heat flux reproduced qualitatively the experimental data of Brundrett & Burroughs (1967). High-fidelity computational studies of heat transfer in square ducts have been quite limited so far. Vázquez & Métais (2002) first studied turbulent flow through a heated square duct by means of large-eddy simulations (LES) at moderate Reynolds number (${{Re}}_b=6000$), considering the case in which one of the walls is hotter than the other three. Accounting for fluid viscosity variation with temperature, they found that turbulent structures near the hot wall become larger than near the other walls, in such a way that wall scaling is satisfied. LES with different wall thermal states, also in the presence of duct rotation, were carried out by Pallares & Davidson (2002). Yang, Chen & Zhu (2009) carried out under-resolved direct numerical simulations (DNS) of turbulent flow in a square duct with uniform volumetric heating, finding good agreement with the temperature profiles measured by Hirota et al. (1997), and with predictions of the heat transfer coefficient resulting from Gnielinski's analogy (Gnielinski 1976). Large-eddy simulations of space-developing turbulent flow through a heated square duct were carried out by Schindler, Younis & Weigand (2019), at bulk Reynolds number up to ${{Re}}_b = 10^4$. As noted by Brundrett & Burroughs (1967), strong correlation of the local wall heat flux with the wall shear stress was found. Sekimoto et al. (2011) carried out DNS of mixed convection in square duct flow ranging from very low to unit Richardson number. They found that modifications of the secondary flows due to buoyancy start at Richardson number ${\approx }0.025$, and the effect of buoyancy becomes dominant at Richardson ${\approx }0.25$ and the eight cross-stream vortices are replaced by two large-scale ones. Turbulence modelling of forced and natural thermal convection is a topic of primary interest for industrial applications, but it has received less attention than its momentum counterpart. Modelling approaches for the turbulent heat flux vector were reviewed by Hanjalić (2002), who pointed out inadequacy of the turbulent Prandtl number concept in complex flows, and advocated the use of second-moment closures or algebraic models for turbulent heat transfer. However, the state of the art for modelling turbulent convection has not advanced much since the time of that review article, and linear eddy diffusivity models relying on use of the turbulent Prandtl number are routinely used.

In this paper, we study heat transfer in fully developed square duct flow with uniform internal heating and isothermal walls, by carrying out DNS at much higher Reynolds numbers than in previous computational studies. Although relatively simple, this set-up includes the main complicating effects involved with non-trivial cross-sectional geometries, and primarily associated with the formation of secondary motions of the second kind (Prandtl 1927; Nikuradse 1930). The assumption herein made of unit molecular Prandtl number (defined as the ratio of the kinematic viscosity to the scalar diffusivity, ${{Pr}}=\nu /\alpha$) is instrumental to scrutinizing more closely differences between the behaviour of the streamwise velocity field and of passively advected scalars. The present study is the continuation of previous efforts (Modesti et al. 2018; Orlandi, Modesti & Pirozzoli 2018; Pirozzoli et al. 2018) targeted to studying turbulent flows in square ducts by means of DNS. In those studies, we found that the intensity of the secondary motions is of the order of a few per cent of the duct bulk velocity, and despite their persistence, they do contribute to the mean duct friction by at most a few per cent. In this study we aim to extend our analysis to the temperature field and to quantify the effect of secondary flows on the wall heat flux. Additionally, we assess the accuracy of the turbulent Prandtl number concept and of constitutive relations for the turbulent heat flux vector.

2. Methodology

The numerical simulation of incompressible turbulent flow in non-circular ducts is a more challenging task for numerical algorithms than the canonical cases of plane channel and pipe flow. The main reason resides in the availability of only one direction of space homogeneity, which prevents the use of efficient inversion procedures for Poisson equations based on double trigonometric expansions (Kim & Moin 1985; Orlandi 2012). This difficulty has been tackled using two-dimensional Poisson solvers based on cyclic reduction (Gavrilakis 1992), using algebraic multigrid methods (Vinuesa et al. 2014; Marin et al. 2016), or fast diagonalization (Pinelli et al. 2010), resulting in a larger computational cost than for channel and pipe flow simulations. In the present work, we use a fourth-order co-located finite-difference solver, previously used for DNS of compressible turbulence, also in the low-Mach-number regime (Pirozzoli & Bernardini 2013; Modesti & Pirozzoli 2016; Bernardini et al. 2021). Here, the convective terms in the Navier–Stokes equations are expanded preliminarily to quasi-skew-symmetric form, in such a way as to preserve discretely total kinetic energy from convection (Pirozzoli 2010). Semi-implicit time stepping is used for time advancement in order to relax the acoustic time step limitation, thus allowing efficient operation at low Mach number, also through the use of the entropy evolution equation rather than the total energy equation (Modesti & Pirozzoli 2018). The streamwise momentum equation is forced in such a way as to maintain a constant mass flow rate (the spatially uniform driving term is hereafter referred to as ${\varPi }$), periodicity is exploited in the streamwise direction, and isothermal no-slip boundary conditions are used at the channel walls. Let $h$ be the duct half-side; the DNS have been carried out for a duct with $[-h:h] \times [-h:h]$ cross-section, and length $6 {\rm \pi}h$.

Six DNS have been carried out at bulk Mach number $M_b=u_b/c_w=0.2$ (where $c_w$ is the speed of sound at the wall temperature) and bulk Reynolds number ${{Re}}_b=4400\text {--}84\,000$ (table 1), and hereafter these are labelled with letters from A to F. In order to quantify the effect of secondary flows we have also carried out simulations at the same bulk Reynolds numbers, wherein secondary flows are suppressed numerically, and which are denoted with the suffix 0. The turbulence Mach number $M_t=u'/{c_w}$ nowhere exceeds $0.01$ for any of the simulations, hence the present DNS may be regarded as representative of genuinely incompressible turbulence. Issues related to the size of the computational box, mesh resolution and statistical convergence were discussed in a previous publication (Pirozzoli et al. 2018), which focused on the analysis of the velocity statistics. In the present study, the Navier–Stokes equations are augmented with the transport equation for a passive scalar field, with molecular diffusivity $\alpha = \nu$, in such a way that the molecular Prandtl number is unity for all simulations. Although the study of passive scalars may be relevant in several fluid dynamics applications, from now on we will always refer to the transported scalar as representative of the temperature field, and the associated fluxes as representative of heat fluxes. Similarly to the streamwise velocity field, the passive scalar equation is also forced with a time-varying, spatially homogeneous forcing term (say $Q$), in such a way that its mean value is maintained in time (Pirozzoli, Bernardini & Orlandi 2016). This forcing method corresponds to adding uniform bulk heating to balance heat losses, although other forcing strategies are possible (e.g. Piller 2005; Abe & Antonia 2017). The most common alternative to uniform bulk heating is enforcement of constant heat flux in time. Both approaches have advantages and disadvantages. Uniform bulk heating is less realistic as it is difficult to attain in experiments (Piller 2005); however, it is more efficient from a computational point of view because of faster convergence towards a statistically stationary state. In the present work, we opt for uniform bulk heating in light of higher computational efficiency in DNS of high-Reynolds-number flows. We expect that this approach may yield a slightly higher heat flux number as compared to the case of strictly constant heat flux (Abe & Antonia 2017; Alcántara-Ávila, Hoyas & Pérez-Quiles 2021). In the following, we use a different notation for the mean temperature in the duct,

(2.1)\begin{equation} \theta_b = \frac 1{A_c} \int_{A_c} (\varTheta-\theta_w)\,\mbox{d} A_c , \end{equation}

and for the bulk temperature,

(2.2a,b) \begin{equation} \theta_m = \frac{1}{u_b A_c} \int_{A_c} U (\varTheta-\theta_w)\,\mbox{d} A_c, \quad u_b = \frac{ 1}{A_c} \int_{A_c} U\,\mbox{d} A_c, \end{equation}

where $\theta _w$ is the wall temperature, and $A_c$ is the cross-sectional area of the duct. We use capital letters to denote flow properties averaged in the homogeneous spatial directions and in time, angle brackets to denote the averaging operator, and lower-case letters to denote fluctuations from the mean. We also exploit geometrical symmetries, and average the flow statistics over the duct quadrant or duct octant, whenever possible. The effect of quadrant averaging was discussed in detail in our previous publication Pirozzoli et al. (2018).

Table 1. Flow parameters for square duct DNS. Box dimensions are $6{\rm \pi} h \times 2h \times 2h$ for all flow cases. Flow cases denoted with the 0 suffix are carried out by suppressing the secondary motions. We denote by ${{Re}}_b = 2 h u_b / \nu$ the bulk Reynolds number, and by ${{Re}}_{\tau }^* = h u_{\tau }^* / \nu$ the friction Reynolds number; $C_f = 2 \tau _w /(\rho u_b^2)$ is the friction coefficient, and ${{Nu}}$ is the Nusselt number; $\Delta x$ is the mesh spacing in the streamwise direction, and $\Delta z$, $\Delta y_w$ are the maximum and minimum mesh spacings in the cross-stream direction, all given in global wall units $\delta _v^*=\nu /u_{\tau }^*$.

The $+$ superscript is here used to denote local wall units, namely quantities made non-dimensional with respect to the local friction velocity $u_{\tau } = (\tau _w/\rho )^{1/2}$ (where $\tau _w= \rho \nu \,{\mbox {d} U}/{\mbox {d} y}\vert _w$ is the local wall shear stress), and the local viscous length scale $\delta _v=\nu /u_{\tau }$. The $*$ superscript is reserved to denote global wall units, with friction velocity based on the perimeter-averaged wall shear stress,

(2.3)\begin{equation} u^*_{\tau} = \left( \frac {\tau^*_w}{\rho} \right)^{1/2} = \left( \frac { h \left\langle \varPi \right\rangle }{2 \rho} \right)^{1/2}, \end{equation}

and viscous length scale $\delta _v^*=\nu /u_{\tau }^*$. Likewise, for normalization of the temperature field, we consider either the local friction temperature,

(2.4)\begin{equation} \theta_{\tau} = \frac {q_w}{u_{\tau}}, \end{equation}

where $q_w = \alpha \,{\mbox {d} {\varTheta }}/{\mbox {d} y}\vert _w$ is the local wall heat flux, or the global friction temperature,

(2.5)\begin{equation} \theta^*_{\tau} = \frac {q^*_w}{u^*_{\tau}} = \frac{ \left\langle {Q} \right\rangle h}{2 u^*_{\tau}}. \end{equation}

Equations (2.3) and (2.5) make it clear that the global friction velocity and temperature are related to the imposed forcing of the streamwise momentum and temperature equations, rather than to the detailed distribution of wall momentum and temperature fluxes. Notably, local and global wall units collapse in the case of canonical wall-bounded flows because of homogeneity, but the two are distinctly different in internal flows in non-circular ducts. DNS studies such as the present one can then be used to distinguish between the two scalings, which is impossible in more canonical set-ups.

3. Temperature statistics

We begin by inspecting the instantaneous streamwise velocity and temperature fields in figure 1, which shows the flow in both the cross-stream and wall-parallel planes. The instantaneous flow structures are not different from those typical of canonical wall-bounded flows, and the near-wall region is populated with temperature and streamwise velocity streaks, resulting from sweep and ejections events visible in the cross-stream plane. Large structures are also visible in the cross-stream plane, which convey high-speed, hot flow from the duct core towards the corners, and which interact with the smaller near-wall flow structures. Return motions of low-speed, cool fluid are also observed to emerge from around the bisector of each side of the square duct. These observations suggest that the secondary motions resulting from space- and time-averaging are not just an artefact, but they are well present in the instantaneous flow realizations. Similarity between streamwise velocity and temperature is apparent, with the main difference that the latter reveals finer structures and sharper fronts owing to the absence of the smoothing action of pressure (Pirozzoli et al. 2016).

Figure 1. Instantaneous streamwise velocity (a) and temperature (b) fields for flow case F. Wall-parallel planes are taken at fifteen wall units from the walls.

Figure 2 depicts the mean temperature (contours) and velocity (isolines) fields in the duct cross-section, along with representative cross-stream velocity vectors. This representation brings to light the presence of a pair of counter-rotating secondary eddies in each quarter of the domain, whose apparent role is bringing high-speed and high-temperature fluid from the duct core towards the corners, to compensate for the deficit. As a result, both the temperature contours and the velocity isolines bend towards the corners, featuring a bulging first noticed by Nikuradse (1930). This effect seems to be non-monotonic with the Reynolds number, being most evident for DNS-A and DNS-E/F. As expected, temperature and velocity contours bear close similarity, being nearly coincident near the wall, and retaining the same shape farther off. This analogy is the result of the formal similarity of the controlling equations, and it is further investigated in figure 3, where we show a scatter plot of mean velocity and temperature. In the near-wall region (say, $U^* \le 10$), perfect matching is recovered (namely, $\varTheta ^* = U^*$), which is not surprising, as at unit Prandtl number, $\varTheta ^+ \equiv U^+ \equiv y^+$ in the wall vicinity. However, near equality of the two distributions in global wall units is less trivial and is a symptom of strong similarity between the two fields. At non-unit Prandtl number, we still expect a linear relationship to be present like $\varTheta ^* = {{Pr}} \, U^*$, since in general $\varTheta ^+ \equiv U^+ \equiv {{Pr}} \, y^+$ in the wall vicinity (Kader 1981). Farther from walls ($U^* \gtrsim 15$, which corresponds to the root of the logarithmic layer in canonical wall-bounded flows), a linear imprinting is still visible, although with larger scatter than in the near-wall region, and with slope shallower than unity. This can be explained by recalling (Pirozzoli et al. 2016) that in canonical flows, the distributions of both velocity and temperature are nearly logarithmic, but with different slope. In particular, currently accepted values of the Kármán constant are $\kappa \approx 0.39$ for the mean velocity, and $\kappa _{\theta } \approx 0.46$ for the mean temperature. As a matter of fact, the scatter plot is well fitted with a linear function with slope $\kappa /\kappa _{\theta } \approx 0.84$. At ${{Pr}} \ne 1$, we thus expect a similar relationship, but with different additive constant, as a result of the different intercept in the log law for the temperature field. This observation is quite interesting in our opinion, as the presence of a (nearly) universal relationship between mean velocity and temperature in a non-trivial flow in a square duct in principle would allow us to reconstruct directly the whole temperature field based on the velocity field.

Figure 2. Mean streamwise temperature (${\varTheta }^*$, flooded contours), mean streamwise velocity (${U}^*$, lines), and mean cross-stream velocity vectors ($V^*, W^*$) for flow cases A (a), B (b), C (c), D (d), E (e), F ( f). For clarity, only a subset of the velocity vectors is shown, and only a quarter of the full domain is shown. The dashed diagonal lines indicate the corner bisector.

Figure 3. Scatter plot of mean temperature ($\varTheta ^*$) versus mean velocity ($U^*$) for DNS data of all flow cases, and corresponding fitting functions for the near-wall region, $\varTheta ^* = U^*$ (dashed line), and for the outer region, $\varTheta ^* = 2.57 + 0.84 U^*$ (solid line). Colours are as in table 1.

Based on the previous observations, it is natural to study and compare the temperature statistics expressed in local wall units ($+$) and in global wall units ($*$). Figure 4 shows the mean temperature profiles as a function of the wall-normal distance up to the corner bisector (white dashed line in figure 2), in local wall units. For reference purposes, the mean temperature profiles obtained from the experimental fitting by Kader (1981) at matching ${{Re}}_{\tau }$ are also reported. Excellent collapse of the locally scaled profiles is recovered near the wall, also including the near-corner region. The distributions become more widespread past $y^+ \approx 10$, with maximum scatter of about $\pm$5 % in the ‘logarithmic’ layer. When the global friction temperature is used for normalization (see figure 5), scatter among the temperature profiles is observed near the wall as a result of the local variation of the wall heat flux. Perhaps unexpectedly, this normalization does yield similarly good universality away from walls, and near coincidence with the pipe temperature profiles, at least at high enough Reynolds numbers. This finding is probably related to the fact that mean temperature transport away from solid walls is controlled by the imposed spatially uniform heat source rather than by the non-uniform wall heat flux. Inspection of figure 5 shows that transition from wall scaling to ‘global’ scaling (which is controlled by the spatially uniform heat source) occurs at a wall distance of about $0.2 h$, which is also the lower limit for the core region in canonical flows (Pope 2000). This observation seems to support the validity of Townsend's outer layer similarity also for the temperature field, in that the outer flow perceives mainly the global heat flux and not its detailed local distribution at the wall. Validity of outer-layer similarity for the temperature field was reported for rough walls (MacDonald, Hutchins & Chung 2019), whereas we are not aware of similar conclusions for the case of more complex geometries because most studies on forced convection are limited to canonical wall-bounded flows, in which the wall heat flux is spatially uniform. These results support findings reported previously for the streamwise velocity (Pirozzoli et al. 2018), and confirm that square duct flow is a convenient testbed for evaluating differences between local and global scaling.

Figure 4. Mean wall-normal temperature profiles scaled in local wall units for flow cases A (a), B (b), C (c), D (d), E (e), F ( f). Profiles are plotted at several distances from the left wall, up to the corner bisector (see figure 2 for reference): $z^*=15$ (diamonds), $z/h=0.1$ (right triangles), $z/h=0.25$ (triangles), $z/h=0.5$ (circles), $z/h=1$ (squares). The dashed lines denote fits of experimental data (Kader 1981) at matching ${{Re}}_\tau ^*$.

Figure 5. Mean wall-normal temperature profiles scaled in global wall units for flow cases A (a), B (b), C (c), D (d), E (e), F ( f). Profiles are plotted at several distances from the left wall, up to the corner bisector (see figure 2 for reference): $z^*=15$ (diamonds), $z/h=0.1$ (right triangles), $z/h=0.25$ (triangles), $z/h=0.5$ (circles), $z/h=1$ (squares). The dashed lines denote fits of experimental data (Kader 1981) at matching ${{Re}}_\tau ^*$.

The fluctuating velocity and temperature fields in the vicinity of the duct corners are shown in figure 6 upon global inner normalizations, as well as the corresponding correlation coefficients. To facilitate comparison across the Reynolds number range, the horizontal coordinate is reported in inner units up to $z^*=150$, and the vertical one in outer units. This allows simultaneous visualization of both the corner and the core region, at all Reynolds numbers. Far from the corner, the distributions of the velocity and temperature variances are nearly unaffected by Reynolds number variation, and feature a buffer-layer peak at a distance of about fifteen wall units from the closer wall. The amplitude of the peak is reduced moving towards the corner, whose effect is felt to within a distance of $O(100)$ global wall units. The $u$–$\theta$ correlation coefficient shows strict similarity of velocity and temperature fluctuations also in the instantaneous sense, with the exception of the close neighbourhood of the corner. Correlation in the outer layer is less than near the wall, but it tends to increase at higher Reynolds number, similar to the case of plane channel flow (Pirozzoli et al. 2016).

Figure 6. Streamwise velocity and temperature fluctuation intensities (in global inner units), and their correlation coefficient, in the proximity of a duct corner for flow cases A (a–c), B (d–f), C (g–i), D (j–l), E (m–o), F (p–r). Note that the horizontal coordinate ($z$) is given in inner units, and the vertical coordinate $y$ is given in outer units and at a logarithmic scale.

A quantitative analysis of the buffer-layer peaks of the velocity and temperature variances is reported in figure 7. Specifically, for each horizontal location ($z$, measured from the corner), we show the peak variances of $u$ and $\theta$ expressed in local wall units. In this respect, we note that $z^+$ can be interpreted as effective local friction Reynolds numbers, based on the distance of the wall point from the corner bisector (see figure 2). It is known that in canonical flows (Marusic et al. 2013), the peak velocity variances increase logarithmically with the friction Reynolds number. Hence it is interesting to verify whether the same trend also occurs in flows involving non-trivial cross-sections, and featuring secondary motions. For that purpose, in both panels of figure 7 we show the logarithmic trend predicted by Marusic et al. (2013) as benchmark. It is noteworthy that away from the corner (say $z^+ \gtrsim 100$), local wall units yield good universality of the distributions, and suggest an increasing trend of the variances with the local friction Reynolds number. However, the growth rate seems to be far less than in canonical flows, possibly violating Townsend's attached-eddy hypothesis. We speculate that this peculiar behaviour may be due to either disruption of the hierarchy of wall-attached eddies because of the action of secondary motions, or perhaps impeded spanwise growth caused by the presence of corners, both effects being absent in planar channel or pipe flow.

Figure 7. Inner-layer peak variances of axial velocity (a), and temperature (b), expressed as functions of the local friction Reynolds number ($z^+$). The dashed lines denote the logarithmic fit for canonical flows (Marusic, Baars & Hutchins 2017), namely $\langle u^2\rangle ^+_{{IP}} = 0.63 \log {{Re}}_{\tau } + 3.8$. Symbols as in table 1.

The distributions of the local friction and heat transfer coefficients along the duct wall are shown in figure 8. As shown previously by Pirozzoli et al. (2018), the distribution of the wall friction coefficients is non-monotonic, featuring a peak at the side bisector, and a secondary peak near the corner, at a distance scaling with global wall units. As from theoretical predictions (Spalart, Garbaruk & Stabnikov 2018), the friction tends to equalize over the duct perimeter at higher ${{Re}}$. The heat transfer coefficient seems to have the same qualitative behaviour, thus corroborating the experimental observations of Brundrett & Burroughs (1967) and Hirota et al. (1997). Quantitative assessment is shown in figure 8(c), where we show the Reynolds analogy factor $2 {{St}} / C_f$, which is expected to be unity at ${{Pr}} = 1$ (White & Majdalani 2006). In fact, whereas deviations are observed at low Reynolds numbers, the distribution corresponding to DNS-F varies between 0.98 and 1.2, right at the corner.

Figure 8. Distributions of local wall friction coefficient $C_f=2 u_{\tau }^2/u_b^2$ (a), local Stanton number ${{St}}=(u_{\tau } \theta _{\tau })/(u_b \theta _m)$ (b), and Reynolds analogy factor (c), along the left part of the lower wall. Symbols as in table 1.

The overall heat transfer performance of the duct is quantified in terms of the global Stanton number,

(3.1)\begin{equation} {{St}}^* = \frac{q_w^*}{u_b \theta_m} = \frac 1{u_b^* \theta_m^*}, \end{equation}

or more frequently in terms of the Nusselt number,

(3.2)\begin{equation} {{Nu}}^* = {{St}}^* \, {{Re}}_b \, {{Pr}} . \end{equation}

In figure 9, we compare the distributions of the heat transfer coefficients (inverse Stanton number and Nusselt number) obtained from the present DNS (squares) with those resulting from DNS of circular pipe flow (circles), from previous LES (diamonds), and from experiments (triangles) in square ducts. We note that the latter two datasets were obtained for ${{Pr}}=0.71$, hence the ${{Nu}}$ data are rescaled by a factor $1/\sqrt {{{Pr}}}$ to compare with the present ones. As a reference, we also report correlations used widely in engineering practice, including that from Gnielinski (1976),

(3.3) \begin{equation} {{Nu}} = \frac{{{Pr}}\,(C_f/32) ( {{Re}}_b - 1000 )}{1 + 12.7 (C_f/32)^{1/2} ( {{Pr}}^{2/3} - 1)}, \end{equation}

and from Kays & Crawford (1993),

(3.4)\begin{equation} {{Nu}} = 0.022 \, {{Re}}_b^{0.8} \, {{Pr}}^{0.5}. \end{equation}

Finally, direct fitting the present DNS data (at ${{Pr}}=1$) with a power-law expression yields

(3.5)\begin{equation} {{Nu}} = 0.0216 \, {{Re}}_b^{0.805}. \end{equation}

Figure 9. Inverse global Stanton number (a) and Nusselt number (b) as functions of bulk Reynolds number, from the present DNS (squares), from DNS of circular pipe flow (circles, Pirozzoli et al. 2021), from previous LES (diamonds, Pallares & Davidson 2002), and from experiments (triangles, Brundrett & Burroughs 1967). The dashed lines denote the correlation (3.4), and the dot-dashed lines denote the correlation (3.3).

The DNS data show good correspondence with the pipe DNS data at matching ${{Re}}_b$, with the exception of the lower Reynolds number case, which supports validity of the hydraulic diameter concept as the relevant length scale for heat transfer data reduction. Agreement with experiments and LES in square ducts is also quite good, on account of experimental uncertainties and modelling errors in LES. It is interesting that differences are levelled off when the popular representation in terms of the Nusselt number is used, as in figure 9(b), hence we believe that the $1/{{St}}$ representation should be used when relatively small differences must be discriminated. Classical correlations seem to suggest systematic difference of up to $5\,\%$ in the prediction of the heat transfer coefficient. This difference may be due partly to inaccuracy of correlations based on old experimental data, or to the fact that those are trained mainly for the ${{Pr}}=0.71$ case, whereas here ${{Pr}}=1$. Discrepancies can also be attributed partly to our heat forcing scheme, in which a spatially uniform heat source is prescribed, which tends to slightly overpredict the heat flux as compared to other approaches in which the total heat flux is maintained strictly constant in time (Abe & Antonia 2017; Alcántara-Ávila et al. 2021). A slight adjustment of the Kays–Crawford power-law formula coefficients as after (3.5) seems to yield very good representation of the DNS data.

An issue that deserves further investigation is why use of the hydraulic diameter yields excellent results, at least in the case of square ducts under consideration. Pirozzoli et al. (2018) showed that universality of friction is connected with near applicability of the logarithmic velocity law in the direction normal to the nearest wall. It is worthwhile verifying whether this is also the case for the heat flux coefficient. As shown previously, the temperature distributions are nearly universal away from walls, even when expressed in global wall units. Hence, approximating the outer layer profiles with the classical log law, namely $U^*=1/{\kappa } \log y^* + C$, $\varTheta ^*=1/\kappa _{\theta } \log y^* + C_{\theta }$, and integrating over the duct cross-section, gives the following expression for the inverse Stanton number velocity:

(3.6)\begin{align} 1 / {{St}}^* &= u_b^* \theta_m^* = \frac{8}{4{{Re}}_\tau^{*2}} \int_0^h\int_0^z\left[\frac{1}{\kappa \kappa_\theta}\log{y^*}^2 + \log{y^*} \left(\frac{C_\theta}{\kappa} + \frac{C}{\kappa_\theta}\right) +CC_\theta\right]{{\rm d}y}^*\,{\rm d} z^* \nonumber\\ &=\frac{1}{\kappa \kappa_\theta}\left(\log^2{{{Re}}_\tau^*} + 3\log{{{Re}}_\tau^*} + \frac{7}{2}\right) + \left(\frac{C_\theta}{\kappa}+\frac{C}{\kappa_\theta}\right) \left(\log{{{Re}}_\tau^*}-\frac{3}{2}\right) + CC_\theta. \end{align}

Equation (3.6) should be compared with the corresponding expression for a circular duct with diameter $D$,

(3.7)\begin{equation} 1 / {{St}} = \frac{1}{\kappa \kappa_\theta}\left(\log^2{{{Re}}_\tau} + 3\log{{{Re}}_\tau} + \frac{7}{2}\right) + \left(\frac{C_\theta}{\kappa}+\frac{C}{\kappa_\theta}\right) \left(\log{{{Re}}_\tau}-\frac{3}{2}\right) + CC_\theta, \end{equation}

where ${{Re}}_{\tau } = D u_{\tau } / (2 \nu )$. The two expressions are identical for $2 h = D$, hence the Reynolds number based on the hydraulic diameter is the same. It is interesting that (3.6) is arrived at basically by neglecting the local wall shear stress and heat flux variation along the duct perimeter, and disregarding the flow deceleration at corners. Apparently, these effects very nearly cancel out upon integration.

4. Contribution of secondary flows to heat transfer

In order to quantify the effects of secondary motions on heat transfer, we derive a generalized version of the Fukagata–Iwamoto–Kasagi (FIK) identity (Fukagata, Iwamoto & Kasagi 2002) for the Stanton number, by following the approach of Modesti et al. (2018). We consider the mean temperature balance equation,

(4.1)\begin{equation} \alpha \nabla^2 \varTheta = \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{q}_C + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}_T - \langle{Q}\rangle, \end{equation}

where $\boldsymbol {q}_C={\varTheta }\,{\boldsymbol {U}}_{yz}$ is associated with mean cross-stream convection (hence, with the secondary motions), $\boldsymbol {q}_T={\langle \theta \boldsymbol {u}_{yz} \rangle }$ is associated with turbulence convection, $\boldsymbol {u}_{yz}=({v},{w})$ is the cross-stream velocity vector, and $\langle {Q}\rangle$ is the mean driving bulk heating.

Equation (4.1) may be interpreted as a Poisson equation for the mean temperature, with source terms $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {q}_T$, $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {q}_C$ and $\langle Q\rangle$ obtained from the DNS dataset. Hence the solution of (4.1) may be cast formally as the superposition of three parts, namely ${\varTheta }={\varTheta }_D+{\varTheta }_T+{\varTheta }_C$, and these temperature fields can be obtained by solving three Poisson equations,

(4.2a–c)\begin{equation} \alpha \nabla^2{\varTheta}_D={-}\langle{Q}\rangle,\quad \alpha \nabla^2{\varTheta}_T=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}_T,\quad \alpha \nabla^2{\varTheta}_C=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}_C, \end{equation}

with homogeneous boundary conditions, where ${\varTheta }_D$, ${\varTheta }_T$ and ${\varTheta }_C$ denote the diffusion, turbulence and convection contributions to the mean temperature field. Accordingly, the mean temperature in the duct may be evaluated as

(4.3a,b)\begin{equation} \theta_b=\frac {\langle{Q}\rangle A}{\alpha}\,{\theta_b}_1+{\theta_b}_T+{\theta_b}_C, \quad {\theta_b}_X = \frac{1}{A} \int_A {\varTheta}_X\,\mathrm{d}A, \end{equation}

where we have introduced the unitary Stokes temperature field $\theta _1$, defined as the solution of $\nabla ^2 \theta _1 = -1/A$, which by construction is a function of only the duct geometry, and which allows us to express the diffusion-associated temperature field as $\varTheta _D= A \, \theta _1 \langle {Q}\rangle / \alpha$. For convenience, in the present analysis we use a modified form of the Stanton number, based on the mean temperature rather than on the bulk temperature, namely

(4.4)\begin{equation} {{St}}_b^* = \frac{q_w^*}{u_b \theta_b}= \frac 1{u_b^* \theta_b^*}, \end{equation}

instead of (3.1), which when substituted into (4.3a,b) yields

(4.5)\begin{equation} {{St}}_b^* = \frac{1}{{{Pr}}\,{{Re}}_P\,{\theta_b}_1}\left(1-\frac{{\theta_b}_T}{\theta_b} - \frac{{\theta_b}_C}{\theta_b} \right)={{St}}_D + {{St}}_T + {{St}}_C, \end{equation}

where ${{Re}}_P=u_b P / \nu$ is the bulk Reynolds number based on the duct perimeter $P$. Equation (4.5) shows clearly that the heat transfer coefficient may be regarded as the sum of contributions due to diffusion, turbulence and convection (labelled as ${St}_X$ in (4.5)), and it yields an expression similar to that from the classical FIK identity for the friction coefficient.

We note that this formalism has additional generality compared to the classical FIK identity and its extensions (Peet & Sagaut 2009; Jelly, Jung & Zaki 2014), and it allows us to isolate the effects of mean convection, turbulence and diffusion terms to the wall shear stress distributions along the duct perimeter, resulting from

(4.6a–c)\begin{equation} {q_w}_D = \alpha \left.\frac{\partial {\varTheta}_D}{\partial y}\right\rvert_w, \quad {q_w}_T = \alpha \left.\frac{\partial {\varTheta}_T}{\partial y}\right\rvert_w, \quad {q_w}_C = \alpha \left.\frac{\partial {\varTheta}_C}{\partial y}\right\rvert_w. \end{equation}

Regarding (4.6a–c), it is important to note that only ${q_w}_D$ has non-zero mean, whereas integration of (4.2a–c) for ${\varTheta }_C$ and ${\varTheta }_T$ shows readily that their integrated contributions vanish, as $\boldsymbol {q}_C$ and $\boldsymbol {q}_T$ are both zero at walls. The contributions to the mean heat transfer coefficients are given in table 2, both in absolute terms and as a fraction of the total. Consistent with physical expectations, the diffusion contribution is observed to decline at increasing Reynolds number with respect to the turbulent term. The contribution of cross-stream convection is found to be roughly constant across the explored Reynolds number range; however, it remains much less than the turbulence contribution. Extrapolating the DNS data, we expect the convection contribution to exceed the diffusion one at high enough Reynolds number. Figure 10 shows the mean temperature fields associated with the diffusion, turbulence and convection terms. For the sake of clarity, each contribution is normalized with the corresponding bulk value, ${\theta _b}_X$, where ${\theta _b}_T$ and ${\theta _b}_C$ are both negative, thus providing an additive effect to the heat transfer (see (4.5)). The diffusion-associated temperature field $\varTheta _V$ arises from the solution of a Poisson equation with uniform right-hand side (the heat source), hence its shape is identical to the case of laminar flow (Shah & London 2014), depending only on the duct cross-sectional geometry. The turbulence-associated temperature field $\varTheta _T$ is everywhere negative and topologically similar to the diffusion-associated field, highlighting the cooling effect of turbulence on the bulk flow, with incurred increase of the heat transfer coefficient. The temperature field $\varTheta _C$ associated with mean convection has a more complex organization. Positive values are found near the duct corners, whereas negative values are found near the duct bisectors, the zero crossings being located halfway in between. This finding is consistent with the intuitive expectation that secondary motions tend to equalize temperature across the duct cross-section, thus quantitatively corroborating claims made by early investigators (Prandtl 1927). Notably, all the distributions shown in figure 10 are not qualitatively affected by Reynolds number variation when scaled with respect to their mean integral value, thus showing that changing the Reynolds number changes only the relative importance of the three terms, as quantified in table 2.

Figure 10. Contributions of diffusion, turbulence and mean convection terms to the mean temperature field, as determined by solving (4.2a–c), for flow cases A (a–c), B (d–f), C (g–i), D (j–l), E (m–o), F (p–r).

Table 2. Contributions of diffusion, turbulence and mean convection terms to the heat transfer coefficient, as from (4.5).

The observations made regarding the organization of the temperature field have a direct impact on the distribution of the local wall heat flux, shown in figure 11. The heat flux distribution induced by diffusion terms is the same for all cases, and nearly parabolic in shape. The turbulence terms have a complex behaviour, exhibiting multiple peaks that change with the Reynolds number. In general, they yield heat flux increase toward the duct bisectors, and slight attenuation at the corners (compare the lines with circles and with dots), thus making the heat flux distributions more non-uniform. Mean convection yields large positive peaks in the corner vicinity, whose distance from the neighbouring walls scales in inner units (Pirozzoli et al. 2018), and negative values around the duct bisector. As a result, the distributions of the heat flux tend to be rather flat, especially at higher Reynolds numbers.

Figure 11. Contributions to mean wall heat flux (shown along half a duct side) for flow cases A (a), B (b), C (c), D (d), E (e), F ( f): diffusive (${q_w}_D(z)$, dotted), turbulent (${q_w}_T(z)$, dashed), diffusive + turbulent (${q_w}_D(z)+{q_w}_T(z)$, circles), mean convection (${q_w}_C(z)$, dash-dotted), and total ($q_w(z)$, solid line).

As an alternative to quantify the effect of secondary flows we have also considered the approach of Modesti et al. (2018), and carried out numerical experiments at the same bulk Reynolds numbers as the baseline cases, whereby the secondary motions are suppressed artificially (flow cases with the 0 suffix in table 1). For that purpose, we force the streamwise-averaged cross-stream velocity components to have zero mean by setting

(4.7a)$$\begin{gather} v(x,y,z,t) \rightarrow v(x,y,z,t) - \overline{v}^x(y,z,t), \end{gather}$$

(4.7b)$$\begin{gather}w(x,y,z,t) \rightarrow w(x,y,z,t) - \overline{w}^x(y,z,t), \end{gather}$$

at each Runge–Kutta sub-step, where $\overline {(.)}^x$ denotes the streamwise averaging operator. Figure 12 shows the temperature contours with superimposed streamwise velocity isolines for flow cases A0–D0. We note that the typical bulging of velocity and temperature towards the corners is absent for these cases, and the core region shows near axial symmetry as in circular pipe flow. The mean velocity and temperature fields are nearly coincident at low Reynolds numbers, whereas differences become more evident as ${{Re}}_b$ increases. Although the mean flow fields in figure 12 obviously differ from the full DNS cases shown in figure 2, the heat transfer is found to be barely affected. In fact, similar to what is reported by Modesti et al. (2018) for frictional drag, suppressing the secondary flows yields even higher heat transfer at low Reynolds numbers, whereas a bare $3\,\%$ decrease of the Nusselt number is found for flow case $D0$; see table 1. These figures are consistent with those suggested by the FIK identity, which actually returns about $3\,\%$ contribution of secondary flows to global heat transfer for flow case D.

Figure 12. Flow cases without secondary motions: mean streamwise velocity (${U}^*$, lines), and mean temperature (${\varTheta }^*$, flooded contours), for flow cases A0 (a), B0 (b), C0 (c), D0 (d). Only a quarter of the full domain is shown. The dashed diagonal lines indicate the corner bisector.

5. Modelling the turbulent heat flux

The numerical solution of the Reynolds-averaged Navier–Stokes equations is a common approach for flows of industrial interest. Solution of the temperature field requires a suitable closure model for the turbulent heat flux vector, $\boldsymbol {q}=\langle {\boldsymbol {u}\theta }\rangle$. Such closures are typically less sophisticated than those developed for momentum flux, and generally rely on validity of the Reynolds analogy. The primary outcome of this school of thought is the popularity retained over the years by the turbulent Prandtl number,

(5.1)\begin{equation} {{Pr}}_t=\frac{\nu_t}{\alpha_t}, \end{equation}

where $\nu _t$ and $\alpha _t$ are the eddy viscosity and eddy diffusivity coefficients, respectively. These coefficients are the heart of the linear eddy-viscosity hypothesis and its heat transfer counterpart,

(5.2a,b)\begin{equation} \tilde{\boldsymbol{a}} ={-}\nu_t \left( \boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{\nabla} \boldsymbol{U}^{\rm T} \right), \quad \tilde{\boldsymbol{q}} ={-}\alpha_t \boldsymbol{\nabla} \varTheta, \end{equation}

where $\tilde {\boldsymbol {a}}$ is the modelled counterpart of the anisotropic Reynolds stress tensor ($\boldsymbol {a}=\boldsymbol {\tau }-2/3 k$), and $\tilde {\boldsymbol {q}}$ is the modelled turbulent heat flux vector, with $\boldsymbol {\tau }=\langle {\boldsymbol u \boldsymbol u}\rangle$ and $k=1/2 \, \mathrm {tr} (\boldsymbol {\tau })$. Most heat transfer models assume a uniform turbulent Prandtl number (${{Pr}}_t \approx 1$) to avoid solving ad hoc transport equations for $\alpha _t$, although this hypothesis is not always accurate (Pirozzoli et al. 2016; Kaller et al. 2019). Even assuming that the approximation ${{Pr}}_t \approx 1$ holds, limitations of the linear eddy-viscosity hypothesis (5.2a,b) are evident, and reliance on the turbulent Prandtl number concept replicates modelling inaccuracies. Pope (1975) extended the eddy-viscosity ansatz by using a tensor polynomial to model $\tilde {\boldsymbol {a}}$. This approach is mathematically sound because it relies on the Cayley–Hamilton theorem. In fact, it turns out that five tensor bases are in most cases sufficient to recover an exact representation of $\tilde {\boldsymbol {a}}$ (Gatski & Speziale 1993; Jongen & Gatski 1998; Modesti 2020). The same mathematical framework was extended to scalar flux modelling, both for heat transfer (So, Jin & Gatski 2004a; Younis, Speziale & Clark 2005; Wang et al. 2007) and for buoyancy (So, Jin & Gatski 2004b). Similar to the original tensor polynomial expansion of Pope (1975), the turbulent heat flux is expanded as

(5.3)\begin{equation} \tilde{\boldsymbol{q}} = \sum_{n=1}^{N}\alpha^{(n)}\boldsymbol{V}^{(n)},\quad N=1, \ldots, 10, \end{equation}

where the coefficients $\alpha ^{(n)}$ are scalar functions, and $\boldsymbol {V}^{(n)}$ is a set of ten vector bases,

(5.4) \begin{equation} \left.\begin{array}{lll@{}} &\boldsymbol{V}^{(1)}= \boldsymbol\nabla\varTheta, &\qquad\quad \boldsymbol{V}^{(6)} = \boldsymbol{\varOmega}\boldsymbol\nabla\varTheta,\\ &\boldsymbol{V}^{(2)}= \boldsymbol{a}\boldsymbol\nabla\varTheta, &\qquad\quad \boldsymbol{V}^{(7)} = \boldsymbol{\varOmega}\boldsymbol{\varOmega}\boldsymbol\nabla\varTheta,\\ &\boldsymbol{V}^{(3)}= \boldsymbol{S}\boldsymbol\nabla\varTheta, &\qquad\quad \boldsymbol{V}^{(8)} = \left(\boldsymbol{S}\boldsymbol{\varOmega}+\boldsymbol{\varOmega}\boldsymbol{S}\right)\boldsymbol\nabla\varTheta,\\ &\boldsymbol{V}^{(4)}= \boldsymbol{a}\boldsymbol{a}\boldsymbol\nabla{\varTheta}, &\qquad\quad \boldsymbol{V}^{(9)} = \left(\boldsymbol{a}\boldsymbol{S}+\boldsymbol{S}\boldsymbol{a}\right)\boldsymbol\nabla\varTheta,\\ &\boldsymbol{V}^{(5)}= \boldsymbol{S}\boldsymbol{S}\boldsymbol\nabla{\varTheta}, &\qquad\quad \boldsymbol{V}^{(10)} = \left(\boldsymbol{\varOmega}\boldsymbol{S}+\boldsymbol{S}\boldsymbol{\varOmega}\right)\boldsymbol\nabla\varTheta, \end{array}\right\} \end{equation}

where $\boldsymbol {S} = 1/2 ( \boldsymbol {\nabla } \boldsymbol {U} + \boldsymbol {\nabla } \boldsymbol {U}^\textrm {T} )$ and $\boldsymbol {\varOmega } = 1/2 ( \boldsymbol {\nabla } \boldsymbol {U} - \boldsymbol {\nabla } \boldsymbol {U}^\textrm {T} )$ are the mean strain-rate and rotation-rate tensors, respectively. As for the Reynolds stress tensor, in most practical flow cases, the ten bases are not independent, and an exact representation of the turbulent heat flux ($\tilde {\boldsymbol {q}}=\boldsymbol {q}$) is recovered for $N=3$. In order to assess the accuracy of different truncations of the polynomial expansion (5.3), we compute the coefficients $\alpha ^{(n)}$ following the approach of Jongen & Gatski (1998) and Modesti (2020): namely, we take the scalar product of (5.3) with each vector basis $\boldsymbol {V}^{(m)}$, thus obtaining a linear system of $N$ equations for the unknown coefficients $\alpha ^{(n)}$:

(5.5)\begin{equation} \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{V}^{(m)} = \sum_{n=1}^{N} \alpha^{(n)}\boldsymbol{V}^{(n)}\boldsymbol{\cdot}\boldsymbol{V}^{(m)},\quad m=1,\ldots, N. \end{equation}

For $N=1$, this approach yields the standard definition of linear eddy diffusivity,

(5.6)\begin{equation} \alpha_t={-}\alpha^{(1)}={-}\frac{\boldsymbol{q}\boldsymbol{\cdot}\boldsymbol\nabla\varTheta}{\boldsymbol\nabla \varTheta\boldsymbol{\cdot}\boldsymbol\nabla\varTheta}, \end{equation}

in analogy with the definition of linear eddy viscosity,

(5.7)\begin{equation} \nu_t={-}\frac{\left\{\boldsymbol{a}\boldsymbol{S}\right\}}{2\left\{\boldsymbol{S}\boldsymbol{S}\right\}}, \end{equation}

where $\{\boldsymbol {AB}\} = A_{ik}B_{ki}$. In the more general case, one can solve (5.5) to find the coefficients $\alpha ^{(n)}$ for different $1\leq N \leq 3$, and then evaluate the modelled heat flux vector $\tilde {\boldsymbol {q}}$ from (5.3).

In order to assess the validity of analogy between turbulent momentum and heat transport, in figure 13 we compare the linear eddy viscosity with the linear eddy diffusivity estimated from DNS flow case F. Visual scrutiny suggests that the assumption ${{Pr}}_t\approx 1$ is not valid everywhere, and in fact contours of $\nu _t$ and $\alpha _t$ do not match. Furthermore, the eddy diffusivity seems to be more affected by the presence of secondary flows as its isolines protrude more markedly towards the corner and tend to be more parallel to the walls as compared to $\nu _t$.

Figure 13. Estimated fields of linear eddy viscosity $\nu _t$ as after (5.7) (a), and linear eddy diffusivity $\alpha _t$ as after (5.6) (b), for flow case F (${{Re}}_\tau =2000$). Only one quadrant of the duct is shown.

The spatial distribution of $\nu _t$ and $\alpha _t$ suggests that the value of the two scalars is controlled by the closest wall, and therefore in figure 14 we report bundles of $\nu _t$ and $\alpha _t$ profiles up to the corner bisector. Besides the obvious increase of $\nu _t$ and $\alpha _t$ with the Reynolds number, the figure shows that the bundles of $\alpha _t$ have more limited scatter than $\nu _t$, confirming that the simple picture in which square duct flow is regarded as the superposition of two independent walls is more accurate for temperature than for momentum transport.

Figure 14. Wall-normal distributions of linear eddy viscosity $\nu _t$ (a), and linear eddy diffusivity $\alpha _t$ (b), at all $z$ locations, up to the corner bisector ($y=z$, see figure 2) for flow cases E (orange lines) and F (black lines). The red dashed lines indicate predictions of the mixing length hypothesis $\kappa y^+$ (a), and $\kappa _\theta y^+$ (b).

This effect may be interpreted as due to the additional effective diffusivity brought about by the non-local action of the pressure gradient term in the streamwise momentum equation, which as previously noted tends to smoothen the velocity field. Likewise, it appears that it also has the effect of enhancing interactions between neighbouring walls. The red dashed lines in figure 14 show the result of the mixing length hypothesis, namely $\nu _t/\nu \approx \kappa y^+$ and $\alpha _t/\alpha \approx \kappa _\theta y^+$, which is revealed to be only partially accurate in the overlap layer, especially for the eddy diffusivity.

Figure 15 shows bundles of the turbulent Prandtl number (see (5.1)) for flow case E and F, as a function of the wall distance. We also report predictions of available refined models for the turbulent Prandtl number, including that developed by Cebeci (1973),

(5.8a,b)\begin{equation} {{Pr}}_t = \frac{k}{\kappa_\theta}\left(\frac{1-\exp{({-}y^+{/}A)}}{1-\exp{({-}y^+{/}B)}}\right),\quad B=\frac{1}{{{Pr}}^{1/2}}\sum_{i=1}^5C_i\left(\log_{10}{{Pr}}\right)^{i-1}, \end{equation}

where $A=26$, $C_1=34.96$, $C_2=28.79$, $C_3=33.95$, $C_4=6.3$, $C_5=-1.186$, and that developed by Kays & Crawford (1993),

(5.9)\begin{equation} {{Pr}}_t=\left[\frac{1}{2{{Pr}}_{t,b}} + \frac{C\,{{Pe}}_t}{\sqrt{{{Pr}}_{t,b}}} - (C{{Pe}}_t)^2\left(1-\exp{\left(-\frac{1}{C\,{{Pe}}_t\,\sqrt{{{Pr}}_{t,b}}}\right)}\right)\right]^{{-}1}, \end{equation}

where ${{Pe}}_t={{Pr}} \, \nu _t/\nu$ is the turbulent Péclet number, ${{Pr}}_{t,b}=0.84$ is the assumed bulk turbulent Prandtl number (namely, away from walls), and $C=0.3$. Figure 15 shows that ${{Pr}}_t$ is far from universal across the wall-normal stations, and large scatter of the bundles is observed near the wall, whereas slightly better behaviour is observed at $y/h \gtrsim 0.5$. We further note that in the outer region, the turbulent Prandtl number attains values, around 0.7, that are significantly lower than in plane channel flow (Pirozzoli et al. 2016). The engineering correlations do not perform well. The model of Cebeci (1973) largely overpredicts ${{Pr}}_t$ in the duct core, and it does not reproduce the correct trend close to the wall. Also, the model of Kays & Crawford (1993) overpredicts the turbulent Prandtl number in the bulk flow, and it does not account accurately for spatial non-uniformities.

Figure 15. Wall-normal distributions of turbulent Prandtl number ${{Pr}}=\nu _t/\alpha _t$, in outer coordinates (a), and inner coordinates (b). Profiles are plotted at all $z$ locations, up to the corner bisector ($y=z$, see figure 2). Flow cases E (${{Re}}_\tau =1000$, orange) and F (${{Re}}_\tau =2000$, black) are compared to (5.8a,b) of Cebeci (1973, red dashed) and (5.9) of Kays & Crawford (1993, blue).

We assess further the validity of the nonlinear eddy diffusivity model (5.3) for different truncations of the vector polynomial basis. Here we focus only on the $q_1$ and $q_2$ components of the heat flux vector because $q_2(y,z)=q_3(z,y)$ after geometrical symmetry. Several nonlinear constitutive relations can be obtained with different combinations of the vector integrity bases (5.4), providing different accuracy for the modelled heat flux vector. In order to quantify the error with respect to DNS, we use the averaged correlation coefficient between $q_{i}$ and $\tilde {q}_{i}$,

(5.10)\begin{equation} \tilde{C}_{i} = \frac{\langle q_{i}\tilde{q}_{i}\rangle-\langle q_{i}\rangle \langle \tilde{q}_{i}\rangle }{\left[\left(\langle q_{i}^2\rangle- \langle q_{i}\rangle^2\right)\left(\langle \tilde{q}_{i}^2\rangle- \langle \tilde{q}_{i}\rangle^2\right)\right]^{1/2}}, \end{equation}

where the angle brackets denote average over the duct cross-section. By construction, $\tilde {C}_i\in [-1,1]$, where $\tilde {C}_i=1$ indicates perfect correlation, and $\tilde {C}_i=-1$ negative correlation.

Table 3 shows the correlation coefficients $\tilde {C}_{1}$ and $\tilde {C}_{2}$ for flow case F, where we consider all possible base combinations that can be obtained using the first five bases, and up to $N=3$. The analysis shows that it is possible to find optimal base subsets providing the maximum accuracy for a given $N$ (boldface values in table 3). The linear eddy diffusivity hypothesis ($N=1$) is able to predict accurately the turbulent heat flux component $\tilde {q}_2$, whereas $\tilde {q}_1=0$ by construction. In the case of $N=2$, we find that the base combination $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)}$ is the one that brings the highest accuracy on $\tilde {q}_1$ and also a very good prediction of $\tilde {q}_2$, although slightly lower than the linear eddy diffusivity hypothesis. The choice $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(4)}$ leads to the same accuracy on $\tilde {q}_2$ as for the linear case, but it increases the error on $\tilde {q}_1$. We further note that $\boldsymbol {V}^{(4)}$ contains the anisotropic Reynolds stress tensor, therefore it is not a practical choice from the modelling point of view. Finally, we find that using three vector bases is sufficient to recover the exact representation of the turbulent heat flux vector if the subset $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)},\boldsymbol {V}^{(5)}$ is used. Notably, these vector bases do not contain the anisotropic Reynolds stress tensor, whose presence would certainly lead to a much lower accuracy in a practical scenario, as it also requires modelling.

Table 3. Correlation coefficient $\tilde {C}_{i}$ (see (5.10)) between the turbulent heat flux of DNS $q_{i}$ and modelled turbulent heat flux $\tilde {q}_{i}$ (see (5.3)) for flow case F and for different combinations of vector bases. The optimal subset for each $N$ is highlighted in boldface.

Figure 16 shows scatter plots of the modelled turbulent heat flux components versus those obtained from DNS, in which the axis bisector indicates ideal model behaviour. Here, we limit ourselves to plotting the results for the linear eddy diffusivity hypothesis and for the optimal base subsets highlighted in table 3. The linear eddy diffusivity hypothesis returns an accurate prediction of $\tilde {q}_2$, but it is not able to predict $\tilde {q}_1$, which reminds us of difficulties with linear eddy-viscosity models in predicting the normal Reynolds stress components. Using two vector bases improves the prediction of $\tilde {q}_1$, but it yields a less accurate approximation of $\tilde {q}_2$, which would negatively impact prediction of the turbulent heat flux. Similar considerations apply to the anisotropic Reynolds stress tensor when two tensor bases are used (Modesti 2020), although less evident than here. Finally, exact representation of the heat flux vector is recovered for $N=3$.

Figure 16. Scatter plots of the modelled turbulent heat flux vector components ($\tilde {q}_1$, $\tilde {q}_2$) versus those obtained from DNS ($q_1$, $q_2$) for flow case F, for an increasing number of vector bases: $N=1$ (a,b), $N=2$, $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)}$ (c,d), and $N=3$, $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)},\boldsymbol {V}^{(5)}$ (e, f). The dashed diagonal line indicates ideal model performance.

Additional information on the accuracy of nonlinear eddy diffusivity models can be extracted from the spatial distribution of $\tilde {q}_1$ and $\tilde {q}_2$ for an increasing number of tensor bases; see figure 17. When linear representation is used, $\tilde {q}_2$ closely resembles the DNS data, whereas $\tilde {q}_1=0$. Using two vector bases improves the prediction of $\tilde {q}_1$, which becomes indistinguishable from DNS data, whereas the modelling accuracy of $\tilde {q}_2$ substantially decreases, especially close to the side wall. For $N=3$, the exact heat flux vector is recovered. This analysis reveals that the linear eddy diffusivity hypothesis has the potential to predict accurately $\langle {v'\theta '}\rangle$, and therefore the overall heat flux. However, this good performance is likely to be hampered by the use of a constant turbulent Prandtl number, which is far from being accurate. Slightly better prediction of $\tilde {q}_2$ can be achieved using a nonlinear eddy diffusivity with at least three vector bases, but improvement in the prediction $\tilde {q}_2$ might not justify the additional modelling complexity.

Figure 17. Modelled turbulent heat flux vector components ($\tilde {q}_1$, $\tilde {q}_2$) in a duct quadrant for flow case F and an increasing number of vector bases: $N=1$ (a,b), $N=2$, $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)}$ (c,d), and $N=3$, $\boldsymbol {V}^{(1)},\boldsymbol {V}^{(3)},\boldsymbol {V}^{(5)}$ (e, f).

6. Conclusions

We have carried out DNS of square duct flow at the unprecedented Reynolds number ${{Re}}_\tau \approx 2000$, wherein Navier–Stokes equations have been augmented with the transport of a passive scalar representing the temperature field. This configuration corresponds to the case of forced convection in which hot fluid is pumped through the duct and cooled at the walls. The choice of using a unity Prandtl number and homogeneous cold walls allows us a direct comparison between heat and momentum transport. In this respect, we find that the streamwise velocity and temperature fields are very similar, both instantaneously and on average. Close to the walls, $U$ and $\varTheta$ are nearly identical because the fluid viscosity and thermal diffusivity dominate, and ${{Pr}}=1$, whereas in the overlap layer they both follow a logarithmic law, although with a different slope. At low Reynolds numbers, the correlation coefficient between $u$ and $\theta$ is less than 1 at the duct core, but it but it approaches the unit value as the flow becomes well mixed, restoring the similarity inherent in the underlying equations.

The Nusselt number of square duct flow is in excellent agreement with pipe flow data at matching Reynolds number, which supports the validity of the hydraulic diameter concept. We explain this good match by using a simple model in which the duct flow can be regarded as the superposition of four concurrent walls, where the flow is controlled by the closest one. This simple cartoon is well supported by the mean temperature profiles in a duct octant, which follow with good accuracy the canonical law-of-the-wall.

We have derived a generalized version of the popular FIK identity for the heat flux, which confirms the common idea that secondary flows increase the heat transfer by mixing hot fluid at the duct core with cold fluid at the corners. We have also carried out numerical simulations whereby secondary flows have been suppressed artificially, which also confirms the minor effect of the secondary currents on the global heat transfer. Although the secondary flows are able to bend the isolines of the mean velocity and temperature, their global effect on heat transfer is rather small and can be estimated to be about $5\,\%$.

The DNS data have been used to establish the accuracy of linear and nonlinear eddy diffusivity closures for the turbulent heat flux. The commonly accepted hypothesis of uniform turbulent Prandtl number is highly inaccurate in a large part of the duct, and at the duct core ${{Pr}}_t\approx 0.7$, which is much lower than what is found in channel and pipe flow (${{Pr}}_t\approx 0.85$). We show that an exact representation of the heat flux vector can be recovered using a vector polynomial integrity basis expansion with at least three bases, although this seems impractical from a modelling point of view. The linear eddy diffusivity formulation is able to predict accurately the wall-normal turbulent heat flux, therefore it can be a much more reliable alternative to the turbulent Prandtl number, provided that an accurate transport equation for the eddy conductivity is available.