of 20
PHYSICAL REVIEW FLUIDS
6
, 024602 (2021)
Direct numerical simulations of a statistically stationary streamwise periodic
boundary layer via the homogenized Navier-Stokes equations
Joseph Ruan
and Guillaume Blanquart
California Institute of Technology, Pasadena, California 91125, USA
(Received 4 February 2020; accepted 6 January 2021; published 5 February 2021)
We demonstrate a method for direct numerical simulations (DNS) of incompressible,
flat-plate, zero pressure gradient, turbulent boundary layers, without the use of auxiliary
simulations or fringe regions, in a streamwise periodic domain via the homogenized
Navier-Stokes equations. This approach is inspired by Spalart’s original (1987) method,
but improves upon his drawbacks while simplifying the implementation. Most simulations
of flat-plate boundary layers require long streamwise domains owing to the slow boundary
layer growth and inflow generation techniques. Instead, we use anticipated self-similarity
to solve the equations in a normalized coordinate system to allow for streamwise period-
icity, similar to Spalart’s original method. The resulting integral values, the skin friction
coefficient and shape factor,
H
12
and
C
f
, are within
±
1% and
±
3% of the empirical fits.
The mean profiles show good agreement with spatially developing DNS and experimental
results for a wide range of Reynolds numbers from Re
δ
∗
=
1460 to 5650. The method
manages to reduce computational costs by an estimated one to two orders of magnitude.
DOI:
10.1103/PhysRevFluids.6.024602
I. INTRODUCTION
Turbulent boundary layers have played a crucial role in numerous engineering applications since
Prandtl introduced the concept of the boundary layer in 1904 [
1
]. Wall-bounded flows contain a
wealth of open questions from the near-wall generation cycle [
2
] to the large scale structures [
3
]
and even the precise nature of the logarithmic layer [
4
]. All of these phenomena more fully manifest
themselves in higher Reynolds number flows. To thoroughly investigate them from a computational
perspective, there is an ever-present need for simulations to match or exceed experimental Reynolds
numbers.
At the moment, pipe and channel flow direct numerical simulations (DNS) have been able to
do so. Lee and Moser (Ref. [
4
]) conducted a channel DNS at a friction Reynolds number of 5186,
and, for comparison, experimental channel measurements by Schultz and Flack (Ref. [
5
]) were
at a friction Reynolds number of 6000. In contrast, boundary layer DNS have only been able to
reach Reynolds numbers of at most Re
δ
∗
=
13 000 [
6
], while modern boundary layer experiments
regularly achieve far higher Reynolds numbers, up to Re
δ
∗
≈
120 000 [
7
]. We define the Reynolds
number for boundary layers by Re
δ
∗
=
u
∞
δ
∗
/ν
, with the displacement thickness
δ
∗
, the free-stream
streamwise velocity
u
∞
, and the kinematic viscosity
ν
. As an example, to verify numerically the
Kármán constant of the logarithmic region, one would need a boundary layer DNS with a Reynolds
number at least one order of magnitude larger than the current state-of-the-art [
8
]. This status quo
has two origins: the need for inflow boundary conditions and the growing nature of the boundary
layer. Previous computational solutions can be classified into two main categories based on their
treatments of both issues.
Method No. 1: Inflow generation.
Recognizing that simulation of transition is computationally
expensive and complex [
9
], and long streamwise domains are required to reach high Reynolds
2469-990X/2021/6(2)/024602(20)
024602-1
©2021 American Physical Society
JOSEPH RUAN AND GUILLAUME BLANQUART
TABLE I. Summary of computational characteristics of previously published numerical approaches.
Statistically
Statistically
Closed
Method
homogeneous (in
x
)
stationary
form
Inflow generation
–
x
–
Temporal DNS
x
–
x
Spalart (1988)
x
x
–
Proposed method
x
x
x
numbers, multiple studies have used inflow generation techniques to bypass transition and initialize
their simulations with much higher Reynolds number inflows. There have been several popular
methods such as synthetic generation methods [
10
], strong recycling methods [
11
], and weak
recycling methods [
12
–
16
]. While these methods have been used to great effect, they still suffer
from multiple computational challenges.
First, current state-of-the-art methods still waste large portions of the computational domain as
a result of inflow generation techniques. The recycling domain alone can take up to 25% of the
domain [
12
]. Synthetic methods can overcome these recycling costs, but Sillero
et al.
[
16
]showed
that all inflow generation methods have an additional “eddy-turnover recovery length” over which
the simulation is heavily influenced by the inflow generation method. The authors found that for
over 25% of their so-called “production” domain, none of the simulation statistics match empirical
results. Overall, these inflow-generation methods underexploit available computational resources.
Second, as the boundary layer grows in space (
d
δ
∗
/
dx
>
0), momentum is displaced away from
the wall, and mass leaves through the top surface of the computational domain (
V
∞
>
0). To enforce
global mass conservation in incompressible DNS, an
apriori
streamwise dependence of
V
∞
must be
imposed over the entire domain [
12
,
15
,
16
]. Unfortunately, the imposition of
V
∞
enforces a particular
boundary layer growth rate. This is easily seen by integrating continuity (
V
∞
=
U
∞
d
δ
∗
dx
). Under these
conditions, one might ask if boundary layer simulations with fixed growth rates are true DNS.
Method No. 2: Streamwise periodicity.
One of the crucial differences between channel and
boundary layer flows from a numerical point of view is the streamwise growth of the boundary
layer. A statistically stationary, periodic method for boundary layer simulations would be ideal
since it would bypass the use of an inflow generation technique. This was proposed by Spalart
[
17
] who used a multiscale decomposition of velocity and subsequently scaled the wall-normal
coordinate to take advantage of the flow’s self-similarity. In the new coordinates, the mean quantities
were only functions of the rescaled wall-normal coordinate, and so Spalart obtained a streamwise
periodic and statistically stationary boundary layer. Unfortunately, this transformation generated
several additional terms in the governing equations that could not be closed without performing
auxiliary simulations of boundary layer flow for lower Reynolds numbers.
Eschewing these unclosed source terms, several studies [
18
–
21
] have nevertheless sought to
impose streamwise periodic boundary conditions, attempting to leverage the boundary layer’s
“quasiparallel” nature. While the numerical solution is now homogeneous in the streamwise direc-
tion (owing to the periodicity), it can never reach a statistically stationary state. These “temporal”
DNS (as opposed to the previously discussed “spatial” DNS) have a limited sample time (
<
1
δ
99
/
u
τ
)
before the boundary layer thickness drifts significantly.
In summary, there is currently no immediate computational framework for the simulation of
boundary layers that is (1) statistically homogeneous in the streamwise direction, (2) statistically
stationary, and (3) fully closed (see Table
I
). The objective of the present work is to propose
a framework to improve on Spalart’s preliminary method such that it does not rely on auxiliary
simulations for closure. The overall concept remains the same: to use periodic boundary conditions
and a scaling enforced by a coordinate transformation to keep the boundary layer statistically
stationary.
024602-2
DIRECT NUMERICAL SIMULATIONS OF A ...
We will detail the principles of the transformation in Sec.
II
and conduct extended domain simu-
lations in Sec.
III
to justify streamwise statistical homogeneity. We highlight the numerical methods
in Sec.
IV
. In Sec.
V
, we present validation and comparisons to DNS data from Refs. [
16
,
22
] and
empirical fits by Refs. [
23
,
24
]. Finally, in Sec.
VI
, we discuss computational savings.
II. ANALYSIS OF STATIONARY BOUNDARY LAYER
The goal of this section is to describe the proposed method of simulating flat-plate turbulent
boundary layers. It begins with a description of the spatial transformation, and a discussion of the
simplifications leading to the final set of equations.
A. Transformation of the Navier-Stokes equations
The flat-plate turbulent boundary layer is analyzed in the Cartesian coordinate system using index
notation such that the velocity components in the Cartesian streamwise (
x
1
), wall-normal (
x
2
), and
spanwise (
x
3
) directions are
u
1
,
u
2
, and
u
3
, respectively. With pressure and density as
P
and
ρ
,
respectively, the incompressible Navier-Stokes equations for mass and momentum conservation are
∂
u
j
∂
x
j
=
0
,
(1)
∂
u
i
∂
t
+
u
j
∂
u
i
∂
x
j
=−
1
ρ
∂
P
∂
x
i
+
ν
∂
2
u
i
∂
x
2
j
.
(2)
We now apply a coordinate transformation from
x
i
to
ξ
i
which rescales the wall-normal coordinate
by a streamwise varying
C
2
function
q
=
q
(
x
1
).
ξ
1
=
x
1
,ξ
2
=
q
0
q
x
2
,ξ
3
=
x
3
,
(3)
where
q
0
=
q
(
x
0
) is a normalization constant that is yet to be determined. Applying this coordinate
transformation directly to the Navier-Stokes equations yields the following set of equations for mass
and momentum conservation.
∂
u
j
∂ξ
j
=
ξ
2
q

q
∂
u
1
∂ξ
2
+
H
c
,
(4)
∂
u
i
∂
t
=−
u
j
∂
u
i
∂ξ
j
−
1
ρ
∂
P
∂ξ
i
+
ν
∂
2
u
i
∂ξ
2
j
+
ξ
2
q

q
u
1
∂
u
i
∂ξ
2
+
H
p
(
u
i
)
+
H
ν
(
u
i
)
,
(5)
where
H
c
=
(
1
−
q
0
q
)
∂
u
2
∂ξ
2
,
(6)
H
p
(
u
i
)
=
δ
1
i
1
ρ
q

q
ξ
2
∂
P
∂ξ
2
+
(
q
−
q
0
q
)(
δ
2
i
1
ρ
∂
P
∂ξ
2
+
u
2
∂
u
i
∂ξ
2
)
,
(7)
H
ν
(
u
i
)
=
ν
[
1
−
(
q
0
q
)
2
+
(
ξ
2
q

q
)
2
]
∂
2
u
i
∂ξ
2
2
+
ν
[
2
(
q

q
)
2
−
q

q
]
ξ
2
∂
u
i
∂ξ
2
−
2
νξ
2
q

q
∂
2
u
i
∂ξ
1
∂ξ
2
,
(8)
where
δ
ij
is the delta-Dirac function,
H
c
is an additional continuity term,
H
p
contains convective and
pressure additional metric terms, and
H
ν
contains the viscous metric terms. The equations as shown
024602-3
JOSEPH RUAN AND GUILLAUME BLANQUART
FIG. 1. Budgets of (a) streamwise momentum and (b) wall-normal momentum equations from DNS data
(Ref. [
16
]). Lines: (solid blue)
|
C
NS

ξ
3
,
t
|
; (solid magenta)
|
V
NS

ξ
3
,
t
|
; (solid green)

P
NS

ξ
3
,
t
; (dashed black)
|
Src

ξ
3
,
t
|
; (dashed green)
|
H
p

ξ
3
,
t
|
; (dashed magenta)
|
H
ν

ξ
3
,
t
|
.
are exact and equivalent to the original Navier-Stokes equations. At the moment,
q
(
x
1
) still requires
a closure equation for the transformed Navier-Stokes equations to be complete. There are several
possible choices to choose from such as the 99% boundary layer thickness
δ
99
, the displacement
thickness
δ
∗
, and the momentum thickness
θ
. Each choice yields a unique and mathematically valid
coordinate transformation.
B.
apriori
analysis
We perform a budget analysis of the streamwise and wall-normal momentum equations (
4
)–(
7
)
and the turbulent kinetic energy equations (
C1
)–(
C4
). This
apriori
analysis is performed using the
DNS data from Ref. [
16
] near Re
θ
0
=
4000.
Any
apriori
analysis of Eqs. (
4
)–(
7
) requires estimates for the function
q
(
x
1
). This function is
here approximated by
θ
(
x
1
), and justification for the estimate will be given in Sec.
II C
. Empirical
fits from [
23
] provide value for
θ

θ
θ
0
at Re
θ
0
=
4000.
We start with the streamwise momentum equation. First, we evaluate all terms at
q
=
q
0
.This
reduces
H
p
to a single term, removes a term from
H
ν
, and completely eliminates
H
c
.Wealso
group the main convective terms,
C
NS
=
u
1
∂
u
1
/∂ξ
1
+
u
2
∂
u
1
/∂ξ
2
+
u
3
∂
u
1
/∂ξ
3
, and the main vis-
cous terms,
V
NS
=
ν
(
∂
2
u
1
/∂ξ
2
1
+
∂
2
u
1
/∂ξ
2
2
+
∂
2
u
1
/∂ξ
2
3
). The source term for this equation is given
by Src
=
q

/
q
ξ
2
u
1
∂
u
1
/∂ξ
2
. Figure
1(a)
shows the budget analysis of the streamwise momentum
equation. All terms have been first averaged over time and spanwise coordinate, represented by
·
ξ
3
,
t
, and then the inner-scaled absolute values of these averages are plotted. Notably, the main
convective terms change sign near the wall to balance the main viscous terms.
As expected, the most dominant terms are the convective and viscous terms,
C
NS
and
V
NS
.
Furthermore, Fig.
1(a)
clearly shows that the convective metric term
ξ
2
u
1
(
q

/
q
)
∂
u
1
/∂ξ
2
is the
most dominant of the additional metric terms. It balances the main convective terms near the
end of the logarithmic region and throughout the wake region (80
<ξ
+
2
<
2000). In contrast,
the viscous metric term
H
ν
is over six orders of magnitude smaller than the streamwise convective
term throughout the entire boundary layer. Similarly, the
H
p
term is at least three orders of
magnitude smaller than the streamwise convective metric term until the very end of the wake region
near the free stream. From an
apriori
perspective, the neglecting of
H
ν
and
H
p
is justified.
One can also apply a similar analysis to the wall-normal momentum equation. In this case,
the convective terms are bundled as
C
NS
=
u
1
∂
u
2
/∂ξ
1
+
u
2
∂
u
2
/∂ξ
2
+
u
3
∂
u
2
/∂ξ
3
, and the main
024602-4
DIRECT NUMERICAL SIMULATIONS OF A ...
FIG. 2. Mean turbulent kinetic energy budget from DNS data (Ref. [
16
]). Lines: (solid blue) turbulent pro-
duction; (red) dissipation; (green) pressure diffusion; (black) turbulent diffusion; (cyan) advection; (magenta)
viscous diffusion; (dashed green)
H
p
contribution; (dashed magenta)
H
ν
contribution; (dashed blue) metric
source term contribution.
viscous terms are collected in
V
NS
=
ν
(
∂
2
u
2
/∂ξ
2
1
+
∂
2
u
2
/∂ξ
2
2
+
∂
2
u
2
/∂ξ
2
3
). The remaining terms
are the mean pressure gradient term
P
NS
=
1
/ρ∂
P
/∂ξ
2
and the source term Src
=
q

/
q
ξ
2
u
1
∂
u
2
/∂ξ
2
.
Figure
1(b)
shows the budget analysis of the wall-normal momentum equation. Again, all terms have
been first averaged over time and spanwise coordinate, and then the inner-scaled absolute values of
these averages are plotted.
In this case, the balance between the pressure and the convective terms dominates the entire
budget. The magnitude of the source term is between that of the convective and viscous terms. The
viscous metric term remains seven orders of magnitude smaller than the pressure and convective
terms throughout the boundary layer and thus can justifiably be neglected in the wall-normal
momentum equation. Near the free stream (
ξ
+
2
>
2000), the source term and the convective term
balance the pressure gradient term.
Finally, one can apply a similar analysis to the turbulent kinetic energy equation and track relative
contributions of the
H
p
,
H
ν
, and the source term. The results are shown in Fig.
2
. The metric source
term contribution balances the turbulent advection term which is at least three orders of magnitude
below the dominant budget terms. Throughout the boundary layer, the contributions of both
H
ν
and
H
p
to the kinetic energy budget remain several orders of magnitude lower than the dominant budget
terms. From an
apriori
perspective, their impact on turbulent intensities is negligible.
C. Simplified equations and closure
We can now make the following two critical assumptions:
(1) The governing equations evaluated at
q
(
x
1
)
=
q
0
are valid for a narrow streamwise domain
centered at
x
1
=
x
0
.
(2) There exists a function
q
(
x
1
) such that ensemble-averaged quantities are both statistically
stationary and statistically homogeneous in the
ξ
1
,ξ
3
directions.
Given both assumptions and with the neglecting of
H
ν
and
H
p
, the governing equations simplify
to
∂
u
j
∂ξ
j
=
ξ
2
q

0
q
0
∂
u
1
∂ξ
2
,
(9)
∂
u
i
∂
t
=−
u
j
∂
u
i
∂ξ
j
−
1
ρ
∂
P
∂ξ
i
+
ν
∂
2
u
i
∂ξ
2
j
+
ξ
2
q

0
q
0
u
1
∂
u
i
∂ξ
2
.
(10)
024602-5
JOSEPH RUAN AND GUILLAUME BLANQUART
These are the final governing equations to be solved via streamwise periodic simulation. Section
III
presents an
a posteriori
analysis justifying both the neglecting of
H
ν
and
H
p
and the streamwise
statistical homogeneity of Eqs. (
9
) and (
10
).
The use of both assumptions and the neglecting of
H
ν
and
H
p
mean that the governing equations
can be more accurately described as homogenized Navier-Stokes equations (HNSE). Consequently,
simulations utilizing this set of equations are still DNS but do not directly solve the NSE. The rest of
the document seeks to compare the solutions of the HNSE to experimental and numerical solutions
to the NSE.
We now seek to generate a closure equation for
q

0
/
q
0
by considering the
ξ
2
integrated continuity
and streamwise momentum equations in conservative form.
∫
∞
0
(
∂
u
1
∂ξ
1
+
∂
u
3
∂ξ
3
)
d
ξ
2
+
u
2
,
∞
=
q

0
q
0
∫
∞
0
(
ξ
2
∂
u
1
∂ξ
1
)
d
ξ
2
,
(11)
∫
∞
0
(
∂ρ
u
1
∂
t
+
∂ρ
u
1
u
j
∂ξ
j
+
∂
P
∂ξ
1
)
d
ξ
2
=
q

0
q
0
∫
∞
0
(
ξ
2
∂ρ
u
1
u
1
∂ξ
1
)
d
ξ
2
+
∫
∞
0
μ
∂
2
u
1
∂ξ
2
k
d
ξ
2
.
(12)
We now ensemble average Eqs. (
11
) and (
12
) and denote ensemble-averaged quantities by
·
.
Applying assumption 2 yields the following equations.
u
2
,
∞
=
q

0
q
0
∫
∞
0
(
u
1
,
∞
−
u
1

)
d
ξ
2
=
q

0
q
0
u
1
,
∞
δ
∗
,
(13)
u
1
,
∞
u
2
,
∞
+
τ
w
ρ
=
q

0
q
0
∫
∞
0
(
u
1
,
∞
2
−
u
1
u
1

)
d
ξ
2
,
(14)
where
τ
w
=
μ

∂
u
1
/∂ξ
2
|
ξ
2
=
0
is the wall shear stress. Some further manipulation gives our final
closure equation.
q

0
q
0
=
τ
w
/ρ
∫
∞
0
(
u
1
,
∞

u
1
−
u
1
u
1

)
d
ξ
2
.
(15)
This expression fully closes the simplified governing equations [Eqs. (
9
) and (
10
)] and is used in
the streamwise periodic numerical simulations of Sec.
IV
.
Before discussing the results, it is interesting to estimate
apriori
the right-hand side of the above
equation. The von Kármán momentum integral equation for a flat-plate boundary layer estimates
the growth rate of the momentum thickness by
θ

=
τ
w
ρ
u
2
1
,
∞
.
(16)
By definition,
u
2
1
,
∞
θ
=
∫
∞
0
(
u
1
,
∞

u
1
−
u
1

u
1

)
d
ξ
2
.
(17)
Assuming that

u
1

2
≈
u
2
1

, the von Kármán momentum integral equation becomes
θ

θ
≈
τ
w
/ρ
∫
∞
0
(
u
1
,
∞

u
1
−
u
1
u
1

)
d
ξ
2
=
q

0
q
0
.
(18)
024602-6
DIRECT NUMERICAL SIMULATIONS OF A ...
TABLE II. DNS parameters for the streamwise nonperiodic turbulent boundary layer simulation cases.
Dataset
L
x
Governing equations
Closure equation
Sample time
δ
99
/
u
τ
BL_Cart
15
δ
99
Eqs. (
1
)and(
2
)
None
15
BL_Full
15
δ
99
Eqs. (
4
)–(
8
)
Empirical value
30
BL_Simp
15
δ
99
Eqs. (
9
)and(
10
)
Empirical value
30
BL_Plus
15
δ
99
Eqs. (
4
)–(
8
)
1.1
×
empirical value
15
BL_Minus
15
δ
99
Eqs. (
9
)and(
10
)
0.9
×
empirical value
15
BL_Per
7
.
5
δ
99
Eqs. (
9
)and(
10
)Eq.(
15
)30
This completes the intuition that
q
(
x
) scales like
θ
(
x
) and furthermore provides an estimate for
q

0
/
q
0
.
q

0
q
0
≈
θ

θ
=
1
2
δ
∗
C
f
H
12
,
(19)
where
C
f
=
τ
w
/
(
1
2
ρ
u
2
1
,
∞
) is the skin-friction coefficient and
H
12
=
δ
∗
/θ
is the shape factor.
Note that the closure equation ensures a statistically stationary flow and consequently the solution
will be specific to a single Reynolds number. This is in direct contrast with recycling and rescaling
methods which solve for a range of Reynolds numbers but also use flow at high Reynolds number
stations as a substitute for a low Reynolds number inflow. The net effect of the closure equation
[Eq. (
15
)] is to allow the current method to avoid unphysical inflows by focusing on a single
Reynolds number.
III. SPATIALLY DEVELOPING SIMULATIONS
We conduct six sets of boundary layer simulations, each solving a different set of governing equa-
tions and boundary conditions, summarized in Table
II
. The results are used to justify assumptions
(1) and (2), and the simplifications made to the governing equations in Sec.
II B
.
A. Simulations and numerical methods
With the exception of BL_Per, all of the cases have streamwise nonperiodic boundaries in an
inflow
/
outflow setup. Case BL_Per corresponds to the most “modified” case: it solves Eqs. (
9
) and
(
10
) with Eq. (
15
) and implements streamwise periodic boundary conditions. Case BL_Simp also
solves Eqs. (
9
) and (
10
) but does not use streamwise periodic boundary conditions. Case BL_Full
solves Eqs. (
4
)–(
8
) and contains all of the previously neglected terms. Cases BL_Full and BL_Simp
use empirical relations for low Reynolds number [
24
] for the closure of
q

0
/
q
0
by approximating
q
≈
θ
. Cases BL_Plus and BL_Minus differ from BL_Simp by using a closure for
q
0
artificially
increased and decreased by 10%, respectively. Finally, case BL_Cart solves the regular Cartesian
Navier-Stokes equations [Eqs. (
1
) and (
2
)].
All of the cases have periodic spanwise directions and nonperiodic wall-normal directions. The
bottom of the domain is treated with a no-slip boundary condition, and the top of the computational
domain is treated with a Neumann boundary condition. Each of the five inflow
/
outflow cases use
planes from case BL_Per as an inflow (at
ξ
1
=
0). All of the streamwise nonperiodic cases use
convective outflow conditions at the streamwise outlet and have mass conservation conducted at
the streamwise outlet. In contrast, case BL_Per has mass conservation conducted at the wall-normal
outlet. The top of the computational domain requires vertical transpiration for all six cases. BL_Cart
imposes a transpiration velocity given by Ref. [
23
], similar to Ref. [
16
]. For the remaining cases,
Eq. (
11
) shows that any closure for
q

0
/
q
0
directly provides a value for
u
2
,
∞
.
All of the cases have the same spanwise length of
L
z
=
2
.
6
δ
99
and wall-normal height of
L
y
=
3
.
4
δ
99
. They all have the same spatial resolution:
x
+
=
9
,
y
+
min
=
0
.
3
,
z
+
=
6. The key
024602-7
JOSEPH RUAN AND GUILLAUME BLANQUART
FIG. 3. Streamwise variation of (a) displacement thickness (
δ
∗
), (b) momentum thickness (
θ
), (c) shape
factor (
H
12
), and (d) skin friction coefficient (
C
f
). Lines: (solid black) BL_Full; (solid green) BL_Simp; (solid
cyan) BL_Minus; (solid magenta) BL_Plus; (solid red) BL_Cart; (dashed black) BL_Per; (dashed blue) scaled
empirical fit [
24
] to match inlet skin friction coefficient.
difference between the streamwise periodic and streamwise nonperiodic cases is the streamwise
domain length. BL_Per has a domain length of
L
x
=
7
.
5
δ
99
, whereas the rest of the cases have
a domain length of
L
x
=
15
δ
99
. It is known from [
25
] that the flow recovers from this particular
inflow technique after
∼
4–5
δ
99
. Accounting for potential outflow effects of at most
∼
2
δ
99
,this
leaves about 8
δ
99
of uncontaminated statistics.
Each set of governing equations is solved using the computational solver NGA [
26
]. The numer-
ical code solves the conservative-variable formulation of the low-Mach Navier-Stokes equations
with staggered finite difference operators and uses a fractional step method to enforce continuity.
The code is run fully second order in space and time.
B. Results
Figures
3(a)
and
3(b)
present the normalized displacement thickness and momentum thickness
averaged in
ξ
3
and in time for the streamwise nonperiodic cases. Naturally, each of the streamwise
nonperiodic simulation cases is affected by the convective outflow condition, most clearly seen in
case BL_Cart. Each of the displacement thickness and momentum thickness plots deviate in slope at
about 1
δ
99
from the streamwise outlet. Because BL_Cart represents a spatially developing boundary
layer, the displacement thickness increases from its original inflow value. The present increase by
024602-8
DIRECT NUMERICAL SIMULATIONS OF A ...
about 20% is expected given that
d
δ
∗
/
dx
=
u
2
,
∞
/
u
1
,
∞
and the imposed value of
u
2
,
∞
/
u
1
,
∞
≈
3
×
10
−
3
. Cases BL_Full and BL_Simp are indistinguishable and show relatively constant values of
δ
∗
and
θ
with fluctuation magnitudes of
±
0
.
3% of the inflow nominal value. The thicknesses of
BL_Plus and BL_Minus show immediate departures from the nominal value, by approximately 2%
of the original inflow value.
References [
9
,
27
] underscore the need for consistency when comparing DNS profiles of bound-
ary layers. For the sake of comparison, integral and global quantities are computed as described by
Ref. [
9
]. Specifically, the shape factor,
H
12
, is evaluated as
H
12
=
∫
δ
99
0
(1
−
u
1
/
u
1
,
∞
)
d
ξ
2
∫
δ
99
0
(
u
1
/
u
1
,
∞
)(1
−
u
1
/
u
1
,
∞
)
d
ξ
2
,
(20)
where, for the remainder of this section, the overbar represents temporal and spanwise averaging.
Similarly, the wall shear stress is evaluated by
τ
w
=
μ∂
u
1
/∂ξ
2
|
ξ
2
=
0
. Figure
3(c)
presents the shape
factor for the nonstreamwise periodic cases. BL_Cart has a shape factor that monotonically drops by
2% from its inflow value, as expected from empirical fits by [
23
] with respect to Reynolds number.
BL_Plus and BL_Minus also exhibit a slowly varying shape factor, changing by approximately
±
0
.
3% from the inflow value. This is in contrast with BL_Full and BL_Simp, whose shape factors
are virtually identical and do not exhibit major mean variations.
Figure
3(d)
presents the skin friction coefficient, averaged in
ξ
3
and in time for the nonstreamwise
periodic cases. Case BL_Cart features a decreasing skin friction coefficient over the domain,
consistent with increasing Reynolds number. The skin friction coefficients of BL_Plus, BL_Minus,
BL_Full, and BL_Simp have fluctuations of about 1% of their expected mean value. Any variations
of the value with streamwise distance are masked by these fluctuations. Again, BL_Simp and
BL_Full are virtually indistinguishable.
Figure
4
shows temporal and spanwise-averaged profiles of
u
1
,
u
1
,
rms
,
u
2
,
rms
, and
−
u

1
u

2
from case
BL_Simp and BL_Full. These profiles are extracted from three streamwise locations: near the inlet
(
ξ
1
=
0) and outlet (
ξ
1
=
13
δ
99
) for case BL_Simp, and in the middle of the domain (
ξ
1
=
7
.
5
δ
99
)
for both cases BL_Simp and BL_Full. The mean streamwise velocity profiles are within
±
0
.
5% of
each other. The streamwise and wall-normal rms collapse within
±
1% of each other. The Reynolds
stress profiles show a strong collapse in both the inner and outer regions.
Overall, these streamwise nonperiodic simulations show that under a rescaling by
q
(
x
), the
resulting flow does not feature immediately observable streamwise inhomogeneities over a sizable
streamwise domain. The neglecting of
H
ν
and
H
p
terms and the use of streamwise periodic
conditions under Eqs. (
9
) and (
10
) are consequently well justified.
IV. NUMERICAL SETUP OF STREAMWISE PERIODIC SIMULATIONS
The present section outlines the simulations conducted in streamwise periodic domains. It
clarifies domain constraints and initial conditions, and describes additional numerical techniques
used during simulation.
A. Simulation cases
We now solve Eqs. (
9
) and (
10
) with streamwise periodic boundary conditions for four different
Reynolds numbers, summarized in Table
III
. Case BL1460 is equivalent to BL_Per. Cases BL2830,
BL3550, and BL5650 were chosen for direct comparison against the DNS and experiments of [
22
].
The domain size, (
L
x
,
L
y
,
L
z
), is determined primarily by the sizes of the largest turbulent struc-
tures. The pressure fluctuations are known to reach the furthest out of the boundary layer to about
2
.
4
δ
99
[
28
], setting the minimum requirement for wall-normal height. We set our domain height
to 18
δ
∗
∼
3
δ
99
to fully capture these fluctuations. Since low-momentum streaks are approximately
0
.
5
δ
99
in width [
15
,
29
], we opt for a spanwise width of 14
δ
∗
∼
2
.
5
δ
99
, which is comparable to the
domain size of Ref. [
22
]. The large-scale motions corresponding to bulges or hairpin packets have a
maximum streamwise length of 3
δ
99
[
30
–
34
]. In contrast, the very large-scale motions have lengths
024602-9
JOSEPH RUAN AND GUILLAUME BLANQUART
FIG. 4. Inner scaled (a) mean streamwise velocity
u
+
1
, (b) streamwise rms (
u
+
1
,
rms
), (c) wall-normal rms
(
u
+
2
,
rms
), and (d) Reynolds shear stress
−
u
+
1
u
+
2
, averaged over time and spanwise direction (
ξ
3
) for case
BL_Simp at different streamwise locations and BL_Full at the middle of the domain. Symbols: (green)
BL_Simp at
ξ
1
=
0 (inlet); (blue) BL_Simp at
ξ
1
=
7
.
5
δ
99
; (red) BL_Simp at
ξ
1
=
13
δ
99
;(
◦
)BL_Fullat
ξ
1
=
7
.
5
δ
99
.
of up to 10
δ
99
in the streamwise direction [
32
–
34
]. Lee and Sung (Ref. [
15
]) have found that these
structures have a mean streamwise length of less than 6
δ
99
and that statistically, over 95% of the
turbulent structures in their DNS had streamwise lengths of
<
6
δ
99
. And so, we opt for a domain of
40
δ
∗
∼
7
δ
99
in streamwise length.
The resolution is chosen so that the smallest turbulent structures can be adequately resolved.
The streamwise and spanwise grids are uniform with
x
+
=
9
,
z
+
=
6, which is comparable to
the resolution parameters of Sillero
et al.
[
16
](
x
+
≈
7,
z
+
≈
4
.
7) and Orlu and Schlatter [
22
]
(
x
+
≈
8
.
5,
z
+
≈
4). The wall-normal domain uses a hyperbolic stretching with eight points in
the viscous sublayer, (
ξ
2
<
5
δ
ν
), with
y
+
min
≈
0
.
3. This is comparable to the wall-normal resolu-
tion of Sillero
et al.
[
16
] who also had eight points in the viscous sublayer at the inlet and that of
Orlu and Schlatter [
22
] who had ten points in the viscous sublayer at their lowest Reynolds number.
To improve accuracy, we opt to use fourth-order finite difference spatial operators. Appendix A
compares the effect of both higher and lower finite difference spatial operators on case BL1460.
Cases BL2830, BL3550, and BL5650 are sampled over a period of 15
δ
99
/
u
τ
,
whereas case
BL1460 is sampled over a period of 30
δ
99
/
u
τ
. BL1460 was run for longer specifically to gather
temporal statistics of the global quantities, e.g., skin friction coefficient and shape factor.
024602-10