PHYSICAL REVIEW FLUIDS
4
, 073905 (2019)
Predicting the response of turbulent channel flow to varying-phase
opposition control: Resolvent analysis as a tool for flow control design
Simon S. Toedtli,
1
,
*
Mitul Luhar,
2
and Beverley J. McKeon
1
1
Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, California 91125, USA
2
Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles,
California 90089, USA
(Received 3 November 2018; published 31 July 2019)
The present study evaluates the capabilities of a low-order flow model based on the
resolvent analysis of McKeon and Sharma [B. J. McKeon and A. S. Sharma,
J. Fluid Mech.
658
,
336
(
2010
)] for the purpose of controller design for drag reduction in wall-bounded
turbulent flows. To this end, we first show that the model is able to approximate the change
in mean wall shear stress, which is commonly used as measure for drag reduction. We also
derive an analytical expression that decomposes the drag reduction in internal flows into
terms that can be predicted directly by the model and terms that allow for quantification
of model error if high-fidelity data are available. We then show by example of varying-
phase opposition control in a low-Reynolds-number turbulent channel flow that the drag
reduction predicted by the resolvent model captures the trend observed in direct numerical
simulation (DNS) over a wide range of controller parameters. The DNS results confirm
the resolvent model prediction that the attainable drag reduction strongly depends on the
relative phase between sensor measurement and actuator response, which raises interesting
flow physics questions for future studies. The good agreement between the resolvent model
and DNS further reveals that resolvent analysis, which at its heart is a linear technique, is
able to approximate the response of the full nonlinear system to control. We also show
that in order to make accurate predictions the model only needs to resolve a small subset
of the DNS wave numbers and that the controlled resolvent modes obey the Reynolds-
number scaling laws of the uncontrolled resolvent operator derived by Moarref
et al.
[R.
Moarref
et al.
,
J. Fluid Mech.
734
,
275
(
2013
)]. As a consequence, our results suggest
that resolvent analysis can provide a suitable flow model to design feedback flow control
schemes for the purpose of drag reduction in incompressible wall-bounded turbulent flows
even at technologically relevant Reynolds numbers.
DOI:
10.1103/PhysRevFluids.4.073905
I. INTRODUCTION
Turbulence may be aesthetically pleasing to the observer’s eye, yet it is undesired in most
engineering applications. For example, turbulent drag accounts for more than 50% of the total drag
force exerted on a large oil tanker [
1
] and thus a reduction of the turbulence level by just a few
percent could lead to tremendous energetic and monetary savings. The potential benefits associated
with targeted manipulations of turbulent flows have spurred enormous interest in the field of flow
control for well over a century and while important advances have been made in various aspects,
equally many questions remain open even at the very fundamental level (see, for example, Ref. [
2
]).
*
stoedtli@caltech.edu
2469-990X/2019/4(7)/073905(23)
073905-1
©2019 American Physical Society
TOEDTLI, LUHAR, AND MCKEON
One important open question of the field is the formulation of tractable and effective flow models
to enable systematic controller design for various control objectives. The present study aims at
contributing towards this question and assesses the potential of a flow model based on the resolvent
analysis of McKeon and Sharma [
3
] for the purpose of controller design. We limit the scope of this
study to the class of incompressible wall-bounded turbulent flows, to the control objective of drag
reduction, and to the mechanism of active feedback flow control. The resolvent model may also be
useful to design controllers for other types of flows, objectives, and control mechanisms, but the
assessment of such applications is beyond the scope of this study.
A. Resolvent analysis
Research efforts over the past 25 years have shown that linear mechanisms play a key role
in the transition to turbulence (among others, Ref. [
4
]) and in the dynamics of fully developed
wall-bounded turbulent flows (for example, Ref. [
5
]). One line of work devoted to these linear
mechanisms in wall-bounded turbulence considers the amplification characteristics of the forced
linear Navier-Stokes operator. Different analysis techniques and types of forcing have been proposed
in the past (see Ref. [
6
] for an overview) and the resemblance of the associated linear responses to
different features of turbulent flows has underscored the importance of linear mechanisms and the
robustness of the Navier-Stokes system to the type of forcing. The foundation for this study is
one particular analysis of the forced linear Navier-Stokes system, namely, the resolvent analysis
of McKeon and Sharma [
3
]. Under the resolvent formulation, the full nonlinear Navier-Stokes
equations (NSE) are rewritten as a linear operator that maps the nonlinear advection term to velocity
and pressure fluctuations. Unique to this approach is that the nonlinear term is viewed as intrinsic
forcing to the remaining linear system and that an efficient representation of the linear operator is
obtained from identifying the directions which are most amplified by the operator in an energy norm
sense. The resulting decomposition of the operator lends itself to low-order modeling and will be
at the heart of the flow model considered here. A clear understanding of the method and results
of this study requires basic familiarity with resolvent analysis and we therefore briefly review its
key elements below by example of a fully developed turbulent channel flow with periodic boundary
conditions in the wall-parallel directions.
Following McKeon and Sharma [
3
], we first Reynolds decompose all flow quantities into a spa-
tiotemporal mean
(
·
) and turbulent fluctuations (
·
)
. For example,
u
(
x
,
y
,
z
,
t
)
=
u
(
y
)
+
u
(
x
,
y
,
z
,
t
),
where
x
,
y
, and
z
denote the streamwise, wall-normal, and spanwise coordinates, respectively,
with corresponding components
u
,
v
, and
w
of the velocity vector
u
, and
t
represents time. The
streamwise mean velocity profile
u
(
y
) is assumed to be known (from an experiment, numerical
simulation, or eddy viscosity model), which allows us to restrict the problem to the turbulent
fluctuations.
The model development progresses by seeking an optimal basis in all three spatial dimensions
and time, where the notion of optimality is defined in an energy norm sense: In each dimension, we
seek that basis whose first
N
basis functions capture more kinetic energy
|
u
|
2
(under an
L
2
-norm)
of the flow than the first
N
functions of any other basis. It can be shown analytically that the Fourier
basis is the optimal one in the statistically homogeneous streamwise and spanwise directions and
in the stationary temporal coordinate [
7
]. This motivates expanding the fluctuating velocity and
pressure (
p
) fields in terms of Fourier modes
u
(
x
,
y
,
z
,
t
)
p
(
x
,
y
,
z
,
t
)
=
∞
l
=−∞
∞
m
=−∞
∞
−∞
ˆ
u
(
l
,
m
,ω,
y
)
ˆ
p
(
l
,
m
,ω,
y
)
e
i
[(2
π/
L
x
)
lx
+
(2
π/
L
z
)
mz
−
ω
t
]
d
ω,
(1)
with coefficients [
ˆ
u
,
ˆ
p
]
T
. Recall that the domain is periodic in
x
and
z
, with period
L
x
and
L
z
,
respectively. The periodicity restricts the streamwise (
k
x
) and spanwise (
k
z
) wave numbers to
integer multiples of the fundamental wave numbers
k
x
=
l
(2
π/
L
x
) and
k
z
=
m
(2
π/
L
z
), with
{
l
,
m
}∈
Z
, and reduces the Fourier transforms in these directions to Fourier series. The temporal
073905-2
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
coordinate is unbounded and therefore
ω
∈
R
. A single Fourier mode can be characterized either
by its streamwise and spanwise indices and temporal frequency [
l
,
m
,ω
]
T
or by its wave-number
vector
k
=
[
k
x
,
k
z
,ω
]
T
, and both notations will be used interchangeably when working in Fourier
domain.
The fluctuating part of the NSE can be written for each mode with
k
=
0
as
⎡
⎢
⎣
ˆ
u
(
k
,
y
)
ˆ
v
(
k
,
y
)
ˆ
w
(
k
,
y
)
ˆ
p
(
k
,
y
)
⎤
⎥
⎦
=
⎡
⎢
⎢
⎢
⎣
−
i
ω
+
ik
x
u
−
Re
τ
d
u
dy
0
ik
x
0
−
i
ω
+
ik
x
u
−
Re
τ
0
d
dy
00
−
i
ω
+
ik
x
u
−
Re
τ
ik
z
ik
x
d
dy
ik
z
0
⎤
⎥
⎥
⎥
⎦
−
1
H
(
k
)
×
⎡
⎢
⎢
⎣
ˆ
f
u
(
k
,
y
)
ˆ
f
v
(
k
,
y
)
ˆ
f
w
(
k
,
y
)
0
⎤
⎥
⎥
⎦
ˆ
f
(
k
,
y
)
,
(2)
where
=
d
2
/
dy
2
−
k
2
x
−
k
2
z
is the Laplacian,
{
ˆ
f
u
,
ˆ
f
v
,
ˆ
f
w
}
denote the Fourier transforms of the
nonlinear advection terms, and Re
τ
=
u
τ
h
/ν
is the Reynolds number based on the channel half-
height
h
, friction velocity
u
τ
=
√
τ
w
/ρ
, and kinematic viscosity of the fluid
ν
(
τ
w
is the mean wall
shear stress and
ρ
is the fluid density). As implied by Re
τ
, all quantities are made dimensionless
with
h
and
u
τ
and, unless stated otherwise, the same nondimensionalization will be employed in
the remainder of this study. The linear operator
H
(
k
) that maps the nonlinear forcing terms to the
velocity and pressure fluctuations at each
k
is called the resolvent operator. Recall that we assumed
the mean velocity profile to be known and therefore
H
(
k
) is fully determined.
The optimal basis in the remaining wall-normal direction is obtained from a singular-value
(Schmidt) decomposition of the resolvent operator, as suggested by McKeon and Sharma [
3
], and
interested readers are referred to Refs. [
8
,
9
] for an in-depth discussion of the technicalities. The
singular-value decomposition returns an ordered basis pair
{
ˆ
ψ
j
,
ˆ
φ
j
}
which can be used to rewrite
Eq. (
2
)as
ˆ
u
(
k
,
y
)
ˆ
p
(
k
,
y
)
=
H
(
k
)
ˆ
f
(
k
,
y
)
=
∞
j
=
1
σ
j
(
k
)
ˆ
ψ
j
(
k
,
y
)
ˆ
φ
j
(
k
,
y
)
,
ˆ
f
(
k
,
y
)
,
(3)
where
σ
1
σ
2
···
>
0 denote the singular values (gains) of
H
(
k
),
ˆ
ψ
j
=
[ˆ
u
j
,
ˆ
v
j
,
ˆ
w
j
,
ˆ
p
j
]
T
and
ˆ
φ
j
=
[
ˆ
f
uj
,
ˆ
f
v
j
,
ˆ
f
w
j
,
0]
T
are, respectively, the left and right singular vectors (wall-normal basis
functions) associated with
σ
j
, and
·
,
·
denotes the
L
2
-inner product in the wall-normal direction for
vector-valued functions. We will refer to the left singular vectors
ˆ
ψ
j
as resolvent modes hereafter.
The basis elements are ranked in order of descending gain
σ
j
and the singular values come in pairs
of equal or at least same order of magnitude values, which is a peculiarity of the channel geometry
(see Ref. [
8
] for details). The paired singular vectors are not unique (any linear combination of them
is also a legitimate singular vector) and an additional constraint is required to ensure uniqueness.
Following [
8
], we therefore impose an additional wall-normal symmetry constraint on the paired
singular vectors, which results in distinct wall-normal symmetries of the resolvent modes: ˆ
u
j
,ˆ
w
j
,
and ˆ
p
j
of the first paired mode are even while ˆ
v
j
is odd in
y
and the symmetries of all components
are inverted for the second mode of the pair.
Previous work has shown that the resolvent operator at each
k
is low rank [
3
,
8
] and owing to
the particular choice of bases the low-rank nature is reflected in expansion (
3
): Typically, the first
two singular values are an order of magnitude larger than the remaining ones and it can therefore be
073905-3
TOEDTLI, LUHAR, AND MCKEON
assumed that the flow is reasonably approximated by just the first two terms of the expansion. Since
almost all singular values are paired and the corresponding singular vectors only differ in their
wall-normal symmetry, it is often possible to simplify further and just consider the first singular
value and vector, which gives a so-called rank-1 approximation of the operator. It will be shown in
Sec.
II B
that a rank-1 approximation is sufficient for the purpose of this study.
While the singular values and vectors can be calculated directly from the resolvent, the nonlinear
forcing
ˆ
f
is not known
apriori
and has to be modeled to make further progress. It is well known
that the resolvent operator is a very selective amplifier [
3
] and it is therefore reasonable to assume
that the exact form of forcing is irrelevant as long as the real flow contains some forcing in the
dominant directions. Note, however, that there is an increasing body of work devoted to modeling
the shape of the forcing to accurately reconstruct the flow field [
10
–
12
]. In the present study we
seek the simplest representation of the forcing that captures the control trends for this system and
therefore we use the so-called rank-1 broadband forcing model
ˆ
φ
1
(
k
)
,
ˆ
f
(
k
)
=
1
∀
k
instead of the
more complex models mentioned above. To guard against a common misconception, we would like
to emphasize that the resolvent formulation of the NSE with broadband forcing is not a linearization.
A linearization would set
ˆ
f
(
k
)
=
0
, while the rank-1 broadband forcing assumption employed here
explicitly retains the nonlinearity and models it as
ˆ
f
(
k
)
=
ˆ
φ
1
(
k
).
The resulting rank-1 broadband forcing approximation of the NSE in resolvent form (
3
)isgiven
by
ˆ
u
(
k
,
y
)
ˆ
p
(
k
,
y
)
≈
σ
1
(
k
)
ˆ
ψ
1
(
k
,
y
)
.
(4)
To facilitate the following discussion, we will refer to Eq. (
4
) as the resolvent model hereafter. We
note that this approach was developed in a series of previous publications and was shown to capture
many key statistical and structural features of wall-bounded turbulent flows [
3
,
8
,
13
]. The subject
of this study is to assess whether the resolvent model is a suitable flow model to design controllers
for turbulent drag reduction in internal flows. The assessment is based on a generalized version of
the well-known opposition control scheme [
14
], which often serves as a benchmark within the flow
control community and which will be briefly reviewed hereafter. It is important to point out that
there is nothing particular about the choice of controller and we expect that the results generalize to
any control scheme with a linear control law.
B. Opposition control and low-order modeling
The idea at the heart of opposition control is to reduce turbulent drag in wall-bounded flows by
detecting and suppressing the sweeps and ejections associated with the coherent structures of the
near-wall cycle. To this end the control scheme measures the wall-normal velocity at a detection
plane located at a distance
y
d
above the wall (
y
w
=±
1) and applies blowing and suction with equal
amplitude but opposite sign at the wall:
v
(
y
w
)
=−
v
(
y
d
).
Previous direct numerical simulation (DNS) studies of turbulent channel [
14
] and pipe [
15
]flow
at Re
τ
≈
180 show that opposition control can reduce drag by up to 25%. The control effectiveness
strongly depends on the sensor location: Drag reduction (DR) is only observed for sensors located
below
y
+
d
23 (the plus superscript denotes normalization with respect to inner units, i.e.,
y
+
=
yu
τ
/ν
) and the maximum DR occurs around
y
+
d
=
15 [
16
]. The range of sensor locations leading to
drag reduction can be increased by reducing the controller gain [
16
] and the maximum DR can be
increased by either using upstream sensor information [
17
] or adding an integral term to the control
law [
18
]. The effectiveness of the control scheme depends on the Reynolds number and decreases
with increasing Re
τ
: For example, the maximum DR drops from approximately 25% at Re
τ
=
180
to 18% at Re
τ
=
1000 [
19
]. The computational cost of DNS currently inhibits studies at higher
Re
τ
, so DNS cannot be used to analyze and design a version of the control law for technologically
relevant Reynolds numbers.
073905-4
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
The prohibitive cost of DNS necessitates the development of low-order flow models for controller
design. Grounded in the insight that linear mechanisms play a key role in wall-bounded turbulence,
one line of work approaches flow control for internal flows from a linear system point of view. Lim
and Kim [
20
] used a singular-value analysis to assess the effect of opposition control on optimal
initial disturbances to the linearized NSE. In agreement with previous DNS results, their analysis
shows that the growth of the disturbances, which is interpreted as a proxy for turbulence intensity,
decreases for a narrow range of sensor locations
y
+
d
and increases if the sensor is located too far away
from the wall. In a similar manner but considering a different control strategy, Duque-Daza
et al.
[
21
] used the linearized NSE to study how streamwise-traveling waves of spanwise wall velocity
affect the energy amplification of initial perturbations that resemble elongated streamwise vortices.
Their results indicate a strong correlation between the change in energy amplification of the initial
perturbation and the drag reduction observed in DNS for traveling waves of various frequencies
and wavelengths. Of particular relevance for the present study is previous work by Luhar
et al.
[
13
], who used the resolvent model (
4
) to analyze opposition control in a turbulent pipe flow. In
contrast to the aforementioned work, the resolvent model does not capture how control affects the
growth of specific initial disturbances, but rather how control alters the response of the flow to
intrinsic forcing from the nonlinear advection terms. Luhar
et al.
[
13
] showed that the resolvent
model can qualitatively reproduce the changes in DR with sensor location and Reynolds number
previously observed in DNS. They further used the resolvent model to analyze a generalized version
of the opposition control law given by ˆ
v
(
k
,
y
w
)
=−
ˆ
A
d
ˆ
v
(
k
,
y
d
), where
ˆ
A
d
∈
C
is the complex
controller gain at each Fourier mode characterized by
k
and the original opposition control scheme
corresponds to
|
ˆ
A
d
|=
1 and
∠
ˆ
A
d
=
0. The resolvent model predicts that the effectiveness of
opposition control in turbulent pipe flow strongly depends on the phase
∠
ˆ
A
d
between the sensor
measurement and the actuation response, but this prediction has not yet been confirmed in the
full nonlinear system. The same resolvent model was subsequently used by Nakashima
et al.
[
22
]
to study the effect of suboptimal control in a turbulent channel flow and building on this work,
Kawagoe
et al.
[
23
] used resolvent analysis to design modified versions of the suboptimal control
law. In this case the agreement between resolvent prediction and DNS was mixed: The study
considered two modified controllers and good agreement between the model and DNS was observed
for one of them, while the two predictions deviated for the other controller. Similar ideas were
recently also applied to other classes of flows. Yeh and Taira [
24
] used the gains of the uncontrolled
resolvent operator as proxy to identify adequate control parameters for separation control over an
airfoil and showed that parameter combinations with high uncontrolled gain correspond to favorable
controller settings. Leclerq
et al.
[
25
] considered and controlled a two-dimensional incompressible
open-cavity flow using structured
H
∞
synthesis on a sequence of linear flow models, where each
flow model of the sequence is based on the resolvent operator about the mean of the previous
controlled flow. Their study shows that the controller is able to successfully suppress the oscillations
in a sequence of five controllers.
The literature results reviewed above suggest that linear models are able to capture at least
aspects of the response of internal flows to control. With regard to resolvent-based models, it
remains unclear to what extent the model is able to capture the response of the nonlinear system,
since the available data set for comparison (compiled in [
13
,
22
,
23
]) is small and shows mixed
agreement. Furthermore, previous resolvent-based control designs in internal flows have been
limited to low Reynolds numbers and it remains to be clarified how the approach can be generalized
to technologically relevant Re
τ
.
These two questions will be the main thrust of this study. Concretely, we consider the generalized
opposition control law with unit gain (
|
ˆ
A
d
|=
1) but nonzero phase (
∠
ˆ
A
d
=
0) and we will refer
to this control scheme as varying-phase opposition control hereafter. The resolvent model will
be used to predict the drag reduction of varying-phase opposition control in a turbulent channel
flow at Re
τ
=
180 and a DNS of the same flow and control scheme will be used to validate the
predictions over a wide range of parameters, which will allow the drawing of stronger conclusions
about the capabilities of the resolvent model. We also formalize an approach for wave-number
073905-5
TOEDTLI, LUHAR, AND MCKEON
space subsampling and outline how scaling laws can be used to make resolvent model predictions
at technologically relevant Reynolds numbers.
II. APPROACH
Two approaches are used to study a turbulent channel flow at Re
τ
=
180: The resolvent model
introduced in Sec.
IA
and a DNS. This section presents the details of the two frameworks,
outlines their commonalities and differences, and describes wave-number space subsampling in
the evaluation of the resolvent model.
A. Direct numerical simulation
The direct numerical simulations of the full nonlinear system are performed with a code
framework developed by Flores and Jimenéz [
26
]. The DNS follows the method proposed by Kim
et al.
[
27
] and solves the NSE in a velocity-vorticity formulation by integrating the wall-normal
vorticity and the Laplacian of the wall-normal velocity in time. The details about the numerical
method can be found in [
26
]; it is only mentioned here that the streamwise and spanwise directions
are discretized using a pseudospectral Fourier-Galerkin method. The DNS therefore naturally
provides access to the Fourier coefficients
ˆ
u
(
k
x
,
k
z
,
t
,
y
), which greatly simplifies the connection
and comparison between the DNS and the resolvent model: Leaving the order reduction in
y
out of
the account, the two formulations differ only in the treatment of the temporal coordinate and one
can map from the resolvent model to the DNS by means of an inverse Fourier transform in time.
The parameters of the DNS are chosen to match the domain size and resolution of previous
numerical simulations such as [
26
,
28
]. The nondimensional size of the computational domain in
the streamwise and spanwise directions is
L
x
=
4
π
and
L
z
=
2
π
, respectively, and
N
x
=
N
z
=
256 Fourier modes are used in these directions. This corresponds to resolutions of
x
+
≈
8
.
8
and
z
+
≈
4
.
4 in terms of Fourier modes before dealiasing at the nominal Re
τ
=
180 of the
uncontrolled flow. A sinusoidal grid with
N
y
=
172 points is used in the wall-normal direction,
which gives a resolution of
y
+
min
≈
0
.
37 at the wall and
y
+
max
≈
3
.
09 at the channel center.
Note that Re
τ
, as well as therefore the resolution in inner units, changes when control is applied.
The resolution increases if the drag is decreased and vice versa and runs with Re
τ
>
245 may
be considered slightly underresolved. Two runs (one for
y
+
d
=
15 and
∠
ˆ
A
d
=−
π/
4 and another
for
y
+
d
=
15 and
∠
ˆ
A
d
=
π/
2) with higher resolution (
N
x
=
512,
N
y
=
272, and
N
z
=
512) were
performed to rule out the possibility of grid effects in the results and it was indeed confirmed that
the results were almost identical.
The flow is driven by a constant mass flux so that the effect of control can be captured in the
change of wall shear stress. All control experiments are started from the same initial condition
and statistics are collected over at least ten eddy turnover times (
h
/
u
τ
) once a statistically steady
state is reached. It is also worth mentioning for later comparison that the cost of a single DNS was
around 120 core hours at the time of writing. The adequacy of the current settings is confirmed
by comparing the results of the present DNS to the literature data of an uncontrolled channel flow
at Re
τ
=
180 by Lee and Moser [
28
]: The mean Reynolds stress profiles of the two DNSs are
shown in Fig.
1(a)
(dash-dotted blue and solid black lines). It is apparent that the two profiles are
in excellent agreement throughout the channel and so are all the other statistical quantities available
for comparison (data not shown).
B. Rank-1 resolvent model
The low-order flow model put to test by DNS is the rank-1 resolvent model (
4
) introduced in
Sec.
IA
. Following a large body of work (e.g., [
8
,
29
,
30
]), we use the semiempirical turbulent eddy
viscosity model proposed by Reynolds and Tiedermann [
29
] with parameters
α
=
25
.
4 (constant in
the van Driest wall law) and
κ
=
0
.
426 (constant in the von Kármán logarithmic law) to represent
the mean velocity profile
u
of the resolvent operator. Such an approach may be extended to
073905-6
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
(a)
(b)
FIG. 1. Uncontrolled Reynolds stress profiles from DNS and the resolvent model: (a)
, DNS data of
Lee and Moser [
28
] (uncontrolled flow);
, present DNS (uncontrolled flow) and (b)
, resolvent model
baseline (uncontrolled flow);
, subsampled model (uncontrolled flow);
, linear fit to model Reynolds
stress profile in the outer flow.
regimes in which experimental or numerical results are not available and is appealing because of
the analytical formulation it provides, which is readily incorporated into resolvent codes. In this
study, we use the mean profile obtained using the eddy viscosity for resolvent analyses of both
uncontrolled and controlled flows, i.e., we do not obtain the mean profile for the controlled flows.
The latter constitutes a further modeling simplification, some implications of which are discussed
below. An energy norm is chosen for the singular-value decomposition of the resolvent and the
approach described by Luhar
et al.
[
9
] is used to compute the singular values and vectors.
As indicated earlier, a rank-1 approximation of the operator is sufficient for the present study:
It will be shown in Sec.
III B
that in order to evaluate the effect of control on the DR, we need to
compute the mean Reynolds stress profile and evaluate an integral of the form
I
=
1
−
1
y
(
u
v
)
dy
.
(5)
For a rank-2 model the contribution of a single Fourier mode to the mean Reynolds stress
u
v
=
u
v
(
k
=
0
) is proportional to
u
v
(
0
,
y
)
∼
σ
2
1
ˆ
u
1
(
k
,
y
)ˆ
v
∗
1
(
k
,
y
)
+
σ
1
σ
2
ˆ
u
1
(
k
,
y
)ˆ
v
∗
2
(
k
,
y
)
+
σ
1
σ
2
ˆ
u
2
(
k
,
y
)ˆ
v
∗
1
(
k
,
y
)
+
σ
2
2
ˆ
u
2
(
k
,
y
)ˆ
v
∗
2
(
k
,
y
)
,
(6)
where
denotes the real part of a complex number and the asterisk superscript denotes its complex
conjugate. Recall from the discussion about the wall-normal symmetry of resolvent modes (Sec.
IA
)
that if the singular values are paired, then ˆ
u
1
and ˆ
v
2
are even, while ˆ
u
2
and ˆ
v
1
are odd in
y
.
Consequently, the cross terms (
σ
1
σ
2
ˆ
u
1
ˆ
v
∗
2
and
σ
1
σ
2
ˆ
u
2
ˆ
v
∗
1
)areevenin
y
. The cross terms integrate
to zero when weighted with
y
and therefore they do not contribute to the integral
I
. Moreover, the
real parts of the two diagonal terms (
σ
2
1
ˆ
u
1
ˆ
v
∗
1
and
σ
2
2
ˆ
u
2
ˆ
v
∗
2
) are equal and odd in
y
. The diagonal terms
therefore do contribute to
I
, but since their real parts are equal there is no difference between a
rank-1 and a rank-2 approximation. Recall that almost all singular values are paired, so the above
statement is true for almost all resolvent modes; it can therefore be expected that the DR (which is
based on integrals of the form
I
) predicted by a rank-1 and a rank-2 model is almost identical and
this has indeed been confirmed in numerical experiments (data not shown). For the purpose of the
present study it thus suffices to consider a rank-1 approximation, as anticipated in Sec.
IA
.
073905-7
TOEDTLI, LUHAR, AND MCKEON
(a)
(b)
FIG. 2. Effect of wave speed and streamwise wave number on the mean Reynolds stress contribution
of individual resolvent modes: (a) effect of wave speed
c
(for
k
x
=
5
.
5,
k
z
=
2, and
c
+
={
8
,
10
,
12
,
14
}
),
with darker colors indicating faster wave speeds, and (b) effect of streamwise wave number
k
x
(for
k
x
=
{
5
,
5
.
5
,
6
,
6
.
5
}
,
k
z
=
2, and
c
+
=
10), with darker colors indicating smaller streamwise wave numbers.
The resolvent model uses the same numerical parameters as the DNS wherever possible, so
model errors due to resolution discrepancies can be ruled out and a meaningful comparison can
be made. This choice has the additional benefit that minimal empirical knowledge is required to
select an appropriate set of model parameters. In this spirit, we resolve the same streamwise and
spanwise wave numbers in the model as in the DNS (after dealiasing):
k
x
=
l
/
2,
l
∈
[
−
84
,
84],
and
k
z
=
m
,
m
∈
[
−
84
,
84]. In practice, only the non-negative wave numbers are evaluated, since
the wave-number symmetries discussed in Sec.
II D
can be used to account for the contribution of
the negative
k
x
and
k
z
. In the wall-normal direction the resolvent operator is discretized on a set of
N
y
=
172 Chebyshev collocation points so that the grid of the model and the DNS are identical up
to a stretching factor described in [
28
].
The only coordinate that requires empirical input for parameter selection is the temporal
frequency: The continuous range of frequencies resolved in the DNS is too large to be handled
reasonably by the numerical framework of the model and furthermore the DNS does not provide a
well-defined sampling rate because the time step changes subject to a Courant-Friedrichs-Lewy
(CFL) condition. However, there is compelling empirical evidence that the temporal frequency
content of each Fourier mode is approximately sparse, so the entire frequency spectrum does not
need to be resolved. As summarized by Bourguignon
et al.
[
31
], the dominant temporal frequencies
at
k
x
=
0 can be parametrized by the wave speed
c
=
ω/
k
x
and fall within the empirical range
10
c
U
cl
, where
U
cl
is the channel centerline velocity. This suggests picking the temporal
frequency vector
ω
at each
k
x
such that at least the range 10
ω/
k
x
U
cl
can be resolved. For
reasons that will become clear shortly, we use the slightly more conservative range 0
ω/
k
x
U
cl
in this study. Note that since the range of wave speeds is kept constant across all
k
, the corresponding
temporal frequencies change for different
k
x
.
An appropriate sampling rate in
c
can be derived from the wall-normal localization of the
resolvent modes: Fig.
2(a)
shows that the mean Reynolds stress contribution of a single resolvent
mode is localized around its critical layer, i.e., around the wall-normal location
y
c
, where its
wave speed equals the local mean velocity
c
≈
u
(
y
c
). Conversely, one can say that the dominant
contribution to the Reynolds stress at a fixed
y
comes from modes whose critical layer
y
c
corresponds
to that
y
. The discretization of the wall-normal coordinate naturally samples the mean velocity
profile and defines the critical wave speeds that are resolved by the grid. From the previous argument
we expect that the dominant Reynolds stress contribution at each grid point
y
i
is given by modes
with wave speed
c
=
u
(
y
i
). This suggests that the sampled wave speeds should correspond to
073905-8
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
the discretized mean velocity profile so that 0
ω
k
x
U
cl
(as mentioned before) and
ω
i
(
k
x
)
=
k
x
u
(
y
i
). Note that empirical knowledge is only required to justify this range of
ω
; no empirical
knowledge is required to evaluate it. The temporal frequency vector is fully determined by the given
mean velocity profile and wall-normal grid.
The exceptions to the above discussion are the streamwise constant modes (
k
x
=
0). For such
modes the wave speed
c
is not well defined and a different, for now empirical, approach is needed to
determine an appropriate frequency vector. The amplification of streamwise constant modes peaks
at
ω
=
0 and drops off very quickly and symmetrically as
ω
moves away from the origin. For
example, the maximum
σ
1
across all
k
z
drops from
σ
1
=
696 at
ω
=
0to
σ
1
=
34 at
ω
=±
0
.
4
and
σ
1
=
2at
ω
=±
3
.
5 and further decreases as
|
ω
|
increases (data not shown in the interest of
brevity). Recall that the Reynolds stress contribution depends on the square of the singular value and
therefore contributions beyond
ω
=±
3
.
5 can be neglected. We thus restrict the temporal frequency
vector of
k
x
=
0to
−
3
.
5
ω
3
.
5 and sample densely around the origin (
ω
=
0
.
01 for
−
0
.
25
ω
0
.
25) and more coarsely away from it (
ω
=
0
.
05 for 0
.
3
|
ω
|
0
.
5 and
ω
=
0
.
25 for
0
.
75
|
ω
|
3
.
5).
It was verified that the results presented in this study do not depend on the particular choice of
ω
: Across all
k
x
the DR obtained for various
∠
ˆ
A
d
varies only minimally if the bandwidth of
ω
is
increased or if the sampling frequency
ω
is changed. For reasons that will become more clear
in Sec.
II C
, we will refer to resolvent model runs with the above parameter settings as the model
baseline. Also, for future reference it is worth mentioning that one model baseline run costs around
160 core hours at the time of writing, i.e., the evaluation of the model baseline and DNS is similarly
expensive at this Re
τ
.
The Reynolds stress profile predicted by the model baseline for an uncontrolled channel flow
at Re
τ
=
180 is shown as a dash-dotted blue line in Fig.
1(b)
. The model baseline resolves all
dynamically relevant spatial and temporal scales and therefore can be considered a resolvent model
prediction at DNS resolution. However, it is apparent from the figure that the model profile does
not match the DNS Reynolds stress: The model overpredicts the peak magnitude [compare the
y
axes of Figs.
1(a)
and
1(b)
] and the modeled profile peaks closer to the wall than the real one
and exhibits a plateau within
−
0
.
75
y
/
h
<
−
0
.
4 instead of a near-linear decrease. Recall that the
model is not constrained to satisfy the mean momentum equation, since we neglected the latter in the
model development. However, it is interesting to note that the modeled Reynolds stress still exhibits
a near-linear decrease in the outer flow, as indicated by the gray line in Fig.
1(b)
. In agreement
with the overpredicted magnitude of the Reynolds stress, the slope is steeper than what one would
expect from an analysis of the mean momentum equation. These observations are in agreement with
previous results from a similar study of a turbulent pipe flow by Luhar
et al.
[
13
]. The only difference
from the pipe flow result is the heavier tail in the outer channel (
y
/
h
>
−
0
.
5), which can be ascribed
to the streamwise constant (
k
x
=
0) modes with zero temporal frequency and low spanwise wave
numbers: These dominant streamwise constant modes peak close to the channel center and have
the largest singular values of all resolved wave numbers. Under the broadband forcing assumption
their contributions have a strong footprint on the integrated Reynolds stress profile and therefore
give rise to the observed heavy tail. The role of the
k
x
=
0 modes within the resolvent framework
remains to be fully understood and an in-depth analysis is beyond the scope of this investigation, but
we would like to point out that the largest amplification occurring for streamwise constant modes
is reminiscent of transient growth, where streamwise constant vortices are well known to be the
optimal perturbations [
32
].
Given that the resolvent model is not constrained to satisfy the mean momentum equation and the
numerous model assumptions, in particular the broadband forcing, the observed deviations between
the model and DNS profile may not be surprising. Nevertheless, in view of recent work on more
complex forcing models that approximate the correct flow statistics [
10
–
12
] one may ask whether
a rank-1 model with broadband forcing is still adequate. However, the key point not to forget in a
control context is this: A useful control-oriented model mainly needs to capture the relation between
inputs and outputs; it need not produce the correct flow statistics of the system [
33
]. In this spirit the
073905-9
TOEDTLI, LUHAR, AND MCKEON
goal should be to find the simplest possible model that is able to approximate the response of the full
nonlinear system to control and in this sense the broadband forcing assumption is preferred over the
aforementioned more complex models. Furthermore, unlike the other models, the broadband forcing
assumption does not require flow data as input, which minimizes the dependence on empirical data.
C. Subsampling in wave-number space
While the model baseline is the appropriate starting point for comparison with DNS, its eval-
uation is computationally too expensive for practical purposes like parameter studies or controller
design at higher Re
τ
. An analysis of the mean Reynolds stress contribution of individual resolvent
modes reveals that their contribution varies smoothly in at least parts of the wave-number space.
This observation suggests that the resolvent model could take advantage of subsampling across
wave-number space to reduce the number of resolved scales compared to the DNS. This approach
has already been used in previous studies (for example, [
8
,
13
]), but without formal justification. The
goal of this section is to demonstrate that subsampling is indeed appropriate.
As can be seen from Fig.
2(a)
, the mean Reynolds stress contribution of a single mode is localized
around the critical layer
y
c
and therefore its wall-normal profile is mainly determined by the wave
speed
c
. The wall-normal profiles for different
c
look very similar and the location of their peak
moves slowly away from the wall as the wave speed is increased. This suggests that not all the
wave speeds resolved in the model baseline are required to capture the wall-normal shape of the
mean Reynolds stress profile, since resolvent modes with similar wave speeds largely overlap in
y
.
Figure
2(b)
further shows the mean Reynolds stress contribution for modes with various
k
x
at fixed
k
z
and
c
. It is again apparent that the wall-normal localization of the modes is fixed by the wave
speed. Furthermore, it can be seen that the wall-normal profiles look very similar for various
k
x
and
that the peak magnitude slowly decreases as
k
x
increases. This observation holds for all sufficiently
large streamwise and spanwise wave numbers (
k
x
5 and
k
z
10) and suggests that the spatial
wave numbers can be subsampled as well as the wave speeds.
Based on these insights, the wave-number space is subsampled as follows to evaluate the
mean Reynolds stress (the contributions of negative frequencies can be obtained from the
symmetries discussed next):
k
x
=
[0
,
0
.
5
,... ,
5
,
7
.
5
,
10
,
20
,
31
,
42],
k
z
=
[0
,
1
,... ,
10
,
12
,
14
,
16
,
20
,
25
,
30
,
40
,
50
,
60
,
70
,
84], and only every other value of
u
(
y
) is used to populate
c
.This
corresponds to a reduction from 85
×
85
×
86 (baseline) to 16
×
22
×
44 (subsampled model) re-
solved wave numbers. The missing wave numbers are linearly interpolated and
ω
is not subsampled
if
k
x
=
0.
The resulting Reynolds stress profile of the uncontrolled flow is shown as a solid black line in
Fig.
1(b)
. The good agreement with the model baseline (dash-dotted blue line) confirms that the
wave-number space can indeed be sampled very sparsely. The subsampled model only resolves
about 2% of the wave numbers of the baseline, which reduces the computational time to 2 core
hours. This is cheap enough to allow, for example, extensive parameter searches for controller
design. Furthermore, the memory requirements and operation count of the subsampled model are
small enough so that it can be run on an off-the-shelf laptop and no high-performance computing
facilities are required to evaluate it.
Motivated by the results of this section, we will use the subsampled model for all resolvent
predictions hereafter and we will make all comparisons between the subsampled resolvent model
and DNS.
D. Wave-number symmetries
Before moving on to the description of the control scheme, we take a closer look at the
wave-number symmetries of the DNS and the resolvent model. The different symmetries present
in the two frameworks will have a bearing on the formulation of the control scheme and, generally
speaking, are an important constraint to keep in mind when linking the resolvent model and DNS.
073905-10
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
The first wave-number symmetry, which is present in both the DNS and the resolvent model, follows
from the properties of the Fourier transform: The Fourier coefficients of the velocity and pressure
fluctuations are Hermitian symmetric, because the underlying physical quantities are real-valued
functions. For example, the Fourier coefficients of the velocity field satisfy
ˆ
u
(
−
k
x
,
−
k
z
,
t
,
y
)
=
ˆ
u
∗
(
k
x
,
k
z
,
t
,
y
) (DNS) and
ˆ
u
(
−
k
,
y
)
=
ˆ
u
∗
(
k
,
y
) (resolvent model). A second symmetry, which is
present in the resolvent model but not in the DNS, follows from the structure of the resolvent
operator
H
(
k
) and the broadband forcing assumption: Evaluating Eq. (
2
)at
̃
k
=
[
k
x
,
−
k
z
,ω
]
T
instead of
k
=
[
k
x
,
k
z
,ω
]
T
only changes the sign of two off-diagonal terms in the resolvent
operator. The sign change can be absorbed into the coefficient vectors by making the substitutions
ˆ
u
†
(
̃
k
,
y
)
=
[ˆ
u
,
ˆ
v
,
−
ˆ
w
,
ˆ
p
]
T
and
ˆ
f
†
(
̃
k
,
y
)
=
[
ˆ
f
u
,
ˆ
f
v
,
−
ˆ
f
w
,
0] so that
ˆ
u
†
(
̃
k
,
y
)
=
H
(
k
)
ˆ
f
†
(
̃
k
,
y
). The
broadband forcing assumption sets
ˆ
f
†
(
̃
k
,
y
)
=
ˆ
f
(
k
,
y
) and thus imposes the additional symmetry
ˆ
u
†
(
̃
k
,
y
)
=
ˆ
u
(
k
,
y
).
As a consequence of these two symmetries, the phase of a resolvent mode has a clear interpreta-
tion in physical domain. To see this, we consider the example of the wall-normal velocity component
and write the rank-1 approximation as ˆ
v
1
(
k
,
y
)
=
ˆ
v
1
(
̃
k
,
y
)
=
A
1
e
i
φ
1
, where
{
A
1
,φ
1
}∈
R
. Evaluating
the inverse Fourier transform (
1
) of only ˆ
v
1
(
k
,
y
), ˆ
v
1
(
̃
k
,
y
), and their complex conjugates gives a
structure
v
RES
k
in the physical domain
v
RES
k
(
x
,
y
,
z
,
t
)
=
4
A
1
cos(
k
x
x
−
ω
t
+
φ
1
) cos(
k
z
z
)
,
(7)
where the factor 4 results from Hermitian symmetry and the trigonometric sum-to-product identity.
Due to the symmetries, the phase
φ
1
only enters the first term and a nonzero
φ
1
can be interpreted as
a streamwise or temporal shift of the structure. Conversely, it is apparent that the spanwise position
of the structure is unaffected by a change in
φ
1
.
With regard to the DNS, the situation is more complex: Even though the DNS Fourier
coefficients are symmetric in
k
z
in a statistical sense, their instantaneous values are unrelated,
i.e.,
ˆ
u
(
k
x
,
−
k
z
,
t
,
y
)
=
ˆ
u
(
k
x
,
k
z
,
t
,
y
), and the missing symmetry in
k
z
complicates the interpre-
tation of the phase in the physical domain. To see this we consider again the wall-normal
velocity component and introduce the notation ˆ
v
(
k
x
,
k
z
,
t
,
y
)
=
A
1
e
i
φ
1
and ˆ
v
(
k
x
,
−
k
z
,
t
,
y
)
=
A
2
e
i
φ
2
,
where
{
A
1
,
A
2
,φ
1
,φ
2
}∈
R
. Evaluating the inverse Fourier transform (
1
) of only ˆ
v
(
k
x
,
k
z
,
t
,
y
),
ˆ
v
(
k
x
,
−
k
z
,
t
,
y
), and their complex conjugates gives again a structure
v
DNS
[
k
x
,
k
z
]
in the physical domain
v
DNS
[
k
x
,
k
z
]
(
x
,
y
,
z
,
t
)
=
4
A
1
cos
k
x
x
+
φ
1
+
φ
2
2
cos
k
z
z
+
φ
1
−
φ
2
2
+
2(
A
2
−
A
1
) cos(
k
x
x
−
k
z
z
+
φ
2
)
,
(8)
where we can assume
A
1
<
A
2
without loss of generality. The first summand, in which
x
and
z
factor
into different terms, is the equivalent of
v
RES
k
. The phases of the Fourier coefficients
φ
1
and
φ
2
have
a clear physical interpretation for this structure: If the mean of the two phases changes, the structure
is shifted in the streamwise direction, and if their difference changes, the structure is shifted in the
spanwise direction. Due to the missing symmetry in
k
z
, an additional term with mixed arguments in
x
and
z
appears which has no counterpart in the resolvent model. For this structure there is no clear
interpretation of the phase either: Since the argument is mixed, a change in
φ
2
can correspond to a
streamwise or spanwise shift.
Note that Eq. (
8
) is a generalization of
v
RES
k
and the latter can be recovered by setting
A
2
=
A
1
and
φ
2
=
φ
1
. Whenever the (broadband forcing) resolvent model and DNS are linked, it should
therefore be kept in mind that the model only represents a subset of the real physics and caution is
required to ensure that the same problem is studied in both frameworks.
III. CONTROL LAW
After characterizing the two frameworks to be compared, we formally introduce the varying-
phase opposition control law and the DR measure in this section.
073905-11
TOEDTLI, LUHAR, AND MCKEON
A. Varying-phase opposition control
As mentioned in Sec.
IB
, Luhar
et al.
[
13
] studied varying-phase opposition control in a pipe
flow by means of the resolvent model. They defined the control law as ˆ
v
(
k
,
y
w
)
=−
ˆ
A
d
ˆ
v
(
k
,
y
d
),
with
ˆ
A
d
=
e
i
φ
for all
k
, and one may think that the same control law can be incorporated directly
in the DNS. However, the control law is a good example of the importance of the wave-number
symmetries discussed in Sec.
II D
and due to the different symmetries present in the DNS and
resolvent model, additional care is required to ensure that both frameworks incorporate the same
physics.
The previous discussion about symmetries shows that a phase shift
∠
ˆ
A
d
=
φ
in the resolvent
model has a clear physical interpretation: Eq. (
7
) implies that multiplying the Fourier coefficient
with a nonzero phase corresponds to a streamwise or temporal shift of the mode. A phase shift
between the actuator and sensor in the Fourier domain can therefore be thought of as a streamwise
or temporal lead or lag of the actuation in physical domain. However, this physical interpretation is
lost when applying the same control law in the DNS, due to the missing symmetry in
k
z
: In addition
to the streamwise shift, the mixed-argument term in Eq. (
8
) introduces a spanwise shift between
the sensor and the actuator, which cannot be captured by the resolvent model. In order to guarantee
a streamwise-only shift in the DNS so that the same physics is studied in both frameworks, the
varying-phase opposition control law has to be modified slightly,
ˆ
v
(
k
,
y
w
)
=−
ˆ
A
RES
d
ˆ
v
(
k
,
y
d
)
,
ˆ
A
RES
d
=
⎧
⎨
⎩
0if
k
x
=
k
z
=
0
1if
k
x
=
0
,
k
z
=
0
e
i
φ
otherwise
,
ˆ
v
(
k
x
,
k
z
,
t
,
y
w
)
=−
ˆ
A
DNS
d
ˆ
v
(
k
x
,
k
z
,
t
−
t
,
y
d
)
,
ˆ
A
DNS
d
=
⎧
⎪
⎪
⎨
⎪
⎪
⎩
0if
k
x
=
k
z
=
0
1if
k
x
=
0
,
k
z
=
0
e
i
φ
if
k
x
=
0
,
k
z
=
0
min[
|
ˆ
v
(
k
x
,
k
z
,
t
−
t
,
y
d
)
|
,
|
ˆ
v
(
k
x
,
−
k
z
,
t
−
t
,
y
d
)
|
]
|
ˆ
v
(
k
x
,
k
z
,
t
−
t
,
y
d
)
|
e
i
φ
otherwise
,
(9)
where the formulation for the resolvent model (RES) is shown at the top and the equivalent
formulation for the DNS at the bottom.
The particular values for
ˆ
A
d
are explained as follows. In order to ensure that control does not
change the net mass flux of the flow we have to set
ˆ
A
d
=
0if
k
x
=
k
z
=
0. No phase shift, i.e.,
ˆ
A
d
=
1, is applied to the remaining streamwise constant modes (
k
x
=
0 and
k
z
=
0), since a nonzero
phase would by definition lead to a spanwise shift (at least in the DNS). The values of
ˆ
A
d
for
k
x
=
0 are chosen such that a streamwise shift is guaranteed: A streamwise shift is ensured for
the spanwise constant modes (
k
x
=
0 and
k
z
=
0) of the DNS and we can therefore set
ˆ
A
DNS
d
=
e
i
φ
for those modes. The last value of
ˆ
A
DNS
d
requires some more explanation. Spanwise shifts can be
prevented if the control input only generates the first term of Eq. (
8
). This can be achieved by
setting
A
2
=
A
1
, i.e., enforcing
|
ˆ
v
(
k
x
,
k
z
,
y
w
,
t
)
|=|
ˆ
v
(
k
x
,
−
k
z
,
y
w
,
t
)
|
. While either amplitude would
be a valid choice, it is desirable from a control perspective to pick the smaller of the two in order to
ensure a minimal controller input. The particular choice of
ˆ
A
d
ensures just that and sets
A
2
=
A
1
=
min[
|
ˆ
v
(
k
x
,
k
z
,
y
d
,
t
)
|
,
|
ˆ
v
(
k
x
,
−
k
z
,
y
d
,
t
)
|
]inEq.(
8
). Note that this choice leaves the mixed-argument
term uncontrolled, but ensures a streamwise-only shift. The value of
ˆ
A
d
for the
k
x
=
0 resolvent
modes is simpler to explain: Recall from Eq. (
7
) that the additional symmetry of the resolvent
modes ensures a streamwise shift, so we can set
ˆ
A
RES
d
=
e
i
φ
for all
k
x
=
0 modes.
Two comments are important to mention before moving on to the DR measure. First, following
[
13
], the potential changes in the mean velocity profile and nonlinear forcing due to control are, for
simplicity, not incorporated in the resolvent model. In other words, the uncontrolled mean velocity
profile and broadband forcing assumption are also used in the controlled cases and the controlled
073905-12
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
and uncontrolled resolvent models differ only in the boundary conditions on
v
. Second, note that
there is a time delay of one time step between the sensor measurement and the actuator response
in the DNS, because the actuation is implemented as a Dirichlet boundary condition. The time
delay is equivalent to an additional frequency-dependent phase
e
i
ω
t
in the Fourier domain, which
is not represented in the resolvent model. The time step
t
varies subject to a CFL condition, but
is of the order 10
−
3
. The phase error introduced to the structures of the near-wall cycle, which are
expected to dominate the control signal and can be characterized by
λ
+
x
≈
1000 and
c
+
≈
10–12
[
13
], is
ω
t
≈
π/
450. This is two orders of magnitude smaller than the phase shifts induced by the
controller and it can therefore be expected that the phase error due to the time delay is negligible.
Simulations with smaller
t
and therefore smaller phase errors were performed to validate this
assumption and it was indeed confirmed that the drag reduction obtained for smaller
t
was nearly
the same.
B. Drag reduction measure
Following previous DNS studies [
14
,
16
], we quantify the DR in terms of change in mean wall
shear stress. For a constant mass flux channel flow, we can use the Fukagata-Iwamoto-Kasagi (FIK)
identity [
34
] (adjusted to the coordinate system and nondimensionalization of the present paper) to
express the change in mean wall shear stress as
τ
=
(
τ
w
)
u
−
(
τ
w
)
c
(
τ
w
)
u
=
1
−
1
y
[(
u
v
)
u
−
(
u
v
)
c
]
dy
Re
b
(Re
2
τ
)
u
+
1
−
1
y
(
u
v
)
u
dy
,
(10)
where we used the friction velocity of the uncontrolled flow (
u
τ
)
u
to nondimensionalize all velocities
(the subscripts
u
and
c
label quantities of the uncontrolled and controlled flows, respectively) and
where Re
b
=
2
U
b
h
/ν
is the Reynolds number based on twice the bulk velocity
U
b
. As explained in
[
34
], the denominator of Eq. (
10
) is a sum of a laminar contribution (first term), fully determined
by the mean velocity profile, and a turbulent contribution (second term), fully determined by the
weighted integral of the Reynolds stresses.
The discussion in Sec.
II B
underpinned that the resolvent model does not satisfy the mean
momentum equation and Fig.
1
in particular demonstrated that the modeled Reynolds stress does
not match the true one. As a consequence, the resolvent model does not predict the correct relative
magnitude of the laminar and turbulent drag contribution and it would therefore not be meaningful
to evaluate the DR according to Eq. (
10
). Instead, we use the mean momentum equation to recast
Eq. (
10
) into a more appropriate form for the purpose of the present study
τ
=
1
−
1
y
[(
u
v
)
u
−
(
u
v
)
c
]
dy
1
−
1
y
(
u
v
)
u
dy
=
T
1
(
y
d
,
∠
ˆ
A
d
)
1
−
3
2
Re
b
Re
2
τ
u
=
T
2
×
1
−
1
y
[(
u
v
)
DNS
u
−
(
u
v
)
DNS
c
]
dy
1
−
1
y
[(
u
v
)
u
−
(
u
v
)
c
]
dy
=
T
3
(
y
d
,
∠
ˆ
A
d
)
1
−
1
y
(
u
v
)
u
dy
1
−
1
y
(
u
v
)
DNS
u
dy
=
T
4
,
(11)
where the superscript DNS labels a quantity that is obtained from DNS and which therefore satisfies
the mean momentum equation. This is in contrast to quantities without a superscript label which,
depending on the context, originate from DNS or the resolvent model. The term
T
1
(
y
d
,
∠
ˆ
A
d
)
represents the turbulent DR predicted by the model under consideration (resolvent or DNS) and
as such depends on the sensor location
y
d
and the phase shift
∠
ˆ
A
d
. The second term
T
2
is the ratio
of the turbulent drag to the total drag in the real uncontrolled flow, which is a constant for fixed Re
b
.
The product
T
1
(
y
d
,
∠
ˆ
A
d
)
T
2
represents the total change in wall shear stress, possibly subject to model
errors represented by the last two terms:
T
3
(
y
d
,
∠
ˆ
A
d
) is the model error in turbulent DR, which is a
073905-13
TOEDTLI, LUHAR, AND MCKEON
−
3
/
4
−
1
/
2
−
1
/
4
0
1
/
4
1
/
2
3
/
4
ˆ
A
d
/π
5
10
15
20
25
y
+
d
−
7
.
00
−
3
.
00
0
.
00
0
.
00
0
.
96
−
20
−
15
−
10
−
5
0
1
(a)
−
3
/
4
−
1
/
2
−
1
/
4
0
1
/
4
1
/
2
3
/
4
ˆ
A
d
/π
5
10
15
20
25
y
+
d
−
7
.
00
−
7
.
00
−
3
.
00
−
3
.
00
0
.
00
0
.
96
−
20
−
15
−
10
−
5
0
1
(b)
FIG. 3. Contour map showing the normalized drag reduction
ξ
as a function of sensor location
y
+
d
and
phase shift
∠
ˆ
A
d
for (a) the prediction of the resolvent model and (b) the DNS data. Note that the color scale is
nonlinear in order to highlight the region of drag reduction. The vertical dashed black line denotes
∠
ˆ
A
d
=
0,
which is closely related to opposition control, and the contour lines of (a) are replotted as dotted blue lines
in (b). Nondimensionalization to inner units is based on (
u
τ
)
u
of the uncontrolled flow and the maximum DR
used to normalize the contour maps is (a) resolvent model max DR
=
5% and (b) DNS max DR
=
21%.
function of the control parameters, and the constant
T
4
is the model error in the uncontrolled profile.
Note that for DNS data
T
3
(
y
d
,
∠
ˆ
A
d
)
=
T
4
=
1, so the usual definition for DR is recovered.
With DNS data we can evaluate all terms of Eq. (
11
), while with the resolvent model we can
evaluate
T
1
(
y
d
,
∠
ˆ
A
d
), which solely depends on model data, and
T
2
, which can be obtained from the
input (Re
τ
)
u
and the corresponding mean profile. The model errors
T
3
(
y
d
,
∠
ˆ
A
d
) and
T
4
cannot be
obtained from the resolvent model alone. In other words, with the resolvent model we can estimate
the change in mean wall shear stress up to a model error
DR
=
T
1
(
y
d
,
∠
ˆ
A
d
)
T
2
=
τ
T
3
(
y
d
,
∠
ˆ
A
d
)
T
4
.
(12)
From here on, the term drag reduction will be used to denote the quantity in Eq. (
12
). When
comparing different control configurations, as will be done in Sec.
IV
, it is advantageous to
normalize the DR by a reference value, for example the maximum DR (max DR for short) obtained
over the parameter space considered
ξ
=
DR
max DR
=
τ
max
τ
T
3
(
y
d
,
max
,
∠
ˆ
A
d
,
max
)
T
3
(
y
d
,
∠
ˆ
A
d
)
.
(13)
This allows the elimination of
T
4
, the model error in the uncontrolled profile.
IV. RESULTS AND DISCUSSION
A. Comparison between the rank-1 resolvent model and DNS
The varying-phase opposition control scheme (
9
) introduced in the preceding section
is used as a test to evaluate the capabilities of the resolvent model for controller de-
sign. To this end, a total of 50 varying-phase opposition control runs covering a pa-
rameter range of five sensor locations
y
+
d
=
[5
,
10
,
15
,
20
,
25] and ten phase shifts
∠
ˆ
A
d
=
[
−
3
π/
4
,
−
π/
2
,
−
3
π/
8
,
−
π/
4
,
−
π/
8
,
0
,π/
8
,π/
4
,π/
2
,
3
π/
4] were performed in the subsam-
pled resolvent model and the DNS. The raw DR data were subsequently interpolated using bilinear
splines and normalized with the respective maximum DR to generate the two contour maps of
073905-14
PREDICTING THE RESPONSE OF TURBULENT CHANNEL ...
ξ
shown in Fig.
3
. The capabilities of the resolvent model will be judged based on the overall
agreement between the two maps.
It should be mentioned that the largest drag increase among all DNS runs was observed for
∠
ˆ
A
d
=
3
π/
4 and
y
+
d
=
[20
,
25]. The Re
τ
of these runs increased by more than a factor of 2 and the
simulations developed numerical instabilities due to insufficient grid resolution before a statistically
steady state was reached. The two runs were not repeated with higher resolution, first because the
DR trend could already be estimated from the partially completed runs and second because the exact
flow statistics of parameter combinations leading to high drag increase are not of particular interest
to this study. Instead, the maximum drag increase observed in all successfully completed runs was
assigned to these two control cases (which is a lower bound for the effective drag increase) and the
saved computational time was used towards additional runs in the parameter region of most interest
(note the higher resolution in
∠
ˆ
A
d
within the range
−
π/
2
∠
ˆ
A
d
π/
4).
The resulting map of
ξ
as a function of
y
d
and
∠
ˆ
A
d
is shown in Fig.
3
, where Fig.
3(a)
is the
resolvent model prediction and Fig.
3(b)
shows the DNS results. Recall from Eq. (
13
) that each map
is normalized by its maximum DR, which corresponds to 5% [resolvent model, Fig.
3(a)
] and 21%
[DNS, Fig.
3(b)
], respectively. Bright shading (positive numbers) represent drag reduction, while
dark colors (negative numbers) indicate drag increase. The solid black lines outline a few selected
contour levels and the contour lines of Fig.
3(a)
are replotted as dotted blue lines in Fig.
3(b)
to
facilitate comparison between resolvent model and DNS. Note that the contour levels of both plots
are identical.
Figure
3
shows that the resolvent model is able to capture the trend observed in DNS over a wide
range of parameters: In both frameworks, the effect of the controller strongly depends on the phase
shift
∠
ˆ
A
d
and generally speaking a small negative shift (e.g.,
∠
ˆ
A
d
=−
π/
4) leads to improved drag
reduction, while a positive phase shift deteriorates the control performance and eventually leads to
drag increase. Furthermore, both frameworks show that for a fixed positive phase shift
∠
ˆ
A
d
>
0
the control performance decreases as
y
d
increases, while for a fixed negative phase shift
−
π/
2
∠
ˆ
A
d
<
0 the control performance initially increases, reaches a maximum, and then decreases as the
sensor moves away from the wall. The qualitative agreement between the resolvent model and DNS
can be made more clear by comparing the overlaid contour lines in Fig.
3(b)
: It is apparent that the
contours of
ξ
of the resolvent model and the DNS collapse for
∠
ˆ
A
d
>
0 over the entire range of
y
d
. This suggests that the model error in turbulent DR,
T
3
(
y
d
,
max
,
∠
ˆ
A
d
,
max
)
/
T
3
(
y
d
,
∠
ˆ
A
d
), is constant
and approximately equal to unity for positive phase shifts. Further note that
T
3
(
y
d
,
max
,
∠
ˆ
A
d
,
max
)is
a constant and therefore
T
3
(
y
d
,
∠
ˆ
A
d
)
≈
const for all positive shifts. The agreement between the
resolvent model and DNS is reduced for negative phase shifts: The contour levels do not overlap for
∠
ˆ
A
d
<
−
π/
4 and the model error in turbulent DR increases for more negative
∠
ˆ
A
d
. However, the
resolvent model is still able to capture the trend of the DNS reasonably well and also the parameter
combination leading to maximum drag reduction is similar in both frameworks and corresponds to
y
+
d
≈
10
−
15 and
∠
ˆ
A
d
≈−
π/
4.
The explanation of the optimal sensor location and phase shift is beyond the scope of the present
work, as is a more complete understanding of the physics of the phase shifts. The multimode nature
of the problem complicates signaling time (ballistic) arguments, which would attempt to explain the
optimal sensor location and phase shift by relating the time the control signal requires to travel from
the wall to the sensor location to the time a dominant flow structure requires to travel the streamwise
shift between the sensor and actuator. Yet it cannot currently be ruled out that control mainly acts
on the near-wall modes and a more detailed analysis of the optimal phase shift and sensor location
is left for future work.
While the resolvent model is able to capture the qualitative trend of the DNS, it underestimates
the absolute DR values observed in the real flow. For example, the resolvent model predicts a
maximum DR of 5% at the aforementioned optimal parameter combination
y
+
d
=
15 and
∠
ˆ
A
d
=
−
π/
4, while a 21% DR is observed in the DNS. In combination with Fig.
3(b)
, which shows that
the resolvent model underpredicts
|
ξ
|
, it follows that the model underpredicts the drag reduction
over the entire parameter range. We can therefore conclude from Eq. (
12
) that the model error terms
073905-15