PHYSICAL REVIEW FLUIDS
8
, 064612 (2023)
Towards real-time reconstruction of velocity fluctuations
in turbulent channel flow
Rahul Arun
,
*
H. Jane Bae
, and Beverley J. McKeon
†
Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, California 91125, USA
(Received 18 January 2023; accepted 1 June 2023; published 27 June 2023)
We develop a framework for efficient streaming reconstructions of turbulent velocity
fluctuations from limited sensor measurements with the goal of enabling real-time ap-
plications. The reconstruction process is simplified by computing linear estimators using
flow statistics from an initial training period and evaluating their performance during a
subsequent testing period with data obtained from direct numerical simulation. We address
cases where (i) no, (ii) limited, and (iii) full-field training data are available using estimators
based on (i) resolvent modes, (ii) resolvent-based estimation, and (iii) spectral proper
orthogonal decomposition modes. During training, we introduce blockwise inversion to ac-
curately and efficiently compute the resolvent operator in an interpretable manner. During
testing, we enable efficient streaming reconstructions by using a temporal sliding discrete
Fourier transform to recursively update Fourier coefficients using incoming measurements.
We use this framework to reconstruct with minimal time delay the turbulent velocity
fluctuations in a minimal channel at Re
τ
≈
186 from sparse planar measurements. We
evaluate reconstruction accuracy in the context of the extent of data required and thereby
identify potential use cases for each estimator. The reconstructions capture large portions
of the dynamics from relatively few measurement planes when the linear estimators are
computed with sufficient fidelity. We also evaluate the efficiency of our reconstructions and
show that the present framework has the potential to help enable real-time reconstructions
of turbulent velocity fluctuations in an analogous experimental setting.
DOI:
10.1103/PhysRevFluids.8.064612
I. INTRODUCTION
A. Real-time estimation and control
Processing measurements generated by turbulent flows in real time is an increasingly important
problem in applications including flow estimation and control. As depicted in Fig.
1
, estimates
can be used to inform control schemes in real-time applications. Flow estimation problems may
be further subdivided into smoothing, filtering, and prediction problems. These problems aim to
inform estimates of past, present, and future states (respectively) using incoming measurements that
are usually incomplete and noisy. Real-time predictions can be useful for active control, but they
are typically more challenging than (smoothing or filtering) reconstructions of turbulent flows due
to multiscale dynamics and sensitivity to initial conditions.
Real-time flow estimation tasks are often practically limited by low-fidelity measurements and
task-related time constraints, e.g., in the context of aviation meteorology [
1
]. Measurements can be
improved by optimizing the design [
2
] and configuration [
3
] of sensors. Computational processing
*
rarun@caltech.edu
†
Present address: Department of Mechanical Engineering, Stanford University, Stanford, California 94305,
USA.
2469-990X/2023/8(6)/064612(42)
064612-1
©2023 American Physical Society
ARUN, BAE, AND McKEON
FIG. 1. High-level block diagram depicting a framework in which a real-time control scheme is informed
by real-time estimates generated by noisy measurements. We focus on the real-time estimation problem (red)
in the context of flow reconstruction.
can be expedited by leveraging data-driven and physics-based methods [
4
] that often reduce the
dimensionality of the problem or introduce a simplifying model. In the context of turbulent jets,
Sasaki
et al.
[
5
] demonstrated that the parabolized stability equations can be used to estimate (in
real time) the pressure at downstream sensors using upstream pressure measurements.
Beyond estimation, real-time control schemes require yet further empirically or heuristically
prescribed simplifications to ensure in-time actuation. Maia
et al.
[
6
] used a reactive control scheme
to experimentally attenuate centerline velocity fluctuations in a forced jet (Re
D
=
50 000) by
generating wave packets that destructively interfere with stochastic disturbances generated at the
nozzle. Abbassi
et al.
[
7
] demonstrated that actuating an array of cross-flowing jets that penetrate
into the log region of a turbulent boundary layer (Re
τ
=
14 400) can influence the resulting skin
friction drag and turbulence intensity.
Previous real-time flow estimation and control schemes demonstrate that, by combining gov-
erning equations with measurements and simplifying models, hybrid methods allow for a tailored
balance between the efficacy and the speed of the scheme. One commonality of these methods is
that their simplifying assumptions are often not well justified from first principles. In the present in-
vestigation, we focus on the problem of real-time flow reconstruction in the context of wall-bounded
turbulence. Even for this relatively simple estimation problem, the task of efficiently and accurately
representing turbulent dynamics from first principles [i.e., the Navier-Stokes equations (NSE)] using
minimal assumptions has remained challenging.
B. Reconstruction techniques
Reduced-order models (ROMs) often provide efficient flow representations to reduce the com-
plexity of high-dimensional turbulent flows, e.g., using proper orthogonal decomposition (POD)
[
8
] and dynamic mode decomposition (DMD) [
9
]. For statistically stationary flows, spectral POD
(SPOD) provides an efficient means of identifying coherent structures that retains a direct relation-
ship to resolvent analysis and DMD [
10
]. Unless otherwise stated, the SPOD modes we refer to
are those of the velocity fluctuations. Truncated SPOD mode estimation (TSME) often provides
an efficient means of flow reconstruction since SPOD modes form an optimal orthogonal basis in
terms of the variance captured by a given subset of modes [
11
]. Data-driven methods like SPOD
are often limited in that they neglect the governing equations and require extensive postprocessing.
However, Ghate
et al.
[
12
] showed that supplementing SPOD-based truncations with physics-based
enrichment using Gabor modes enables representations of flows over a broad range of scales.
Moreover, SPOD is amenable to a streaming formulation [
13
] and convolution-based strategies
for time-domain analysis [
14
], thereby reducing its computational requirements.
Equation-based frameworks have been developed in conjunction with data-driven methods to
efficiently estimate turbulent flows. Techniques from control theory are often used for flow esti-
mation in the time domain (e.g., Kalman filters [
15
–
17
]) and the frequency domain (e.g., Wiener
filters [
18
]). Techniques like linear stochastic estimation (LSE) [
19
] and its spectral variant [
20
]
are rooted in conditional estimation. LSE-based techniques are often used with data-driven [
21
,
22
]
064612-2
TOWARDS REAL-TIME RECONSTRUCTION OF VELOCITY ...
and equation-based [
23
,
24
] models to augment flow reconstructions. Other methods use simplified
governing models to produce forward and backward estimates that augment reconstructions from
low-resolution and multiresolution temporal measurements [
25
,
26
]. One central challenge in these
estimation techniques is addressing the nonlinearity of the governing equations.
Resolvent analysis [
27
,
28
] provides a powerful framework for addressing the nonlinearity of
the NSE and identifying energetic linear amplification mechanisms using minimal assumptions. In
a manner related to the work of Farrell and Ioannou [
29
], the turbulent fluctuation dynamics are
recast as a linear system forced by their nonlinearity in the resolvent framework. The resolvent
operator often admits a low-rank representation, and it can be constructed using a base flow profile,
which can be modeled or learned from data, e.g., via data assimilation [
30
,
31
].
One flavor of flow reconstruction using resolvent analysis is truncated response mode estimation
(TRME). TRME revolves around the singular value decomposition (SVD) of the resolvent operator,
which produces orthonormal bases for the forcing and response that are related through correspond-
ing gains or singular values. Modal truncation is performed based on the (often) low-rank nature of
the resolvent operator, which occurs when the gains associated with a small number of leading
modes dominate those of the remaining (suboptimal) modes. Moarref
et al.
[
32
] demonstrated
that streamwise energy amplification in each half of a turbulent channel (Re
τ
=
2003) is well
captured by a rank-1 approximation of the resolvent operator. Moarref
et al.
[
33
] then used convex
optimization techniques to compute resolvent mode weights, which encode dynamical information,
using a similar low-rank approximation to capture the velocity spectra in a turbulent channel. In
similar fashions, TRME was used to approximate energetic structures in flow around a cylinder
[
31
] and in lid-driven cavity flow (Re
=
1200) with a known (2D) mean flow [
34
]. In all of these
cases, the resolvent modes are weighted using spatially isolated velocity measurements in either
Fourier or physical space. In such setups, it is important that the measurements at least partially
capture the most energetic regions of the flow. For example, Symon [
30
] showed that tailoring probe
locations to capture energetic regions of both low-frequency wake modes and high-frequency shear
layer modes significantly improves reconstruction accuracy. Beneddine
et al.
[
35
] found a similar
result by separately capturing low-frequency and high-frequency regions in a backward-facing
step flow configuration. Recently, a quasisteady resolvent analysis has been used to reconstruct
high-frequency fluctuations using the phase of a lower-frequency periodic motion [
36
]. Resolvent
truncations have also recently been used to design
H
2
-optimal estimators and controllers [
37
] and
select optimal sensor and actuator locations [
38
] in cylinder flows.
While TRME provides an efficient and informative approximation, the assumption of a
fixed-rank [
39
] truncation may not be appropriate across all frequencies. In contrast to TRME,
resolvent-based estimation (RBE) revolves around the SVD of an input-output operator that mod-
ifies the resolvent operator to reflect sensor locations. Since flow sensors are usually limited and
sparsely spaced, RBE has shown promise in producing reconstructions of both flow statistics [
40
]
and dynamics [
41
] from limited data without truncating the resolvent operator
apriori
. RBE [
40
]
relates response measurement statistics to observable forcing statistics via an input-output operator
and (typically) models the unobservable forcing as uncorrelated with the observable forcing. The
estimated forcing statistics then yield estimates of the response statistics via the resolvent operator.
Building from this approach, Martini
et al.
[
41
] derived optimal, noncausal transfer functions
relating the measurements to the full response when the second-order forcing statistics and resolvent
operator are known. Further, Morra
et al.
[
42
] showed that the statistics in a turbulent channel can
be reasonably modeled using a low-rank SPOD approximation of the forcing. These studies reveal
that the tradeoff between accurately capturing the forcing color and retaining a sparse measurement
configuration can be somewhat alleviated by modeling the structure of the nonlinearity. For example,
in turbulent channel flows, including an eddy-viscosity model [
17
] can improve the accuracy of RBE
techniques for the response [
43
] and its statistics [
40
,
42
,
44
] with minimal input data.
We primarily focus on noncausal forms of linear estimators, including that derived for optimal
resolvent-based estimation (ORBE) [
41
,
43
]. More recently, a causal resolvent-based formulation
for estimation and control [
45
] was developed via the Wiener-Hopf formalism. However, the
064612-3
ARUN, BAE, AND McKEON
numerical machinery for forming this causal estimator is typically more complex than that of the
noncausal estimator due to an additional term in the estimation problem. By instead applying the
noncausal formulation to the real-time reconstruction problem, we avoid these complex methods and
retain a relatively simple framework. Noncausal estimation techniques (including ORBE) are not
typically amenable to real-time applications [
43
], but we circumvent this limitation by continuously
updating the estimated temporal Fourier coefficients over a small window at every time step. Further,
using this sliding window enables efficient recursive techniques for updating the coefficients [
46
].
The implementation of this technique, which has linear algorithmic complexity, differs from the
streaming Fourier sums technique of Schmidt and Towne [
13
], which has quadratic algorithmic
complexity but a smaller memory footprint.
Nonlinear estimation techniques, especially those incorporating machine learning techniques,
can also reconstruct flows with competitive accuracy. For example, Fukami
et al.
[
47
] reported
relative errors of less than 15% for spatial superresolution of a channel flow using a coarse input with
six elements across each half of the channel. Although they achieve significant data compression,
the coarsened input for estimation, obtained via averaging, is not conducive to typical experimental
measurements. By contrast, our investigation focuses on practical considerations for experimental
implementation, including time delays, measurement configurations, scalability, and uncertainty
quantification. Working in the linear estimation setting is further advantageous in that it is amenable
to a wealth of estimation and control theory for linear systems. Further, in contrast to many
machine learning techniques, the present linear estimation techniques are interpretable in terms
of their optimality and through the basis and flow statistics imposed by the linear model. Finally,
whereas machine learning techniques can take significantly longer (e.g., multiple days [
47
]) to train
estimators on optimized hardware, we introduce estimators that can be constructed from training
data within an hour on a laptop.
C. Contributions
The goal of this study is to construct a method for reconstruction of turbulent velocity fluctuations
from sparse measurements, adding to the vast flow estimation literature. For this, we develop a new
framework to enable efficient streaming capabilities with the potential for real-time implementation.
In Sec.
II
we detail the governing equations and formulate the statistical framework, based on the
generalized Wiener filter, that we use to reconstruct the fluctuations. We also detail the application
of methods inspired by TRME, ORBE, and TSME to this more general framework. In Sec.
III
,
we introduce efficient methods that enable streaming flow reconstructions with minimal time delay
and improved scalability compared to standard techniques. When forming linear estimators from
training data, we introduce blockwise inversion to efficiently compute the resolvent operator.
When reconstructing the flow from testing data, we apply a recursive sliding discrete Fourier
transform (SDFT) to efficiently update the inputs to the estimators in a streaming fashion. In
Sec.
IV
we evaluate the performance of the reconstructions using each method. We specifically
address (i) the fidelity of the training data and the testing data, (ii) the validity of the models
underlying the estimators, (iii) the potential use cases for each reconstruction method, and (iv)
the reconstruction efficiency in the context of real-time applications. In doing so, we demonstrate
that practical reconstruction speeds can be achieved using modest computational resources while
retaining competitive accuracies with other linear estimation techniques. Finally, we summarize
our reconstruction methods, the performance of the estimators, and promising future prospects and
applications in Sec.
V
.
II. BACKGROUND METHODOLOGY
A. Governing equations
We focus on turbulent channel flow between two parallel walls, for which we denote the stream-
wise, wall-normal, and spanwise directions by
x
,
y
, and
z
and the corresponding velocities by
U
,
V
,
064612-4
TOWARDS REAL-TIME RECONSTRUCTION OF VELOCITY ...
and
W
, respectively. This flow is governed by the nondimensional, incompressible Navier-Stokes
equations (NSE)
∂
t
U
+
(
U
·
∇
)
U
+
∇
P
=
1
Re
τ
∇
2
U
,
∇
·
U
=
0
,
(1)
where
U
=
[
U
,
V
,
W
]
T
is the velocity vector,
P
is pressure, and
t
denotes time. For channel
flow, Re
τ
=
u
τ
h
/ν
is defined in terms of the channel half-height
h
, the friction velocity
u
τ
, and
the kinematic viscosity
ν
. Lengths, velocities, and pressures are normalized by
h
,
u
τ
, and
ρ
u
2
τ
,
respectively, where
ρ
is the fluid density.
We combine the velocity and pressure into a single state variable,
Q
=
[
U
,
V
,
W
,
P
]
T
,
whose evolution is governed by the NSE. The state fluctuations are expressed by (Reynolds)
decomposing the state vector as
Q
=
Q
+
q
, where
Q
and
q
represent the base flow and the
fluctuations, respectively. Here we specify the base flow as the mean state profile,
Q
(
y
)
=
Q
(
y
)
=
[
U
(
y
)
,
V
(
y
)
,
W
(
y
)
,
P
(
y
)]
T
, where
V
(
y
)
=
W
(
y
)
=
0 for channel flow [
48
]. The corresponding fluc-
tuations are Fourier transformed in time and in the spatially homogeneous and periodic directions
(
x
and
z
here), giving
ˆ
q
(
y
,
k
)
=
ˆ
u
(
y
,
k
)
ˆ
p
(
y
,
k
)
=
F
{
q
(
x
,
t
)
}=
1
8
π
3
∞
−∞
q
(
x
,
t
)e
i(
ω
t
−
k
x
x
−
k
z
z
)
dt dxdz
,
(2)
where i
=
√
−
1,
x
=
[
x
,
y
,
z
]
T
and each triplet,
k
=
[
k
x
,
k
z
,ω
]
T
, contains the streamwise wave num-
bers, spanwise wave numbers, and (angular) temporal frequencies, respectively. The transformed
differential operators are given by
ˆ
∇
=
[i
k
x
,∂
y
,
i
k
z
]
T
and
ˆ
∇
2
=
∂
yy
−
k
2
, where
k
2
=
k
2
x
+
k
2
z
.
Applying the Reynolds decomposition and the Fourier transform defined in (
2
), (
1
) can be used
to express the state fluctuation dynamics for each triplet via a linear subsystem that is externally
forced by the nonlinear advection terms. Correspondingly, each linear subsystem may be written as
L
(
k
)
ˆ
q
(
y
,
k
)
=
L
B
(
k
)
ˆ
∇
ˆ
∇
T
0
ˆ
u
(
y
,
k
)
ˆ
p
(
y
,
k
)
=
ˆ
f
(
y
,
k
)
0
,
(3)
where
k
=
[0
,
0
,
0]
T
.Here
L
B
represents the “basic” terms in the unsteady momentum equations,
which are linear in the
velocity
fluctuations, and
ˆ
f
(
y
,
k
)
=
F
{−
(
u
·
∇
)
u
+
(
u
·
∇
)
u
}
represents the
nonlinear forcing. While the subsystem for each triplet is
externally
forced by velocity fluctuations
satisfying
k
=
k
, the system is
internally
forced by its nonlinearity. In Sec.
II B
we discretize
this state-space formulation of the Navier-Stokes system and express the corresponding resolvent
operator in discrete form.
B. Discretized linear system
When formed directly from the continuous governing equations, the linear subsystems in (
3
)
are infinite-dimensional. However, by discretizing the equations onto a grid, they assume a finite-
dimensional representation based on the spatial resolution in the wall-normal direction. As depicted
in Fig.
2
, the direct numerical simulation (DNS) data in the present case are staggered [
49
]ona
Cartesian grid. For a wall-normal domain discretized into
N
y
cells, the discrete representations of
ˆ
u
and
ˆ
q
are column vectors of sizes
N
u
=
3
N
y
+
1 and
N
q
=
N
u
+
N
y
, respectively. The numerical
quadrature (i.e., wall-normal integration) weights associated with center and edge quantities are
y
j
and 0
.
5(
y
j
−
1
+
y
j
), respectively.
Since state variables are Fourier transformed in
x
and
z
, we accurately reference the cell-center
values for ˆ
u
and ˆ
w
by shifting their Fourier transform origins by half a cell width in
x
and
z
,
respectively. To retain a direct analogy with the (discrete) governing equations used in simulating
the flow, we maintain ˆ
v
and the mean shear,
∂
y
U
, at the wall-normal edges,
y
e
j
and
y
e
j
+
1
, of each cell.
Correspondingly, we represent the mean shear using a bidiagonal matrix,
∂
y
U
, that acts to average
064612-5
ARUN, BAE, AND McKEON
FIG. 2. Staggered state variable locations on an arbitrary cell in the computational domain. This diagram
is not drawn to scale since
x
=
z
=
y
j
,and
y
j
varies in the wall-normal direction. Velocities reside on
the faces of the cell and pressure resides at the cell center. The superscripts
·
c
and
·
e
represent the (wall-normal)
centers and edges of the cell, respectively.
the contributions of ˆ
v
∂
y
U
at the cell edges to the cell centers. The diagonal matrices,
U
c
,
e
, represent
the mean profile at the cell centers and edges, respectively.
In discrete form, the linear operator,
L
, is a matrix,
L
∈
C
N
q
×
N
q
,oftheform
L
(
k
)
=
⎡
⎢
⎢
⎢
⎢
⎣
L
c
(
k
)
∂
y
U
0
i
k
∗
x
I
c
0
L
e
(
k
)
0
∂
y
00
L
c
(
k
)i
k
∗
z
I
c
i
k
∗
x
I
c
∂
y
i
k
∗
z
I
c
0
⎤
⎥
⎥
⎥
⎥
⎦
=
L
B
ˆ
∇
ˆ
∇
T
0
,
(4)
where the
0
matrices are conformable and
L
B
is the basic subblock, analogous to
L
B
in (
3
). The
diagonal subblocks take the form
L
c
,
e
=−
i
ω
I
c
,
e
+
i
k
∗
x
U
c
,
e
−
1
Re
τ
ˆ
∇
2
c
,
e
,
(5)
where
ˆ
∇
2
c
,
e
and
I
c
,
e
are the Laplacian and identity matrices of appropriate size. The discrete
representations of
∂
y
and
∂
yy
are computed via central finite differences with modified boundary
terms to enforce the relevant boundary conditions at the walls. In each case, boundary values can
be determined by applying central finite differences on a ghost grid that symmetrically extends the
center grid about the wall (edge) locations. Further,
k
∗
x
=
2
x
sin
k
x
x
2
and
k
∗
z
=
2
z
sin
k
z
z
2
(6)
are the modified wave numbers corresponding to central finite differences in the streamwise and
spanwise directions, respectively. The differences between
k
∗
x
and
k
x
induce differences in the
advection and Fourier phase speeds which grow as
|
k
x
|
increases.
Correspondingly, the spatially discretized form of (
3
) is given by
L
(
k
)
S
(
k
)
0
0
I
c
ˆ
q
(
k
)
=
BS
(
k
)
ˆ
f
(
k
)
,
(7)
where we have assumed that the streamwise and spanwise components of
ˆ
q
and
ˆ
f
are unshifted.
Here
S
is a diagonal matrix that shifts the origin for staggered flow quantities (in
x
and
z
) to the cell
064612-6