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

INFLUENCE OF WALL SHEAR STRESS ON THE SECONDARY FLOW INSQUARE DUCTS

Alexander Doehring
Technical University of Munich
TUM School of Engineering and Design
Chair of Aerodynamics and Fluid Mechanics
Boltzmannstr. 15, 85748 Garching, Germany
alex.doehring@tum.deThomas Kaller and Steffen J. Schmidt and Nikolaus A. Adams
Technical University of Munich
TUM School of Engineering and Design
Chair of Aerodynamics and Fluid Mechanics
Boltzmannstr. 15, 85748 Garching, Germany
thomas.kaller@tum.de
steffen.schmidt@tum.denikolaus.adams@tum.de

ABSTRACT

at one of the four walls.
Turbulent flows through rectangular ducts are affected by turbulence induced secondary flow, also known as Prandtl’s flowof the second kind. The flow is characterized by a vortex system consisting of a counter rotating vortex pair in each corner.
According to Salinas Vázquez & Métais (2002) the secondary
flow of the second kind is relatively weak with a strength ofapproximately 1 − 3% of the bulk flow, but it increases the
mixing by transporting high momentum fluid from the duct
core to the corners. The secondary flow is induced due to theanisotropy of the Reynolds stress tensor (Demuren & Rodi
(1984)). Thus, turbulence models using the isotropic eddy
viscosity ansatz fail. Demuren & Rodi (1984) analysed thegeneration of Prandtl’s flow of second kind using the streamwise vorticity equation. They determined that gradients of theanisotropic Reynolds stress tensor act as source term for thevorticity production. The influence of the duct geometry on
corner vortex system has been investigated by Vinuesa et al.
(2014), performing DNS of different aspect ratio duct flows.
They observed an array of secondary vortices along the long
sidewalls. Vidal et al. (2017) investigated the influence ofround corners in duct flows on the secondary flow. Unexpectedly, rounding the corner did not reduce the secondary flow,
but the magnitude of the cross-flow was similar compared to
sharp corners ducts.
Modelling difficulties were shown by Vane & Lele (2015) performing wall-modelled LES with an equilibrium wall-modelfor different duct flows. Their results were improved by introducing the Buleev length scale and adjusting the coupling
point, but the wall shear stress distribution was not matching
the DNS results. For this reason, they obtained a shifted vortex
system, as it is strongly affected by the wall shear stress distribution. Based on our well-resolved LES we are investigating
the role of the wall shear stress for the secondary flow.

We investigated turbulent duct flows with a square crosssection using well-resolved large-eddy simulations (LES). A
physically consistent subgrid-scale turbulence model is incorporated, that is based on the Adaptive Local Deconvolution
Method (ALDM) for implicit LES. The wall shear stress is
artificially modified at one of the four walls to mimic different surface properties. A direct numerical simulation (DNS)
of a symmetrical duct flow is used as reference to asses theinfluence of a modified wall shear stress. The modification
results in an asymmetrical distribution of the secondary flow
source terms, affecting the momentum distribution. Further,
the anisotropy of the Reynolds stress tensor, which induces thesecondary flow vortices is considerably affected by the wall
shear stress change.

Introduction
Turbulent flow within a straight duct with a rectangular cross section is of great importance for many engineering
applications. This includes for instance ventilation systems,
cooling channels in rocket engines, nuclear reactors and heatexchangers. A thorough understanding of the turbulent flow isnecessary in order to predict the cooling efficiency, the pressure drop and the system lifetime. These properties are highly
affected by the specific manufacturing technology that is employed, since it determines the wall surface properties. Especially additive manufacturing processes provide the possibility
for new design options of a product, but introduce rougher surfaces. The roughness parameters are affected by the production parameters as for instance the build direction. In order toinvestigate the effect of an asymmetric wall roughness variation on the square duct flow field, we performed well-resolved
LES simulations with an artificially modified wall shear stress1

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

Y2

Z2 2h

Z1

y
z

Figure 2: Intended mean wall shear distribution for dif-ferent cases. The perimeter-averaged wall shear stressτw∗ of LES P is used for normalization. The lines correspond to LES P (), LES T+ (), LES T), LES L+ (), LES L- (
).
(

Y1
2h

Figure 1: Cross section of the computational domain.

The walls are labelled and the modified lower y-wall(Y1) is indicated green. The stars are locations which
are used for further evaluation with z/h = −0.9, z/h =−0.5, z/h = 0 from left to rigth.

Table 1: Summary of resolution parameters.

Numerical ModelLES simulations are performed solving the threedimensional, fully compressible Navier-Stokes Equations.
The ideal gas law is used to close the system and the viscosity is modelled by the Sutherland law. For time advancement
an explicit second-order low-storage four-stage Runge-Kutta
method with enhanced stability region is applied (Schmidt
et al. (2006)). The governing equations are spatially discretized on a block structured, curvilinear grid using the finitevolume method. Viscous fluxes are discretized by a linear
second-order centered scheme. Convective fluxes are determined using a compact four cell stencil by Egerer et al. (2016).
The switch between the two schemes is made by sensor functionals. The convective flux calculation is based on the Adaptive Local Deconvolution Method (ALDM), which is a nonlinear finite volume method, that provides a physically consistent subgrid-scale turbulence model for implicit LES, see
Hickel et al. (2006, 2014).

h+ymin, Y1 × h+ymin, Y2

1.554215 × 1.130634

h+ymax, Y1 × h+ymax, Y2

31.114517 × 22.634663

h+zmin, Z1 × h+zmin, Z21.122118 × 1.119152

h+zmax, Z1 × h+zmax, Z2

22.464188 × 22.464188

with the friction velocity u2τ = (τw /ρw ) and the wall shear
stress τw = (µ∂ u/∂ y)∣w . Here, the wall shear stress is averaged over each wall and the resolution is based on the wholecell size, but the flow variables are evaluated at the cell center.
Therefore, the effective minimum wall distance is ∆ymin /2.
The values in the table represent the highest values obtained
from all five simulations. The domain is shown in figure 1.
The walls are labelled and the modified wall Y1 is highlighted
in green. A forcing source term in the momentum and energy
equations ensures a constant bulk velocity and the respective
bulk Reynolds number. The additional body force is based on
Brun et al. (2008). First, a LES of case C is performed without modifications labelled LES P (Pirozzoli), that serves as areference. Four more LES are performed with a modified wall
shear stress τw,mod (Y 1) = C(y) × τw,orig . Two LES simulations
contain a wall shear stress which is augmented (LES T+) anddecreased (LES T-) by a factor of two, compared to LES P.
Thus, C(y) = 2 and C(y) = 0.5. In the remaining two LES thewall shear stress is linearly increasing (LES L+) and decreasing (LES L-). The intention is to obtain a linear distribution
with a minimum and maximum value at the center, which corresponds to simulation LES T+ and LES T-. The targeted wall
shear stresses for the performed LES are shown in figure 2.
Note, that the actual profiles will decrease towards z/h = ±1
due to the influence of the sidewalls.

Setup
We are investigating a generic duct configuration, focusing on the influence of an asymmetric wall shear stress modulation on secondary flow. For this purpose, case C from the
DNS study by Pirozzoli et al. (2018) with a bulk Reynolds
number of Reb = 17800 and a friction Reynolds number Reτ =519 is used. The wall shear stress of the lower wall in y direction is artificially changed by multiplying the obtained shear
stress within the simulation by a constant factor. This adjustment is used to mimic a variation of surface properties. Periodic boundary conditions are used in streamwise direction.
The walls are isothermal and a no slip boundary condition isapplied. The bulk pressure is set to pb = 1.013 25 bar. The
duct geometry is 10h × 2h × 2h in streamwise, wall-normal
and spanwise direction with the duct half height h, respectively. We use a hyperbolic mesh distribution in wall-normal
directions and a uniform grid in the streamwise direction.
All simulations have been performed on the same grid with
Nx = 117, Ny = 88, Nz = 88 cells. Table 1 provides the obtained
resolution with respect to wall units ∆y+ = h+y = ∆yρw uτ /µw ,Results
At first we evaluate the mean flow field by averaging
in streamwise direction and in time after reaching a quasistationary state. Mean values are indicated with an overbar
and fluctuations with a prime. The mean velocity profiles arepresented in figure 3. The profiles are extracted at wall Y12

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

Figure 3: Scaled mean velocity u+ over wall units y+ at wall Y1. The black lines indicate LES P and the blue lines themodified wall with a) LES T+, b) LES T-, c) LES L+, d) LES L-. Different line styles correspond to different wall), z/h = −0.5 () and z/h = −0.90 (). The law of the wall () ispositions including z/h = 0 (included with u+ = 1/0.41 ⋅ lny+ + 5.2.

Figure 4: Mean wall shear stress at wall Y1 in a) and wall Y2 in b). The profiles are normalized with the perimeter-averaged

wall shear stress τw∗ of each individual LES. The lines correspond to LES P (LES L+ (), LES L- ().
at three different positions, z/h = 0, z/h = −0.5 and z/h = −0.9.
The law of the wall velocity profile consists of the viscous sublayer u+ = y+ and logarithmic region u+ = 1/κ ⋅ ln y+ + B, with
κ = 0.41 and B = 5.2. Based on LES P it is observed, that at
the center of wall Y1 the normalized velocity profile follows
the law of the wall. Approaching the side wall results in a flattening of the profile. At z/h = −0.5 the profile is identical withthe profile at the channel center until y+ ≈ 200, with a flattening in the outer layer due to secondary flow. This can also be
observed for the constant value cases LES T+ and LES T- in
figure 3a,b). For z/h = −0.9 the flattening affects the whole
logarithmic layer. Both cases with an increased wall shear
stress (LES T+, LES L+) feature a downward shift of the velocity profile as it can be observed in rough channel flows, cf.
Thakkar et al. (2018). In general, an approach of the wall shear
stress towards the value of LES P results in an approach of thevelocity profile towards the law of the wall superimposed withthe local secondary flow effect.

), LES T+ (), LES T- (),

wall shear stress (τw∗ ) of each individual LES which is usedhere for normalization. For instance, τw∗ of LES T+ is higher
due to the modification compared to τw∗ of LES P resulting ina maximum value of ≈ 1.6 and not 2. Although the secondary
flow is not strong, it is capable to alter the wall shear stress distribution. This can be seen in terms of the dent at z/h ≈ −0.7
leading to a peak value close to the corner and at z/h ≈ −0.5.
This is also observed in DNS studies as for example in Pirozzoli et al. (2018) and Zhang et al. (2015).
Due to the one-sided modification of the wall shear√ stress,the mean cross-sectional velocity distribution uc f = v2 + w2
becomes asymmetrical compared to the symmetrical LES P,
compare figures 5a)-d) with 5e). When the wall shear stressis reduced the vortices at the modified wall Y1 are getting
smaller, weaker and are pushed towards the corners, whereas
the vortices at the opposite wall Y2 are expanding into thelower duct half, see 5b). The vortices at the sidewalls Z1 andZ2 are getting stronger and bigger. The opposite effect occurs
for an increase of the wall shear stress at Y1. A linear C(y)distribution intensifies the previous observations made at Y1
and Y2. Interestingly, the intensification is not observed forthe vortices at the sidewalls, since they are weaker comparedto the constant value modification. This coincides with the

The obtained wall shear stress at wall Y1 and Y2 is presented in figure 4. As already mentioned, the stress is decreasing towards the sidewalls at z/h = ±1. Furthermore, the obtained absolute stress values do not match the intended ones,
cf. figure 2. This is a consequence of the perimeter-averaged

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022√
v2 + w2 . a) LES T+, b) LES T-, c) LES L+, d) LES Land e) LES P. The counter-rotating corner vortex system is visualized by velocity vectors.
Figure 5: Mean cross-sectional velocity distribution with uc f =4

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

Figure 6: Mean secondary flow v and w normalized with the bulk velocity ub at Y1 for z/h = −0.5. The lines correspond to

LES P (), LES T+ (et al. (2018) ( ).

), LES T- (), LES L+ (), LES L- (), and DNS data from Pirozzoli

Figure 7: Distribution of diagonal components a) and off-diagonal components b) of the Reynolds stress anisotropy tensor

at z/h = −0.5. The lines correspond to LES P (), LES T+ (), LES T- ().

Kronecker delta. Figure 7 presents the components of theanisotropy tensor influenced by the wall shear stress. Patelet al. (2016) demonstrated that any change in viscosity or density has only a marginal effect on the wall normal component
a22 . This also holds for changes in the wall shear stress up to
the buffer layer at approximately y+ = 10. An increase of τw
results in a higher anisotropy value for the streamwise component a11 and a decrease for the tangential component a33
compared to the unmodified LES. The opposite is the case fora decrease of the wall shear stress. The ratio between the turbulent kinetic energy and turbulent shear stress indicates themomentum transfer in figure 7b). Using a higher wall shear
stress as in case LES T+ causes a higher momentum transfer.
Further, we apply the anisotropy-invariant map (AIM) in figure 8 (Emory & Iaccarino (2014)). Here, the Lumley triangle
defines the borders of all possible states and the 2nd invariant(I2 = ai j a ji /2) and the 3rd invariant (I3 = ai j a jn ani /3) of the
Reynolds stress tensor are used as axis. A higher anisotropy
for an increased τw is also observed in AIM, since LES T+
starts closer and gets closer to the one-component (1-C) limit.
The opposite is the case for a decrease of the wall shear stress.
In addition, a kink is present for LES + starting at the 2component edge close to the 1-C limit, moving along the edge
towards the 2-component axis-symmetric limit(2-C) and turning again towards 1-C, see zoom in figure 8a). This phenomenon is also observed close to the sidewall at z/h = −0.9 infigure 8b) but not for LES T- and LES P.

wall shear stress values of the linear modified LES simulations
close to the sidewalls, that are approaching LES P.
The effects of the wall shear stress modification on themean secondary flow velocities v and w are shown in figure 6
at z/h = −0.5. DNS data from Pirozzoli et al. (2018) are included to assess the numerical accuracy of LES P. The maximum of the wall normal velocity v close to Y1 refers of the two
smaller corner vortices pushing fluid upwards and the following minimum refers to the two larger corner vortices pushingfluid downwards. The tangential velocity w transports fluid
from the duct corner to the midplane. LES P is in good agreement with the DNS data with a slight overestimation of the vortex intensity close to the wall. The modification of τw mainly
changes the intensity of w whereas the peak position remains
fixed and for v both the intensity and peak position are altered.
These observations can be similarly made for asymmetricalwall heating, which modifies the wall shear stress as well as
the turbulence intensities due to viscosity modulation. Kaller
et al. (2019) showed that a heated wall for water is weakening
the intensity of the counter rotating vortices. An τw increase
moves the vortices closer into the corner and a τw decreaseaway from the corner, see 6a). This is similarly shown in thecontour plots of figure 5. The observed shrinkage of the vortex at Y1 from LES P over LES T- to LES L- in figure 5 isreflected in the reduction of the maximum value for w.
A study by Patel et al. (2016) demonstrated the influence of density and viscosity gradients on turbulence using theanisotropy tensor, which is defined asConclusion
ai j =

u′i u′j
2k

−

δi j
,
3

We have performed well-resolved LES simulations of aturbulent straight duct flow with a rectangular cross section. A
physically consistent subgrid-scale turbulence model based on
the Adaptive Local Deconvolution Method (ALDM) for implicit LES was incorporated. The wall shear stress τw was arti-(1)with the turbulent kinetic energy k = u′l u′l and δi j the512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 8: Reynolds stress anisotropy invariant map evaluated at z/h = 0 a) and z/h = 0.9 b). The lines correspond to LES P(), LES T+ (), LES T- (states are defined by the Lumley triangle.

). The lines start at the modified wall Y1 indicated with dots and the limiting

ficially modified at one of the four walls in order to investigate
the effect of an asymmetric wall roughness variation on thesquare duct flow field. A LES simulation based on the DNS
study by Pirozzoli et al. (2018) was used to assess the influence of the modified wall shear stress. Four modifications havebeen considered, consisting of two with a constant decrease
and increase of τw and two with a linear decrease and increaseof τw . The analysis of the law of the wall showed a flatteningof the velocity profile while approaching the sidewalls. Basedon the cross-sectional velocity distribution we observed thatthe vortices at the modified wall are getting smaller, weaker
and are pushed towards the corners whereas the vortices at the
opposite wall are expanding into the lower duct half. Usinga linear profile intensifies the observations made at the modified wall but only a moderate strengthening or weakening
was present at the sidewalls. Furthermore, the modificationof τw has mainly led to a change in intensity of the secondary
flow velocity w whereas the peak position remain fixed. Forv both the intensity and peak position have been altered. The
anisotropy of the Reynolds stress tensor was investigated bymeans the anisotropy invariant map and the tensor components. A higher wall shear stress causes a higher anisotropy
which consequently led to an increase in the strength of thesecondary flow since the are induced by the anisotropy of the
Reynolds stress tensor.

Emory, M. & Iaccarino, G. 2014 Visualizing turbulenceanisotropy in the spatial domain with componentality contours. Center for Turbulence Research Annual Research
Briefs pp. 123–138.
Hickel, S., Adams, N. A. & Domaradzki, J. A. 2006 An adaptive local deconvolution method for implicit LES. Journalof Computational Physics 213, 413–436.
Hickel, S., Egerer, C. P. & Larsson, J. 2014 Subgrid-scalemodeling for implicit large eddy simulation of compressible flows and shock-turbulence interaction. Physics of Fluids 26.
Kaller, T., Pasquariello, V., Hickel, S. & Adams, N. A. 2019
Turbulent flow through a high aspect ratio cooling duct withasymmetric wall heating. Journal of Fluid Mechanics 860,
258–299.
Patel, A., Boersma, B. J. & Pecnik, R. 2016 The influence ofnear-wall density and viscosity gradients on turbulence inchannel flows. Journal of Fluid Mechanics 809, 793–820.
Pirozzoli, S., Modesti, D., Orlandi, P. & Grasso, F. 2018 Turbulence and secondary motions in square duct flow. Journalof Fluid Mechanics 840, 631–655.
Salinas Vázquez, M. & Métais, O. 2002 Large-eddy simulation of the turbulent flow through a heated square duct.
Journal of Fluid Mechanics 453, 201–238.
Schmidt, S. J., Sezal, I. H. & Schnerr, G. H 2006 Compressible simulation of high-speed hydrodynamics with phasechange. In European Conference on Computational Fluid
Dynamics. TU Delft, The Netherlands.
Thakkar, M., Busse, A. & Sandham, N. D. 2018 Direct numerical simulation of turbulent channel flow over a surrogate
for nikuradse-type roughness. Journal of Fluid Mechanics
837.
Vane, Z. P. & Lele, S. K. 2015 Prediction of turbulent secondary flows in ducts using equilibrium wall-modeled les.
In 53rd AIAA Aerospace Sciences Meeting, p. 1274.
Vidal, A., Vinuesa, R., Schlatter, P. & Nagib, H. M. 2017 Influence of corner geometry on the secondary flow in turbulent
square ducts. International Journal of Heat and Fluid Flow67, 69–78.
Vinuesa, Ricardo, Noorani, Azad, Lozano-Durán, Adrián,
Khoury, George K El, Schlatter, Philipp, Fischer, Paul F &Nagib, Hassan M 2014 Aspect ratio effects in turbulent ductflows studied through direct numerical simulation. Journalof Turbulence 15 (10), 677–706.
Zhang, H., Trias, F. X., Gorobets, A., Tan, Y. & Oliva, A. 2015
Direct numerical simulation of a fully developed turbulent
square duct flow up to reτ= 1200. International Journal of
Heat and Fluid Flow 54, 258–267.

Acknowledgments
The authors gratefully acknowledge the financial support provided by the German Research Foundation (Deutsche
Forschungsgemeinschaft-DFG) within the framework of theSonderforschungsbereich Transregio 40, SFB-TRR40 (Technological foundations for the design of thermally and mechanically highly loaded components of future space transportation
systems). Computational resources have been provided by the
Leibniz Supercomputing Centre Munich (LRZ).

REFERENCES
Brun, C., Boiarciuc, M. P., Haberkorn, M. & Comte, P. 2008
Large eddy simulation of compressible channel flow. Theoretical and Computational Fluid Dynamics 22, 189 – 212.
Demuren, A. O. & Rodi, W. 1984 Calculation of turbulencedriven secondary motion in non-circular ducts. Journal of
Fluid Mechanics 140, 189–222.
Egerer, C. P., Schmidt, S. J., Hickel, S. & Adams, N. A. 2016
Efficient implicit LES method for the simulation of turbulent cavitating flows. Journal of Computational Physics
316, 453–469.

6