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

ASSESSMENT OF THREE WALL MODELS IN A SPATIALLY DEVELOPING
THREE-DIMENSIONAL BOUNDARY LAYER IN A BENT SQUARE DUCT

Xiaohan Hu, Imran Hayat, & George Ilhwan Park
Department of Mechanical Engineering and Applied Mechanics
University of Pennsylvania
Philadelphia, Pennsylvania 19104, USA
gipark@seas.upenn.edu

ABSTRACT

number is moderately high, with Reθ ranging between 4100
and 8500. Upstream of the bend, the flow along the centerline is a canonical 2D zero pressure gradient (ZPG) flat-plateboundary layer. Mean flow three-dimensionality was generated in the bend region approximately between x′ = 1626 mm
and x′ = 2224 mm due to the cross-stream pressure gradientinduced by the bend. Downstream of the bend, the 3DTBL
gradually returned to a 2DTBL owing to the vanished spanwise pressure gradient. The experimental study focused on theboundary layer along the local centerline where the streamwisepressure gradient is small.

We investigate and compare the predictive capabilitiesof three common wall models in a pressure-driven threedimensional turbulent boundary layer (3DTBL) developing on
the floor of a bent square duct. The wall models employed area simple equilibrium stress (ODE) model, an integral nonequilibrium wall model, and a PDE nonequilibrium model basedon unsteady RANS. The mean flow statistics from WMLES
are compared with the experiment of Schwarz & Bradshaw(J. Fluid Mech. (1994), vol. 272, pp. 183–210). Although
the difference in the wall-stress magnitudes predicted by thethree wall models is negligible, the two nonequilibrium wallmodels are shown to produce a more accurate prediction of
the wall-stress direction due to the three-dimensionality incorporated in their formulation, with the PDE nonequilibrium
wall model being more accurate among the two. However,the LES solution away from the wall is agnostic to the typeof wall model used, resulting in nearly identical predictionsof the mean and turbulent statistics in the outer layer for all
the wall models. This is explained by the vorticity dynamics
and the inviscid skewing mechanism of generating the mean
three-dimensionality.

Fig. 1 shows the test section in the experiment, which
consisted of a square duct (D × D = 0.762m × 0.762m) witha total curved length of L = 3.748m. The computational domain is identical to the experimental test section. The present
study was conducted at two grid resolutions: a coarse meshwith 8 million control volumes and a fine mesh with 38 million control volumes. The local boundary layer contains approximately 16∼23 and 32∼45 cells across its thickness in thecoarse and fine computational meshes, respectively.

Flow configuration

Inflow characterization

The reference configuration for the present study is theexperimental setup of Schwarz & Bradshaw (1994). In thisexperiment, a spatially developing turbulent boundary layer
grows along the floor of a square duct with a 30◦ bend (Fig. 1).
The boundary layer on the floor was very thin compared tothe duct height, with δ99 /D ranging between 0.026 and 0.07
throughout the test section, where D is the width (or height)of the square duct. The flow was far from being fully developed. The primary region of interest in the experiment was
the centerline region of the duct and it was expected that thesecondary flow near the corner regions would have negligible
influence on the centerline.
As in the reference study of Schwarz & Bradshaw (1994),
the global Cartesian coordinate system is denoted by (x, y, z),whereas a curvilinear coordinate system aligned with the local duct centerline is denoted by (x′ , y′ , z′ ). x′ and z′ represent
the local streamwise and cross-stream directions respectively,
and y = y′ represents the wall-normal distance from the floorof the duct. In the experiment, the boundary layer on the floor
was tripped using a trip wire at the duct inlet located at x′ = 0.
Boundary layers on the other three walls of the duct were nottripped (Schwarz, private communication, 2019). Reynolds

The experiment reports flow statistics at the 22 locations
shown in Fig. 1 along the duct centerline, with the first measurement location being far downstream of the test section inlet (at x′ = 826 mm). Instead of tripping the boundary layer
with a wire as done in the experiment, we employ a synthetic
turbulence generator based on a digital filter approach (Klein
et al., 2003) for approximating the inflow boundary conditionat x′ = 0 mm. This approach requires iterative guesses on thelength of the development region (if any) to be appended upstream of the nominal trip location in the experiment (x′ = 0mm), and the state of the inflow to be prescribed at the new
inlet location. It should be noted that the goal here is to reproduce the 2DTBL upstream of the bend reasonably well, whichthen acts as the inflow for the 3DTBL within the bend, rather
than to exactly match the flow conditions at the test section inlet. After iterating on several inflow conditions, we found thatprescribing a flat-plate turbulent boundary layer at Reθ = 2560
(Schlatter et al., 2010) at the inlet (x′ = 0 mm) reproduces theboundary layer statistics well at the first measurement location(station 0: x′ = 826 mm).Fig. 2 shows reasonable agreement
between the simulation and the experiment in terms of the distributions of the boundary layer and momentum thicknesses.

1

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

Figure 1. A Schematic of the floor of the duct (reproduction from figure 1 of (Schwarz & Bradshaw, 1994)). The measurement
locations in the experiment are marked as numbers 0-21 along the duct centerline. Two coordinate systems are employed. (x, y, z) is
a fixed coordinate system with the origin located at the inlet. (x′ , y′ , z′ ) is a curvilinear coordinate system aligned with the local duct
centerline (measurements in mm).

10000

0.06

8000
0.04

6000
4000

0.02

2000
0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

Figure 2. Centerline distributions of (a) boundary layer thickness and (b) momentum thickness (coarse mesh). Symbols, experiment;
red dash-dotted line, equilibrium wall model; blue solid line, PDE nonequilibrium wall model; green dashed line, integral nonequilibrium wall model. Black vertical dashed lines denote the start and end of the bend region.

Flow solver and SGS / near-wall modelingthough never completely recovering to zero. Consistent withthe respective complexity of the three wall models, within thebend region where the 3D effects are most prominent, the
PDE NEQWM predicts the flow turning angle most accurately
among the three wall models, followed by integral NEQWM,
and then EQWM. Fig. 3(b) shows the near-wall flow direction
predicted by different wall models through the select surface
streamlines calculated from the mean wall shear-stress vector.
Although all three wall models predict the flow deviation fromthe local centerline, the magnitude of deviation is not predicted
evenly across the different wall models. This is because thetotal flow turning is an accumulative effect of the local flowchange depicted in Fig. 3(a), and the area under the curve in
Fig. 3(a) can be interpreted as an approximation of the nearwall total flow turning angle. Thus, the NEQWM with thehighest prediction of local flow angle in the bend shows thelargest deviation from the centerline downstream of the bend,
even though the local flow angle downstream of the bend isalmost equal among the three wall models.

The simulations were performed with CharLES, an
unstructured cell-centered finite-volume compressible LES
solver developed at Cascade Technologies, Inc. The solver
employs an explicit third-order Runge-Kutta (RK3) scheme
for time advancement and a second-order central scheme forspatial discretization. The Vreman model (Vreman, 2004) isused to close the SGS stress and heat flux. The three wall models considered in the present study are: an equilibrium stressmodel (EQWM) in the form of ordinary differential equations(ODE) (Kawai & Larsson (2012)), an integral nonequilibrium wall model (integral NEQWM) that solves the verticallyintegrated Navier-Stokes equations with assumed mean velocity profiles for the wall stress, (Yang et al. (2015)), and a PDE
nonequilibrium wall model (PDE NEQWM) that retains thecomplexity of the full Navier-Stokes equations (Park & Moin
(2014)).

RESULTS

Fig. 4 shows the mean-velocity magnitude profiles for thethree wall models, the no-slip simulation and the experiment atseveral locations along the centerline, including upstream of,within, and downstream of the bend. The no-slip LES with nowall-flux modeling, predicts the mean velocity poorly, highlighting the need for wall modeling for the present flow withthe coarse mesh resolution. Here, a higher momentum is imparted to the boundary layer as a consequence of the under-

To highlight the characteristics of the 3DTBL, we will
focus on the results characterizing the mean flow threedimensionality here. The variation of the surface flow direction relative to the freestream direction is shown in Fig. 3(a).
The flow turning angle is essentially zero in the upstream,
grows significantly within the bend, reaches a maximum near
the end of the bend, and decays downstream of the bend,

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

Figure 3. (a) Centerline distribution of surface flow turning-angles with respect to the freestream (γw is the wall shear stress direction,
γ∞ is the freestream direction). (b) Streamlines of wall shear stress. Squares: experiment; red line: equilibrium wall model; blue line:
PDE nonequilibrium wall model; green line: integral nonequilibrium wall model. Solid line: coarse grid resolution; dashed line: finegrid resolution

0.1
0.08
0.06
0.04
0.02
0

0

1

2

3

4

5

Figure 4. Profiles of the mean-velocity magnitude at 5 measurement locations (stations 4, 8, 12, 16, and 20, from left to right).
Station 4 is upstream of the bend; station 8 is within the bend; stations 12, 16 and 20 are downstream of the bend. Red dash-dotted
line, EQWM; blue solid line, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black circle,
experiment. Profiles are shifted along the abscissa by 1.

predicted wall shear force. A significant improvement is observed in the predicted mean-velocity profiles with wall modeling. As the predicted skin-friction coefficient magnitude isalmost identical from the three wall models (not shown here),
we observe very little difference in the mean velocity profiles
across the wall models.
To characterize the three-dimensionality of this flow, next
we look at the flow direction which is defined by the angle
between the mean velocity vector and the freestream velocity
vector. Fig. 5 shows the variation in flow direction along thewall-normal direction. The flow direction varies strongly, withthe strongest mean flow three-dimensionality observed at thewall, which becomes weaker away from the wall, as evident
from the diminishing crossflow away from the wall. A consistent difference of approximately 3 degrees is observed for all
three WMLES compared with the experiment, indicating thatthe difference in the wall-model outputs (the wall-shear force
direction observed in Fig. 3(b)) has very little impact on the
LES solutions away from the wall.
In this flow, the mean three-dimensionality in the outerlayer is created by the inviscid skewing mechanism, where thestreamwise vorticity is produced by the reorientation of the

spanwise vorticity. We can further analyze the characteristics of such flow from the perspective of the vorticity transport equation. Equation 1 is the transport equation for themean streamwise vorticity (Bradshaw, 1987), where the termson the right hand side of this equation represent the contributions from vortex stretching, vortex tilting, Reynolds stress,and viscous effect, in the shown order. The ”inviscid skewing”
mechanism generates the mean three-dimensionality by reorienting the mean spanwise vorticity as the streamwise vorticity.
If we assume that the ”inviscid skewing” mechanism is thedominant mechanism for the mean three-dimensionality, thethird term on the right hand side would contribute the most to
the generation of mean streamwise vorticity. Keeping only thedominant terms in this equation leads to equation 2, where thesecond equality in this equation is obtained by assuming the
− ∂W
= 0).
mean wall-normal vorticity to be zero (Ωy = ∂U
∂z
∂x
∂ Ωx
∂ Ωx
∂ Ωx
∂U
∂U
∂U
+V
+W
= Ωx
+ Ωy
+ Ωz
(1)
∂x
∂y
∂z
∂x
∂y
∂z
 2

∂
∂2
∂2
+
− 2 (−v′ w′ ) +(v′2 − w′2 ) + ν∇2 Ωx
2∂ y∂ z∂y∂z

U3

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

25
20
15
10
5
0
-5

0

0.01

0.02

0.03

0.04

0.05

0.06

Figure 5. Mean flow direction relative to the local freestream as a function of wall distance. Red dash-dotted line, EQWM; blue solidline, PDE NEQWM; green dashed line, integral NEQWM; magenta dotted line, no-slip LES; black solid line, experiment. Symbols areused to differentiate stations along the duct floor centerline only (Square, station 0; triangle, station 6; diamond, station 10). Lines withthe same symbols denote the results at the same stations.

Integrating equation 7 along the streamline and substituting thedefinitions of vorticity, we get,
U

∂U
∂W
∂ Ωx
= Ωz
= Ωz
∂x
∂z
∂x

(2)


Un
Us
.
= 2γe 1 −Ue
Ue

Reorganizing equation 2 by dividing UΩz gives equation 3.
Assuming that Ωz and U vary little in the streamwise direction,
we can re-express this equation as equation 4.
1 ∂W
1 ∂ Ωx=Ωz ∂ x
U ∂x

(3)

∂ Ωx
∂ W
=
∂ x Ωz∂x U(4)When the ”inviscid skewing” mechanism is the dominant generation mechanism of the mean streamwise vorticity (meanthree-dimensionality), Un plotted against Us under the local
freestream coordinates will form a straight line. Fig. 6 shows
the Johnston triangular plot for the current flow in bent duct
and a temporally developing shear-driven 3D channel flowfrom Lozano-Durán et al. (2020). It is observed that the mean
velocities in the outer layer from the duct flow satisfy the SWH
formula well, whereas they deviate from the SWH relation in
the shear-driven case. This indicates the ”inviscid skewing”
mechanism is the dominant contribution to the mean threedimensionality in the outer part of the boundary layer in theduct flow. While in the transient 3D channel flow, ”inviscid
skewing” does not exist, i.e. the vortex tilting term is zero,as expected. The slope in the SWH relation represents thefreestream turning angle with respect to the upstream flow, andthe freestream slope in the triangular plot (Fig. 6(a)) therefore
increases toward the downstream direction.
In addition to the LES solution, Fig. 7 also shows thewall-model solutions (i.e. the mean velocities solved within
the wall-models) in the triangular plot. It shows that different wall models have different capabilities in predicting theskewed mean-velocity profiles. The wall model solutions areplotted from the origin to the LES matching locations. The
EQWM, due to its unidirectional-flow assumption, cannot describe skewed mean-velocity profiles. It shows up as a straightline starting from the origin in the triangular plot. For PDE
NEQWM and integral NEQWM, the profiles appear as curvedlines in the triangular plot, which means the flow directionchanges with the wall distance. Thus, the two NEQWM areable to capture the skewed mean-velocity profiles. At station 8
where the three dimensionality is strongest, the PDE NEQWM
is able to express a richer wall-normal dependent skewness

Equation 4 is referred to as the SWH relation in the global coordinate system. Equation 2 in the local freestream coordinate
system is given by equation 5 Horlock & Lakshminarayana(1973):

Us

∂Us
∂Us
∂Us
∂ Ωs Us Ωn
+
= Ωs
+ Ωn
+ Ωb
,
∂s
R
∂s
∂n
∂b

(5)

here, (s, n, b) represents the streamwise, spanwise and wall
normal directions in the local freestream coordinate system,ssrespectively. Since Ωn = ∂Uand Ωb = −( ∂U+ URs ), we get,
∂b
∂n
∂
∂s



Ωs
Us


=−

2Ωn
Us R

(6)

Using the relation Rdθ = ds, equation 6 can be re-written as
equation 7, which is the SWH formula in the local freestream
coordinate system.
dΩs= −2dθ .
Ωn

(8)

(7)

4

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

0.2

1

0.2

0.8

0.8

0.6

0.1

0.6

0.1

0.4

0.4

0.2

0
0.5

0.75

1

0.2

0
0.5

0

0.75

1

0

Figure 6. Johnston triangular plot (a) WMLES (EQWM) of the bent square duct: square, station 0; circle, station 4; triangle, station 6;
cross, station 8; diamond, station 10; star, station 12. (b) DNS of the shear-driven 3DTBL from the transient channel flow at Reτ = 546
at t + = 192 (Lozano-Durán et al., 2020). Red straight line, the SWH formula Eq. (8). Color bar denotes the wall distance normalized
by the local boundary layer thickness.
1

0.2

1

0.2

0.8

0.8

0.6

0.1

0.6

0.1

0.4

0.4

0.2

0

0

0.25

0.5

0.75

1

0.2

0

0

0

0.25

0.5

0.75

1

0

Figure 7. Johnston triangular plot of the WM and LES solutions. The outer layer LES solution, which predicted almost identicallywith different wall models, is colored by the wall distance normalized by the local boundary layer thickness. (a) Crossflow developing
stage: square, station 0; circle, station 4; upward-pointing triangle, station 6; downward-pointing triangle, station 8; diamond, station
10; star, station 12. (b) Crossflow decaying stage: square, station 12; circle, station 14; upward-pointing triangle, station 16; downwardpointing triangle, station 18; diamond, station 20; star, station 21. Red solid straight lines from the bottom right corner are given by
the SWH formula Eq.(8). Red dash-dotted lines from the origin represent the wall-model solution: EQWM; blue solid line: PDE
NEQWM; green dashed line: integral NEQWM.

merical or large eddy simulations. J. Comput. Phys. 186,
652–665.
Lozano-Durán, A., Giometto, M. G., Park, G. I. & Moin, P.
2020 Non-equilibrium three-dimensional boundary layers
at moderate reynolds numbers. J. Fluid Mech. 883, A20.
Park, G. I. & Moin, P. 2014 An improved dynamic nonequilibrium wall-model for large eddy simulation. Phys.
Fluids 26, 015108.
Schlatter, P., Li, Q., Brethouwer, G., Johansson, A. V. & Henningson, D. S. 2010 Simulations of spatially evolving turbulent boundary layers up to reθ = 4300. Int. J. Heat. Fluid.
Fl. 31, 251–261.
Schwarz, W. R. & Bradshaw, P. 1994 Turbulence structuralchanges for a three-dimensional turbulent boundary layer
in a 30 ◦ bend. J. Fluid Mech. 272, 183–210.
Vreman, A. W. 2004 An eddy-viscosity subgrid-scale modelfor turbulent shear flow: Algebraic theory and applications.
Phys. Fluids 16, 3670.
Yang, X. I. A., Sadique, J., Mittal, R. & Meneveau, C. 2015
Integral wall model for large eddy simulations of wallbounded turbulent flows. Phys. Fluids 27, 025112.

than the integral NEQWM. During the crossflow developing
stage (Fig. 7(a)), the difference between the two NEQWM solutions and the EQWM solution gradually grows. Downstreamof the bend, where the crossflow starts to decay (Fig. 7(b)),
the triangular plots for the three wall model solutions progressively collapse onto each other, until they become almost identical toward the end of the duct, where the flow is essentially
unidirectional.

REFERENCES
Bradshaw, P. 1987 Turbulent secondary flows. Ann. Rev. Fluid
Mech. 19, 53–74.
Horlock, J. H. & Lakshminarayana, B. 1973 Secondary flows:theory, experiment, and application in turbomachinery aerodynamics. Ann. Rev. Fluid Mech. 5, 247–280.
Kawai, S. & Larsson, J. 2012 Wall-modeling in large eddy
simulation: Length scales, grid resolution, and accuracy.
Phys. Fluids 24, 015105.
Klein, M., Sadiki, A. & Janicka, J. 2003 A digital filter based
generation of inflow data for spatially developing direct nu-

5