of 23
PHYSICAL REVIEW FLUIDS
4
, 084606 (2019)
Effective forcing for direct numerical simulations of the shear
layer of turbulent free shear flows
Chandru Dhandapani,
*
Kyupaeck Jeff Rah, and Guillaume Blanquart
California Institute of Technology, Pasadena, California 91125, USA
(Received 20 December 2017; published 28 August 2019)
A numerically efficient configuration to simulate turbulent flows is to use triply periodic
domains, with numerical forcing techniques to sustain turbulence. Previous homogeneous
shear turbulence simulations considered only idealized homogeneous shear flows and
not the statistically stationary shear turbulence observed in practical free shear flows. In
contrast, the current study mathematically derives the complete forcing technique from
the large scales of the turbulent free shear flows. Different statistically stationary free
shear flows are considered in this study, namely, a nearly homogeneous shear turbulent
flow, turbulent mixing layer, a turbulent planar jet, and a turbulent round jet. The
simulations are performed on triply periodic, statistically homogeneous cubic domains in
the vicinity of the shear layer in the self-similar region. An
apriori
analysis is performed
to calculate the effects of the different forcing terms and to predict the expected turbulence
quantities. The forcing technique is then used to perform direct numerical simulations at
different Reynolds numbers. Numerical results for the different cases are discussed and
compared with results from experiments and other simulations of free shear turbulent
flows. Anisotropy is observed both in the components of velocity and vorticity, with
stronger Reynolds number dependence in the anisotropy of vorticity. Energy spectra
obtained from the present homogeneous shear turbulence agree well with the spectra from
temporally evolving shear layers. The results also highlight the effects of the additional
forcing terms that were neglected in previous studies and the role of shear convection
and the associated splitting errors in the unbounded evolution of previous numerical
simulations.
DOI:
10.1103/PhysRevFluids.4.084606
I. INTRODUCTION
Turbulent free shear flows are found in a multitude of industrial applications and in nature, and
their analysis gives a lot of insight into turbulence and its structure. However, owing to the range
of scales and the stochastic and unsteady nature of turbulence, simulating such flows has proven
to be quite challenging. Various configurations have been used to simulate turbulent flows using
direct numerical simulations (DNS), which are resolved down to the smallest turbulent lengths
scales.
The most obvious configuration is to use the entire domain to solve the spatially evolving flow
[
1
3
]. In this configuration, the turbulence statistics reach a stationary state after a transient period,
and hence the results are ultimately independent from the initial conditions. Unfortunately, the
overall flow field depends strongly on the boundary conditions. These simulations typically include
a near-field region where the turbulence is not fully developed, and so this configuration is not
*
cdhandap@caltech.edu
2469-990X/2019/4(8)/084606(23)
084606-1
©2019 American Physical Society
DHANDAPANI, RAH, AND BLANQUART
computationally efficient. Since the entire flow field needs to be solved, these simulations are usually
performed at lower Reynolds numbers to reduce the computational costs.
Another configuration is to consider temporally evolving turbulent flows. A perfect example
is the mixing layer simulation by Rogers and Moser [
4
], which introduces homogeneity in the
streamwise direction. The extra periodic direction increases the computational efficiency, but
at the expense of physics. It also aids in calculating the energy spectra, which can be used to observe
the different scales of turbulence. Unfortunately, the statistics never reach a stationary state, and
hence the results still depend heavily on the initial conditions.
The third configuration is that of a triply periodic domain where the turbulence statistics are
homogeneous and there are no boundary conditions to implement because of the triple periodicity.
A good example is the numerical simulation of isotropic turbulence by Orszag and Patterson [
5
].
The computational efficiency is much higher as the flow is fully turbulent throughout. However,
in the absence of mean shear, the turbulent kinetic energy in the domain decays over time due to
viscous dissipation [
6
]. Hence, to keep the turbulent kinetic energy stationary over time, the missing
mean shear needs to be emulated through a method of forcing the turbulence. Turbulence in the
past has been forced by different techniques, including spectral forcing [
7
10
] and linear forcing
[
11
13
]. While these techniques are required to generate turbulence in the domain, these numerical
forcing techniques have been mostly arbitrary, and do not capture the physics of the large scale flows
accurately.
Recently, Rah
et al.
[
14
] combined the numerical tractability of the third configuration with
the physical accuracy of the first configuration. They used a triply periodic computational domain,
with the forcing calculated from the flow physics at a small region at the centerline of a turbulent
round jet, and forcing turbulence in a mathematically consistent way. The current study extends this
work by considering a small region in the self-similar shear layer of multiple statistically stationary
free-shear flows and using a triply periodic computational domain to simulate this shear-dominated
flow.
Several homogeneous shear turbulence (HST) simulations have been performed in the literature
using similar techniques [
15
21
]. While a shear production term was included in each study,
considering an idealized homogeneous shear flow, the forcing terms were not derived for practical
turbulent flows, and ultimately lacked the mathematical background to be compared to realistic
turbulent flows. Most of these simulations use shear periodic boundary conditions, but simulate
idealized homogeneous shear flow and the turbulence statistics are not stationary [
16
20
]. Some
simulations include a wall boundary in the cross-stream direction, and are not homogeneous [
21
].
By simulating only the velocity fluctuations, the simulations of the current study aim to obtain
the elusive combination of homogeneity and statistical stationarity in homogeneous shear turbulence
calculations.
In Sec.
II
, the forcing technique will be mathematically derived, and then calculated, from
locations in the shear layers of four different turbulent flows shown in Fig.
1
.An
apriori
analysis
will be performed for the forcing technique in Sec.
III
to observe the effects of the different forcing
terms. Section
IV
describes the simulations, and contains a discussion of the numerical results and
comparison with experiments and other simulations. Section
V
includes additional simulations
including linear diagonal terms, nonlinear terms, and mean advection terms. Section
VI
makes
concluding remarks about the observations from the study.
II. MATHEMATICAL DERIVATION
We start by reviewing Lundgren’s mathematical approach, which uses a Reynolds decomposition
to identify the effects of the large turbulent scales on the small scales. Then, four canonical flows are
considered (see Fig.
1
), and the forcing matrix is calculated for each of them. A forcing technique,
common to the self-similar shear layer of these four flows, is discussed after.
084606-2
EFFECTIVE FORCING FOR DIRECT NUMERICAL ...
FIG. 1. Different turbulent free shear flows considered for the current study with the computational domain
chosen (red cube): (a) nearly homogeneous shear turbulence (NHST), (b) mixing layer (ML), (c) planar jet (PJ),
(d) round jet (RJ).
A. Methodology: Review of Lundgren’s approach
First, we consider the Navier-Stokes (NS) equations for the velocity field
u
for a fluid flow with
constant density
ρ
, where
p
is the pressure and
ν
is the kinematic viscosity,
u
t
+
u
·
u
=−
1
ρ
p
+
ν
2
u
.
(1)
For any turbulent flow phenomenon, the instantaneous velocity field can be decomposed into mean
and fluctuating velocity fields (i.e., Reynolds decomposition),
u
=
u
+
u

, where
·
represents the
ensemble average. Transport equations for the fluctuations are obtained by calculating the difference
between the NS equations for the full velocity field and the transport equations for the mean velocity
field, namely,
NS
(
u
+
u

)
NS
(
u
+
u

)
.
(2)
This leads to
u

t
+
(
u
+
u

)
·
u

=−
1
ρ
p

+
ν
2
u

+
·
u

u

u

·
u
.
(3)
The extra terms in the transport equations for the fluctuating velocity, when compared with
Eq. (
1
), are the mean-flow advection term (
u
·
u

), the divergence of the Reynolds stress term
(
·
u

u

), and the production term (
u

·
u
). Lundgren focused on the production term as the
only contributor to turbulent kinetic energy production [
11
], but this is not the case as will be seen
later in section C. The major contribution to the turbulent kinetic energy comes from the production
084606-3
DHANDAPANI, RAH, AND BLANQUART
term, which is rewritten as a forcing term
A
·
u

,
u

t
+
u

·
u

=−
1
ρ
p

+
ν
2
u

+
A
·
u

,
(4)
where
A
is the forcing matrix, given by
A
=−
u
. The source term is linear in
u

, forces velocity
along all scales, and keeps the turbulent kinetic energy from decaying due to viscous dissipation.
Lundgren [
11
] further assumed that the forcing matrix,
A
is a diagonal matrix that generates
isotropic turbulence,
A
Lundgren
=
A
00
0
A
0
00
A
.
(5)
This isotropic forcing term was implemented as
A
u

, where
A
is an arbitrary forcing constant,
calculated based on the required turbulent Reynolds number [
13
]. In practice, the forcing matrix
depends on the gradients of the mean velocity.
B. Mean velocity gradients
Different free shear flows are considered in this study, namely a nearly homogeneous shear
turbulence (NHST) flow, a turbulent mixing layer (ML), a turbulent planar jet (PJ), and a turbulent
round jet (RJ). The mean velocity gradients can be calculated from the mean velocity profiles
obtained from experiments, for each free shear flow. Once again, the intent is to perform simulations
on triply periodic, statistically homogeneous cubic domains in the vicinity of the shear layer in their
respective self-similar region as shown in Fig.
1
.
1. Nearly homogeneous shear turbulence
For a homogeneous shear turbulence flow, the mean flow is in the streamwise direction (
x
). The
freestream velocity is constant along
x
and varies linearly in
y
, away from the walls located at
y
=−
h
/
2 and
y
=
h
/
2. The mean streamwise velocity at the center of the wind tunnel,
y
=
0is
U
C
.
Far downstream, the quantities are self-similar and are homogeneous in the
y
direction away from
the walls. However, the integral length scale

increases linearly with
x
[
22
], and consequently, the
Reynolds stresses and the velocity fluctuation magnitudes increase with
x
, hence the name “nearly”
homogeneous shear turbulence. Equation (
3
) for the HST flow in the center of the wind tunnel
becomes
u

t
+
u

·
u

=−
1
ρ
p

+
ν
2
u

u
x
y
u

y
e
x
u
x
u

x
+
u

x
u

x
x
e
x
+
u

x
u

y
x
e
y
.
(6)
Most simulations of homogeneous shear turbulence use periodic boundary conditions in the
x
direction without rescaling the velocity, and choose to neglect the divergence of the Reynolds
stress terms. The forcing matrix for NHST at
y
=
0is
A
NHST
=−
0
u
x
y
0
000
000
=
B
NHST
010
000
000
.
(7)
The only element of the forcing matrix is due to the shear strain rate
u
x
y
, and the matrix is normalized
by that quantity.
2. Mixing layer
For a spatial mixing layer, the mean flow is primarily in the streamwise direction (
x
). The
freestream velocity is constant along
x
, and is 0 for
y
→+∞
and
U
S
for
y
→−∞
. The center of the
084606-4
EFFECTIVE FORCING FOR DIRECT NUMERICAL ...
shear layer is at
y
1
/
2
, where the mean streamwise velocity is
U
S
/
2. Far downstream, the mixing layer
quantities are self-similar and are only a function of the similarity variable,
η
[
y
y
1
/
2
(
x
)]
(
x
).
The mixing layer thickness
δ
increases linearly with
x
, and
y
1
/
2
is linear in
x
[
23
]. There is no mean
flow in the spanwise direction (
u
z
=
0) and the flow is statistically homogeneous in the spanwise
direction (
u
x
z
=
0 and
u
y
z
=
0).
Equation (
3
) at the center of the mixing layer becomes,
u

t
+
u

·
u

=−
1
ρ
p

+
ν
2
u

u
x
x
u

x
e
x
u
x
y
u

y
e
x
u
y
x
u

x
e
y
u
y
y
u

y
e
y
u
x
u

x
u
y
u

y
+
u

x
u

x
x
e
x
+
u

x
u

y
x
e
y
.
(8)
The forcing matrix for the spatial mixing layer at
y
=
y
1
/
2
is
A
ML
=−
u
x
x
u
x
y
0
u
y
x
u
y
y
0
000
=
B
ML
0
.
035
1
0
0
.
001 0
.
035 0
000
,
(9)
calculated from the mean velocity profile given by Lumley [
24
]. The largest element of the forcing
matrix is due to the shear strain rate
u
x
y
, and the matrix is normalized by
B
ML
=−
u
x
y
(
y
1
/
2
)
=
1
.
022
U
S
δ
.
3. Planar jet
In a planar jet, the mean flow is primarily in the streamwise direction (
x
), and the centerline mean
velocity at the jet axis,
U
o
(
x
), decays along
x
as 1
/
x
[
1
,
25
27
]. The mean velocities are self-similar
far from the jet exit, and when normalized by the centerline velocity, are only functions of the
similarity variable,
η
y
/
y
1
/
2
(
x
), where
y
1
/
2
is the half-width of the jet defined by
u
x
[
x
,
y
1
/
2
(
x
)]
=
U
o
(
x
)
/
2.
The jet has no mean flow in the spanwise coordinate (
z
), and no mean gradients along
z
.The
forcing matrix for the planar jet in the middle of the shear layer at
y
=
y
1
/
2
is calculated from mean
velocity profiles given by Bradbury [
25
],
A
PJ
=−
u
x
x
u
x
y
0
u
y
x
u
y
y
0
000
=
B
PJ
0
.
071
1
0
0
.
007 0
.
071 0
000
,
(10)
where
B
PJ
=
0
.
730
U
o
y
1
/
2
. Once again, the largest contribution to the forcing matrix comes from the
off-diagonal shear strain term. The forcing matrix is comparable to the mixing layer forcing matrix
in Eq. (
9
).
4. Round jet
For a round jet, Eq. (
3
) is rewritten in cylindrical coordinates for simplicity. The mean flow
is primarily in the streamwise direction (
x
), and the mean centerline velocity
U
o
(
x
) has a 1
/
x
dependence [
2
,
28
31
]. We recall the flow is self-similar and the jet quantities, when normalized
by the centerline velocity, are only functions of the similarity variable
η
r
/
r
1
/
2
(
x
), where
r
1
/
2
is
the half-width of the jet.
There is no mean flow in the azimuthal direction (
θ
), and no mean gradients along
θ
. Hence, the
forcing matrix for the round jet in the middle of the shear layer at
r
=
r
1
/
2
as shown in Fig.
1(d)
is
084606-5
DHANDAPANI, RAH, AND BLANQUART
calculated from mean velocity profiles taken from Schlichting [
32
],
A
RJ
=−
u
x
x
u
x
r
0
u
r
x
u
r
r
0
00
u
r
r
=
B
RJ
0
.
014
1
0
0
.
001 0
.
037
0
00
0
.
023
,
(11)
where
B
RJ
=
0
.
586
U
o
r
o
1
/
2
. Once again, the largest element in the matrix is the off-diagonal shear strain
u
x
r
. The matrix is comparable to the velocity gradient matrix for planar jets from Eq. (
10
).
C. Additional source terms
In addition to the mean velocity gradients, there are source terms that arise from enforcing
periodic boundary conditions in the simulation domain [
14
]. The velocity fluctuations are appropri-
ately normalized to ensure that their second order statistics are homogeneous, so periodic boundary
conditions can be used, and these normalizations result in additional source terms. Although these
source terms can be calculated for any of the canonical flows mentioned before, the turbulent round
jet case is considered for the following calculations, as it has been researched in literature in greater
detail.
1. Periodicity in x
As mentioned earlier, in a round jet, the centerline velocity
U
o
(
x
) decreases with
x
as 1
/
x
. Since
the velocity fluctuations are proportional to the centerline velocity, they also decay along
x
as 1
/
x
.
Under these conditions, the flow is not statistically homogeneous in the
x
direction, and it would
be inappropriate to assume periodic boundaries. To lift this limitation, the velocity fluctuations are
rescaled by the 1
/
x
dependence as
u

x
=
u
(
x
)
x
x
o
x
,
u

y
=
u
(
x
)
y
x
o
x
,
u

z
=
u
(
x
)
z
x
o
x
,
(12)
where
u
(
x
)
is the velocity fluctuation that is statistically homogeneous along
x
in the vicinity of
x
=
x
o
. This rescaling produces extra elements in the forcing matrix from the
u
·
u

term. At
x
=
x
o
, the forcing matrix due to the periodicity correction in
x
is given by
A
x
=
u
x
x
o
00
0
u
x
x
o
0
00
u
x
x
o
.
(13)
2. Periodicity in r
The simulation assumes periodicity along
r
as well, but the velocity fluctuations depend on the
radial distance. To maintain statistical homogeneity along
r
, the velocity fluctuations are rescaled
by their individual
r
dependencies,
u
(
x
)
x
=
u
(
r
)
x
f
(
η
)
,
u
(
x
)
r
=
u
(
r
)
r
g
(
η
)
,
u
(
x
)
θ
=
u
(
r
)
θ
h
(
η
)
,
(14)
where
u
(
r
)
is the velocity fluctuation that is statistically homogeneous along
r
in the vicinity of
r
=
r
o
1
/
2
=
r
1
/
2
(
x
o
). The forcing matrix due to
u
·
u

applied on Eq. (
14
) is then
A
r
=
u
r
S
u
x
r
o
1
/
2
C
1
00
0
C
2
0
00
C
3
,
(15)
084606-6
EFFECTIVE FORCING FOR DIRECT NUMERICAL ...
where
S
=
dr
1
/
2
/
dx
is the spreading rate,
C
1
=−
df
d
η
(1)
,
C
2
=−
dg
d
η
(1), and
C
3
=−
dh
d
η
(1). From the
velocity fluctuations profiles from Hussein
et al.
[
31
], we have at
r
=
r
1
/
2
,
u
x
=
0
.
5
U
o
,
u
r
=
0
.
014
U
o
,
S
=
0
.
0935
,
C
1
=
0
.
517
,
C
2
=
0
.
398
,
C
3
=
0
.
345
.
(16)
3. Continuity
The original continuity equation for
u

is
u

x
x
+
1
r
(
ru

r
)
r
+
1
r
u

θ
∂θ
=
0
.
(17)
After the normalization in
x
and
r
for periodicity [Eqs. (
12
) and (
14
)], the continuity equation for
u
(
r
)
becomes
u
(
r
)
x
x
+
1
r
(
ru
(
r
)
r
)
r
+
1
r
u
(
r
)
θ
∂θ
=
(1
C
1
)
u
(
r
)
x
x
o
+
C
2
u
(
r
)
r
r
1
/
2
.
(18)
The continuity equation for
u
(
r
)
has two extra terms. While it is possible to solve the NS equations
with additional terms in the continuity equation, it is preferable to have no source terms. That is why
u
(
r
)
is rewritten in terms of
u

, under the conditions that
u

=
u
(
r
)
at
{
x
o
,
r
1
/
2
}
and
u

is divergence
free:
u
(
r
)
x
=
u

x
exp [(1
C
1
)(
x
/
x
o
1)]
,
u
(
r
)
r
=
u

r
exp

C
2
r
/
r
o
1
/
2
1
,
u
(
r
)
θ
=
u

θ
.
(19)
The forcing matrix due to the continuity correction is
A
C
=
u
x
x
o
(1
C
1
)00
0
S
u
x
u
r
r
o
1
/
2
C
2
0
000
.
(20)
The complete transformation from the original velocity fluctuation
u

to the statistically homo-
geneous, divergence-free velocity fluctuation
u

is given by
u

x
=
u

x
x
o
x
f
(
η
)exp[(1
C
1
)(
x
/
x
o
1)]
,
u

r
=
u

r
x
o
x
g
(
η
)exp

C
2
r
/
r
o
1
/
2
1
,
u

θ
=
u

θ
x
o
x
h
(
η
)
,
(21)
and the transport equation for
u

at
{
x
o
,
r
1
/
2
}
is calculated as
u

t
+
u

·
u

=−
1
ρ
p

+
ν
2
u

+
A
RJ
·
u

u
·
u

+
·
u

u

+
C
1
u
r
r
o
1
/
2
u

x
+
C
1
r
o
1
/
2
(
u

x
u

r
u

x
u

r
)
e
x
+

u
x
x
o
u

x
+
1
x
o
(
u

x
u

r
u

x
u

r
)

e
r
+

u
x
(1
C
3
)
x
o
u

x
+
1
C
3
x
o
(
u

x
u

θ
u

x
u

θ
)

e
θ
+
u
r
r
o
1
/
2
u

x
+
C
3
r
o
1
/
2
(
u

r
u

θ
u

r
u

θ
)
e
θ
+
visc
,
(22)
with the gradients of the normal stress in
·
u

u

being exactly zero, as
u

is homogeneous in
magnitude. Theoretically, the gradient of the Reynolds shear stress would still exist. However, at
084606-7
DHANDAPANI, RAH, AND BLANQUART
r
=
r
1
/
2
, the correlation coeffecient,
ρ
xr
=
u

x
u

r
/
(
u

x
u

x
u

r
u

r
)
1
/
2
is near constant [
33
], and its gradient
is near zero. The additional viscous terms are negligibly small for highly turbulent flows.
The final forcing matrix is calculated as a sum of all of the contributions from Eqs. (
11
), (
13
),
(
15
), and (
20
), and is given by
A
F
=
A
RJ
+
A
x
+
A
r
+
A
C

B
RJ
0
.
039
1
0
0
.
001 0
.
117
0
000
.
038
.
(23)
It is clear that the final forcing matrix is very close to the matrix from Eq. (
11
), with less than 6%
difference compared to the largest element. The periodicity in
x
and
r
, and the continuity correction
do not have significant contributions in the shear layer of a round jet, whereas it had significant
effects at the jet axis [
14
].
4. Nonlinear terms
All the source terms in Eq. (
23
) are linear in
u

; but the transformation from
u

to the statistically
homogeneous and divergence-free
u

in Eq. (
21
) gives rise to some nonlinear source terms owing
to the term
u

·
u

, as seen in Eq. (
22
). These nonlinear source terms can be written as
A
NL
·
u

A
NL
·
u

, where
A
NL
is given by
A
NL
=
C
1
r
o
1
/
2
u

r
00
0
1
C
2
x
0
u

x
0
00
1
C
3
x
0
u

x
+
C
3
r
o
1
/
2
u

r
.
(24)
These terms have similar magnitudes to the linear source terms from Eqs. (
13
) and (
15
), as
u

x
2
/
u
x

0
.
48 and
u

r
2
/
u
x

0
.
36.
D. Summary
The simulation considers the forcing matrix calculated at
{
x
,
r
}={
x
o
,
r
o
1
/
2
,
0
}
, and hence the
r
-
θ
direction in the jet coordinates can be replaced by
y
and
z
in the cartesian coordinate system of
the DNS. The velocity solved for in the simulation correspond to values at the half-width of the jet,
{
u

x
,
u

r
,
u

θ
}
(
x
o
,
r
o
1
/
2
)
={
u

x
,
u

y
,
u

z
}
. For simplicity,
u

would be represented as
u

henceforth.
Some key aspects of this derived forcing must be emphasized. First, the forcing is a direct result
of the physics of the free shear turbulent flows considered; the forcing term is not arbitrary, and is
derived mathematically from the large scales of the mean flow. Second, the forcing is not isotropic,
which is consistent with results from experiments of free shear flows, where
u

2
x
>
u

2
y
[
25
,
31
].
Third, the forcing in this case is not purely from the diagonal terms as suggested by Lundgren’s
isotropic turbulence, but rather dominated by an off-diagonal shear term.
Comparing with other Homogeneous Shear Turbulence (HST) simulations, where the only
production term is
Bu

y
ˆ
e
x
, there are additional linear forcing terms on the diagonal due to mean
velocity gradients, renormalizations to maintain periodicity in the
x
and
y
/
r
directions, and
continuity corrections. In addition to the linear diagonal forcing terms, there are also additional
forcing terms that are nonlinear in
u

. Finally, the mean advection term is calculated as
u
·
u

=
By
u

x
, which has been included in past simulations. To avoid confusion with the shear strain (i.e.,
energy production) term, this term is referred to as shear convection.
III. A PRIORI ANALYSIS
Multiple source terms have been computed in the previous section. Their effect on the turbulence
quantities can be estimated using an
apriori
analysis. Once the most dominant source terms have
084606-8
EFFECTIVE FORCING FOR DIRECT NUMERICAL ...
been selected, the relationship between the source terms and other turbulence quantities can be
established.
A. Contribution of source terms
The effect of all the source terms on the turbulence can be observed from the effects on the
turbulent kinetic energy,
k
=
1
2
u

2
x
+
u

2
y
+
u

2
z
(
·
represents ensemble average). The transport
equation for the turbulent kinetic energy can be obtained from the velocity fluctuations transport
equation as
dk
dt
=

u

i
u

i
t

.
(25)
The turbulent kinetic energy equation for the simulation including all the additional linear and
nonlinear terms and mean advection terms is given by
dk
dt
=−
ε
+
P
+
P
diag
+
P
NL
+
P
conv
,
(26)
where
ε
=
2
ν
s
ij
s
ij
is the energy dissipation rate. All other terms vanish under statistical ho-
mogeneity. The contribution by each of the terms to turbulent kinetic energy production can be
calculated and compared with the most dominant shear term
P
=
B
u

x
u

y
. The contribution from
the diagonal terms is calculated as
P
diag
P
=
0
.
039
B
u

x
u

x
+
0
.
089
B
u

y
u

y
+
0
.
038
B
u

z
u

z
B
u

x
u

y
.
(27)
Using Reynolds stress values from the round jet results from Hussein
et al.
[
31
],
P
diag
P
=
0
.
117. The
contribution from the nonlinear terms can also be calculated as
P
NL
P
=
0
.
82
B
u

x
u

x
u

y
+
0
.
10
B
u

x
u

y
u

y
B
u

x
u

y
U
o
+
0
.
10
B
u

x
u

z
u

z
+
0
.
59
B
u

y
u

z
u

z
B
u

x
u

y
U
o
.
(28)
Using velocity triple correlation values from the round jet results from Hussein
et al.
[
31
],
P
NL
P
=
0
.
209. The contribution from the shear convection term is computed as
P
conv
P
=
B

y
u

x
x
u

x
+
y
u

y
x
u

y
+
y
u

z
x
u

z

B
u

x
u

y
=
By
k
x
B
u

x
u

y
.
(29)
Because of statistical homogeneity in the
x
direction,
P
conv
P
=
0. In other words, the advection by the
mean term does not contribute to kinetic energy production, as mentioned earlier. Hence, the shear
convection terms are not included in the current simulation. Further analyses and justifications are
provided in Secs.
III C
and
VB
.
In summary, the shear strain is the most dominant term, contributing to 75% of the production of
turbulent kinetic energy. The linear terms in the diagonal of the forcing matrix and the nonlinear
terms contribute to 9% and 16% of the production, respectively. Similar results are obtained
for mixing layers and planar jets. The off-diagonal shear strain element is at least one order of
magnitude larger than the other elements in the matrix and is the major driving force for turbulence
production in these aptly named free shear flows, accounting for at least 75% of the turbulent kinetic
energy production.
In the current study, for a triply periodic simulation of HST, it is a good approximation to use the
off-diagonal shear strain, B, as the only forcing term, with the forcing matrix given by
A
HST
=
0
B
0
000
000
,
(30)
084606-9
DHANDAPANI, RAH, AND BLANQUART
where
B
can be chosen based on the simulation parameters and the desired turbulent Reynolds
number. While this term does not inject any external energy, it represents the injection of energy
into the velocity fluctuations by the mean flow, hence it is an “effective forcing term” in the spirit of
Lundgren’s approach, and is henceforth referred to as a forcing term for simplicity. This is similar to
conventional simulations of HST, where the off-diagonal shear strain term is the only mechanism for
turbulence production [
18
,
20
]. Those studies do not include any of the linear diagonal and nonlinear
forcing terms, but they include the shear convection term that does not contribute to turbulent kinetic
energy.
B. Stationary-state analysis
The entire mathematical framework presented in the previous section relies on the assumption
that the velocity field can be decomposed into mean and fluctuating quantities. Then, the simulations
in the current study solve for the velocity fluctuations. By construction, these velocity fluctuations
represent the fluctuations of the flow field in a small region of a
statistically stationary
turbulent
flow. Hence, the fluctuating quantities and their related statistics
must
reach a statistically stationary
state. This applies to turbulent kinetic energy, dissipation rate, Reynolds stress, and so on.
Before performing the HST simulations, the target Reynolds number of the simulation needs
to be decided, so that the required grid resolution can be evaluated to fully resolve down to the
smallest turbulent scales. The relationship between the forcing constant
B
and the Reynolds number
needs to be established, to calculate the required shear strain,
B
. The expected eddy turnover time
is also calculated from the turbulent kinetic energy and the energy dissipation rate, to determine the
total simulation time. These expected turbulence quantities are estimated from the stationary state
of these simulations.
The turbulent kinetic energy equation for this HST forcing, assuming spatial homogeneity, is
dk
dt
=−
ε
+
B
u

x
u

y
.
(31)
At statistically stationary state, the energy dissipation rate is
ε
=
B
u

x
u

y
.
(32)
This should be compared to the stationary state with the isotropic forcing [
13
],
ε
=
2
Ak
.
(33)
The cross correlation in Eq. (
32
) can be written in terms of the turbulent kinetic energy,
u

x
u

y
=
β
k
, where
β
is a nondimensional parameter.
The integral length scale,

, is defined as

=
u
3
rms
ε
=
2
u
rms
3
β
B
,
(34)
where rms is root-mean square and
u
rms
=

2
k
3
=
3
2
β
B
.
(35)
The Taylor microscale,
λ
, is calculated as
λ
=

15
ν
ε
u
rms
.
(36)
The expected Taylor microscale Reynolds number for HST is calculated as
Re
o
λ
=

45
β
B

2
2
ν
.
(37)
084606-10
EFFECTIVE FORCING FOR DIRECT NUMERICAL ...
FIG. 2. (a) Integral length scale normalized by the domain width, (b) Reynolds shear stress
u

x
u

y
normalized by turbulent kinetic energy, for the four DNS, and (c) comparison of Reynolds number dependence
of Reynolds shear stress
u

x
u

y
with other studies. Dashed lines corresponds to the averaged value obtained
from all simulations in the current study.
For isotropic turbulence simulations, it was given by Carroll and Blanquart [
13
],
Re
o
λ
=

45
A

2
ν


9
AL
2
5
ν
,
(38)
as
/
L

0
.
2 for isotropic turbulence in a triply periodic box domain [
12
,
13
], where L is the domain
width. As will be shown from numerical results in Figs.
2(a)
and
2(b)
,
β

0
.
4 and
/
L

0
.
28 for
HST. So, given the same domain width and viscosity, DNS of HST can be performed with the
same Reynolds number as DNS of homogeneous isotropic turbulence, using the forcing constant
B

3
.
2
A
.
The expected values for turbulent kinetic energy,
k
o
, and energy dissipation rate,
ε
o
, can be
calculated as
k
o
=
3
2
u
2
rms
=
27
8
β
2
B
2

2
(39)
and
ε
o
=
u
3
rms

=
27
8
β
3
B
3

2
.
(40)
The expected eddy turnover time
τ
o
is given by
τ
o
=
k
o
ε
o
=
1
β
B

25
32
A
,
(41)
which is slightly higher than for the isotropic case, where
τ
o
=
1
2
A
[
13
].
C. Shear convection
The proposed HST simulation has a key difference from most simulations of shear turbulence
[
17
20
]; it does not include the shear convection term
By
u

x
. The shear convection term requires
either a remeshing scheme after every few iterations [
34
] or implementing shear periodicity along
the
y
direction to avoid boundary discontinuities [
18
20
]. It is often accomplished by using operator
splitting [
16
,
20
], which may introduce further errors in the computational solution.
As mentioned earlier, the shear convection term does not contribute to turbulent kinetic energy
production (see Sec.
III A
), as
P
conv
=
By
u

x
·
u

=
0 due to spatial homogeneity. That being said,
it may still impact the turbulent flow. To quantify this impact, we evaluate the shear strain produced
by the advection term and compare it to the existing shear strain due to the turbulence. The shear
strain due to the convection term can be calculated as
u
x
y
=
B
and compared against the existing
084606-11
DHANDAPANI, RAH, AND BLANQUART
TABLE I. Simulation parameters of the different cases.
No Re
o
λ
N
3
L
ν
B
Forcing matrix Re
λ
u

x
/
u
rms
u

y
/
u
rms
u

z
/
u
rms
u

x
u

y
/
k
13664
3
2
π
0.159
7.33
A
HST
32
±
6 1.26 0.90 0.77
0.42
2 54 128
3
2
π
0.159
16.5
A
HST
52
±
9 1.23 0.93 0.78
0.41
3 80 192
3
2
π
0.159
37.1
A
HST
80
±
13 1.23 0.92 0.78
0.39
3a 80 192
3
2
π
0.159
37.1
A
F
80
±
15 1.21 0.94 0.79
0.39
3b 80 192
3
2
π
0.159
37.1
A
F
+
A
NL
80
±
13 1.20 0.95 0.80
0.38
3c 80 192
3
2
π
(1) 0
.
159
(2)
0
.
0159
37.1
A
HST
121
±
20 1.23 0.93 0.79
0.37
3d 80 192
3
2
π
0.1431
37.1
A
HST
85
±
10 1.23 0.94 0.78
0.39
4 128 384
3
0.126
1.5e-5
2.77
A
HST
135
±
23 1.22 0.93 0.80
0.38
shear strain due to the turbulence,
u

x
y
. Since
u

x
y
=
0, the second-order statistics are compared as
B
2

u

x
y
2

=
15
ν
2
B
2
15
ν
2

u

x
y
2


15
ν
B
2
2
ε
=
50
β
2
Re
2
λ
,
(42)
with the isotropic assumption that
ε

15
ν
2
(
u

x
y
)
2
.ForRe
λ
=
100, the ratio is 0.031. Hence, the
impact of the shear convection term is small, and decreases with increasing Reynolds number.
Thus, the shear convection term is omitted for true spatial homogeneity and numerical efficiency.
Its impact will be discussed in Sec.
VB
.
IV. NUMERICAL RESULTS
A. Simulation
Direct numerical simulations of homogeneous shear turbulence are performed in a triply periodic
box domain that is statistically homogeneous in all three directions. Simulations are performed with
a domain width of
L
=
2
π
, and various Reynolds numbers Re
λ
.
The simulations are performed using NGA [
35
], a semi-implicit velocity solver with an energy-
conserving finite difference scheme on a standard staggered grid. The code solves the Navier-Stokes
equations with the derived source term from Eq. (
30
) for constant density, temperature, and
viscosity.
The initial velocity fields are generated randomly, using the method suggested by Eswaran and
Pope [
7
]. These velocity fields conform to a specified Passot-Pouquet energy spectrum [
36
] and are
divergence free, as is required for constant density flows.
Multiple simulations are performed at different expected values of Re
λ
. The simulation pa-
rameters for the four different cases are tabulated in Table
I
. Cases 1 and 2 were performed to
investigate low Reynolds number effects, if any. Cases 3 and 4 are chosen so the Reynolds number
is comparable to simulations and experiments, published in literature (see Table
II
for a full list of
experimental and full domain DNS studies). More precisely, case 3 has a similar Reynolds number
to cases 1 [
4
], 4 [
1
], and 8 [
2
] from Table
II
; case 4 has a Reynolds number close to cases 5 [
26
] and
9[
3
] in Table
II
.
The simulations were performed for a total of 50 eddy turnover times, during which the
simulations were stable. The average values for the numerical results were calculated in the range,
10
τ
o
to 50
τ
o
.
084606-12