12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

DIRECT NUMERICAL SIMULATION OF FORCED THERMALCONVECTION IN SQUARE DUCTSDavide Modesti
Faculty of Aerospace Engineering
Delft University of Technology
Kluyverweg 2, 2629 HS Delft, The Netherlands
d.modesti@tudelft.nl

Sergio Pirozzoli
Dipartimento di Ingegneria Meccanica e Aerospaziale
Sapienza Università di Roma
via Eudossiana 18, 00184 Roma, Italia
sergio.pirozzoli@uniroma1.it

ABSTRACT

is duct side length and also duct hydraulic diameter, ub is thebulk velocity, and ν the fluid kinematic viscosity) 33000 and
67000. Close similarity between the wall shear stress and theheat flux distributions was shown (with quoted discrepancy of±2%), which the authors connected with similar mixing actionof the secondary currents on momentum and temperature. As
a result, they observed that the ratio of the average friction andheat-transfer coefficients for a square duct is approximately thesame as for a circular pipe. Measurements of the wall-normal
temperature profiles highlighted close universality when thelocal wall heat flux is used for normalization, and the presenceof a sizeable logarithmic layer, with extrapolated value of thescalar Kármán constant of κθ ≈ 0.51. Those findings werequalitatively 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 Kármán constant was κθ ≈ 0.46, hence more similar to the values generally quoted for circular pipe flow, namely κθ ≈ 0.47 (Kader &Yaglom, 1972).
Early computational studies of heat transfer in squareducts mostly relied on the use of RANS models (Launder &Ying, 1973), but heat transfer coefficient showed general underprediction by about 10%. High-fidelity computational studies of heat transfer in square ducts have been quite limited
so far. Vázquez & Métais (2002) first studied turbulent flow
through a heated square duct by means of large-eddy simulation (LES) at bulk Reynolds number Reb = 6000, considering
the case in which one of the walls is hotter than the other three.
Accounting for fluid viscosity variation with temperature, theyfound that turbulent structures near the hot wall become larger
than near the other walls, in such a way that wall scaling issatisfied.
In this paper we study heat transfer in fully developedsquare duct flow with uniform internal heating and isothermal
walls, by carrying out DNS at unit molecular Prandtl number(defined as the ratio of the kinematic viscosity to the scalardiffusivity, Pr = ν/α), and at much higher Reynolds number than in previous computational studies, The present study

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 Reτ ≈ 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. We find good agreement between
the Nusselt number of square duct and circular pipe flow when
the Reynolds number based on the hydraulic diameter is used,thus corroborating the common engineering practice. Popularengineering correlations for the heat transfer reveal deviations
up to 5% with respect to DNS data, which are nicely fitted bya power law.

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 havebeen carried out in the past, and a variety of analytical andsemi-empirical prediction tools have been developed, which
are extensively reported in classical textbooks (Kays & Crawford, 1993). Most studies have been carried out for the canonical case of ducts with circular cross section, whereas muchless is known about the case of ducts with more complex geometry, which also have great practical relevance, for instancein water draining or ventilation systems, nuclear reactors, heatexchangers, 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 ofthe duct (Kays & Crawford, 1993). However, the large scatterin experimental data makes it difficult to quantify the actual
accuracy of semi-empirical prediction formulas.
Heat transfer in square ducts was first studied experimentally by Brundrett & Burroughs (1967), who considered air
flow at bulk Reynolds number Reb = Hub /ν (where H = 2h1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Table 1. Flow parameters for square duct DNS. Box dimensions are 6πh × 2h × 2h for all flow cases. Reb = 2hub /ν is thebulk Reynolds number, Re∗τ = hu∗τ /ν is the friction Reynolds
number, and Nu is the Nusselt number. Nx × Ny × Nz are themesh points in the respective coordinate directions.

Case

Reb

Re∗τ

Nu

Nx

Ny

Nz

A

4410

150

18.1

512

128

128

B

7000

227

26.9

640

144

144

C

12000

365

41.9

768

208

208

D

17800

519

57.6

1024

256

256

E

40000

1055

109.0

2048

512

512

F

84000

2041

197.0

4096

1024

1024

streamwise velocity is populated with streaks, resulting fromsweep and ejections visible in the cross-stream plane. In the
cross-stream plane we note bulges of high-speed and hot flow
spanning two opposite walls, and interacting with the smaller
structures close to the wall. We immediately note a strong similarity between streamwise velocity and temperature, with themain difference that the latter reveals much finer structures dueto the absence of the smoothing action of the pressure gradient.
The picture that emerges from the instantaneous snapshots issimilar to the one of other canonical wall flows, instead, themean flow field is characterized by the presence of secondary
flows in the cross-stream plane.
Figure 2 depicts the mean temperature (contours) and velocity (iso-lines) fields in the duct cross section, along withrepresentative 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, which are relatively depleted with those. As a result, both the temperature contours
and the velocity isolines bend towards the corners, featuring abulging 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 andvelocity contours bear close similarity, being nearly coincident
near the wall, and retaining the same shape farther off.
Figure 3 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 global wall units (basedon the perimeter-averaged wall-shear stress). For reference
purposes, the mean temperature profiles obtained with the experimental fitting by Kader (1981) at matching Reτ are also
reported. 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 yields good
universality away from walls, and near coincidence with thepipe 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 thenonuniform wall heat flux. Inspection of figure 3 shows thattransition to the ’global’ scaling (which is controlled by thespatially uniform heat source) occurs at a wall distance ofabout 0.2h, which is also the lower limit for the core region
in canonical flows (Pope, 2000). 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 andglobal scaling.
The global heat transfer performance of the duct is quantified in terms of the Stanton number, St∗ = 1/(u∗b θm∗ ), whereub and θm are the bulk flow velocity and bulk temperature temperature,is the continuation of previous efforts (Pirozzoli et al., 2018;
Modesti et al., 2018) targeted to studying turbulent flows insquare ducts by means of DNS.

1METHODOLOGY

In the present work we use a fourth-order co-locatedfinite-difference solver, previously used for DNS of compressible turbulence (Bernardini et al., 2021). Here, 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 (Modesti & Pirozzoli, 2018). The
streamwise momentum equation is forced in such a way as tomaintain a constant mass flow rate (the spatially uniform driving term is hereafter referred to as Π), periodicity is exploitedin the streamwise direction, and isothermal no-slip boundary
conditions are used at the channel walls. Six DNS have been
carried out at bulk Mach number Mb = ub /cw = 0.2 (where
cw is the speed of sound at the wall temperature), and bulk
Reynolds number Reb = 4400 − 84000 (table 1), and hereafter labeled with letters from A to F. Navier–Stokes equations
are augmented with the transport equation for a passive scalar
field, with molecular diffusivity α = ν, in such a way that themolecular Prandtl number is unity for all simulations. Similarly to the streamwise velocity field, the passive scalar equation is also forced with a time-varying, spatially homogeneous
forcing term Q, in such a way that its mean value is maintainedin time.
The + superscript is used to denote quantities madenondimensional with respect to the local wall friction, namely
with uτ = (τw /ρ)1/2 (τw = νdU/dy|w is the local wall shear
stress), and δv = ν/uτ , and the ∗ superscript is used to denotequantities made nondimensional with respect to the perimeteraveraged friction, τw∗ = h hΠi /2, namely u∗τ = (τw∗ /ρ)1/2 , andδv∗ = ν/u∗τ . Likewise, for normalization of the temperaturefield we will consider either the local friction temperatureθτ = qw /uτ , where qw = αdΘ/dy|w is the local wall heat flux,
or the global friction temperature θτ∗ .

θm =

1
ub Ac

Z
Ac

U(Θ − θw )dAc ,ub =1
UdAc .
A c AcZ(1)

Additionally, we also quantify the heat transfer using theNusselt number, Nu∗ = St∗ Reb Pr.
In figure 4 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

RESULTS
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. Instantaneous flow structure are not different from the ones typical of
canonical wall flows, and the near wall of both temperature and212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
2h = D, hence provided the Reynolds number based on the hydraulic diameter is the same. It is interesting that equation (2)is basically arrived at by neglecting the local wall shear stressand heat flux variation along the duct perimeter, and disregarding the flow deceleration at corners. Apparently, these effects
very nearly cancel out upon integration.

ducts. We note that the latter two datasets were obtained forPr√= 0.71, hence the Nu data have been rescaled by a factor
1/ Pr to compare with the present ones. As a reference, we
also report correlations widely used in the engineering practice, including that by Gnielinski (1976), from Kays & Crawford (1993) and direct fitting the present DNS data with apower-law expression which yields, Nu = 0.0216 Reb0.805 . The
DNS data show good correspondence with the pipe DNS data
at matching Reb , 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 squareducts 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 interms of the Nusselt number is used, as in figure 4(b), hence
we believe that the 1/St representation should be used whenrelatively small differences must be discriminated. Classicalcorrelations 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 basedon old experimental data, or to the fact that those are mainly
trained for the Pr = 0.71 case, whereas here Pr = 1. Additionally, discrepancies can be attributed to our heat forcing
scheme, in which a spatially uniform heat source is prescribed,
which tends to slightly overpredict the heat flux as comparedto other approached in which the total heat flux is kept strictly
constant in time (Abe & Antonia, 2017). Slight adjustment of
the Kays–Crawford power-law formula coefficients seems toyield very good representation of the DNS data.
An issue which deserves further investigation is why useof 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 withnear 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/κ log y∗ + C, Θ∗ = 1/κθ log y∗ +Cθ , and integrating over the duct cross section, the following
expression for the inverse Stanton number results

Z hZ z8
1
1/St =
log y∗ 2
κκ
4Re∗2
0
0
θ
τ



C
∗ Cθ
+ log y+CCθ dy∗ dz∗ =
+
κ
κθ


1
7
log2 Re∗τ + 3 log Re∗τ +
κκθ
2



Cθ
C
3
+
+
log Re∗τ −+CCθ .
κκθ
2

CONCLUSIONS AND FUTURE WORK
We have carried out DNS of square duct flow at theunprecedented Reynolds number of Reτ ≈ 2000, wherein
Navier–Stokes equations have been augmented with the transport of a passive scalar representing the temperature field. Thisconfiguration corresponds to the case of forced convection inwhich hot fluid is pumped through the duct and cooled at thewalls. In this respect, we find that the streamwise velocity andtemperature fields are very similar, both instantaneously andon average.
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. An extensive analysis of this DNS
dataset is available in our recent publication (Modesti & Pirozzoli, 2022), where we also discuss the temperature fluctuations
and the accuracy of eddy viscosity models for the turbulent
heat transfer. Future effort will be devote to study the flow inasymmetrically heated square and rectangular ducts, which is
a common configuration in heat exchangers.

∗(2)

Equation (2) should be compared with the corresponding
expression for a circular duct with diameter D,17
1/St =
log2 Reτ + 3 log Reτ +
κκθ
2



Cθ
C
3
+
+
log Reτ −+CCθ ,
κκθ2(3)

where Reτ = Duτ /(2ν). The two expressions are identical for3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
REFERENCES

heat transfer in ducts of square cross-section. Proc. Instn. Mech. Engrs. 187 (1), 455–461.
Modesti, D. & Pirozzoli, S. 2018 An efficient semi-implicit
solver for direct numerical simulation of compressible flows
at all speeds. J. Sci. Comput. 75 (1), 308–331.
Modesti, D. & Pirozzoli, S. 2022 Direct numerical simulationof forced thermal convection in square ducts up to Reτ ≈2000. J. Fluid Mech. 941, A16.
Modesti, D., Pirozzoli, S., Orlandi, P. & Grasso, F. 2018 On
the role of secondary motions in turbulent square duct flow.
J. Fluid Mech. 847.
Nikuradse, J. 1930 Turbulente strömung in nichtkreisförmigen rohren. Ing. Arch. 1, 306–332.
Pallares, J. & Davidson, L. 2002 Large-eddy simulations of
turbulent heat transfer in stationary and rotating squareducts. Phys. Fluids 14 (8), 2804–2816.
Pirozzoli, S., Modesti, D., Orlandi, P. & Grasso, F. 2018 Turbulence and secondary motions in square duct flow. J. Fluid
Mech. 840, 631–655.
Pirozzoli, S., Romero, J., Fatica, M., Verzicco, R. & Orlandi,
P. 2021 One-point statistics for turbulent pipe flow up toReτ ≈ 6000. J. Fluid Mech. 926, A28.
Pope, S.B. 2000 Turbulent flows. Cambridge University Press.
Vázquez, M.S. & Métais, O. 2002 Large-eddy simulation of
the turbulent flow through a heated square duct. J. Fluid
Mech. 453, 201–238.

Abe, H. & Antonia, R.A. 2017 Relationship between the heattransfer law and the scalar dissipation function in a turbulent
channel flow. J. Fluid Mech. 830, 300–325.
Bernardini, M., Modesti, D., Salvadore, F. & Pirozzoli, S.
2021 STREAmS: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows.
Comput. Phys. Commun. 263, 107906.
Brundrett, E. & Burroughs, P.R. 1967 The temperature innerlaw and heat transfer for turbulent air flow in a verticalsquare duct. Int. J. Heat Mass Transfer 10, 1133–1142.
Gnielinski, V. 1976 New equations for heat and mass transfer
in turbulent pipe and channel flow. Int. Chem. Eng. 16, 359–
367.
Hirota, M., Fujita, H., Yokosawa, H., Nakai, H. &Itoh, H. 1997 Turbulent heat transfer in a square duct.
Int. J. Heat Fluid Flow 18, 170–180.
Kader, B.A. 1981 Temperature and concentration profiles infully turbulent boundary layers. Int. J. Heat Mass Transf.
24 (9), 1541–1544.
Kader, B.A. & Yaglom, A.M. 1972 Heat and mass transfer
laws for fully turbulent wall flows. Int. J. Heat Mass Trans.
15 (12), 2329–2351.
Kays, W.M. & Crawford, M.E. 1993 Convective heat and masstransfer, 3rd edn.
Launder, B.E. & Ying, W.M. 1973 Prediction of flow and412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

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

(a)

(b)

(c)

(e)

(f )

y/h

1

0.5

0

(d )

y/h

1

24
20
16
12
8
4
0

0.5

0
0

0.5

z/h

1

0

0.5

1

z/h

0

0.5

U ∗ , Θ∗1z/hFigure 2. Mean streamwise temperature (Θ∗ , 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. Onlya quarter of the full domain is shown. The dashed diagonal lines indicate the corner bisector.

5Θ+

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

(a)

(b)

(c)

25

25

25

20

20

20

15

15

15

10

10

10

5

5

5

0 0
10

10

1

2

3

10

10

0 0
10

10

1

2

10

3

10

0 0
10

(e)

(f )

25

25

20

20

20

15

15

15

10

10

10

5

5

5

Θ+

(d )
25

0 0
10

101

102

y

103

0 0
10

101

102

+

y

103

0 0
10

10

1

2

10

1

2

10

+

10

10

y

3

3

10

+

Figure 3. 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 Re∗τ .

600

1000

Nu*

1/St*

500
400

100

300
200

10000

10

100000Reb

10000

100000Reb

Figure 4. Inverse 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 & Davidson, 2002), and fromexperiments (Brundrett & Burroughs, 1967, triangles). The dashed line denotes the correlation by Kays & Crawford (1993), and the
dot-dashed lines the correlation Gnielinski (1976).

6