Model-Based Estimation of Vortex Shedding in Unsteady Cylinder Wakes

TL;DR

This paper develops a harmonic decomposition-based model and two estimation methods, LE and LTE, to accurately estimate vortex shedding in cylinder wakes from single sensor data at different Reynolds numbers.

Contribution

It introduces LTE, a novel nonlinear estimation method that outperforms LE in reconstructing vortex shedding dynamics from limited sensor data.

Findings

LTE outperforms LE in flow field reconstruction.

Higher harmonic motions are subordinate to the fundamental frequency.

The model effectively captures vortex shedding dynamics at specified Reynolds numbers.

Abstract

This paper considers single-sensor estimation of vortex shedding in cylinder wakes at $Re=100$ in simulations and at $Re=1036$ in experiments. A model based on harmonic decomposition is developed to capture the periodic dynamics of vortex shedding. Two model-based methods are proposed to estimate time-resolved flow fields. First, Linear Estimation (LE) which implements a Kalman Filter to estimate the flow. Second, Linear-Trigonometric Estimation (LTE), which utilizes the same Kalman Filter together with a nonlinear relationship between harmonics of the vortex shedding frequency. LTE shows good performance and outperforms LE regarding the reconstruction of vortex shedding. Physically this suggests that, at the Reynolds numbers considered, the higher harmonic motions in the cylinder wake are slave to the fundamental frequency.

Figures

Equations

$$u_i$$

$$U$$

$$\tilde{U}$$

$$\begin{array}{rl}U& =\left[\begin{array}{ccc}{u}_1(t_1)& \dots & u_1(t_M)\\ ⋮& \ddots & ⋮\\ u_I(t_1)& \dots & u_I(t_M)\\ v_1(t_1)& \dots & v_1(t_M)\\ ⋮& \ddots & ⋮\\ v_I(t_1)& \dots & v_I(t_M)\end{array}\right],\\ \Phi & =\left[\begin{array}{ccc}{c}_{1,1}& \dots & c_{1,(2N)}\\ ⋮& \ddots & ⋮\\ c_{2I,1}& \dots & c_{2I,(2N)}\end{array}\right]\\ & =\left[\begin{array}{ccc}{\varphi }_1& \dots & {\varphi }_{2N}\end{array}\right],\end{array}\begin{array}{rl}\tilde{U}& =\left[\begin{array}{ccc}{\tilde{u}}_1(t_1)& \dots & {\tilde{u}}_1(t_M)\\ ⋮& \ddots & ⋮\\ {\tilde{u}}_I(t_1)& \dots & {\tilde{u}}_I(t_M)\\ {\tilde{v}}_1(t_1)& \dots & {\tilde{v}}_1(t_M)\\ ⋮& \ddots & ⋮\\ {\tilde{v}}_I(t_1)& \dots & {\tilde{v}}_I(t_M)\end{array}\right],\\ X& =\left[\begin{array}{ccc}\cos (\omega t_1)& \dots & \cos (\omega t_M)\\ \sin (\omega t_1)& \dots & \sin (\omega t_M)\\ ⋮& \ddots & ⋮\\ \cos (N\omega t_1)& \dots & \cos (N\omega t_M)\\ \sin (N\omega t_1)& \dots & \sin (N\omega t_M)\end{array}\right]\\ & =\left[\begin{array}{ccc}x(t_1)& \dots & x(t_M)\end{array}\right],\end{array}$$

$$\hat{\Phi }$$

$$\dot{x}(t)$$

$$A=\left[\begin{array}{ccc}{A}_{{\omega }_1}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & A_{{\omega }_N}\end{array}\right],A_{{\omega }_n}=\left[\begin{array}{cc}0& -n\omega \\ n\omega & 0\end{array}\right],$$

$$\dot{x}(t)$$

$$\tilde{U}(t)$$

$$\dot{x}(t)$$

$$\tilde{U}(t)$$

$$u_m(t)$$

$$\dot{\hat{x}}(t)$$

$$\hat{\tilde{U}}(t)$$

$$\left[\begin{array}{c}{x}_{2n-1}\\ x_{2n}\end{array}\right]$$

$$X(t)=({\hat{\Phi }}^T\hat{\Phi }{)}^{-1}{\hat{\Phi }}^T\tilde{U}.$$

$$\gamma (f)={\int }_y{\int }_x||\mathcal{F}\{\hat{\tilde{U}}(x,y,t)-\tilde{U}(x,y,t)\}|{|}_{2}dxdy,$$

Full text

I Introduction

The ability to accurately estimate flow fields is needed for a diverse range of fluid mechanics applications. These include flow prediction, drag reduction, and separation control. Generally, flow estimation aims to reconstruct the spatial and temporal evolution of the flow of interest, using limited time-resolved measurements. Flows are commonly estimated using statistical tools or dynamic models. Perhaps the most widely used statistical tool in flow estimation is Linear Stochastic Estimation (LSE) [1, 2]. The popularity of LSE has led to a number of modifications, including POD-LSE [3], Spectral LSE [4, 5] and Multi-Time-Delay LSE [6].

Alternatively one can estimate flows using dynamic models. Model-based estimation requires a dynamic model to describe the evolution of the flow with time, and an estimator to estimate the flow using limited measurements. Generally, system modelling can be either physics-driven or data-driven. Physics-driven approaches derive models based on the Navier-Stokes equations, including the Galerkin-POD method [7] and a recently developed resolvent-based method [8, 9, 10]. Alternatively, models can be built via data-driven approaches, including Dynamic Mode Decomposition (DMD) [11], Artificial Neural Networks (ANN) [12], Sparse Identification of Nonlinear Dynamics (SINDy) [13], and the Eigensystem Realization Algorithm (ERA) [14, 15].

These various estimation methods have been applied across a broad range of flows, including cavity flows [16, 17, 9]; wall-bounded flows [10]; jets [18]; and the flow over a flat plate [19], an airfoil [20] and a cylinder [21, 13]. This paper focuses particularly on the estimation problem in cylinder wakes, for which the flow first becomes unsteady (i.e. vortex shedding appears) at low Reynolds numbers ($Re>49$) [22]. Previous work has mostly focused on the estimation of the cylinder wake in the laminar vortex shedding regime ($49<Re<188.5$ [23]), and has done so mostly in simulations [21, 13]. The aim of this work is to develop estimation methods that i) can be applied in both simulations and experiments; ii) can be applied at higher Reynolds numbers where turbulent structures also appear.

The paper is organized as follows. Section II introduces the computational cylinder flow at $Re=100$ and the experimental cylinder flow at $Re=1036$. The harmonic-based system modelling procedure is introduced in Section III. The two estimation methods, Linear Estimation and Linear-Trigonometric Estimation, are first introduced in Section IV, and then applied to the $Re=100$ case and the $Re=1036$ case in Section V. Conclusions are drawn in Section VI.

II Computational and Experimental Setup

This section introduces the two setups considered in this paper: the unsteady cylinder wake at $Re=100$ in simulations and the unsteady cylinder wake at $Re=1036$ in an experiment.

II.1 Computational Setup

The cylinder wake at a Reynolds number of $Re=100$ is solved using Direct Numerical Simulation (DNS) in two dimensions using FEniCS [24]. The simulation follows the numerical scheme used in [25]. The Reynolds number is defined as $Re=\frac{UD}{\nu}$, where $U$ is the freestream velocity, $D$ is the cylinder diameter and $\nu$ is the kinematic viscosity. The spatial discretization is performed using a Taylor-Hood finite element scheme. A first-order difference scheme is used for time integration. The final mesh contains $5.46\times 10^{4}$ triangles, with a smallest radial increment of 0.01 near the cylinder. Once the flow reaches the limit-cycle state, 5400 instantaneous velocity fields are recorded (with a time step between them of $\Delta t=0.01s$). Together these constitute approximately nine shedding cycles. These are divided into a modelling set (600 time steps) for system modelling and a reference set (4800 time steps) for estimation.

II.2 Experimental Setup

The cylinder wake at a Reynolds number of $Re=1036$ is conducted in the water tunnel at the University of Melbourne. A schematic of the experimental setup is shown in figure 1. The freestream velocity is $U_{f}=41.5$ mm/s. The tunnel has an open top, from which the cylinder is suspended in water. The test cylinder is a sealed aluminum pipe with outer diameter $D=25$ mm and a thickness of $3$ mm. The total length of the cylinder is $495$ mm giving an aspect ratio of $19.8$.

The velocity fields are obtained using a Time-Resolved Particle Image Velocimetry (TR-PIV) setup. Flow snapshots are recorded at 50Hz for 6001 time steps, with a spatial resolution of $1920\times 1080$ pixels. This frame rate ensures that at least 50 snapshots are available for one vortex shedding period. A 60-watt continuous $532$nm laser source is placed at the bottom of the tunnel, providing a laser sheet with $3$ mm thickness perpendicular to the cylinder axis. The tracer particles used are silver-coated hollow glass spheres ($10$ $\mu$m in diameter) from DANTEC DYNAMICS. The field of interest (shown in Figure 1) covers the vortex formation region and a downstream region containing approximately two vortices simultaneously. An image background subtraction method [26] is used to remove the background patterns consistently appearing in the raw images. Abnormal measurements are filtered out using a $3\times 3$ median number filter. To obtain the velocity fields, a multi-grid interrogation method is used. In the first stage, the window size is $64\times 64$ with an overlap of $0.5$. In the second stage, a smaller window size of $32\times 32$ pixels is used with the same overlap. Finally, the 6000-time-step velocity fields are divided into two equally sized sets for system modelling and estimation.

III Harmonic-Based System Modelling

III.1 Harmonics in cylinder flows

Periodicity is a key feature of vortex shedding in cylinder flows [22]. The normalized vortex shedding frequency, known as the Strouhal Number, has been widely used to characterize the periodic features of vortex shedding [27, 28, 29]. However, the vortex shedding frequency is not the whole story, particularly when it comes to flow reconstruction and modelling. Figure 2(a-c) show the temporal signals of the streamwise velocity fluctuations at different locations in the cylinder wake at $Re=100$. The velocity signals all repeat themselves at the same frequency (vortex shedding frequency), but they have different shapes. These signals are therefore periodic but not sinusoidal. Therefore, from a flow-reconstruction point of view, knowing only the dominant frequency might not be sufficient.

The key to characterizing the shape of a non-sinusoidal periodic signal is its higher harmonics. Figure 2(d-f) show the velocity signals in the frequency domain. Clear peaks at the harmonics of the vortex shedding frequency can be found. Intuitively, the fundamental frequency (vortex shedding frequency) determines the period of a velocity signal, and its harmonics determine the signal shape within one period. Therefore, the vortex shedding frequency and its harmonics all need to be considered when reconstructing the vortex shedding patterns.

Gathering information about the higher harmonics is difficult because velocity signals are dominated by the vortex shedding frequency. Measurement noise in experimental cases makes it even more difficult to extract higher harmonics. Therefore in this work, a phase averaging method, based on Proper Orthogonal Decomposition (POD), is used to filter out the energy of sensor noise and any aperiodic structures. POD captures the most energetic structures in the flow, which makes it suitable for identifying the vortex shedding phase in cylinder flows [30]. The POD-based method proceeds in two parts. First, we compute the phases (in time) of the two most energetic POD modes. Second, we phase-average the flow using these phases. The phase-averaged flow is then used to build the low-order model of the flow in the following subsections.

III.2 Harmonic Decomposition

The strong harmonic peaks suggest the possibility of finding a low-order representation of the fluid system based on these harmonics. Here we use Harmonic Decomposition, a decomposition method used in Harmonic Balance methods[31], to decompose the flow based on the vortex shedding frequency and its harmonics. Using the same notation introduced in [32], the phase-averaged velocity field is decomposed into two parts: i) a time-independent mean flow, and ii) a sum of harmonic modes. This decomposition of the velocity field gives:

[TABLE]

where $\omega$ is the known fundamental harmonic (i.e. vortex shedding frequency). The terms $c_{i,2n-1}$ and $c_{i,2n}$ contain the time-independent coefficients of the $nth$ harmonic at grid point $i$; and $N$ is the number of harmonics retained. To further extend (1) to the entire domain, one can represent the evolution of the velocity fluctuations of the entire field in terms of a finite number of harmonics and their corresponding modeshapes $\Phi$, given as:

[TABLE]

where $U$ contains the two-dimensional velocity fields as column vectors. The term $\bar{U}$ corresponds to the mean velocity field, and $\tilde{U}$ corresponds to the fluctuations about the mean. The term $\Phi$ represents row vectors of modeshapes containing the coefficients $c$. The term $X$ contains all harmonics. These terms are given as:

[TABLE]

where $I$ is the number of grid points and $M$ is the number of time steps. $u_{i}(t_{m})$ and $v_{i}(t_{m})$ represent the streamwise and transverse velocity at time $t=t_{m}$, respectively. Each harmonic has two modeshapes ($\phi_{2n-1}$ and $\phi_{2n}$), corresponding to $\cos(n\omega t)$ and $\sin(n\omega t)$, respectively. $x(t_{m})$ contains the time dynamics of all modes at time $t=t_{m}$.

A more common method to find modeshapes based on frequencies is Dynamic Mode Decomposition (DMD) [11]. DMD finds active frequencies in the flow and generates corresponding modeshapes. DMD presents two challenges for the estimation problem of interest in this work. First, DMD can give rise to multiple modes with similar temporal frequencies—particularly when the flow exhibits turbulent features. Second, it can be challenging to rank DMD modes based on their relative importance. For these reasons we choose to use harmonic decomposition in this work, since only one mode is obtained for each harmonic, and ranking the modes is straightforward.

To summarize, harmonic decomposition decomposes the evolution of the velocity fields $U(x,y,t)$ into time-independent modeshapes $\Phi(x,y)$ and a time-dependent dynamic term $X(t)$. It significantly reduces the number of time-dependent variables from twice the number of grid points to twice the number of harmonics chosen. Less time-dependent variables allow a reduced-order model to be formed in a data-driven way in the following subsection.

III.3 Data-Driven Model Building

The model building procedure consists of two parts: i) finding modeshapes; and ii) forming a dynamic model.

Finding modeshapes: A least-square-error method is used to fit the harmonic modeshapes $\Phi$. The method requires knowledge of the phase-averaged velocity fields and the vortex shedding frequency, which can come either from simulations or from experimental measurements. The method gives

[TABLE]

where $\hat{\Phi}$ are the fitted modeshapes and $\tilde{U}_{pa}$ are the velocity fluctuations of the phase-averaged flow. Direct measurements of the velocity field $\tilde{U}$ may be used to replace $\tilde{U}_{pa}$ in (4) if measurement noise is negligible and any turbulence is weak.

Forming a dynamic model: The term $X(t)$, which represents sinusoidal signals with multiple temporal frequencies, is modelled as a marginally stable linear oscillator with multiple harmonic frequencies. The model is of the form

[TABLE]

where $x$ is the state of the model. The state matrix $A$ is

[TABLE]

where $\omega_{n}$ indicates the $n$th harmonic, and $N$ is the number of harmonics considered in the model.

State-space model: Combining the fitted modeshapes $\hat{\Phi}$ and the linear model of the time dynamics $x$, a state-space model of the velocity fluctuation is formed,

[TABLE]

where $\tilde{U}(t)$ represents one column of $\tilde{U}$ at time $t$. The modeshape matrix $\hat{\Phi}$ is regarded as an output matrix in the model. The state-space model is then used to design an estimator in the following section.

IV Estimation Methods

In this section we introduce two model-based estimation methods: Linear Estimation (LE) and Linear-Trigonometric Estimation (LTE). The model described in Section III is used for both methods.

IV.1 Linear Estimation

The goal of Linear Estimation is to estimate the entire velocity field based on limited velocity sensor measurements (one sensor in the current case). LE includes two steps: i) augment the linear model to include system disturbances, sensor measurement and measurement noise and ii) design a dynamic estimator based on the augmented linear model. Once the dynamic estimator is obtained, the velocity fields can be estimated by feeding the sensor measurement to the dynamic estimator.

In the first step, we augment the model described in Section III (Equation (7)). The augmented model now takes the form:

[TABLE]

where $u_{m}(t)$ is the single sensor measurement, and $C_{m}$ is the corresponding output matrix, which contains a certain row of the modeshape matrix $\hat{\Phi}$. This row contains the values of the modeshapes for all the harmonics at the sensor location. The term $w(t)$ represents disturbances to the system, and $s(t)$ represents measurement noise.

In the second step, we use the augmented model (Equation (8)) to design a Kalman Filter [33], which is an optimal linear dynamic estimator in the sense of minimizing the amplification of the estimation error from unknown disturbances $w(t)$ and measurement noise $n(t)$. The Kalman Filter uses knowledge of the sensor measurement $u_{m}(t)$ to estimate the full state and then the entire velocity field:

[TABLE]

where $\hat{x}(t)$ is the estimate of the full state $x(t)$; $\hat{\tilde{U}}(t)$ is the estimate of the velocity field $\tilde{U}(t)$; and $\hat{u}_{m}(t)$ is the estimate of the sensor measurement. The Kalman gain $L$ is obtained by solving a Riccati equation, which ensures that the error $e(t)=\hat{x}(t)-x(t)$ converges to zero.

Essentially, LE builds linear transfer functions from the sensor measurement to both the fundamental harmonic and higher harmonics. However, it has been pointed out in [34] that linear mechanisms are not dominant at higher harmonics in cylinder flows. Therefore, two questions may be asked. First, can LE estimate the higher harmonics, for which linear mechanisms are not dominant? We will discuss this in more detail in Section V. Second, can we estimate the higher harmonics nonlinearly? To answer the second question, we will introduce Linear-Trigonometric Estimation in the following.

IV.2 Linear-Trigonometric Estimation

The idea behind Linear-Trigonometric Estimation is to include the nonlinear trigonometric relationships among harmonics into Linear Estimation. LTE uses the same method as LE to linearly estimate the dynamics of the fundamental harmonic. The difference is that LTE then nonlinearly estimates the time dynamics of the higher harmonics from the estimate of the fundamental harmonic. We can do this because the term $x(t)$ only contains disturbed sinusoidal signals at different harmonics, and all the information about their phases and amplitudes are averaged and stored in modeshapes $\Phi$. By introducing the nonlinear operator $\mathscr{N}$, the output equation (Equation (9)) is augmented as,

[TABLE]

where $\hat{x}_{1}(t)$ and $\hat{x}_{2}(t)$ are the estimated states containing the cosine and sine of the fundamental harmonic respectively; and $N$ is the number of harmonics retained for estimation. The term $\mathscr{N}(\begin{bmatrix}\hat{x}_{1}(t)&\hat{x}_{2}(t)\end{bmatrix}^{T})$ contains the cosines and sines of the higher harmonics, where $\mathscr{N}$ is a recursive nonlinear operator based on Trigonometric Addition Formulas (TAF), given as,

[TABLE]

where $x_{2n-1}$ and $x_{2n}$ are the states containing the time dynamics of the $n$th harmonic, and $x_{2n-1}$ and $x_{2n}$ are the dynamics of the $n$th harmonic. Based on Equation (11), any harmonic can be expressed as a nonlinear function of the fundamental harmonic. A schematic is shown in Figure 3 to compare the work flows of both LE and LTE.

V Results and Discussions

The performance of LE and LTE are compared in two cases: the $Re=100$ case in DNS and the $Re=1036$ case in an experiment.

V.1 Cylinder flow at $Re=100$

For the $Re=100$ case, the Strouhal Number is $f=0.167$ (obtained from the power spectral density of the velocity), which is consistent with existing literature [29]. The transverse (vertical) velocity is chosen as the input signal because the magnitude of its fluctuation is greater than that of the streamwise velocity. This provides a higher signal-to-noise ratio. The sensor is placed at $x_{0}=2.90D,y_{0}=0.02D$ because this is the location where the perturbation energy is maximized.

First, we examine the reconstruction of the flow fields using the harmonic decomposition (Equation 3). LE estimates the time dynamic term $X(t)$ in Harmonic Decomposition, which means a perfect estimator can only give estimates as good as the reconstruction from Harmonic Decomposition. The reference set of the velocity fields is first decomposed into harmonic modeshapes to obtain the time dynamic term,

[TABLE]

Then the flow field is reconstructed based on Equation (3). Generally, the more harmonics included in $X(t)$, the higher percentage of energy can be captured. Figure 4 shows the percentage of energy captured as more harmonics are included. We observe that higher harmonics contain significantly less energy. Therefore, a third-order harmonic decomposition (first three frequencies) is used in both LE and LTE, capturing 99.9% of the total perturbation energy.

We now look at the performance of LE and LTE for the first three harmonics, comparing to the harmonic decomposition. Figure 5 compares the time dynamics term $X(t)$ decomposed from the DNS field (solid lines) to the estimate from LE and LTE (dashed line) for the first three harmonics. We observe that the estimates match the DNS curves after a short transient at the fundamental harmonic (Figure 5(a,b)) for both LE and LTE. This suggests that the dynamic estimator (Kalman filter) works well at the fundamental harmonic using the linear model and the measurements available.

However, clear discrepancies are observed for higher harmonics in LE (Figure 5(c,e)). One possible explanation for the poor performance of LE at higher harmonics is: The components of higher harmonics in the measurement signal are significantly weaker than the fundamental harmonic, making it difficult to obtain information about the higher harmonics directly from measurements.

Compared to LE, LTE shows better performance for the higher harmonics (Figure 5(d,f)). What differs is that LTE estimates the time dynamics of the higher harmonics using the fundamental harmonic and the nonlinear trigonometric relationship among harmonics. Therefore, good estimates from LTE indicate that the time dynamics of the higher harmonics are dependent on the fundamental harmonic. This suggests that the flow features related to the higher harmonics can be regarded as ”slave” to the fundamental harmonic. We will discuss this in more detail in Section V.4.

Now we compare the estimates of flow fields from LE and LTE to the harmonic decomposition and the DNS fields. Figure 6 (movies provided) compares instantaneous vorticity fields from DNS, flow reconstruction (3 harmonics), LE and LTE. Overall good agreement can be found between the reconstructed field (Figure 6(b)) and the DNS field (Figure 6(a)). This is consistent with the reconstructed flow capturing 99.9% of the perturbation energy. Comparing the estimates of LE (Figure 6(c)) and LTE (Figure 6(d)) to the DNS field (Figure 6(a)), two observations are made. First, the estimator reasonably estimates the streamwise positions of the vortices and their convective motion for both methods. This is consistent with the observation from Figure 5(a)(b), where both methods estimate the fundamental harmonic reasonably well. This is because the fundamental harmonic (vortex shedding frequency) determines the period of the flow patterns, which in turn determines the convective velocity of the vortex street. Second, clear discrepancies can be seen in the LE estimates regarding the shapes and the inclination angles of the vortices. This again is consistent with the previous observation of poor estimates of LE at higher harmonics (Figure 5(c)(e)). This is because, as discussed in Section III.1, higher harmonics determine the shapes of the velocity fluctuation within one period, which in large shapes the vortices.

To quantitatively compare the performance of the two estimation methods, we first introduce a measure of the estimation error $\gamma$:

[TABLE]

where $\mathscr{F}$ represents Fast Fourier Transformation and $||\cdot||_{2}$ represents the two norm of a complex number. $\hat{\tilde{U}}$ and $\tilde{U}$ are the estimate and the true value of the velocity fluctuation respectively. The quantity $\gamma(f)$ is the 2-norm of the spectrum of the estimation error integrated over the entire domain. It represents the distribution of the estimation error (integrated over the entire physical domain) in the frequency domain. Figure 7 compares $\gamma(\omega)$ for LE and LTE. Two observations are made. First, the two methods perform equally well at the fundamental harmonic, which suggests no extra error is introduced (for LTE) at the fundamental harmonic when we non-linearly estimate the time dynamics of the higher harmonics. Second, LTE performs significantly better than LE at the second and the third harmonics.

V.2 Cylinder flow at $Re=1036$

We now look at the $Re=1036$ experimental case. The Strouhal Number is measured as 0.21, which is consistent with the literature [29]. Following the same configuration as the $Re=100$ case, we use the PIV measurements at a single location as the sensor measurement, measuring the transverse velocity fluctuation at $x_{0}=2.8D,y_{0}=0$ , where the perturbation energy is the largest. Then we use LE and LTE to estimate the time dynamics of the harmonics, and further reconstruct the entire flow field.

Estimation in this case is more challenging than in the Re 100 case for three reasons. First, at this Reynolds number, periodic three-dimensional structures and aperiodic small turbulent structures appear [22]. This makes it more challenging to obtain the vortex shedding patterns from the raw PIV measurement. Second, disturbances and measurement noise introduce extra uncertainties. Third, the vortex shedding peak observed in the velocity spectrum becomes more broad-band with increasing Reynolds number [35], which means the vortex shedding frequency is prone to jitter at this higher Reynolds number. Therefore, it is more difficult to approximate the shedding mechanisms with a linear model. In addition to the harmonic decomposition, we also tried using Linear Stochastic Estimation and its variants (SLSE [5] and mtd-LSE [6]), but none gave a reliable estimate of the flow field.

We first look at reconstruction using Harmonic Decomposition. Figure 8 shows the energy captured in the reconstructed flow as the number of harmonics retained is increased. A clear drop of the gradient can be found at the third harmonic, which means harmonics higher than the third harmonic contain significantly less energy. Therefore, only the first three harmonics are considered in the estimation, capturing $72.4\%$ of the total energy. The effect of the truncation order will be discussed in more detail later.

We first compare the estimate of the time dynamics to harmonic decomposition (Equation 12), as shown in Figure 9. One observation is that the time dynamics computed directly from PIV measurement (thick solid lines) show both amplitude variation and phase shifting, which does not occur at $Re=100$. This indicates a disturbance to the linear model. Such disturbance could be the effect of turbulent structures. The other observation is that with good performance at the fundamental harmonic, LE fails at higher harmonics, whereas LTE still estimates the time dynamics at higher harmonics reasonably well. It again suggests strong nonlinear connections between the time dynamics of the fundamental harmonic and higher harmonics.

Now let us consider the instantaneous flow fields in physical space. Figure 10 shows instantaneous vorticity fluctuation fields from PIV, Harmonic Decomposition (3 harmonics), LE and LTE. It can be observed that Harmonic Decomposition (Figure 10(b)) captures the large-scale vortex street reasonably well, but small turbulent structures are not captured. This explains why the first three harmonics only capture $72.4\%$ of the energy for this Reynolds number whereas it captured $99.9\%$ in the $Re=100$ case. The position of the large vortical structures are estimated well by both LE (Figure 10(c)) and LTE (Figure 10(d)), but LTE better estimates both the overall shape and inclination angle of the structures. This is consistent with the observation in the Re 100 case. Furthermore, even though discrepancies exist compared to the PIV field, the LTE estimate agrees well with the field from Harmonic Decomposition. This suggests that the estimation error is mostly contributed by the system modelling procedure rather than by the estimation of that model. Better performance might therefore be achieved with improved system modelling.

A quantitative comparison of the estimator performance is shown in Figure 11. The lower $\gamma$ values at higher harmonics for LTE indicates better performance than LE at these harmonics, which is consistent with our previous observations.

V.3 Effect of the Truncation Order

Now we look at the influence of the number of harmonics on the estimation performance. Figure 12 compares the estimate from LTE using different numbers of harmonics in both the $Re=100$ case and the $Re=1036$ case. Figure 12c,d show the estimate of a instantaneous vorticity field using the fundamental harmonic alone in the two cases. The flow patterns are symmetric with respect to the centre line and only the convective motion of the vortical structures is captured. The vorticity fields reconstructed using the first three harmonics, shown in Figure 12e,f, are able to capture the asymmtric flow features, such as vortex tilting and stretching. This indicates that higher harmonics should be considered if one intends to estimate cylinder wakes with correct asymmetric shedding patterns. As discussed in Section V.1 and Section V.2, the first three harmonics are sufficient to reconstruct the shedding patterns reasonably well. Including more harmonics may only improve the estimation accuracy slightly if at all. Figure 12g,h together demonstrate that using the first 7 harmonics has an insignificant effect on the estimator performance.

V.4 Nonlinear interaction among harmonics

In both the $Re=100$ case and the $Re=1036$ case, LTE provide reasonable estimates of the higher harmonics using only the time dynamics of the fundamental harmonic. This indicates the higher harmonics are in some sense slave to the fundamental harmonic. This is consistent with the findings in a recent paper by [34]. They found that the mode at the second harmonic can be generated by forcing a Navier-Stokes-based resolvent model with the triadic interactions generated by the fundamental harmonic. The forced resolvent mode agrees with the DMD mode reasonably well at the same frequency. From a Navier-Stokes-based point of view, their results support our suggestion of higher harmonics being ”slave” to the fundamental harmonic. Another recent attempt of utilizing the nonlinear interaction among harmonics is by Meliga [36], who predicts the growth rates and dominant frequencies of cavity flows by considering the nonlinear interaction among harmonics as a feedback (via formation of Reynolds stresses) to a self-consistent model.

VI Conclusions

Single-sensor based estimation of the cylinder wake has been considered for two cases: in simulations at $Re=100$ and in experiments at $Re=1036$. The estimation is based on a harmonic-based model, which is first used to demonstrate the effects of the higher harmonics on the vortex street patterns. The model is then used in Linear Estimation, which estimates the dynamics of the fundamental harmonic reasonably well in both cases. By further including the nonlinear relationship between harmonics, Linear-Trigonometric Estimation estimates the dynamics of the higher harmonics based on the fundamental harmonic, and then reconstructs the flow evolution with reasonable accuracy in both cases. By examining the overall estimation error in the frequency domain, it is found that the improvement of LTE over LE occurs at the higher harmonics. The fact that LTE yields better performance—particularly for the higher harmonics—indicates that, for the flow under consideration, the dynamics of the higher harmonics are in some sense slave to the dynamics of the fundamental harmonic.

We expect reasonable performance for Linear-Trigonometric Estimation for cylinder wakes with Reynolds numbers higher than $1036$. This is because we expect LTE to perform well when the flow contains energetic structures that are periodic in time and coherent in space. Although small-scale turbulent structures become increasingly important in cylinder wakes as Reynolds number increases, it is also known that periodic structures are still significant and energetically important even at $Re=1.4\times 10^{5}$ [37]. Not only in cylinder wakes, LTE also has the potential to be applied in other types of flows in which the dominant flow structure is periodic, such as cavity flows.