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

Davide Modesti1<EMAIL_ADDRESS>Sergio Pirozzoli2 1 Faculty of Aerospace
Engineering, Delft University of Technology, Kluyverweg 2, 2629 HS Delft, The
Netherlands 2 Dipartimento di Ingegneria Meccanica e Aerospaziale, Sapienza
Università di Roma, via Eudossiana 18, 00184 Roma, Italia

###### Abstract

We carry out direct numerical simulation (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 $\Rey_{\tau}\approx 2000$ are
considered which are much higher than in previous studies, and which yield
clear scale separation between inner- and outer-layer dynamics. Close
similarity between the behavior 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, however 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 rigorously show here. Temperature and velocity
fluctuations also have similar behavior, 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 modeling improvement may be in
principle obtained by retaining at least the three terms in the vector
polynomial integrity basis expansion.

###### keywords:

##  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
extensively reported in classical textbooks (Kays and Crawford, 1993; Rohsenow
et al., 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 and Crawford, 1993; White and
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 and Sekulib,
1998).

Heat transfer in square ducts was first studied experimentally by Brundrett
and Burroughs (1967), who considered air flow at bulk Reynolds number
$\Rey_{b}=Hu_{b}/\nu$ (where $H=2h$ is duct side length and also duct
hydraulic diameter, $u_{b}$ is the bulk velocity, and $\nu$ the fluid
kinematic viscosity) 33000 and 67000. Close similarity between the wall shear
stress and the heat flux distributions was shown (with quoted discrepancy of
$\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 of
$\kappa_{\theta}\approx 0.51$. Those findings were qualitatively supported
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 and
Yaglom, 1972).

Early computational studies of heat transfer in square ducts mostly relied on
the use of RANS models. For instance, Launder and 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 qualitatively
reproduced the experimental data of Brundrett and Burroughs (1967). High-
fidelity computational studies of heat transfer in square ducts have been
quite limited so far. Vázquez and Métais (2002) first studied turbulent flow
through a heated square duct by means of large-eddy simulation (LES) at
moderate Reynolds number ($\Rey_{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 and Davidson (2002). Yang et al.
(2009) carried out under-resolved DNS (direct numerical simulation) 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 et
al. (2019), at bulk Reynolds number up to $\Rey_{b}=10^{4}$. As noted by
Brundrett and 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 is 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 number than in previous computational studies. Although
relatively simple, this setup 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,
$\Pran=\nu/\alpha$), is instrumental to more closely scrutinizing differences
between the behavior of the streamwise velocity field and of passively
advected scalars. The present study is the continuation of previous efforts
(Pirozzoli et al., 2018; Modesti et al., 2018; Orlandi 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 in the order of a few
percent of the duct bulk velocity, and despite their persistence they do
contribute to the mean duct friction by at most a few percent. 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

Case | $\Rey_{b}$ | $\Rey_{\tau}^{*}$ | $C_{f}\times 10^{3}$ | Nu | $N_{x}$ | $N_{y}$ | $N_{z}$ | $\Delta x^{*}$ | $\Delta z^{*}$ | $\Delta y_{w}^{*}$ | Symbol
---|---|---|---|---|---|---|---|---|---|---|---
A | $4410$ | $150$ | $9.26$ | $18.1$ | $512$ | $128$ | $128$ | $5.6$ | $3.0$ | $0.55$ | $\bigcirc$
B | $7000$ | $227$ | $8.41$ | $26.9$ | $640$ | $144$ | $144$ | $6.6$ | $4.8$ | $0.51$ | $\bigtriangledown$
C | $12000$ | $365$ | $7.40$ | $41.9$ | $768$ | $208$ | $208$ | $9.0$ | $5.3$ | $0.69$ | $\bigtriangleup$
D | $17800$ | $519$ | $6.80$ | $57.6$ | $1024$ | $256$ | $256$ | $9.5$ | $6.3$ | $0.53$ | $\triangleright$
E | $40000$ | $1055$ | $5.57$ | $109.0$ | $2048$ | $512$ | $512$ | $9.6$ | $6.4$ | $0.60$ | $\square$
F | $84000$ | $2041$ | $4.72$ | $197.0$ | $4096$ | $1024$ | $1024$ | $9.3$ | $6.0$ | $0.67$ | $\Diamond$
A0 | $4410$ | $154$ | $9.66$ | $21.1$ | $512$ | $128$ | $128$ | $5.7$ | $3.3$ | $0.54$ | -
B0 | $7000$ | $231$ | $8.58$ | $30.0$ | $640$ | $144$ | $144$ | $6.8$ | $4.7$ | $0.67$ | -
C0 | $12000$ | $367$ | $7.48$ | $44.5$ | $768$ | $204$ | $204$ | $9.0$ | $5.4$ | $0.71$ | -
D0 | $17800$ | $515$ | $6.62$ | $59.6$ | $1024$ | $256$ | $256$ | $9.5$ | $6.0$ | $0.79$ | -
Table 1: Flow parameters for square duct DNS. Box dimensions are $6\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. $\Rey_{b}=2hu_{b}/\nu$ is
the bulk Reynolds number, and $\Rey_{\tau}^{*}=hu_{\tau}^{*}/\nu$ is 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 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 and 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 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 and Bernardini, 2013; Modesti and
Pirozzoli, 2016; Bernardini et al., 2021). Here, the convective terms in the
Navier–Stokes equations are preliminarily expanded to quasi-skew-symmetric
form, in such a way to discretely preserve 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 and
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 ${\Pi}$), 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\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 $\Rey_{b}=4400-84000$ (table 1), and hereafter labeled 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 numerically suppressed, and which are denoted with the
suffix 0. The turbulence Mach number $M_{t}=u^{\prime}/{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 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 et al., 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 and 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 slightly higher heat flux number as compared to the case of strictly
constant heat flux (Abe and Antonia, 2017; Alcántara-Ávila et al., 2021). In
the following we use a different notation for the mean temperature in the
duct,

$\theta_{b}=\frac{1}{A_{c}}\int_{A_{c}}(\Theta-\theta_{w})\mbox{{d}}A_{c},$
(1)

and for the bulk temperature,

$\theta_{m}=\frac{1}{u_{b}A_{c}}\int_{A_{c}}U(\Theta-\theta_{w})\mbox{{d}}A_{c},\quad
u_{b}=\frac{1}{A}_{c}\int_{A_{c}}U\mbox{{d}}A_{c},$ (2)

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).

The $+$ superscript is here used to denote local wall units, namely quantities
made nondimensional 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}|_{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,

$u^{*}_{\tau}=\left(\frac{\tau^{*}_{w}}{\rho}\right)^{1/2}=\left(\frac{h\left<\Pi\right>}{2\rho}\right)^{1/2},$
(3)

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

$\theta_{\tau}=\frac{q_{w}}{u_{\tau}},$ (4)

where $q_{w}=\alpha{\mbox{{d}}{\Theta}}/{\mbox{{d}}y}|_{w}$ is the local wall
heat flux, or the global friction temperature,

$\theta^{*}_{\tau}=\frac{q^{*}_{w}}{u^{*}_{\tau}}=\frac{\left<{Q}\right>h}{2u^{*}_{\tau}}.$
(5)

Equations (3) and (5) make it clear that the global friction velocity and
temperature are related with the imposed forcing of the streamwise momentum
and temperature equations, rather than with 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 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

(a)  (b)

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.

We begin by inspecting the instantaneous streamwise velocity and temperature
fields in figure 1, which shows the flow both in 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 artifact, 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 absence of the smoothing action
of pressure (Pirozzoli et al., 2016).

Figure 2: Mean streamwise temperature (${\Theta}^{*}$, 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 are shown. Only a quarter of
the full domain is shown. The dashed diagonal lines indicate the corner
bisector. Figure 3: Scatter plot of mean temperature ($\Theta^{*}$) versus
mean velocity ($U^{*}$) for DNS data of all flow cases, and corresponding
fitting functions for the near-wall region, $\Theta^{*}=U^{*}$ (dashed line),
and for the outer region $\Theta^{*}=2.57+0.84\,U^{*}$ (solid line). Colors as
in table 1.

Figure 2 depicts the mean temperature (contours) and velocity (iso-lines)
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^{*}\leq 10$), perfect matching is recovered (namely, $\Theta^{*}=U^{*}$),
which is not surprising, as at unit Prandtl number $\Theta^{+}\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 $\Theta^{*}=\Pran\,U^{*}$, since in
general $\Theta^{+}\equiv U^{+}\equiv\Pran\,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, however 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 $\Pran\neq 1$ we thus expect a
similar relationship, however 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 as a square duct would allow in principle to directly reconstruct the
whole temperature field based on the velocity field.

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
fit of experimental data (Kader, 1981), at matching $\Rey_{\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
fit of experimental data (Kader, 1981), at matching $\Rey_{\tau}^{*}$.

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 the experimental fitting by Kader (1981) at matching $\Rey_{\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 number. 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 nonuniform 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.2h$,
which is also the lower limit for the core region in canonical flows (Pope,
2000). This observation seems to support validity of Townsend’s outer layer
similarity also for the temperature field, in that the outer flow mainly
perceives 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 et al., 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
previously reported 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 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), (b), (c); B (d), (e), (f); C (g), (h), (j);
D (k), (l), (m); E (n), (o), (p); F (q), (r), (a). Note that the horizontal
coordinate ($z$) is given in inner units, and the vertical coordinate $y$ is
given in outer units and in a logarithmic scale.

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 instantaneous sense, with the exception of
the close neighborhood 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 7: Inner-layer peak variances of axial velocity (a), and temperature
(b), expressed as a function of the local friction Reynolds number ($z^{+}$).
The dashed lines denote the logarithmic fit for canonical flows (Marusic et
al., 2017), namely $\langle
u^{2}\rangle^{+}_{{}_{IP}}=0.63\log\Rey_{\tau}+3.8$. Symbols as in table 1.

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 all 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 of Townsend’s
attached-eddy hypothesis. We speculate that this peculiar behavior may be
either due to disruption of the hierarchy of wall-attached eddies due to the
action of secondary motions, or perhaps to impeded spanwise growth caused by
the presence of corners, both effects being absent in planar channel or pipe
flow.

Figure 8: Distributions of local wall friction coefficient
$C_{f}=2u_{\tau}^{2}/u_{b}^{2}$ (a), local Stanton number
$\mbox{{St}}=(u_{\tau}\theta_{\tau})/(u_{b}\theta_{m})$ (b), and Reynolds
analogy factor (c) along left part of lower wall. 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 previously shown by Pirozzoli et al.
(2018), the distributions of the wall friction coefficients is non-monotone,
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 et al., 2018), the friction tends to equalize over the duct perimeter
at higher $\Rey$. The heat transfer coefficient seems to have the same
qualitative behavior, thus corroborating the experimental observations of
Brundrett and Burroughs (1967); Hirota et al. (1997). Quantitative assessment
is shown in figure 8(c), where we show the Reynolds analogy factor,
$2\mbox{{St}}/C_{f}$, which is expected to be unity at $\Pran=1$ (White and
Majdalani, 2006). In fact, whereas deviations are observed at low Reynolds
number, the distribution corresponding to DNS-F varies between 0.98 and 1.2,
right at the corner.

Figure 9:  Inverse global Stanton number (a) and Nusselt number (b) as a
function 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 and Davidson, 2002), and from experiments (Brundrett and
Burroughs, 1967, triangles). The dashed line denotes the correlation (9), and
the dot-dashed lines the correlation (8).

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

$\mbox{{St}}^{*}=\frac{q_{w}^{*}}{u_{b}\theta_{m}}=\frac{1}{u_{b}^{*}\theta_{m}^{*}},$
(6)

or more frequently in terms of the Nusselt number,

$\mbox{{Nu}}^{*}=\mbox{{St}}^{*}\,\Rey_{b}\,\Pran.$ (7)

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 $\Pran=0.71$, hence the Nu
data are rescaled by a factor $1/\sqrt{\Pran}$ to compare with the present
ones. As a reference, we also report correlations widely used in the
engineering practice, including that by Gnielinski (1976),

$\mbox{{Nu}}=\frac{\Pran
C_{f}/32\left(\Rey_{b}-1000\right)}{1+12.7(C_{f}/32)^{1/2}\left(\Pran^{2/3}-1\right)},$
(8)

and from Kays and Crawford (1993),

$\mbox{{Nu}}=0.022\,\Rey_{b}^{0.8}\,\Pran^{0.5}.$ (9)

Last, direct fitting the present DNS data (at $\Pran=1$) with a power-law
expression yields

$\mbox{{Nu}}=0.0216\,\Rey_{b}^{0.805}.$ (10)

The DNS data show good correspondence with the pipe DNS data at matching
$\Rey_{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
modeling 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/\mbox{{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 partly due
to inaccuracy of correlations based on old experimental data, or to the fact
that those are mainly trained for the $\Pran=0.71$ case, whereas here
$\Pran=1$. Discrepancies can also be partly attributed 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 approached in which
the total heat flux is maintained strictly constant in time (Abe and Antonia,
2017; Alcántara-Ávila et al., 2021). Slight adjustment of the Kays–Crawford
power-law formula coefficients as after equation (10) seems to yield very good
representation of the DNS data.

An issue which 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 it is
also the case of the heat flux coefficient. As previously shown, 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$,
$\Theta^{*}=1/\kappa_{\theta}\log y^{*}+C_{\theta}$, and integrating over the
duct cross section, the following expression for the inverse Stanton number
velocity results

$\begin{split}1/\mbox{{St}}^{*}=u_{b}^{*}\theta_{m}^{*}=\frac{8}{4\Rey_{\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]dy^{*}dz^{*}=\\\
\frac{1}{\kappa\kappa_{\theta}}\left(\log^{2}{\Rey_{\tau}^{*}}+3\log{\Rey_{\tau}^{*}}+\frac{7}{2}\right)+\left(\frac{C_{\theta}}{\kappa}+\frac{C}{\kappa_{\theta}}\right)\left(\log{\Rey_{\tau}^{*}}-\frac{3}{2}\right)+CC_{\theta}.\end{split}$
(11)

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

$1/\mbox{{St}}=\frac{1}{\kappa\kappa_{\theta}}\left(\log^{2}{\Rey_{\tau}}+3\log{\Rey_{\tau}}+\frac{7}{2}\right)+\left(\frac{C_{\theta}}{\kappa}+\frac{C}{\kappa_{\theta}}\right)\left(\log{\Rey_{\tau}}-\frac{3}{2}\right)+CC_{\theta},$
(12)

where $\Rey_{\tau}=Du_{\tau}/(2\nu)$. The two expressions are identical for
$2h=D$, hence provided the Reynolds number based on the hydraulic diameter is
the same. It is interesting that equation (11) is basically arrived at 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 et al., 2002) for the Stanton number, by following the approach of
Modesti et al. (2018). We consider the mean temperature balance equation,

$\alpha\nabla^{2}\Theta=\nabla\cdot\boldsymbol{q}_{C}+\nabla\cdot\boldsymbol{q}_{T}-\langle{Q}\rangle,$
(13)

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

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

$\alpha\nabla^{2}{\Theta}_{D}=-\langle{Q}\rangle,\quad\alpha\nabla^{2}{\Theta}_{T}=\nabla\cdot\boldsymbol{q}_{T},\quad\alpha\nabla^{2}{\Theta}_{C}=\nabla\cdot\boldsymbol{q}_{C},$
(14)

with homogeneous boundary conditions, where ${\Theta}_{D}$, ${\Theta}_{T}$,
and ${\Theta}_{C}$ denote the diffusive, turbulent, and convective
contributions to the mean temperature field. The mean temperature in the duct
may accordingly be evaluated as

$\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}{\Theta}_{X}\mathrm{d}A,$
(15)

where we have introduced the unitary Stokes temperature field $\theta_{1}$,
defined as solution of $\nabla^{2}\theta_{1}=-1/A$, which by construction is
only a function of the duct geometry, and which allows us to express the
diffusion-associated temperature field as
$\Theta_{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

$\mbox{{St}}_{b}^{*}=\frac{q_{w}^{*}}{u_{b}\theta_{b}}=\frac{1}{u_{b}^{*}\theta_{b}^{*}},$
(16)

instead of (6), which when substituted into equation (15) yields

$\mbox{{St}}_{b}^{*}=\frac{1}{\Pran\Rey_{P}{\theta_{b}}_{1}}\left(1-\frac{{\theta_{b}}_{T}}{\theta_{b}}-\frac{{\theta_{b}}_{C}}{\theta_{b}}\right)=\mbox{{St}}_{D}+\mbox{{St}}_{T}+\mbox{{St}}_{C}\quad,$
(17)

where $\Rey_{P}=u_{b}P/\nu$ is the bulk Reynolds number based on the duct
perimeter $P$. Equation (17) clearly shows that the heat transfer coefficient
may be regarded as the sum of contributions due to diffusivity, turbulence,
and convection (labeled as ${St}_{X}$ in (17)), and it yields a similar
expression to classical FIK identity for the friction coefficient.

Case | $\mbox{{St}}^{*}_{V}\times 10^{3}$ | $\mbox{{St}}^{*}_{T}\times 10^{3}$ | $\mbox{{St}}^{*}_{C}\times 10^{3}$ | $\mbox{{St}}^{*}\times 10^{3}$ | $\mbox{{St}}^{*}_{V}/\mbox{{St}}^{*}(\%)$ | $\mbox{{St}}^{*}_{T}/\mbox{{St}}^{*}(\%)$ | $\mbox{{St}}^{*}_{C}/\mbox{{St}}^{*}(\%)$
---|---|---|---|---|---|---|---
A | 0.095 | 0.144 | 0.030 | 0.268 | 35.3 | 53.7 | 11.0
B | 0.060 | 0.176 | 0.009 | 0.245 | 24.4 | 72.0 | 3.6
C | 0.035 | 0.172 | 0.009 | 0.216 | 16.2 | 79.8 | 4.0
D | 0.024 | 0.170 | 0.007 | 0.201 | 11.7 | 84.7 | 3.6
E | 0.010 | 0.149 | 0.008 | 0.167 | 6.2 | 89.2 | 4.5
F | 0.005 | 0.123 | 0.010 | 0.138 | 3.6 | 89.4 | 7.0
Table 2:  Contributions of diffusion, turbulence and mean convection terms to
the heat transfer coefficient, as from equation (17).  Figure 10: Contribution
of diffusion, turbulence and mean convection terms to the mean temperature
field, as determined by solving equation (14), for flow cases A (a), (b), (c);
B (d), (e), (f); C (g), (h), (j); D (k), (l), (m); E (n), (o), (p); F (q),
(r), (a).

We note that this formalism has additional generality than the classical FIK
identity and its extensions (Peet and Sagaut, 2009; Jelly et al., 2014), and
it allows to isolate the effects of mean convection, turbulent and diffusion
terms to the wall shear stress distributions along the duct perimeter,
resulting from

${q_{w}}_{D}=\alpha\left.\frac{\partial{\Theta}_{D}}{\partial
y}\right\rvert_{w},\quad{q_{w}}_{T}=\alpha\left.\frac{\partial{\Theta}_{T}}{\partial
y}\right\rvert_{w},\quad{q_{w}}_{C}=\alpha\left.\frac{\partial{\Theta}_{C}}{\partial
y}\right\rvert_{w}.$ (18)

Regarding equation (18) it is important to note that only ${q_{w}}_{D}$ has
non-zero mean, whereas integration of equation (14) for ${\Theta}_{C}$ and
${\Theta}_{T}$ readily shows that their integrated contributions vanish as
$\boldsymbol{q}_{C}$, $\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 convective
contribution to exceed the diffusion one at high enough Reynolds number.
Figure 10 shows the mean temperature fields associated with the diffusion,
turbulent and convective 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 equation (17)). The diffusion-
associated temperature field $\Theta_{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 and London, 2014), only
depending on the duct cross-sectional geometry. The turbulence-associated
temperature field $\Theta_{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 $\Theta_{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 half-way 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 change of the Reynolds number only changes the relative
importance of the three terms, as quantified in table 2.

Figure 11: Contributions to mean wall heat flux (shown along half 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).

The observations made regarding the organization of the temperature field have
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
behavior, exhibiting multiple peaks which 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 nonuniform. Mean convection
yields large positive peaks in the corner vicinity, whose distance from the
neighboring 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 number.

Figure 12: Flow cases without secondary motions: mean streamwise velocity
(${U}^{*}$, lines), and mean temperature (${\Theta}^{*}$, 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.

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 artificially suppressed (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

$\displaystyle v(x,y,z,t)$ $\displaystyle\rightarrow
v(x,y,z,t)-\overline{v}^{x}(y,z,t)$ (19a) $\displaystyle w(x,y,z,t)$
$\displaystyle\rightarrow w(x,y,z,t)-\overline{w}^{x}(y,z,t),$ (19b)

at each Runge-Kutta sub-step, where $\overline{(.)}^{x}$ denotes the
streamwise averaging operator. Figure 12 shows the temperature contours with
superposed 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 number, whereas differences become more evident as $\Rey_{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 reported by Modesti et al. (2018) for
frictional drag, suppressing the secondary flows yields even higher hear
transfer at low Reynolds number, whereas 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.

##  5 Modeling the turbulent heat flux

The numerical solution of the Reynolds Averaged Navier–Stokes (RANS) equations
is a common approach for flows of industrial interest. Solution of the
temperature field requires a suitable closure models 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 relying 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,

$\Pran_{t}=\frac{\nu_{t}}{\alpha_{t}},$ (20)

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,

$\widetilde{\boldsymbol{a}}=-\nu_{t}\left(\nabla\boldsymbol{U}+\nabla\boldsymbol{U}^{T}\right),\quad\widetilde{\boldsymbol{q}}=-\alpha_{t}\nabla\Theta,$
(21)

where $\widetilde{\boldsymbol{a}}$ is the modelled counterpart of the
anisotropic Reynolds stress tensor ($\boldsymbol{a}=\boldsymbol{\tau}-2/3k$),
and $\widetilde{\boldsymbol{q}}$ the modelled turbulent heat flux vector, with
$\boldsymbol{\tau}=\langle{\boldsymbol{u}\boldsymbol{u}}\rangle$,
$k=1/2\,\mathrm{tr}(\boldsymbol{\tau})$. Most heat transfer models assume
uniform turbulent Prandtl number ($\Pran_{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 $\Pran_{t}\approx 1$ holds, limitations of the linear
eddy-viscosity hypothesis (21) are evident, and reliance on the turbulent
Prandtl number concept replicates modeling inaccuracies. Pope (1975) extended
the eddy-viscosity ansatz by using a tensor polynomial to model
$\widetilde{\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 exact representation of
$\widetilde{\boldsymbol{a}}$ (Gatski and Speziale, 1993; Jongen and Gatski,
1998; Modesti, 2020). The same mathematical framework was extended to scalar
flux modelling, both for heat transfer (So et al., 2004a; Younis et al., 2005;
Wang et al., 2007) and for buoyancy (So et al., 2004b). Similar to the
original tensor polynomial expansion of Pope (1975), the turbulent heat flux
is expanded as

$\widetilde{\mathbf{q}}=\sum_{n=1}^{N}\alpha^{(n)}\mathbf{V}^{(n)},\quad
N=1,\dots,10$ (22)

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

$\displaystyle\begin{aligned} \mathbf{V}^{(1)}&=\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(2)}&=\mathbf{a}\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(3)}&=\mathbf{S}\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(4)}&=\mathbf{a}\mathbf{a}\boldsymbol{\nabla}{\Theta}\\\
\mathbf{V}^{(5)}&=\mathbf{S}\mathbf{S}\boldsymbol{\nabla}{\Theta}\\\
\end{aligned}$ $\displaystyle\begin{aligned}
\mathbf{V}^{(6)}&=\mathbf{\Omega}\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(7)}&=\mathbf{\Omega}\mathbf{\Omega}\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(8)}&=\left(\mathbf{S}\mathbf{\Omega}+\mathbf{\Omega}\mathbf{S}\right)\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(9)}&=\left(\mathbf{a}\mathbf{S}+\mathbf{S}\mathbf{a}\right)\boldsymbol{\nabla}\Theta\\\
\mathbf{V}^{(10)}&=\left(\mathbf{\Omega}\mathbf{S}+\mathbf{S}\mathbf{\Omega}\right)\boldsymbol{\nabla}\Theta,\end{aligned}$
(23)

where
$\boldsymbol{S}=1/2\left(\nabla\boldsymbol{U}+\nabla\boldsymbol{U}^{T}\right)$,
$\boldsymbol{\Omega}=1/2\left(\nabla\boldsymbol{U}-\nabla\boldsymbol{U}^{T}\right)$,
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 exact representation of the turbulent heat flux
($\widetilde{\boldsymbol{q}}=\boldsymbol{q}$) is recovered for $N=3$. In order
to assess the accuracy of different truncations of the polynomial expansion
(22) we compute the coefficients $\alpha^{(n)}$ following the approach of
Jongen and Gatski (1998); Modesti (2020), namely we take the scalar product of
(22) with each vector basis $\mathbf{V}^{(m)}$, thus obtaining a linear system
of $N$ equations for the unknown coefficients $\alpha^{(n)}$,

$\mathbf{q}\cdot\mathbf{V}^{(m)}=\sum_{n=1}^{N}\alpha^{(n)}\mathbf{V}^{(n)}\cdot\mathbf{V}^{(m)},\quad
m=1,\dots,N.$ (24)

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

$\alpha_{t}=-\alpha^{(1)}=-\frac{\mathbf{q}\cdot\boldsymbol{\nabla}\Theta}{\boldsymbol{\nabla}\Theta\cdot\boldsymbol{\nabla}\Theta},$
(25)

in analogy with the definition of linear eddy viscosity,

$\nu_{t}=-\frac{\left\\{\mathbf{a}\mathbf{S}\right\\}}{2\left\\{\mathbf{S}\mathbf{S}\right\\}},$
(26)

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

Figure 13: Estimated fields of linear eddy viscosity $\nu_{t}$ as after
equation (26) (a), and linear eddy diffusivity $\alpha_{t}$ as after equation
(25) (b), for flow case F ($\Rey_{\tau}=2000$). Only one quadrant of the duct
is shown.

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 $\Pran_{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 iso-lines more markedly protrude towards the corner and tend to be more
parallel to the walls as compared to $\nu_{t}$.

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 prediction of the mixing length
hypothesis $\kappa y^{+}$ (a) and $\kappa_{\theta}y^{+}$ (b).

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. 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 has also the effect of enhancing interactions
between neighboring 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 reveals to be only
partially accurate in the overlap layer, especially for the eddy diffusivity.

Figure 15: Wall-normal distributions of turbulent Prandtl number
$\Pran=\nu_{t}/\alpha_{t}$, in outer coordinates (a) and in inner coordinates
(b). Profiles are plotted at all $z$ locations, up to the corner bisector
($y=z$, see figure 2). Flow cases E ($\Rey_{\tau}=1000$, orange) and F
($\Rey_{\tau}=2000$, black) are compared to equation (27) of Cebeci (1973, red
dashed) and equation (28) of Kays and Crawford (1993, blue).

Figure 15 shows bundles of the turbulent Prandtl number (20) 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),

$\Pran_{t}=\frac{k}{\kappa_{\theta}}\left(\frac{1-\exp{(-y^{+}/A)}}{1-\exp{(-y^{+}/B)}}\right),\quad
B=\frac{1}{\Pran^{(1/2)}}\sum_{i=1}^{5}C_{i}\left(\log_{10}\Pran\right)^{(i-1)},$
(27)

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 and Crawford (1993),

$\Pran_{t}=\left[\frac{1}{2\Pran_{t,b}}+\frac{C\Pen_{t}}{\sqrt{\Pran_{t,b}}}-(C\Pen_{t})^{2}\left(1-\exp{\left(-\frac{1}{C\Pen_{t}\sqrt{\Pran_{t,b}}}\right)}\right)\right]^{-1},$
(28)

where $\Pen_{t}=\Pran\,\nu_{t}/\nu$ is the turbulent Peclet number,
$\Pran_{t,b}=0.84$ is the assumed bulk turbulent Prandtl number (namely, away
from walls), and $C=0.3$. Figure 15 shows that $\Pran_{t}$ is far from
universal across the wall-normal stations, and large scatter of the bundles is
observed near the wall, whereas slightly better behavior 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, which 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
$\Pran_{t}$ in the duct core, and it does not reproduce the correct trend
close to the wall. Also the model of Kays and Crawford (1993) overpredicts the
turbulent Prandtl number in the bulk flow, and it does not account accurately
for spatial non-uniformities.

We further assess the validity of the nonlinear eddy diffusivity model (22)
for different truncations of the vector polynomial basis. Here we only focus
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 (23), 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 $\widetilde{q}_{i}$,

$\widetilde{C}_{i}=\frac{\langle q_{i}\widetilde{q}_{i}\rangle-\langle
q_{i}\rangle\langle\widetilde{q}_{i}\rangle}{\left[\left(\langle
q_{i}^{2}\rangle-\langle
q_{i}\rangle^{2}\right)\left(\langle\widetilde{q}_{i}^{2}\rangle-\langle\widetilde{q}_{i}\rangle^{2}\right)\right]^{1/2}},$
(29)

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

N | $V^{(n)}$ | $\widetilde{C}_{1}$ | $\widetilde{C}_{2}$
---|---|---|---
1 | 1 | - | 0.997
2 | 1,2 | 0.538 | 0.973
1,3 | 1.000 | 0.988
1,4 | 0.871 | 0.997
1,5 | 0.896 | 0.801
3 | 1,2,3 | 0.999 | 1.000
1,2,4 | 0.538 | 0.973
1,2,5 | 0.988 | 1.000
1,3,4 | 1.000 | 0.988
1,3,5 | 1.000 | 1.000
1,4,5 | 0.892 | 0.801
Table 3:  Correlation coefficient $\widetilde{C}_{i}$ (29) between the
turbulent heat flux of DNS $q_{i}$ and modelled turbulent heat flux
$\widetilde{q}_{i}$, equation (22), for flow case F and for different
combinations of vector bases. The optimal subset for each $N$ is highlighted
in boldface.

Table 3 shows the correlation coefficients $\widetilde{C}_{1}$ and
$\widetilde{C}_{2}$ for flow case F, where we consider all possible bases
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 bases subsets providing
the maximum accuracy for a given $N$ (boldface values in table 3). The linear
eddy diffusivity hypothesis ($N=1$) is able to accurately predict the
turbulent heat flux component $\widetilde{q}_{2}$, whereas
$\widetilde{q}_{1}=0$ by construction. In the case of $N=2$ we find that the
bases combination $\mathbf{V}^{(1)},\mathbf{V}^{(3)}$ is the one which brings
the highest accuracy on $\widetilde{q}_{1}$ and also a very good prediction of
$\widetilde{q}_{2}$, although slightly lower than the linear eddy diffusivity
hypothesis. The choice $\mathbf{V}^{(1)},\mathbf{V}^{(4)}$ leads to the same
accuracy on $\widetilde{q}_{2}$ as for the linear case, but it increases the
error on $\widetilde{q}_{1}$. We further note that $\mathbf{V}^{(4)}$ contains
the anisotropic Reynolds stress tensor, and 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
$\mathbf{V}^{(1)},\mathbf{V}^{(3)},\mathbf{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.

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

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 behavior. Here we limit ourself to plotting the results for the linear
eddy diffusivity hypothesis and for the optimal bases subsets highlighted in
table 3. The linear eddy diffusivity hypothesis returns accurate prediction of
$\widetilde{q}_{2}$, but it is not able to predict $\widetilde{q}_{1}$, which
reminds us of difficulties of linear eddy-viscosity models in predicting the
normal Reynolds stress components. Using two vector bases improves the
prediction of $\widetilde{q}_{1}$, but it yields less accurate approximation
of $\widetilde{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 17: Modelled turbulent heat flux vector components
($\widetilde{q}_{1}$, $\widetilde{q}_{2}$) in a duct quadrant for flow case F
and increasing number of vector bases: $N=1$ (a, b), $N=2$
$\mathbf{V}^{(1)},\mathbf{V}^{(3)}$ (c, d), $N=3$
$\mathbf{V}^{(1)},\mathbf{V}^{(3)},\mathbf{V}^{(5)}$ (e, f).

Additional information on the accuracy of nonlinear eddy diffusivity models
can be extracted from the spatial distribution of $\widetilde{q}_{1}$ and
$\widetilde{q}_{2}$ for increasing number of tensor bases, figure 17. When
linear representation is used $\widetilde{q}_{2}$ closely resembles the DNS
data, whereas $\widetilde{q}_{1}=0$. Using two vector bases improves the
prediction of $\widetilde{q}_{1}$, which becomes indistinguishable from DNS
data, whereas the modeling accuracy of $\widetilde{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 accurately predict
$\langle{v^{\prime}\theta^{\prime}}\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 $\widetilde{q}_{2}$ can be achieved using a nonlinear
eddy diffusivity with at least three vector bases, but improvement in the
prediction $\widetilde{q}_{2}$ might not justify the additional modelling
complexity.

##  6 Conclusions

We have carried out DNS of square duct flow at the unprecedented Reynolds
number of $\Rey_{\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 $\Theta$ are nearly
identical because the fluid viscosity and thermal diffusivity dominate, and
$\Pran=1$, whereas in the overlap layer they both follow a logarithmic law,
although with a different slope. At low Reynolds number the correlation
coefficient between $u$ and $\theta$ is less than one 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 support 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 artificially suppressed, which also confirm 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 large
part of the duct, and at the duct core $\Pran_{t}\approx 0.7$, which is much
lower than what found in channel and pipe flow ($\Pran_{t}\approx 0.85$). We
show that 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 accurately predict the wall-normal
turbulent heat flux, and 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.

DNS data are available at http://doi.org/10.4121/19221657 and at
http://newton.dima.uniroma1.it/database/

Declaration of Interests. The authors report no conflict of interest.

Acknowledgements
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, and the DECI Research Infrastructure Beskow at
PDC, Stockholm, Sweden.

## References

  * Abe and Antonia [2017] H. Abe and R.A. Antonia.  Relationship between the heat transfer law and the scalar dissipation function in a turbulent channel flow.  _J. Fluid Mech._ , 830:300–325, 2017.
  * Alcántara-Ávila et al. [2021] F. Alcántara-Ávila, S. Hoyas, and M.J. Pérez-Quiles.  Direct numerical simulation of thermal channel flow for Re$\tau$= 5000 and Pr=0.71.  _J. Fluid Mech._ , 916, 2021.
  * Bernardini et al. [2021] M. Bernardini, D. Modesti, F. Salvadore, and S. Pirozzoli.  STREAmS: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows.  _Comput. Phys. Commun._ , 263:107906, 2021.
  * Brundrett and Burroughs [1967] E. Brundrett and P.R. Burroughs.  The temperature inner-law and heat transfer for turbulent air flow in a vertical square duct.  _Int. J. Heat Mass Transfer_ , 10:1133–1142, 1967.
  * Cebeci [1973] T. Cebeci.  A model for eddy conductivity and turbulent Prandtl number.  _Int. J. Heat Mass Transf._ , 95:227–234, 1973.
  * Fukagata et al. [2002] K. Fukagata, K. Iwamoto, and N. Kasagi.  Contribution of Reynolds stress distribution to the skin friction in wall-bounded flows.  _Phys. Fluids_ , 14:L73–L76, 2002.
  * Gatski and Speziale [1993] T.B. Gatski and C.G. Speziale.  On explicit algebraic stress models for complex turbulent flows.  _J. Fluid Mech._ , 254:59–78, 1993.
  * Gavrilakis [1992] S. Gavrilakis.  Numerical simulation of low-Reynolds-number turbulent flow through a straight square duct.  _J. Fluid Mech._ , 244:101–129, 1992.
  * Gnielinski [1976] V. Gnielinski.  New equations for heat and mass transfer in turbulent pipe and channel flow.  _Int. Chem. Eng._ , 16:359–367, 1976.
  * Hanjalić [2002] K. Hanjalić.  One-point closure models for buoyancy-driven turbulent flows.  _Annu. Rev. Fluid Mech._ , 34(1):321–347, 2002\.
  * Hirota et al. [1997] M. Hirota, H. Fujita, H. Yokosawa, H. Nakai, and H. Itoh.  Turbulent heat transfer in a square duct.  _Int. J. Heat Fluid Flow_ , 18:170–180, 1997.
  * Jelly et al. [2014] T.O. Jelly, S.Y. Jung, and T.A. Zaki.  Turbulence and skin friction modification in channel flow with streamwise-aligned superhydrophobic surface texture.  _Phys. Fluids_ , 26:095102, 2014.
  * Jongen and Gatski [1998] T. Jongen and T.B. Gatski.  General explicit algebraic stress relations and best approximation for three-dimensional flows.  _Int. J. Eng. Sci._ , 36(7-8):739–763, 1998.
  * Kader [1981] B.A. Kader.  Temperature and concentration profiles in fully turbulent boundary layers.  _Int. J. Heat Mass Transf._ , 24(9):1541–1544, 1981.
  * Kader and Yaglom [1972] B.A. Kader and A.M. Yaglom.  Heat and mass transfer laws for fully turbulent wall flows.  _Int. J. Heat Mass Trans._ , 15(12):2329–2351, 1972.
  * Kaller et al. [2019] T. Kaller, V. Pasquariello, S. Hickel, and N.A. Adams.  Turbulent flow through a high aspect ratio cooling duct with asymmetric wall heating.  _J. Fluid Mech._ , 860:258–299, 2019.
  * Kays and Crawford [1993] W.M. Kays and M.E. Crawford.  Convective heat and mass transfer, 3rd edn., 1993.
  * Kim and Moin [1985] J. Kim and P. Moin.  Application of a fractional-step method to incompressible Navier-Stokes equations.  _J. Comput. Phys._ , 59:308–323, 1985.
  * Launder and Ying [1973] B.E. Launder and W.M. Ying.  Prediction of flow and heat transfer in ducts of square cross-section.  _Proc. Instn. Mech. Engrs._ , 187(1):455–461, 1973.
  * MacDonald et al. [2019] M. MacDonald, N. Hutchins, and D. Chung.  Roughness effects in turbulent forced convection.  _J. Fluid Mech._ , 861:138–162, 2019.
  * Marin et al. [2016] O. Marin, R. Vinuesa, A.V. Obabko, and P. Schlatter.  Characterization of the secondary flow in hexagonal ducts.  _Phys. Fluids_ , 28(12):125101, 2016.
  * Marusic et al. [2013] I. Marusic, J.P. Monty, M. Hultmark, and A.J. Smits.  On the logarithmic region in wall turbulence.  _J. Fluid Mech._ , 716:R3, 2013.
  * Marusic et al. [2017] I Marusic, WJ Baars, and N Hutchins.  Scaling of the streamwise turbulence intensity in the context of inner-outer interactions in wall turbulence.  _Phys. Rev. Fluids_ , 2:100502, 2017.
  * Modesti [2020] D. Modesti.  A priori tests of eddy viscosity models in square duct flow.  _Theor. Comput. Fluid Dyn._ , 34(5):713–734, 2020\.
  * Modesti and Pirozzoli [2016] D. Modesti and S. Pirozzoli.  Reynolds and Mach number effects in compressible turbulent channel flow.  _Int. J. Heat Fluid Flow_ , 59:33–49, 2016.
  * Modesti and Pirozzoli [2018] D. Modesti and S. Pirozzoli.  An efficient semi-implicit solver for direct numerical simulation of compressible flows at all speeds.  _J. Sci. Comput._ , 75(1):308–331, 2018.
  * Modesti et al. [2018] D. Modesti, S. Pirozzoli, P. Orlandi, and F. Grasso.  On the role of secondary motions in turbulent square duct flow.  _J. Fluid Mech._ , 847, 2018.
  * Nikuradse [1930] J. Nikuradse.  Turbulente strömung in nicht-kreisförmigen rohren.  _Ing. Arch._ , 1:306–332, 1930.
  * Orlandi [2012] P. Orlandi.  _Fluid flow phenomena: a numerical toolkit_ , volume 55.  Springer Science & Business Media, 2012.
  * Orlandi et al. [2018] P. Orlandi, D. Modesti, and S. Pirozzoli.  DNS of turbulent flows in ducts with complex shape.  _Flow, Turbul. Combust._ , 100(4):1063–1079, 2018\.
  * Pallares and Davidson [2002] J. Pallares and L. Davidson.  Large-eddy simulations of turbulent heat transfer in stationary and rotating square ducts.  _Phys. Fluids_ , 14(8):2804–2816, 2002.
  * Peet and Sagaut [2009] Y. Peet and P. Sagaut.  Theoretical prediction of turbulent skin friction on geometrically complex surfaces.  _Phys. Fluids_ , 21:105105, 2009.
  * Piller [2005] M. Piller.  Direct numerical simulation of turbulent forced convection in a pipe.  _Int. J. Numer. Methods Fluids_ , 49(6):583–602, 2005.
  * Pinelli et al. [2010] A. Pinelli, M. Uhlmann, A. Sekimoto, and G. Kawahara.  Reynolds number dependence of mean flow structure in square duct turbulence.  _J. Fluid Mech._ , 644:107–122, 2010.
  * Pirozzoli [2010] S. Pirozzoli.  Generalized conservative approximations of split convective derivative operators.  _J. Comput. Phys._ , 229(19):7180–7190, 2010\.
  * Pirozzoli and Bernardini [2013] S. Pirozzoli and M. Bernardini.  Probing high-Reynolds-number effects in numerical boundary layers.  _Phys. Fluids_ , 25:021704, 2013.
  * Pirozzoli et al. [2016] S. Pirozzoli, M. Bernardini, and P. Orlandi.  Passive scalars in turbulent channel flow at high Reynolds number.  _J. Fluid Mech._ , 788:614–639, 2016.
  * Pirozzoli et al. [2018] S. Pirozzoli, D. Modesti, P. Orlandi, and F. Grasso.  Turbulence and secondary motions in square duct flow.  _J. Fluid Mech._ , 840:631–655, 2018.
  * Pirozzoli et al. [2021] S. Pirozzoli, J. Romero, M. Fatica, R. Verzicco, and P. Orlandi.  One-point statistics for turbulent pipe flow up to Re${}_{\tau}\approx 6000$.  _J. Fluid Mech._ , 926:A28, 2021.
  * Pope [1975] S.B. Pope.  A more general effective-viscosity hypothesis.  _J. Fluid Mech._ , 72(2):331–340, 1975.
  * Pope [2000] S.B. Pope.  _Turbulent flows_.  Cambridge University Press, 2000.
  * Prandtl [1927] L. Prandtl.  Über die ausgebildete turbulenz.  In _Verh. 2nd Int. Kong. für Tech. Mech._ , 1927.
  * Rohsenow et al. [1998] W.M. Rohsenow, J.P. Hartnett, and Y.I. Cho, editors.  _Handbook of heat transfer_ , volume 3.  McGraw-Hill New York, 1998.
  * Schindler et al. [2019] A Schindler, BA Younis, and B Weigand.  Large-eddy simulations of turbulent flow through a heated square duct.  _Int. J. Therm. Sci._ , 135:302–318, 2019.
  * Sekimoto et al. [2011] A. Sekimoto, G. Kawahara, K. Sekiyama, M. Uhlmann, and A. Pinelli.  Turbulence-and buoyancy-driven secondary flow in a horizontal square duct heated from below.  _Phys. fluids_ , 23(7):075103, 2011.
  * Shah and Sekulib [1998] R K. Shah and D.R. Sekulib.  Heat exchangers.  In W.M. Rohsenow, J.R. Hartnett, and I.C. Young, editors, _Handbook of heat transfer_ , chapter 17, pages 17.76–17.80. McGraw-Hill, New York, USA, 1998.
  * Shah and London [2014] R.k. Shah and A.L. London.  _Laminar flow forced convection in ducts: a source book for compact heat exchanger analytical data_.  Academic Press, 2014.
  * So et al. [2004a] R.M.C. So, L.H. Jin, and T.B. Gatski.  An explicit algebraic Reynolds stress and heat flux model for incompressible turbulence: Part I Non-isothermal flow.  _Theor. Comput. Fluid Dyn._ , 17(5):351–376, 2004a.
  * So et al. [2004b] R.M.C. So, L.H. Jin, and T.B. Gatski.  An explicit algebraic reynolds stress and heat flux model for incompressible turbulence: Part II Buoyant flow.  _Theor. Comput. Fluid Dyn._ , 17(5):377–406, 2004b.
  * Spalart et al. [2018] P.R. Spalart, A. Garbaruk, and A. Stabnikov.  On the skin friction due to turbulence in ducts of various shapes.  _J. Fluid Mech._ , 838:369–378, 2018.
  * Vázquez and Métais [2002] M.S. Vázquez and O. Métais.  Large-eddy simulation of the turbulent flow through a heated square duct.  _J. Fluid Mech._ , 453:201–238, 2002.
  * Vinuesa et al. [2014] R. Vinuesa, A. Noorani, A. Lozano-Durán, G.K.E. Khoury, P. Schlatter, P.F. Fischer, and H.M. Nagib.  Aspect ratio effects in turbulent duct flows studied through direct numerical simulation.  _J. Turbulence_ , 15(10):677–706, 2014.
  * Wang et al. [2007] B-C Wang, E. Yee, J. Yin, and D.J. Bergstrom.  A general dynamic linear tensor-diffusivity subgrid-scale heat flux model for large-eddy simulation of turbulent thermal flows.  _Numer. Heat Tr. B-Fund._ , 51(3):205–227, 2007\.
  * White and Majdalani [2006] F.M. White and J. Majdalani.  _Viscous fluid flow_.  McGraw-Hill, New York, 2006.
  * Yang et al. [2009] H. Yang, T. Chen, and Z. Zhu.  Numerical study of forced turbulent heat convection in a straight square duct.  _Int. J. Heat Mass Transf._ , 52(13-14):3128–3136, 2009.
  * Younis et al. [2005] B.A. Younis, C.G. Speziale, and T.T. Clark.  A rational model for the turbulent scalar fluxes.  _Proc. R. Soc. A: Math., Phys. Eng. Sci._ , 461(2054):575–594, 2005.