Resolvent-based estimation of space-time flow statistics

TL;DR

This paper introduces a resolvent-based method to estimate space-time flow statistics from limited data, capturing the physics of turbulent flows for improved modeling and control.

Contribution

It develops a novel approach that infers nonlinear forcing statistics to predict unknown flow statistics without assuming low-rank structures, suitable for complex turbulent flows.

Findings

Successfully predicts flow statistics in Ginzburg-Landau model

Accurately estimates turbulent channel flow data

Demonstrates applicability to high-rank turbulent flows

Abstract

We develop a method to estimate space-time flow statistics from a limited set of known data. While previous work has focused on modeling spatial or temporal statistics independently, space-time statistics carry fundamental information about the physics and coherent motions of the flow and provide a starting point for low-order modeling and flow control efforts. The method is derived using a statistical interpretation of resolvent analysis. The central idea of our approach is to use known data to infer the statistics of the nonlinear terms that constitute a forcing on the linearized Navier-Stokes equations, which in turn imply values for the remaining unknown flow statistics through application of the resolvent operator. Rather than making an a priori rank-1 assumption, our method allows the known input data to select the most relevant portions of the resolvent operator for describing the data, making it well-suited for high-rank turbulent flows. We demonstrate the predictive capabilities of the method using two examples: the Ginzburg-Landau equation, which serves as a convenient model for a convectively unstable flow, and a turbulent channel flow at low Reynolds number.

Equations

$$\mathbf{y}=\mathcal{C}\mathbf{q},$$

$$\mathrm{given}\text{mathsfbi}{C}_{\mathbf{y}\mathbf{y}}(\mathbf{x},{\mathbf{x}}^{\prime },\tau )=$$

$$\mathrm{estimate}\text{mathsfbi}{C}_{\mathbf{q}\mathbf{q}}(\mathbf{x},{\mathbf{x}}^{\prime },\tau )=$$

$$\text{mathsfbi}S(\mathbf{x},{\mathbf{x}}^{\prime },\omega )=\underset{-\infty }{\overset{\infty }{\int }}\text{mathsfbi}C(\mathbf{x},{\mathbf{x}}^{\prime },\tau )e^{i\omega \tau }d\tau ,$$

$$\mathrm{given}\text{mathsfbi}{S}_{\mathbf{y}\mathbf{y}}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )=$$

$$\mathrm{estimate}\text{mathsfbi}{S}_{\mathbf{q}\mathbf{q}}(\mathbf{x},{\mathbf{x}}^{\prime },\omega )=$$

$$\mathcal{G}\frac{\partial \mathbf{q}}{\partial t}=\mathcal{F}\left(\mathbf{q}\right),$$

$$\mathbf{q}\left(\mathbf{x},t\right)=\bar{\mathbf{q}}\left(\mathbf{x}\right)+{\mathbf{q}}^{\prime }\left(\mathbf{x},t\right),$$

$$\mathcal{G}\frac{\partial {\mathbf{q}}^{\prime }}{\partial t}-\mathcal{A}\left(\bar{\mathbf{q}}\right){\mathbf{q}}^{\prime }=\mathbf{f}\left(\bar{\mathbf{q}},{\mathbf{q}}^{\prime }\right),$$

$$\mathcal{A}\left(\bar{\mathbf{q}}\right)=\frac{\partial \mathcal{F}}{\partial \mathbf{q}}\left(\bar{\mathbf{q}}\right)$$

$${\mathbf{y}}^{\prime }=\mathcal{C}{\mathbf{q}}^{\prime }.$$

$$\hat{\mathbf{y}}=$$

$$\hat{\mathbf{q}}=$$

$${\mathcal{R}}_{\mathbf{y}}(\mathbf{x},\omega )=$$

$${\mathcal{R}}_{\mathbf{q}}(\mathbf{x},\omega )=$$

$$\text{mathsfbi}{S}_{\mathbf{y}\mathbf{y}}=$$

$$\text{mathsfbi}{S}_{\mathbf{q}\mathbf{q}}=$$

$${\mathcal{R}}_{\mathbf{y}}=$$

$$=$$

$$\mathsfbi{S}_{\bm{f}\bm{f}}=\left[\bm{V}_{1}\,\,\,\bm{V}_{2}\right]\left[\begin{array}[]{cc}\bm{E}_{11}&\bm{E}_{12}\\ \bm{E}_{21}&\bm{E}_{22}\end{array}\right]\left[\bm{V}_{1}\,\,\,\bm{V}_{2}\right]^{*},$$

$$\text{mathsfbi}{S}_{\mathbf{y}\mathbf{y}}={\mathbf{U}}_{\mathbf{y}}{\mathbf{\Sigma }}_1{\mathbf{E}}_{11}{\mathbf{\Sigma }}_1{\mathbf{U}}_{\mathbf{y}}^{\ast }.$$

$${\mathbf{E}}_{11}={\mathbf{\Sigma }}_1^{-1}{\mathbf{U}}_{\mathbf{y}}^{\ast }\text{mathsfbi}{S}_{\mathbf{y}\mathbf{y}}{\mathbf{U}}_{\mathbf{y}}{\mathbf{\Sigma }}_1^{-1}.$$

$$\mathsfbi{S}_{\bm{f}\bm{f}}\approx\left[\bm{V}_{1}\,\,\,\bm{V}_{2}\right]\left[\begin{array}[]{cc}\bm{E}_{11}&0\\ 0&0\end{array}\right]\left[\bm{V}_{1}\,\,\,\bm{V}_{2}\right]^{*}=\bm{V}_{1}\bm{E}_{11}\bm{V}_{1}^{*}.$$

$$\text{mathsfbi}{S}_{\mathbf{q}\mathbf{q}}\approx {\mathcal{R}}_{\mathbf{q}}{\mathbf{V}}_1{\mathbf{E}}_{11}{\mathbf{V}}_1^{\ast }{\mathcal{R}}_{\mathbf{q}}^{\ast }.$$

$$\text{mathsfbi}C(\mathbf{x},{\mathbf{x}}^{\prime },\tau )=\underset{-\infty }{\overset{\infty }{\int }}\text{mathsfbi}S(\mathbf{x},{\mathbf{x}}^{\prime },\omega )e^{-i\omega \tau }d\omega ,$$

$$\text{mathsfbi}{S}_{\mathbf{q}\mathbf{q}}\approx \text{mathsfbi}{U}_{\mathbf{q}}{\mathbf{\Sigma }}_{\mathbf{q}}\text{mathsfbi}{S}_{\beta \beta }{\mathbf{\Sigma }}_{\mathbf{q}}\text{mathsfbi}{U}_{\mathbf{q}}^{\ast },$$

$$\text{mathsfbi}{S}_{\beta \beta }={\mathbf{V}}_{\mathbf{q}}^{\ast }{\mathbf{V}}_1{\mathbf{E}}_{11}{\mathbf{V}}_1^{\ast }{\mathbf{V}}_{\mathbf{q}}$$

$$\mathcal{A}=-\nu \frac{\partial }{\partial x}+\gamma \frac{{\partial }^2}{\partial x^2}+\mu (x).$$

$$\mu (x)=({\mu }_0-c_{\mu }^2)+\frac{{\mu }_2}{2}{x}^2$$

$$g(x,x^{\prime })=\frac{1}{\sqrt{2\pi }{\sigma }_f}\exp \left[-\frac{1}{2}{\left(\frac{x-x^{\prime }}{{\sigma }_f}\right)}^2\right]\exp \left[\mathrm{i}2\pi \frac{x-x^{\prime }}{{\lambda }_f}\right],$$

$$\text{mathsfbi}{C}_{\mathbf{q}\mathbf{q}}({\mathbf{x}}_1,{\mathbf{x}}_2,-\tau )=\text{mathsfbi}{C}_{\mathbf{q}\mathbf{q}}^{\ast }({\mathbf{x}}_2,{\mathbf{x}}_1,\tau ).$$

$$\text{mathsfbi}{S}_{\mathbf{q}\mathbf{q}}(y,y^{\prime };k_x,k_z,\omega )=\frac{1}{N_b}\sum _{j=1}^{N_b}{\hat{\mathbf{q}}}_j(y;k_x,k_z,\omega ){\hat{\mathbf{q}}}_j^{\ast }(y^{\prime };k_x,k_z,\omega ).$$

$$\frac{\partial \mathbf{u}}{\partial t}+\bar{\mathbf{u}}\cdot \nabla \mathbf{u}+\mathbf{u}\cdot \nabla \bar{\mathbf{u}}+\nabla p-\frac{1}{R{e}_{\tau }}\nabla \cdot \left[\frac{{\nu }_T}{\nu }\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right]=$$

Full text

1 Introduction

Practical limitations in both experiments and simulations can lead to partial knowledge of flow statistics. For example, an array of probes in an experiment provides information at a limited number of spatial locations and for a single flow quantity, e.g., velocity from hot-wires or pressure from microphones. Similarly, particle image velocimetry might provide velocity data, but not thermodynamic quantities, in a limited field of view. In simulations, one may wish to know flow statistics in a region that is not adequately resolved by the computational grid, such as unresolved near-wall regions or locations outside of the computational domain.

The ability to use available data to estimate statistics of flow quantities that are not directly accessible would be useful in each of these situations. For example, such a method could enable recovery of full-field statistics from a discrete set of measurements, estimates of the statistics of one variable from measurement of another, or estimates for a region outside of the field of measurement or computational domain.

Two methods for estimating unkown flow statistics from a limited set of known entries have recently been developed. Beneddine et al. (2016) proposed a method for estimating unknown power spectral densities (PSDs) using knowledge of the mean flow field and power spectra at a few locations. This is accomplished using a least-squares fit at each frequency between the known power spectra and the leading singular response mode obtained from the resolvent operator (McKeon & Sharma, 2010), which is derived from the linearized Navier-Stokes equations. This strategy explicitly assumes that the spectral content at frequencies of interest is dominated by the leading resolvent mode, and the method performs well when the matching points are located in regions where this hypothesis is valid. Specifically, excellent PSD estimates were obtained for the flow over a backward-facing step (Beneddine et al., 2016) and an initially laminar jet (Beneddine et al., 2017).

Zare et al. (2017) developed a method that uses arbitrary known entries in the spatial covariance tensor to estimate the remaining unknown entries. Their approach is also based on linearized flow equations and entails solving a convex optimization problem that determines a matrix controlling the structure and statistics of the associated nonlinear forcing terms. The optimization problem is subject to two constraints on the estimated covariance tensor: it must reproduce the known entries and obey a Lyapunov equation that relates the forcing and flow statistics. The constrained optimization problem is computationally demanding and requires a customized algorithm (Zare et al., 2015, 2017).

The objective of the present paper is to build on these previous methods to estimate unknown two-point space-time flow statistics. Both the PSDs (one-point temporal statistics) and spatial covariances (two-point spatial statistics) are subsets of two-point space-time correlations, so our approach represents a generalization of these previous methods. This is an important step since two-point space-time statistics contain additional, fundamental information about the flow. In particular, they carry information about coherent motions within the flow, and can even be used to define the concept of a coherent structure (Towne et al., 2018). Moreover, space-time correlations can be used to obtain real-time estimates of the flow state via convolution with a time-varying input signal (Sasaki et al., 2017).

The method developed in this paper borrows ideas from both of the previously mentioned methods. Like Beneddine et al. (2016), our method is built upon the resovent formalism of McKeon & Sharma (2010). The resolvent operator is derived from the Navier-Stokes equations linearized about the turbulent mean flow and constitutes a transfer function in the frequency domain between terms that are nonlinear and linear with respect to fluctuations to the mean. Whereas Beneddine et al. (2016) constructed their model using only the first singular mode of the resolvent operator (obtained via singular value decomposition), our model relaxes this a priori assumption and allows the known data to self-select the relevent portion of the resolvent operator. This makes the method more applicable to turbulent flows, in which the leading resolvent mode may account for only a modest fraction of the total flow energy (Schmidt et al., 2018), and allows us to extend the method to cross-spectra in addition to power spectra.

Our approach also follows the underlying strategy employed by Zare et al. (2017) of using the known data to infer the statistics of the unknown nonlinear terms that act as a forcing on the linearized equations. Their method assumes the nonlinear terms to have the same spatial correlation at all frequencies, in contrast to recent findings (Rosenberg et al., 2016; Towne et al., 2017). Our method relaxes this assumption and allows different spatial correlations for each frequency. Moreover, our frequency-domain formulation is algorithmically simple, requiring only basic linear algebra manipulations and avoiding Lyapunov equations and the need for external optimization routines.

The objective and capabilities of our method are fundamentally different from those of the classical method of linear stochastic estimation (Adrian, 1994; Bonnet et al., 1994) and related approaches that have recently been investigated (e.g., Encinar et al., 2018). In these methods, cross-correlations between input quantities and output quantities of interest must be known a priori and are used to estimate instantaneous values or conditional averages for the quantities of interest. In contrast, our method assumes no knowledge of the output statistics (or its cross correlation with input quantities), and instead uses input data alone to estimate space-time statistics of the output quantities of interest. Accordingly, our method could in fact be used to obtain an estimate of the statistics required to perform linear stochastic estimation.

The remainder of this paper is organized as follows. The method is derived and described in § 2 and demonstrated in § 3 using two examples: a simple model problem given by the Ginzburg-Landau equation and a turbulent channel flow. Finally, § 4 summarizes the paper and discusses further improvements and applications of the method.

2 Method

Our method for estimating space-time flow statistics from limited measurements is developed in this section. After precisely defining the objective, we develop our approach to the problem and provide some alternative interpretations of the method, which help to elucidate its properties.

2.1 Objective

Consider a state vector of flow variables $\bm{q}(\bm{x},t)$ that describe a flow, e.g., velocities and thermodynamic variables. The independent variables $\bm{x}$ and $t$ represent the spatial dimensions of the problem and time, respectively. Now suppose that the two-point space-time statistics are known for a reduced set of variables

[TABLE]

where the linear operator $\mathcal{C}(\bm{x})$ selects any desired subset or linear combination of $\bm{q}$. The problem objective can now be precisely stated in terms of two-point space-time correlation tensors:

[TABLE]

Here, $E\left\{\cdot\right\}$ is the expectation operator over time and the asterisk superscript indicates a Hermitian transpose.

Using the relationship between space-time correlation tensors and the cross-spectral density (CSD) tensors

[TABLE]

this objective can be equivalently stated in the frequency domain for statistically stationary flows:

[TABLE]

where $\hat{\bm{y}}(\bm{x},\omega)$ and $\hat{\bm{q}}(\bm{x},\omega)$ are the temporal Fourier transforms of $\bm{y}$ and $\bm{q}$, respectively, and the expectation is now taken over realizations of the flow (Bendat & Piersol, 1990).

2.2 Approach

Our approach to this problem relies on the resolvent operator obtained from the linearized flow equations and its connection with the remaining nonlinear terms (McKeon & Sharma, 2010). Begin with nonlinear flow equations of the form

[TABLE]

where ${\mathcal{G}}$ and $\bm{\mathcal{F}}$ are linear and nonlinear operators, respectively. Both compressible and incompressible Navier-Stokes equations can be cast in this form, and ${\mathcal{G}}$ is singular in the incompressible case to account for algebraic divergence-free condition. Alternatively, the incompressible equations can be written with a non-singular ${\mathcal{G}}$ by projecting into a divergence-free basis to eliminate the continuity equation (Meseguer & Trefethen, 2003). Additional transport equations can also be included.

Applying the Reynolds decomposition

[TABLE]

where $\bar{\bm{q}}\left(\bm{x}\right)$ is the mean (time-averaged) flow, to (2.5) and isolating the terms that are linear in $\bm{q}^{\prime}$ yields an equation of the form

[TABLE]

where

[TABLE]

is the linearized Navier-Stokes operator and $\bm{f}$ contains the remaining nonlinear terms. Similarly, (2.1) becomes

[TABLE]

In the frequency domain, (2.7) and (2.9) can be manipulated to give

[TABLE]

where

[TABLE]

are resolvent operators associated with $\hat{\bm{y}}$ and $\hat{\bm{q}}$, respectively.

Using (2.4) and (2.10), the CSD tensors can be written in terms of these resolvent operators as

[TABLE]

where $\mathsfbi{S}_{\bm{f}\bm{f}}(\bm{x},\bm{x}^{\prime},\omega)=E\{\hat{\bm{f}}(\bm{x},\omega)\hat{\bm{f}}^{*}(\bm{x}^{\prime},\omega)\}$ is the CSD tensor of the nonlinear term $\bm{f}$ (Semeraro et al., 2016; Towne et al., 2016, 2018). We emphasize that no approximation has been made to this point; (2.12) is an exact expression of the Navier-Stokes equations.

To obtain an approximation of the desired statistics $\mathsfbi{S}_{\bm{q}\bm{q}}$, we use the known statistics $\mathsfbi{S}_{\bm{y}\bm{y}}$ to estimate $\mathsfbi{S}_{\bm{f}\bm{f}}$. The salient question then becomes: how much can we learn about $\mathsfbi{S}_{\bm{f}\bm{f}}$ from $\mathsfbi{S}_{\bm{y}\bm{y}}$? An answer is provided by examining the singular value decomposition (SVD)

[TABLE]

The columns of the matrices $\bm{V}_{\bm{y}}$ and $\bm{U}_{\bm{y}}$ correspond to input and output modes that form orthonormal bases for $\hat{\bm{f}}$ and $\hat{\bm{y}}$, respectively. The rectangular matrix $\bm{\Sigma}_{\bm{y}}$ determines the gain of each of the input modes to the output. Since the rank of $\mathcal{R}_{\bm{y}}$ can be no greater than the number of entries in $\bm{y}$, i.e., the number of locations/quantities for which the statistics are known, many of the input modes have no impact on the output. Accordingly, the SVD can be written in the form of (2.13c), where the diagonal $\bm{\Sigma}_{1}$ contains the non-zero singular values and the blocks $\bm{V}_{1}$ and $\bm{V}_{2}$ contain input modes that have non-zero and zero gain, respectively. It is important to note that these resolvent modes are different from those usually studied, which are given by the SVD of $\mathcal{R}_{\bm{q}}$ (e.g., McKeon & Sharma, 2010; Schmidt et al., 2018).

The distinction between input modes that do or do not impact the output can be used to isolate the part of $\mathsfbi{S}_{\bm{f}\bm{f}}$ that can be educed from knowledge of $\mathsfbi{S}_{\bm{y}\bm{y}}$. Since $\bm{V}_{\bm{y}}$ provides a complete basis for $\hat{\bm{f}}$, $\mathsfbi{S}_{\bm{f}\bm{f}}$ can be expanded as

[TABLE]

where the matrices $\bm{E}_{ij}$ represent correlation between expansion coefficients associated with each input mode (see Towne et al., 2018). Inserting this expansion into (2.12a) and using (2.13c) to simplify the expression gives rise to the equation

[TABLE]

This means that only the part of $\mathsfbi{S}_{\bm{f}\bm{f}}$ associated with $\bm{E}_{11}$ impacts the observed statistics $\mathsfbi{S}_{\bm{y}\bm{y}}$; the remaining $\bm{E}_{ij}$ terms have no impact and are thus unobservable from these known data. Consequently, $\bm{E}_{11}$ contains all of the information about $\mathsfbi{S}_{\bm{f}\bm{f}}$ that can be inferred from $\mathsfbi{S}_{\bm{y}\bm{y}}$. Using the orthonormality of $\bm{U}_{\bm{y}}$, (2.15) gives

[TABLE]

The remaining terms $\bm{E}_{22}$ and $\bm{E}_{12}=\bm{E}_{21}^{*}$ (this equality is required to make $\mathsfbi{S}_{\bm{f}\bm{f}}$ Hermetian) can be arbitrarily chosen without impacting $\mathsfbi{S}_{\bm{y}\bm{y}}$, but these terms will impact $\mathsfbi{S}_{\bm{q}\bm{q}}$ and therefore must be modeled. The simplest choice, and the one used in the remainder of this paper, is to set these unknown terms to zero, leading to the approximation

[TABLE]

We show in Appendix A that this choice is identical to taking the least-squares approximation of $\mathsfbi{S}_{\bm{f}\bm{f}}$, which can be obtained by applying the pseudo-inverse of $\mathcal{R}_{y}$ and its complex conjugate to the left and right sides of (2.12a), respectively. Therefore, this approximation corresponds to choosing the smallest forcing (in an appropriate norm) that reproduces the known flow statistics.

Inserting (2.17) into (2.12b) gives the corresponding approximation of the desired flow statistics

[TABLE]

By construction, the known statistics used as input are exactly recovered, ensuring that the approximation converges in the limit of full knowledge of the flow statistics. Other approximations can be obtained by choosing the unknown $\bm{E}_{ij}$ terms differently; a few possibilities are discussed in § 4. The estimated space-time correlation tensor can be recovered from the estimated CSD via the inverse Fourier transform

[TABLE]

The method can also be understood in terms of a resolvent-mode expansion of $\mathsfbi{S}_{\bm{q}\bm{q}}$. The standard resolvent modes associated with the linearized flow equations are defined by the SVD $\mathcal{R}_{\bm{q}}=\,\bm{U}_{\bm{q}}\bm{\Sigma}_{\bm{q}}\bm{V}_{\bm{q}}^{*}$. Equation (2.18) can then be written as

[TABLE]

where

[TABLE]

is the CSD of the expansion coefficients in a resolvent-mode expansion of $\hat{\bm{q}}$ (Towne et al., 2018). In general, $\mathsfbi{S}_{\beta\beta}$ can project onto any of the resolvent output modes in $\mathsfbi{U}_{\bm{q}}$. Thus, the known statistics $\mathsfbi{S}_{\bm{y}\bm{y}}$, through their influence on $\bm{E}_{11}$, determine which resolvent modes participate in the estimate of $\mathsfbi{S}_{\bm{q}\bm{q}}$. This can be contrasted with the rank-1 model described earlier, in which only the leading resolvent mode is allowed to contribute.

3 Examples

In this section, our method is demonstrated and analyzed using two example problems: the complex Ginzburg-Landau equation and a turbulent channel flow.

3.1 Ginzburg-Landau equation

The Ginzburg-Landau equation has been used by several previous authors (e.g., Hunt & Crighton, 1991; Bagheri et al., 2009; Chen & Rowley, 2011; Towne et al., 2018) as a convenient one-dimensional model that mimics key properties of the linearized Navier-Stokes operator for real flows, such as a turbulent jet (Schmidt et al., 2018). The linearized operator takes the form

[TABLE]

Several variants of the function $\mu(x)$ have been used in the literature; here the quadratic form

[TABLE]

is adopted (Hunt & Crighton, 1991; Bagheri et al., 2009; Chen & Rowley, 2011). All of the parameters in equations (3.1) - (3.2) are set to the values used by Towne et al. (2018). With these parameters, the leading singular value of $\mathcal{R}_{\bm{q}}$ at its peak frequency is 10 times larger than the second singular value, which is a typical value for real flows. Following Bagheri et al. (2009), the equations are discretized with a pseudo-spectral approach using $N=220$ Hermite polynomials.

The discretized equations are stochastically excited in the time domain using forcing terms with prescribed statistics identical to those used by (Towne et al., 2018). In particular, the forcing is generated by convolving band-limited white noise with a kernel of the form

[TABLE]

where $\sigma_{f}$ is the standard deviation of the envelope and $\lambda_{f}$ is the wavelength of the filter. This leads to a forcing that is white-in-time up to the cut-off frequency but that has non-zero spatial correlation in the form of (3.3) but with $\sigma_{f}$ replaced with $\sqrt{2}\sigma_{f}$. This form of the forcing statistics is qualitatively similar to those of the nonlinear forcing terms in real flows, such as a turbulent jet (Towne et al., 2017). We use $\sigma_{f}=4$ and $\lambda_{f}=20$.

Although these forcing statistics are prescribed in this model problem and therefore known, this knowledge is not made available to the estimation procedure. The equations are integrated using a fourth-order embedded Runge-Kutta method (Shampine & Reichelt, 1997), and a total of $10000$ snapshots of the solution are collected with spacing $\Delta t=0.5$, leading to a Nyquist frequency of $\omega_{Nyquist}=2\pi$. The CSD of the solution is computed from these data using Welch’s (1967) method.

For the majority of the following analysis, $\bm{y}$ is defined to correspond to data obtained from three probes located at $x=-10$, 0 and 10. Other choices are considered in § 3.1.4.

3.1.1 Power spectra

The PSD is contained in the diagonal entries of $\mathsfbi{S}_{\bm{q}\bm{q}}$. The true power-spectral density for the Ginzburg-Landau model problem is shown as a function of $\omega$ and $x$ in figure 1(a). A single peak is observed at $\omega\approx-0.2$ and $x\approx 5$, and the amplitude remains above 1% of the peak over a range of about $-0.75<\omega<2$ and $-20<x<15$. The dashed lines show the $x$ locations where the data is taken as known, and the estimation procedure will attempt to reconstruct the PSD elsewhere.

The approximation of the PSD obtained using these three probes is shown in figure 1(b). By construction, the approximation is exact at the probe locations. The peak is well captured and the agreement is good in high-energy regions. In the lower-energy regions, the PSD is under-predicted away from the probe locations. This is a consequence of neglecting the undetermined portions of the forcing. It is likely that additional improvements could be obtained by modeling these undetermined portions of the forcing, as discussed in § 4.

3.1.2 Cross-spectra

The CSD estimates are evaluated next. Figure 2 compares the real part of the true and modeled CSD at eight frequencies, which are listed in the caption. The contour levels are the same for the true and estimated data at each frequency and range from the minimum to maximum values of the true CSD. The circles indicate the locations where the CSD is known and the remaining values are to be estimated. The first six frequencies (panels a-l) fall within the high-energy region observed in figure 1. In these cases, the estimates are accurate and track the length scales and shape of the CSD as a function of frequency. The final two frequencies (panels m-p) fall in low-energy regions. The basic trends in the length scales and shape are still captured, but the estimates are not quantitatively accurate.

3.1.3 Space-time correlations

The space-time correlation tensor $\mathsfbi{C}_{\bm{q}\bm{q}}$ can be recovered from the cross-spectral density $\mathsfbi{S}_{\bm{q}\bm{q}}$ using the inverse Fourier transform of (2.3). As an example, figure 3 shows the true (solid lines) and estimated (dashed lines) correlations as a function of time lag $\tau$ for three spatial locations, $x=-5$, 0 and 5. These locations correspond to a low-energy region, a probe position and the energy peak, respectively. Each curve has been scaled by the maximum value of the corresponding true correlation.

The one-point autocorrelation for each point is shown in figure 3(a-c). The amplitude of the autocorrelation for the low-energy point at $x=-5$ (panel a) is significantly under-predicted, but the correlation length scale is well captured. The estimated autocorrelation at $x=0$ (panel b) is exact since this point corresponds to one of the probe locations. The autocorrelation at $x=5$ (panel c) is accurately estimated apart from a small under-prediction of the peak, which corresponds to an under-prediction of the variance.

The cross-correlations between these three points are shown in figure 3(d-f). The estimates are quite good in all cases, including those involving the low-energy point (panels d and f) and two unknown points (panel f). It is interesting that the cross-correlations involving the low-energy point are more accurate than the autocorrelation at this point. The agreement for the cross-correlation between the known and high-energy points (panel e) is almost perfect.

The spatial distribution of the true and estimated cross-correlation tensors at fixed values of the time lag $\tau$ is shown in figure 4. The plotted time lag values range from $\tau=0$ to 10; negative values need not be considered due to the symmetry

[TABLE]

The contour levels are the same in each panel and range from zero to the maximum values of the true correlation at $\tau=0$. Again, the circles indicate the locations where the correlations are known, and the remaining values are to be estimated.

The spatial correlation tensor is obtained for $\tau=0$ and is shown in figure 4(a). As already observed in figure 3, the amplitudes of the correlations at zero time lag are slightly underpredicted, but the spatial shape and overall amplitude are well captured. As the time lag $\tau$ is increased, the estimates faithfully track the changing shape of the true correlations up to at least $\tau=10$, by which point the magnitudes of the correlations are small.

3.1.4 Impact of probe location and comparisons with the rank-1 model

In this section, comparisons are made between the new method described in this paper and the rank-1 method of Beneddine et al. (2016) that was discussed in § 1. Particular attention is given to the impact of the probe location(s) on the accuracy of the estimates provided by these two methods.

The PSD, which is the target quantity of the rank-1 model, is considered first. figure 5 compares the true PSD (top row) to the estimates from the rank-1 method (second row) and the new model (third row) for four different sets of probe locations (columns).

First, a single probe is places at $x=0$. In this case, the methods provide similar estimates, but the peak amplitude is slightly better predicted by the new method while higher frequencies are better captured by the rank-1 method. Adding two more probes at $x=\pm 10$ (second column) has little impact on the rank-1 estimate. In contrast, the new method is able to use this additional information to improve its estimate, particularly in regard to the shape of the moderately energetic region surrounding the peak and in the vicinity of the new probes.

Next, a single probe is placed at the low-energy location $x=-10$, well away from the peak. At this point, the underlying assumption of the rank-1 model – that the solution is dominated by the leading resolvent mode – is false. This is representative of the situation that will be encountered in real turbulent flows. Because of this, the rank-1 method leads to large over-predictions of the PSD. In contrast, the new method yields a moderate under-prediction of the PSD, which can be attributed to neglecting the unobservable portions of $\mathsfbi{S}_{\bm{f}\bm{f}}$.

In the final case (fourth column), the flow statistics are known in a continuous region, $x\leq 0$, rather than at an isolated set of points. These results can be compared to the first case in which statistics were known only at the boundary of this region, $x=0$. Including the additional data for $x<0$ leads to worse results for the rank-1 model (notice the significant over-prediction of the peak). This is an undesirable property; it means that poorly placed probes (where the rank-1 assumption is invalid) can obscure the information provided by well-placed probes (where the assumption is valid). More generally, this is a manifestation of the fact that the rank-1 method does not necessarily converge with increasing input information, even in the limit of complete knowledge of the flow statistics (Towne et al., 2018), in contrast to the new method. In the current example, it is clear that the additional information for $x<0$ does not degrade the estimate of the new method as it did for the rank-1 model, but instead leads to small improvements in the estimated PSD.

Next, comparisons are made between the the CSD estimates provided by the two methods for the same four sets of probes. While the rank-1 method was not specifically designed to estimate cross-spectra, its form nevertheless implies values for these two-point statistics, and their accuracy is important if the method is to be used for time-domain modeling, as in Beneddine et al. (2017). Comparisons are made for the frequency of maximum gain, $\omega=-0.6$, where the rank-1 model is expected to be most appropriate.

The overall conclusions regarding the CSD estimates are similar to those just discussed for the PSD. The two methods yield equivalent results for the single probe at $x=0$. The estimates for the new model are improved by adding two more probes at $x=\pm 10$, while they have little impact on the rank-1 results. Using a single probe at the low-energy point $x=-10$ leads to reasonable estimates using the new method but large errors for the rank-1 method. Finally, using data from the region $x\leq 0$ degrades the rank-1 estimates compared to using only the boundary point $x=0$ but improves the estimates obtained using the new model. In this case, the CSD estimates from the new model are indistinguishable from the true CSD, even though the known region does not contain the peak PSD.

3.2 Turbulent channel flow

3.2.1 Flow parameters, simulation and data processing

Next, we apply the new method to an incompressible turbulent channel flow at friction Reynolds number $Re_{\tau}=187$, defined in terms of the friction velocity $U_{\tau}$ and the kinematic viscosity $\nu$. Wall units, denoted by $+$ superscripts, are also defined in terms of $U_{\tau}$ and $\nu$. The flow is computed via direct numerical simulation (DNS) of the incompressible Navier-Stokes equations in a domain of size $x/h\times y/h\times z/h\in[0,2\pi]\times[0,2]\times[0,\pi]$, where $x$, $y$, and $z$ are the streamwise, wall-normal, and spanwise dimensions and $h$ is the channel half-width. The periodic directions $x$ and $z$ are discretized using 64 Fourier modes in each direction, and the wall-normal direction $y$ is discretized using 129 Chebyshev polynomials. The flow is driven by impossing a constant mass flux in the streamwise direction. The equations are advanced in time using a variable time step third-order Runge-Kutta integrator with a CFL number of 0.5. To facilitate post processing, the data is interpolated in time to 10000 evenly spaced time instances with $\Delta t^{+}=1.5$. The mean streamwise velocity is shown in figure 7(a).

The simulation data are used to compute the cross-spectral density tensor $\mathsfbi{S}_{\bm{q}\bm{q}}$, where $\bm{q}=\left[u,v,w\right]^{T}$ and $u$, $v$, and $w$ are the streamwise, wall-normal, and spanwise velocities, respectively. Since the flow is periodic in $x$ and $z$, the cross-spectral density is a function of wavenumber in these directions, i.e., $\mathsfbi{S}_{\bm{q}\bm{q}}=\mathsfbi{S}_{\bm{q}\bm{q}}(y,y^{\prime};k_{x},k_{z},\omega)$. The cross-spectral density is estimated using Welch’s (1967) method. The flow data are divided into overlapping blocks each containing $N_{fft}$ instantaneous snapshots of the flow. A discrete Fourier transform in $x$, $z$, and $t$ is applied to each block, leading to Fourier modes of the form $\hat{\bm{q}}_{j}(y;k_{x},k_{z},\omega)$ for $j=1,2,\dots,N_{b}$, where $N_{b}$ is the total number of blocks. Then, the cross-spectral density is estimated as

[TABLE]

Finally, the estimated cross-spectra are further averaged according to the symmetries described by Sirovich (1987), which ensures that the estimated cross-spectra are symmetric with respect to reflection across the channel center line and to $180$ degree rotation about the $x$-axis. We use blocks containing $N_{fft}=256$ instantaneous snapshots with $75\%$ overlap, leading to $N_{b}=156$ blocks, and we have verified that our results are insensitive to these choices.

3.2.2 Linearized Navier-Stokes equations

The resolvent operators required for the model are obtained from the incompressible Navier-Stokes equations

[TABLE]

where $\bm{u}=\left[u,v,w\right]^{T}$ is a vector of velocity disturbances, $\bar{\bm{u}}=\left[\bar{u},0,0\right]$ is the mean velocity, and $p$ is the pressure disturbance. Following previous work (Reynolds & Hussain, 1972; del Álamo & Jiménez, 2006; Illingworth et al., 2018), we have included an eddy viscosity model in the form of the total viscosity function $\nu_{T}(y)$. Details of our formulation are consistent with those of Illingworth et al. (2018) and can be found there.

Since the linearized equations are homogeneous in $x$ and $z$, we can apply Fourier transforms in these directions and obtain an equation for each $(k_{x},k_{z})$ wavenumber pair in the form of equation (2.7) with

[TABLE]

$\Gamma=\mathrm{diag}\left(\left[1,1,1,0\right]\right)$, and $\bm{q}=\left[u,v,w,p\right]^{T}$. The matrices in equation (3.7) are provided in Appendix B. The wall-normal direction $y$ is discretized using 201 Chebyshev polynomials, and no-slip boundary conditions are applied at the walls.

We choose the known quantity $\bm{y}$ to correspond to the three velocity components at $y/h=0.2$, which in inner units corresponds to $y^{+}=37$. This is the same $y/h$ value considered by Illingworth et al. (2018) in their recent Kalman filter study, although the $y^{+}$ value is different due to differing Reynolds numbers. This location is relevant to the application of LES wall modeling, in which one would use data along such a surface to approximate the near-wall flow and/or shear stress.

To visualize the results, we will focus on the velocity energy spectra, which are obtained from the cross-spectral density tensor as

[TABLE]

3.2.3 Root-mean-squared velocities

We begin by examining the root-mean-squared (RMS) velocity fluctuations, which are obtained by integrating $\mathsfbi{E}_{\bm{q}\bm{q}}(y;k_{x},k_{z},\omega)$ in $k_{x}$, $k_{z}$, and $\omega$ and taking the square root. The true RMS velocity fluctuations computed from the DNS data and those obtained from the model are compared in figure 7(b) as a function of $y^{+}$. The RMS values are accurately estimated for all three velocity components in the near-wall region, specifically for $y^{+}\lesssim 45$ ($y/h\lesssim 0.25$). The streamwise velocity estimates are especially accurate, while slightly larger discrepancies are observed for the wall normal velocity. Notably, the model accurately captures both the location and magnitude of the $u_{RMS}$ peak. For larger values of $y^{+}$, the RMS values quickly fall below the DNS values. Results in the following section show that this under prediction is due to missing energy at small scales which do not have a footprint at the probe location. This missing energy could potentially be recovered by appropriate modeling of the $\bm{E}_{ij}$ terms that have been set to zero, as discussed in § 4.

3.2.4 Energy spectra

Figures 8, 9, and 10 show the energy spectrum for each velocity component as a function of $y^{+}$ and $k_{x}^{+}$, $k_{z}^{+}$, and $\omega^{+}$, respectively. In each case, the energy has been integrated over the other two Fourier variables. The energies have been premultiplied by the appropriate wavenumber or frequency to account for the logarithmic axes. The contour levels are logarithmically spaced and span five orders of magnitude, with the highest level equal to the maximum value of the DNS streamwise velocity spectrum. The same levels are used in all subplots so that magnitudes can be directly compared. The true spectra computed from the DNS data appear in the top row of each figure, and the corresponding model estimates appear in the second row.

In all cases, the model accurately captures the energy distribution of all three velocity components for $y^{+}\lesssim 45$, except at the highest wavenumbers and frequencies. The amplitudes and locations of the energy peaks in $(y^{+},k_{x}^{+},k_{z}^{+},\omega^{+})$ space are captured by the model. On the other hand, the model under-predicts the energy at all wavenumbers and frequencies for higher values of $y^{+}$, which is consistent with the under-prediction of the RMS values observed in figure 7. The highest wavenumbers and frequencies are correctly predicted only near the position of the known input data at $y^{+}=37$ (horizontal dashed lines in the figures).

Figure 11 shows the energy of each velocity component as a function of the streamwise wavenumber $k_{x}$ and frequency $\omega$ at the wall-normal position $y/h=0.05$ ($y^{+}=9$). Here, we use linear axes so that the phase velocity $c_{p}=\omega/k_{x}$ can be easily visualized. The contour levels are logarithmically spaced between the maximum value of the streamwise spectrum from DNS and span five orders of magnitude. The same levels are used in all subplots.

Beginning with the DNS spectra shown in the top row, we see that the spectra are dominated by a band of energy that is approximately linear in $k_{x}-\omega$ space. The slope of this line corresponds to the phase velocity of the most energetic disturbances. At this wall-normal location, the dominant phase velocity is $c_{p}/U_{\tau}\approx 11$, which is shown as a dashed line in each subplot. The model accurately predicts the dominant phase velocity, and the main errors are observed primarily at phase velocities that are significantly different from the dominant one.

Figure 12 shows the $k_{x}-\omega$ spectra for the streamwise velocity at five different wall-normal locations: $y^{+}=2,9,19,56,94$ ($y/h=0.01,0.05,0.1,0.3,0.5$). It is clear that the model correctly captures the changes in phase velocity as a function of wall-normal position, even when the absolute energy levels are under predicted far from the wall.

3.2.5 Cross-spectra

In addition to the energy spectra considered so far, the model also provides predictions for cross-spectra. An example is shown in figure 13. The CSD is plotted as a function of $y^{+}$ and $\left(y^{+}\right)^{\prime}$ for the the streamwise and spanwise wavenumbers $\lambda^{+}_{x}=700$ and $\lambda^{+}_{z}=100$, respectively, and the phase speed $c_{p}=10$, which are typical values at which coherent structures are expected to appear (Sharma & McKeon, 2013). The model uses the input data indicated by the black circles, and accurately reproduces the CSD for all three velocity components.

3.2.6 Autocorrelations

Next, we consider the space-time correlations

[TABLE]

where the expectation is taken over all $x$, $z$, and $t$. These correlations can be recovered from the cross-spectra discussed so far by taking inverse Fourier transforms,

[TABLE]

We will focus on the autocorrelations

[TABLE]

As an example, we examine the autocorrelations as a function of the streamwise and temporal lag variables $\delta x$ and $\delta t$, respectively, at a fixed wall-normal location $y/h=0.05$ ($y^{+}=9$). Figure 14 shows the autocorrelation of each velocity component as a function of $\delta x$ and $\delta t$, i.e., the space-time autocorrelations along the streamwise direction. The contour levels are defined in the same way as in the previous three figures. The inverse slope of the band of high correlation in each plot provides a measure of the convection velocity of disturbances. At this wall-normal location, the convection velocity is approximately $11U_{\tau}$, which is consistent with the phase velocity shown in figure 11 as well as the observations of Kim & Hussain (1993). The convection velocity is accurately approximated by the model for all three velocity components. The correlation magnitudes are also well approximated aside from a moderate under prediction of the peak wall-normal velocity correlations.

4 Conclusions

Building on the work of Beneddine et al. (2016) and Zare et al. (2017), this paper introduces a method for estimating space-time flow statistics from a limited set of known values. The method is based on the resolvent methodology developed my McKeon & Sharma (2010) and the statisitical interpretation of this theory proposed by Towne et al. (2018). The central idea of our approach is to (partially) infer the non-linear term of the linearized Navier-Stokes equations using limited flow data. This non-linear term is then used as a forcing acting on the resolvent operator to reconstruct unknown statistics of the flow.

The method has been demonstrated using two examples problems. First, we applied it to the complex Ginzburg-Landau equation, which serves as a convient model of a convectively unstable flow. Using input data from three probe locations, the method provides good estimates of the unknown power spectra, cross-spectra, and space-time correlations within the energetic regions of $\omega$-$x$ or $\tau$-$x$ space. Comparisons are then made with the rank-1 model proposed by Beneddine et al. (2016). The two methods give similar results when the probes are placed at locations dominated by a single resolvent mode, but the new method gives superior results when the probes are placed at locations that violate this underlying assumption of the rank-1 model. The improved behavior in this case is important for turbulent flows, which cannot in general be described by a single resolvent mode. Furthermore, the estimates provided by the new method improve with the addition of more known input data.

Second, we applied our method to a turbulent channel at friction Reynolds number $Re_{\tau}=187$. Using data exclusively from the wall-normal location $y/h=0.2$ ($y^{+}=37$), the method provides good estimates of the velocity energy spectra and autocorrelations for $y\lesssim 0.25$ ($y^{+}\lesssim 45$). The energies and autocorrelations are under-predicted further away from the wall due to missing energy at small scales that requires additional modeling to capture. The success of the model in the near-wall region using knowledge of the interior flow suggests that it could be useful for designing new wall models for large-eddy simulation that are capable of capturing fluctuations of wall quantities such as shear stress and heat transfer and near-wall velocities that play an important role, for example, in particle laden flows.

Additional work is required to understand these observations, further assess the impact of the location of the known input data, and determine whether the results described in this paper will extend to higher Reynolds numbers and other types of turbulent flows. The properties and performance of the method should also be directly compared to other approaches that use the linearized flow equations as the basis for flow estimation, including the recent Kalman-filter-based approach described by Illingworth et al. (2018).

The method itself could also be further improved by modeling the portions of the forcing cross-spectral density that can not be observed using the known data. In the current formulation, these terms are simply set to zero, and there exist several possible alternatives. One is to assume that the unobserved forcing is uncorrelated with the observed part and with itself, leading to the approximation

[TABLE]

An appropriate value for the scalar amplitude $a$ could be determined from the amplitudes of the known $\mathsfbi{E}_{11}$ terms.

Another possibility is to choose the unobservable terms by insisting that the estimated $\mathsfbi{S}_{\bm{f}\bm{f}}$ projects exclusively onto the first $n$ singular modes of $\mathcal{R}_{q}$. This possibility is similar to a suggestion made by Beneddine et al. (2017), except here the expansion coefficients are treated as statistical quantities rather than complex scalars. As shown by Towne et al. (2018), this statistical treatment removes a fundamental accuracy restriction imposed by treating the expansion coefficients as deterministic scalars and allows for a convergent approximation.

Acknowledgments

A.T. thanks Eduardo Martini, André Cavalieri and Peter Jordan for fruitful discussion on this work. A.L.D. was partially supported by NASA under Grant NNX15AU93A and by ONR under Grant N00014-16-S-BA10.

Appendix A Least-squares approximation of forcing

Let us compute the approximation of $\bm{S}_{\bm{f}\bm{f}}$ we would obtain using the pseudo-inverse of $\mathcal{R}_{{y}}$, which can be written in terms of the SVD (2.13) as

[TABLE]

where $\bm{\Sigma}_{\bm{y}}^{-1}=\left[\bm{\Sigma}_{1}^{-1}\,\,\,\bm{0}\right]^{T}$. Applying the pseudo-inverse and its complex conjugate to the left- and right-hand-sides of (2.12a), respectively, gives the pseudo-inverse approximation of $\bm{S}_{\bm{f}\bm{f}}$:

[TABLE]

with $\bm{E}_{11}$ given by (2.16). Equation (A.2g) is identical to (2.17), which shows that setting the unknown $\bm{E}_{ij}$ terms to zero is equivalent to a least-squares, pseudo-inverse approximation of $\mathsfbi{S}_{\bm{f}\bm{f}}$.

Appendix B Linearized incompressible Navier-Stokes operators

The matrices defining the linearized incompressible Navier-Stokes equations in equations (2.7) and (3.7) are:

[TABLE]

[TABLE]

and $\mathsfbi{A}_{xx}=\mathsfbi{A}_{yy}=\mathsfbi{A}_{zz}=-\frac{1}{Re}\mathsfbi{\Gamma}$. We have defined $\nu_{T}^{\prime}=\frac{1}{Re}\frac{\partial\nu_{T}}{\partial y}$.