Experimental Detection of Organised Motion in Complex Flows with Modified Spectral Proper Orthogonal Decomposition

Abstract

Spectral proper orthogonal decomposition (SPOD) has seen renewed interest in recent years due to its unique ability to decouple organised motion at different timescales from large datasets with limited available information. This paper investigated the unsteady components of the flow field within a simplified turbine centre frame (TCF) model by applying SPOD to experimental, time-resolved flow speed data captured by particle image velocimetry (PIV). It was observed that conventional methods failed to capture the two significant active bands in the power spectrum predicted by preliminary hot wire anemometry measurements. Therefore, a modification to the SPOD procedure, which employs subsampling of the time sequence recorded in the experiment to artificially lower the PIV data acquisition frequency, was developed and successfully deployed to analyse the TCF flow field. The two dynamically active bands were identified in the power spectra, resulting in a closer match to the preceding analyses. Within these bands, SPOD’s ability to capture spatial coherence was leveraged to detect several plausible coherent, fluctuating structures in two perpendicular planes. A partial three-dimensional reconstruction of the flow phenomena suggested that both bands were associated with a distinct mode of organised motion, each contributing a significant percentage of the system’s total fluctuating energy.

1. Introduction

In recent years, the improved availability of complex hardware and computational capacity has made applying more demanding data-driven methods in experimental fluid mechanics more feasible. It is a competitive and rapidly progressing space in which a number of different methods vie for the limelight. The space-only formulation of proper orthogonal decomposition (POD) [1] has long remained the default option, benefiting from denser datasets due to improvements in high-speed imaging and simulation, which allowed expanding its domain of possible applications to more complex turbulent flows and heat transfer problems [2,3]. Other popular methods include dynamic mode decomposition (DMD) [4] and its derivatives [5], wavelet analysis [6,7], and independent component analysis (ICA) [8,9]. This work’s focus lied in spectral proper orthogonal decomposition (SPOD), which has recently enjoyed more attention due to its unique capability to optimally represent spatiotemporal coherence in flow data [10,11,12]. While the method is proven to be a valuable tool mainly in rather idealised scenarios, its applications as an investigative aid in exploring the dynamics of complex systems remain largely unexplored in the literature.

One such system, which, due to technical relevance, warrants closer examination, is the turbine centre frame (TCF)—a single module connecting the low- and high-pressure turbines of twin spool gas turbine engines, which forms the gas duct between the two sections, while also providing internal channels for hydraulics, the transfer of structural loads, and cooling air exchange between the hub and the exterior of the engine. Contemporary TCFs incorporate aggressive geometries to accommodate the rapidly diverging gas duct, complicating the design of the struts that house the service channels. Consequently, the flow of cooling air through the structure, which maintains safe operating temperatures of the TCF itself and is also supplied to the neighbouring turbine disks, is profoundly complex: large pressure gradients and narrow channels produce high local flow speed spikes, which result in a convoluted series of unsteady flow structures spanning the TCF’s multiple internal cavities.

This paper documented the use of SPOD to process flow data acquired in a high-speed particle image velocimetry (HSPIV) experiment on a simplified TCF model. The objective was two-fold: first, to assess the utility of SPOD as a tool to detect yet-unknown coherent structures in a complex fluid system, and second, if the first step proves a success, to obtain a result with the potential for immediate technical application. The dynamics analysis can, for instance, be used to gauge the performance of a numerical solver used to design a TCF interior in an engineering context.

Section 2 and Section 3 provide a brief overview of the theoretical framework of POD and the fluids’ experimental setup, respectively. Section 4 details the necessary preliminary steps to contextualise the modal decomposition results. This includes a discussion of the topology of the mean flow through the model, as well as insights into the fluid system’s dynamics obtained through hot wire anemometry. The results of an analysis of the experimental data with conventional SPOD are presented in Section 5; as the base, the SPOD algorithm was found to perform poorly with the available data. Section 6 proposes and validates a modification, which involves the use of subsampling to lower the effective data acquisition frequency of the input. The results obtained with the updated method are presented thereafter.

2. Spectral Proper Orthogonal Decomposition

Modal decomposition produces a set of deterministic modes that capture dynamically significant patterns in the energy distribution of stochastic data. One such technique, known as spectral proper orthogonal decomposition, was the primary analytic instrument employed in this work.

2.1. Background

POD is a family of decomposition algorithms that contains several distinct techniques. Its inception lies in the work of Lumley [13,14], who first proposed the use of the Karhunen–Loève (K-L) theorem as a tool for analysing turbulent flows. Though its initial formulation was deliberately agnostic towards the nature of the kernel, the field was long dominated by an exclusively spatial form of POD popularised by Sirovich [1], which produces spatially orthogonal modes and, therefore, only represents spatial coherence. Our focus lied in spectral proper orthogonal decomposition, which is geared specifically towards the optimal representation of spatiotemporal coherence. The concept was coined by Picard and Delville [15] and recently applied widely by Towne et al. [10], Schmidt et al. [11], Schmidt and Colonius [16], and others. Note that this algorithm is distinct from the procedure by the same name proposed by Sieber et al. [17].

The fundamental principle of POD lies in computing a set of deterministic functions $\mathit{\varphi }\left(\mathit{z}\right)$ that are the optimal representation of a zero-mean stochastic process $\left\{\mathit{q}\right(\mathit{z};\xi \left)\right\}$ with respect to the inner product $\langle ·,·\rangle$ of the Hilbert space $\mathcal{H}$ that each realisation of $\mathit{q}$ inhabits; $\mathit{z}\in \mathsf{\Omega }$ denotes a set of independent variables, whereas $\xi$ parameterises the probability space. The zero-mean property of $\mathit{q}$ can be expressed as $\mathrm{E}\left\{\mathit{q}\right(\mathit{z};\xi \left)\right\}=0$, with $\mathrm{E}\left\{·\right\}$ being the expectation operator within the probability space.

The desired qualities of $\mathit{\varphi }\left(\mathit{z}\right)\in \mathcal{H}$ can be expressed through the maximisation problem:

$$\lambda =\frac{\mathrm{E}\left\{{|\langle \mathit{q}(\mathit{z};\xi ),\mathit{\varphi }(\mathit{z})\rangle |}^2\right\}}{\langle \mathit{\varphi }\left(\mathit{z}\right),\mathit{\varphi }\left(\mathit{z}\right)\rangle }.$$

(1)

(1) can be intuitively contextualised as an effort to optimise $\mathit{\varphi }$ such that the expectation value of the projection of the stochastic process onto a deterministic modal function is maximal. The expression is normalised to decouple $\lambda$ from the magnitude of $\mathit{\varphi }$. The desired functions can be shown to be the eigenfunctions of the Fredholm integral equation:

$${\int }_{\mathsf{\Omega }}\mathbf{R}(\mathit{z},{\mathit{z}}^{\prime })\mathit{\varphi }\left({\mathit{z}}^{\prime }\right)\mathrm{d}{\mathit{z}}^{\prime }=\lambda \mathit{\varphi }\left(\mathit{z}\right)$$

(2)

whose kernel is the two-point correlation tensor across $\mathsf{\Omega }$:

$$\mathbf{R}(\mathit{z},{\mathit{z}}^{\prime })=\mathrm{E}\left\{\mathit{q}(\mathit{z};\xi ){\mathit{q}}^{*}({\mathit{z}}^{\prime };\xi )\right\}.$$

(3)

The K-L theorem states that (2) is solved by a infinite set of functions $\{{\mathit{\varphi }}_k,{\lambda }_k\},k\in \mathbb{N}$, which are mutually orthonormal (i.e., $\langle {\mathit{\varphi }}_j,{\mathit{\varphi }}_k\rangle ={\delta }_{jk}$) and can be ranked by the magnitude of their eigenvalue, ${\lambda }_1\ge {\lambda }_2\ge \dots \ge 0$. The eigenfunctions form an energy-optimal (proper) basis for each realisation of $\mathit{q}$.

Space-only POD, as typically employed in statistical fluid mechanics, defines $\mathit{q}(\mathit{x},t)$ as a finite set of fields of scalar (e.g., pressure, vorticity) or vector (e.g., velocity) properties, i.e., snapshots of a fluid flow. These can be understood as independent realisations of a stochastic process. The domain $\mathsf{\Omega }$ of the statistical ensemble, therefore, consists of a bounded control volume $V\subseteq {\mathbb{R}}^n,\mathit{x}\in V$, and a temporal component. The definition of the latter is quite elastic [18], though for most purposes, the temporal domain constitutes discrete, uniformly sampled times $t_j$.

In contrast, the objective that informs SPOD is to generalise the decomposition to enable the detection of temporally coherent structures within the flow. This requires a different definition of the stochastic ensemble $\mathit{q}$: instead of viewing individual snapshots as a realisation, we considered a series of correlated snapshots. The inner product associated with $\mathcal{H}$ is defined as

$${\langle {\mathit{q}}_1,{\mathit{q}}_2\rangle }_{\mathit{x},t}=\underset{-\infty }{\overset{\infty }{\int }}\underset{V}{\int }{\mathit{q}}_2^{*}(\mathit{x},t)\mathbf{W}\left(\mathit{x}\right){\mathit{q}}_1(\mathit{x},t)\mathrm{d}\mathit{x}\mathrm{d}t,$$

(4)

where $\mathbf{W}$ is the weight tensor, which further specifies the properties of the norm. Under these parameters, the kernel introduced in (3) becomes

$$\mathbf{R}(\mathit{x},{\mathit{x}}^{\prime },t,t^{\prime })=\mathrm{E}\left\{\mathit{q}(\mathit{x},t){\mathit{q}}^{*}({\mathit{x}}^{\prime },t^{\prime })\right\}\mathbf{W}.$$

(5)

However, the eigenmodes produced by (5) have some inconvenient properties; POD only considers statistically stationary flows, which are unbound in time and, therefore, have infinite energy in the space–time norm. One solution relies on the ergodicity assumption, instead formulating a cross-correlation tensor, which only considers differences between points in time:

$$\mathbf{R}(\mathit{x},{\mathit{x}}^{\prime },t,t^{\prime })\to \mathbf{R}(\mathit{x},{\mathit{x}}^{\prime },\tau ),$$

(6)

with $\tau =t-t^{\prime }$. The cross-spectral density tensor $\mathbf{S}$, which serves as the kernel of the SPOD, is the Fourier transform of $\mathbf{R}$:

$$\mathbf{S}(\mathit{x},{\mathit{x}}^{\prime },f)=\underset{-\infty }{\overset{\infty }{\int }}\mathbf{R}(\mathit{x},{\mathit{x}}^{\prime },\tau )e^{-\mathrm{i}2\pi f\tau }\mathrm{d}\tau$$

(7)

We can construct the spectral eigenvalue problem:

$${\int }_V\mathbf{S}(\mathit{x},{\mathit{x}}^{\prime },f)\mathbf{W}\left({\mathit{x}}^{\prime }\right){\mathit{\psi }}_k({\mathit{x}}^{\prime },f)\mathrm{d}\mathit{x}=\lambda \left(f\right){\mathit{\psi }}_k(\mathit{x},f),$$

(8)

to produce modes $\{{\mathit{\psi }}_k(\mathit{x},f),{\lambda }_k\left(f\right)\}$, which both satisfy the space–time eigenvalue problem once transformed back into the time domain and maintain properties equivalent to space-only POD modes. While these spectral eigenfunctions are frequency-dependent, they are spatially orthonormal (${\langle {\mathit{\psi }}_j,{\mathit{\psi }}_k\rangle }_{\mathit{x}}={\delta }_{jk}$) and serve as an optimal basis for Fourier modes of the flow.

In essence, spectral POD may be applied in fluid mechanics to decompose a statistically stationary flow into a series of characteristic eddies [14], which oscillate at a single frequency and optimally capture a space–time energy distribution.

2.2. Computational Algorithms for SPOD

In practice, the successful application of POD relies on algorithms that efficiently perform modal decomposition from large amounts of discrete flow data, whether of experimental or numerical origin. A thoroughly documented spectral scheme termed batch SPOD was presented by Schmidt and Towne [19]. For reproducing a functional SPOD scheme, this work also benefited from the open SPOD repository generously provided by Schmidt and the California Institute of Technology [20]. What follows briefly presents the core data processing methods used in this work.

As discussed in Section 2.1, the individual realisations forming the basis of the statistical analysis must be time-resolved. While the definitive method of obtaining such data is to run the same experiment multiple times, this approach is typically impractical, especially when working with pre-existing sets of data generated for purposes other than SPOD. Additionally, more complex experiments may be impossible to replicate with sufficient precision to be considered the same process. An alternative solution is Welch’s method [21], a spectral estimation technique in which a single long recording of an ergodic signal is separated into blocks whose spectra are then evaluated statistically. Crucially, it can be demonstrated that the spectral estimates computed with this method converge for an increasing number of snapshots $n_t$, unlike those directly evaluated from the time series. Batch SPOD transplants this idea onto continuous recordings of flow data.

Consider a set of $n_t$ 2D flow data snapshots $\mathit{q}(\mathit{x},t_j)$ with $N=n_x{n}_y$ discrete, equidistant points in the control volume V, where each instance is a vector field representing the instantaneous velocity distribution of the flow. For the purposes of POD, the exact nature of the spatial grid is immaterial; it suffices to sequentially index each vector in the field as ${\mathit{u}}_j^k$, i.e., the velocity at the jth position and the kth time step. We first reorganised the dataset: the x- and y-components of the velocity vectors were separated and each assembled into an $N\times n_t$ matrix. These were then vertically concatenated to produce $\mathbf{U}$, an array containing all input data for the decomposition. $\mathbf{U}$ was then sectioned into $n_{\mathrm{blk}}$ blocks ${\mathbf{U}}_m$ of length $n_{\mathrm{FFT}}$. These blocks may overlap to a certain extent to decrease the variance of the power spectral density (PSD) estimate, though in practice, overlaps larger than 50% yield greatly diminishing returns as the independence of the realisations is increasingly compromised [16,21].

The next step applies a row-wise discrete Fourier transform (DFT) to each block, specifically the DFT algorithm known as the fast Fourier transform (FFT). This produces a set of matrices ${\widehat{\mathbf{U}}}_m\in {\mathbb{C}}^{n_f\times N}$, $n_f$ being the number of frequencies resolved by the FFT, which now contain the transformed data sampled at discrete frequencies. SPOD solves the eigenvalue problem (8) one frequency at a time, by rearranging the rows of ${\widehat{\mathbf{U}}}_m$ into a different set of $n_f$ data matrices ${\mathbf{X}}_f\in {\mathbb{C}}^{N\times n_{\mathrm{blk}}}$, which groups the data from each block by frequency. Given this, the cross-spectral density tensor is

$${\mathbf{S}}_f={\mathbf{X}}_f{\phantom{(}}^{*}\mathbf{W}{\mathbf{X}}_f.$$

(9)

(9) produces a small $n_{\mathrm{blk}}\times n_{\mathrm{blk}}$ tensor, for which the eigenvalue problem is computationally inexpensive to solve. Consequently, modes ${\left({\mathit{\theta }}_m\right)}_f$ obtained in the eigendecomposition:

$${\mathbf{S}}_f={\left(\mathbf{\Theta }\mathbf{\Lambda }{\mathbf{\Theta }}^{-1}\right)}_f$$

(10)

must be transformed before interpretation in the physical flow domain:

$${\mathbf{\Psi }}_f={\left(\mathbf{X}\mathbf{\Theta }{\mathbf{\Lambda }}^{-1/2}\right)}_f={\left[\begin{array}{cccc}{\mathit{\psi }}_1& {\mathit{\psi }}_2& \cdots & {\mathit{\psi }}_{n_{\mathrm{blk}}}\end{array}\right]}_f$$

(11)

After assembly into 2D fields, the modes ${\left({\mathit{\psi }}_m\right)}_f$, represent characteristic eddies of the flow that evolve coherently in time, i.e., oscillate at a single discrete frequency f.

3. Experimental Approach

This section introduces the rig on which the experiments were conducted and briefly details the equipment used for recording PIV and other experimental data.

3.1. TCF Model

The data used in this work were recorded on the TCF test rig assembly, located at the University of Stuttgart’s Institute for Aerospace Thermodynamics. The setup is described in detail by Köhler et al. [22]; we limit ourselves to briefly introducing the relevant elements.

The centrepiece of the rig was a machined PERSPEX® (acrylic glass) model of a TCF ring segment (illustrated in Figure 1), the flow through which was examined by PIV. Its geometry was scaled and simplified to improve optical access to the interior: curved surfaces, such as the vane and hot gas duct walls, were generally approximated with flat panels to avoid unnecessarily complex reflections. The block on the air inflow side mimicked a segment of the hollow outer ring, which was connected to the hub cavity of the inner ring on the opposite side of the hot gas duct by a symbolic load-carrying strut. The strut was encased in a fairing, under which the cooling air passed between the cavities.

Air for the experiment was taken from the environment, pumped into a plenum chamber at a constant, adjustable mass flow rate and homogeneously seeded with droplets of di-ethyl-hexyl-sebacate, which acted as the tracer particles for PIV. The air entered the model through an elliptical inlet gap with a hydraulic diameter of $d_{\mathrm{h}}=26.9\mathrm{mm}$, forming an annular jet, passed through the model, and exited through five outlets, three of which were situated around the tip cavity and two around the hub cavity. The flow conditions in the model were adjusted by modulating the mass flow rate through the air supply duct until the desired Reynolds number at the elliptical inflow gap was obtained.

The control volume of the experiment was the TCF model’s interior (hereafter referred to as the fluid domain), shown in Figure 2. Its maximum dimensions were $400\mathrm{mm}\times 500\mathrm{mm}\times 430\mathrm{mm}$. PIV was performed through the transparent walls of the model, with the laser and camera installed outside at different positions. As 2D PIV only captures the velocity distribution within a planar section of the control volume, measurements were confined to a single plane at a time. The two relevant planes for this work were Y00, the symmetry plane, and Z15, both highlighted in Figure 2; the latter derives its name from being the XY plane situated at 15% along the z-length of the fluid domain.

Figure 3 shows the Y00 and Z15 slices in more detail. As the camera field of view (FOV) was limited, making it impossible to capture the entire domain in a single run, the bottom of the hub cavity (highlighted in grey for each plane) was chosen as the focus area.

While measuring in Y00 is best suited to creating an overview of the dynamics, the evaluation of Z15 allows for better understanding of the spatial characteristics of the flow phenomena. Using the information from both angles, the modal behaviour can be partially interpreted in 3D, which aids in the understanding of three-dimensional flow structures.

3.2. Measuring Equipment

The PIV method employed in this work was a high-speed 2D approach, in which the laser and camera were triggered to produce double images at a prescribed frequency. These images were interpreted by software, which matched the $32\times 32$ pixel interrogation windows in two successive frames to produce vector fields representing the instantaneous velocity distribution in the flow. The specifics of PIV are beyond the scope of this work; while the choice of operating parameters affects the evaluation and will be discussed in Section 5.1, we must accept that the raw images can be reliably processed into instantaneous vector fields if the parameters are well chosen. For more technical insight, the statistical background and application of PIV are explored thoroughly by Raffel et al. [23].

The laser employed in the PIV was an INNOSLAB IS200-2-L single-cavity green short-pulse Tm:YLF unit manufactured by EdgeWave GmbH. The laser pulse frequency can be adjusted up to a maximum of $10\mathrm{kHz}$. The images were recorded with a Photron FASTCAM SA-X2, which was triggered at the data acquisition frequency $f_{\mathrm{s}}$ to produce two consecutive images. This resulted in a maximum total of 5000 velocity snapshots per measurement. The commercial software LaVision DaVis 10.2 was used for post-processing.

Aside from PIV, a hot wire anemometer was deployed inside the fluid domain to record 1D velocities. The hardware used was an SVMtec GmbH SVM-4CTA-DAQ anemometry system coupled with HWP10/90 10 $\mathsf{\mu }$m gold-plated wolfram wire probes. The probe calibration data suggested an error of roughly 1% in the detected velocity in the relevant range.

Due to its low bandwidth and storage requirements, the 1D signal can be written at extremely high sampling frequencies for a long time, allowing highly resolved PSDs to be obtained. The probe was inserted into the hub cavity at the base of the exhaust funnel; measurements were conducted at different depths of insertion, as indexed in Figure 3a, as well as different Reynolds numbers. The positions denoted by 1 and 2 were at $x_{\mathrm{HW}}=[60,100]\mathrm{mm}$, respectively, as measured from the inner wall of the cavity.

4. Initial Analysis

While high-speed PIV is by far the most-powerful available tool for experimental analysis of fluid systems, the necessary equipment is expensive and, therefore, in high demand. Gathering information about the behaviour of the system in advance, using more accessible methods, helps ensure the effective use of the equipment later on. At the same time, obtaining similar results through functionally diverse avenues of analysis serves to validate the involved methods.

In this section, several available sources of data were invoked to narrow down the spatial locations and spectral range of the characteristic structures within the TCF model. This information can be used to set the operational parameters of the PIV, such as the image acquisition frequency, as well as hone the post-processing algorithms used to evaluate the recordings.

4.1. Mean Flow

The baseline for understanding the dynamics of any fluid system is examining the time-averaged conditions. Fortunately, low-speed PIV data recorded at $f_{\mathrm{s}}=8\mathrm{Hz}$ are abundant from previous investigations of the same model [22]. While these recordings proved unhelpful in exploring the time evolution of the system, they were useful for gaining an overview of the mean flow through the fluid domain. The time-averaged flow speed for $\mathrm{Re}=\mathrm{37,000}$ is presented in Figure 3a.

Köhler et al. [22] provide a comprehensive overview of the topology of the mean flow field. The most-dominant structure in the fluid domain was a jet that remains attached to the strut approximately halfway along the fairing, before separating, crossing the hub cavity diagonally and impacting the inner wall in the bottom left corner, where a stagnation area developed. At the same time, it featured sufficiently low velocities to be examined by PIV. The anemometer probe locations and PIV windows detailed in Figure 3 were defined to most-effectively capture it.

The mean flow in the transverse Z15 plane, as computed from PIV, is depicted in Figure 3b. Unsurprisingly, the topology was largely symmetrical along the y-axis. It was dominated by the footprint of the hub cavity jet, whose elliptical shape was owed to the orientation of the plane; from comparisons with Figure 3a, we can infer that it intersected the jet at an angle of roughly ${45}^{\circ }$. Consequently, the section contained the vertical component of the jet’s velocity vector. The most-energetic region was found at $x/D=2$, very close to the bottom wall; in fact, the jet already experienced significant deflection this far downstream, as evidenced by the stagnation point and higher speed outwash regions along the boundary of the camera window. Towards the top, the shape split in two, where the flow separated around the central strut. The lowest flow speeds were found in the shear layer on either side of the jet; beyond it, the field showed air symmetrically moving slightly inwards at $x/D=1.2$ before being deflected straight down and finally outward by $x/D=1.8$. This suggested the presence of larger vortical structures that extended beyond the camera FOV.

4.2. Hot Wire Anemometry

As laid out in Section 3.2, a hot wire anemometer was used in the early stages of the work to gather some information about the dynamics of the internal flow field with a method that is well-understood, reliable, and simple to deploy. The main objectives for this step were to identify dynamically active ranges in the PSD and establish that they are robust towards changes in the Reynolds number. These data were used later as a direct reference to validate the PSDs obtained through more sophisticated post-processing (see Section 5.2 and Section 6.2).

The data were recorded at every position and mass flow rate with a sampling frequency of $f_{\mathrm{s}}=5\mathrm{kHz}$ for $T_{\mathrm{HW}}=60\mathrm{s}$. The three flow regimes that were investigated were ${\mathrm{Re}}_j=[37,55.5,74]\times {10}^3$. We interpreted the data by performing a Fourier transform and examining the resulting PSDs. Keeping with established practises in signal processing [21], the PSD was constructed as the average of $n_{\mathrm{blk}}$ blocks of length $n_{\mathrm{FFT}}$; a block size of $n_{\mathrm{FFT}}=2^{15}$ was selected as a suitable compromise between conclusive averaging and sufficiently informative spectral resolution.

The results of this evaluation are presented in Figure 4. The PSDs at different flow regimes ${\mathrm{Re}}_j$ were grouped by position $x_{\mathrm{HW},k}$; the raw spectra are each shown in grey along with a smoothed spectrum, depicted in red, obtained by applying a Savitzky–Golay finite impulse response filter (SGF) of order three with a frame length of 31 data points. The SGF is based on low-order polynomial fits and the linear least-squares method and is commonly used to smooth 1D data in signal processing and experimental applications such as spectroscopy [24,25]. The filter is a convenient tool to identify trends in an otherwise noisy spectrum, which establishes an overview of the flow dynamics; at the same time, the spectral clarity in low frequency bands remains superior to simply evaluating the FFT with smaller block sizes. The frequency axis was nondimensionalised as the Strouhal number $\mathrm{Sr}=fD/U_0$ to ease the comparison of the data at different Reynolds numbers, where $U_0$ is the nominal inlet flow speed.

As expected, the absolute amount of energy in the flow rose with increasing Reynolds number regardless of position. As our primary concern was the shape of the spectrum, the results at both positions were normalised to the highest detected amplitude for the sake of clarity.

The initial analysis of the results was encouraging. At ${\mathrm{Re}}_1$, both Figure 4a,b feature a pair of distinct peaks at roughly $\mathrm{Sr}=0.014$ and $0.07$. A third peak emerged after an increase in Reynolds number to ${\mathrm{Re}}_2$, though the general structure of the spectrum was preserved: for $x_{\mathrm{HW},1}$, we observed peaks at $0.012$ and $0.149$, as well as a broad elevation between $0.056$ and $0.093$. At ${\mathrm{Re}}_3$, the same pattern appeared at $0.012$, $0.07$–$0.091$ and $0.119$. A very similar structure was observed at $x_{\mathrm{HW},2}$, where, for ${\mathrm{Re}}_2$, the three peaks were registered at $0.013$, $0.07$–$0.084$, and $0.131$ before shifting to $0.012$, $0.053$–$0.084$, and $0.116$ after an increase in Reynolds number to ${\mathrm{Re}}_3$.

In summary, we concluded that, for $\mathrm{Re}\ge \mathrm{37,000}$, the system exhibited strong organised unsteady behaviour in the $\mathrm{Sr}=0$ to $0.2$ range, which can be clearly separated from broadband turbulence in the spectra derived from hot wire anemometry. In this range, the general structure of the PSD was invariant to changes in Reynolds number when presented in nondimensional space, which suggested stable oscillating behaviour. Peaks at higher Strouhal numbers notably exhibited bigger shifts. The spectra did not differ significantly based on the position of the probe, suggesting that the dominant motion was sufficiently extensive to cover both examined locations. For future reference, the two frequency bands in which the reference signal was most active at ${\mathrm{Re}}_1$ will be referred to as Band 1 between $\mathrm{Sr}=0.01$ and $0.03$ and Band 2 between $\mathrm{Sr}=0.07$ and $0.1$.

5. Conventional SPOD

Section 4.2 presents us with strong evidence for stable, coherent structures in the fluid system, which was essential to justify further analysis. The combined anemometry data also provided an estimate of the dynamically active range of the system, which we captured with HSPIV.

5.1. Data Acquisition

A PIV experiment can be tailored to optimally resolve any system by adjusting several process parameters. The first barrier is the successful evaluation of the images recorded through PIV, which provides the instantaneous velocity fields that form the basis for all further post-processing. The experiments showed that the maximum particle displacement, beyond which the DaVis cross-correlation algorithm failed to reliably converge, was reached if the inflow Reynolds number meaningfully exceeded ${\mathrm{Re}}_1=\mathrm{37,000}$. This flow regime was, thus, chosen for the HSPIV measurements, as the analysis of anemometry data in Section 4.2 suggested that strong fluctuations were already present in the system at this point, and the behaviour showed significant commonality with higher Reynolds number regimes.

However, while a high pulse frequency was necessary to produce the velocity fields, single-pulse PIV at $f_{\mathrm{s}},\max =10\mathrm{kHz}$ would fill the camera’s internal storage after a mere second of operation. Such a brief run would include too few cycles of the organised motion in the system’s dynamically active range to be statistically meaningful. We, thus, opted for simulated double-pulse operation to produce snapshots at $f_{\mathrm{s}}=500\mathrm{Hz}$, which resulted in a total of 5000 snapshots per measurement run, produced over a period of $T_{\mathrm{s}}=10\mathrm{s}$. This approach maintained a time step of $\mathsf{\Delta }t={10}^{-5}\mathrm{s}$ between frames, which was necessary to correlate the snapshots.

With this setup, four independent runs were performed in the Y00 plane, referred to as Y00${}_{\mathrm{I}}$ through Y00${}_{\mathrm{IV}}$, as well as a run in the Z15 plane. An additional Y00 recording was performed with different parameters, that is a sampling frequency of $f_{\mathrm{s}}=125\mathrm{Hz}$. The rationale for this decision was that the resulting cutoff frequency $f_{\mathrm{N}}=62.5\mathrm{Hz}$ was the lowest possible choice of $f_{\mathrm{N}}$ that did not infringe upon the dynamically relevant range. The immediately apparent drawback of this approach is that the laser must still fire continuously at a rate of $10\mathrm{kHz}$; the energy input would inevitably damage the TCF model if it were subjected to it for too long, which limited the maximum measurement period to $T_{\mathrm{s}}=20\mathrm{s}$. This setup can, therefore, only produce 2500 snapshots. It is henceforth referred to as Y00${}_{125}$.

5.2. Evaluation

As with the anemometry data, we begin our discussion of the SPOD results with spectral analysis. The individual measurements were characterised primarily by their leading modes; the energy in the remaining modes was distributed near-uniformly across the frequency spectrum and trailed the leading mode’s amplitude by roughly an order of magnitude. We, thus, focused on the leading modes when discussing parallels and differences between measurements.

Figure 5 shows the leading modes of the runs introduced in Section 5.1, as well as the Strouhal number bands obtained from evaluating the hot wire data at ${\mathrm{Re}}_1$ in which the reference exhibited increased activity. The filtered hot wire spectrum at $x_{\mathrm{HW},1}$ and ${\mathrm{Re}}_1$ was included as a prediction for the shape of the PSD distribution. For the sake of simpler qualitative analysis, the spectra were normalised by the tallest detected peak in each measurement. Additionally, in working with the PIV data, PSD information below $\mathrm{Sr}=0.01$ was consistently chaotic and unrepeatable. This likely was a result of a very small number of full cycles being included in the recording; consequently, these data were excluded from the following analysis. Note that, due to the complexity of the signal, each individual measurement was unlikely to produce the exact desired result. Instead, we chose to evaluate the superposition of $n=5$ measurements, as with $n\to \infty$, we expected the cumulative spectrum to approach the reference.

The takeaways from Figure 5 are mixed. While it is encouraging that the dominant peaks in all PIV-derived spectra occurred between $\mathrm{Sr}=0.01$ and $\mathrm{Sr}=0.1$, which matched the reference, they failed to replicate the twin peak clusters around $\mathrm{Sr}=0.02$ and $\mathrm{Sr}=0.08$ present in the hot wire spectra. Instead, the tallest peaks of several measurements, including Y00${}_{\mathrm{II}}$, Z15, and notably, Y00${}_{125}$, were located between $\mathrm{Sr}=0.03$ and $\mathrm{Sr}=0.06$, a range in which the reference displayed much lower amplitudes. Evidently, simply lowering $f_{\mathrm{s}}$ did not immediately yield favourable results.

6. S2POD

The results presented above failed to capture the dynamics identified in the reference. The options for continued analysis with conventional SPOD were unpromising: an increase in data acquisition frequency $f_{\mathrm{s}}$ would inevitably lead to even poorer spectral resolution, while the physical barrier of a maximum recording time $T_{\mathrm{s}}$ caused measurements with decreased $f_{\mathrm{s}}$ to produce less data overall. Faced with this dilemma, we propose a novel preprocessing method termed spectral subsampling POD (S2POD).

6.1. Establishing Ideas

Section 2.2 introduced the idea of the subdivision of the data into blocks, then processed separately using Welch’s method. To reduce the effective sampling frequency $f_{\mathrm{s}}^{*}$, we created the blocks ${\mathbf{U}}_m$ by subsampling the original sequence: assembling every other frame into a block created a realisation of the process over the same length of physical time, but halved the sampling frequency. We can increase the subsampling factor $C_{\mathrm{sub}}$ further until we obtain enough blocks to compute a statistically satisfactory number of modes with $f_{\mathrm{s}}^{*}=f_{\mathrm{s}}/C_{\mathrm{sub}}$.

Figure 6 demonstrates the subsampling method for a fictional scalar value at a single point in the control volume, $\left|{\mathit{u}}_j\right|$. The progression through time corresponded to the column index k. For illustrative purposes, a low base sampling frequency was chosen.

Aside from the sample rate, the four sequences created by this method varied from the original in several ways. We observed a shift in the absolute beginning and end times as a result of the frames being discarded, though absolute time was irrelevant, as most further computation occurred in the frequency domain. A distortion of the signal was avoided, as we assumed that the base sampling frequency was much higher than the most-energetic components of the flow; $f_{\mathrm{s}}^{*}$ can be lowered significantly before they begin approaching the Nyquist limit $f_{\mathrm{N}}$.

This approach offers certain advantages. Selecting parameters for traditional SPOD involves a compromise between the number of modes $n_{\mathrm{blk}}$ and the frequency resolution, characterised by $n_{\mathrm{FFT}}$, while the subsampling approach always yields the optimal frequency resolution for a given set of data $\mathsf{\Delta }f=f_{\mathrm{s}}/n_t$. The decision was instead made between a larger number of modes $n_{\mathrm{blk}}=C_{\mathrm{sub}}$ and a broader frequency domain:

$$f_{\mathrm{N}}=\frac{f_{\mathrm{s}}}{2C_{\mathrm{sub}}}$$

(12)

As this method is primarily geared towards applications where the sampling rate is so high that frequency resolution in the dynamically active range is poor, the top end of the spectrum rarely contains relevant information. Consequently, we traded a lower $f_{\mathrm{N}}$ for improved statistical convergence. However, it must be noted that realisations created through subsampling were strongly correlated, and the variation of the system was, therefore, largely captured by the leading mode.

A simple numerical simulation demonstrated the effectiveness of the procedure. Using an SMAC Navier–Stokes solver [26], we computed snapshots of the flow around a cylinder with the boundary conditions tuned to produce a Kármán vortex street at a cylinder with a Reynolds number of $\mathrm{Re}=100$. The complete dataset, trimmed to only contain fully developed flow, comprised 10,000 vector fields with $\mathsf{\Delta }t=0.02\mathrm{s}$, i.e., sampled with $f_{\mathrm{s}}=500\mathrm{Hz}$.

Without subsampling, as seen in Figure 7a, the spectrum contained a single broad peak at $f=0.2441\mathrm{Hz}$—even with the balance shifted heavily towards the frequency resolution, as $n_{\mathrm{FFT}}$ was so large that the data yielded only three modes. The results were drastically more favourable once subsampling was applied: Figure 7b shows a a clear peak at $f=0.25\mathrm{Hz}$, as well as a series of progressively weaker harmonics at $f=0.25n\mathrm{Hz}$, $n=2,3,\cdots$, which was consistent with the results obtained by other analyses of the dynamics of the Kármán vortex street [27]. Simultaneously, $n_{\mathrm{blk}}=10$ and low amplitudes for the lesser modes showed the reliable statistical convergence of the process.

Of course, as a consequence of the method, $f_{\mathrm{N}}$ dropped from $250\mathrm{Hz}$ to just $25\mathrm{Hz}$. However, as displayed in the PSD graphs, even a much lower limit still allowed us to draw informative conclusions about the dynamics of the flow.

The plausibility of the result can be verified by estimating the main frequency of the vortex shedding: Schlichting and Gersten [28] suggest that, for flows around a cylinder, the Strouhal number $\mathrm{Sr}=fd/u_{\infty }$ is correlated with the Reynolds number, which was known to be $\mathrm{Re}=100$ in the simulation. Thus, with $\mathrm{Sr}\approx 0.15$, as well as cylinder diameter $d=0.5625\mathrm{m}$ and $u_{\infty }=1\mathrm{m}/\mathrm{s}$, we obtained $f\approx 0.2667\mathrm{Hz}$ for the main shedding frequency. This aligned with the peaks observed in Figure 7.

Another issue to be addressed is the aliasing error, which inevitably contaminated the downsampled signal. Signal components oscillating at frequencies exceeding the new effective maximum frequency $f_{\mathrm{N}}^{*}$ can no longer be resolved by the Fourier transform. If left unheeded, this would result in a higher overall noise background and artificially inflated eigenvalues on secondary modes. A distortion of the downsampled signal was prevented by applying a lowpass filter at $f_{\mathrm{N}}^{*}$ before breaking the data up into blocks, discarding the energy contribution of the high-frequency components. This basic level of post-processing was applied whenever subsampled data are discussed within this work.

6.2. Results

Having demonstrated the effectiveness of S2POD on a predictable, simulated system, we shall attempt to improve upon the analysis presented in Section 5. The PSDs computed through S2POD with $C_{\mathrm{sub}}=8$ for the $500\mathrm{Hz}$ HSPIV measurements are showcased in Figure 8. Y00${}_{125}$ is not included; any downsampling with $C_{\mathrm{sub}}>2$ would reduce $f_{\mathrm{N}}$ below the acceptable threshold.

Figure 8 significantly improves upon Figure 5. We can identify several key differences. A brief glance reveals that the noise background was substantially lower, even at low frequencies. This was owed primarily to the improved spectral resolution, allowing for a greater dynamic range. The energy distribution of the superimposed S2POD solutions matched the hot wire reference spectrum very well. The highest observed peaks of each measurement, with the exception of Y00${}_{\mathrm{I}}$, fell within a $\mathsf{\Delta }\mathrm{Sr}=\pm 0.01$ window of the dominant peaks in the reference: Y00${}_{\mathrm{II}}$ and Z15 matched the identified bands very well at $\mathrm{Sr}=0.072$ and $\mathrm{Sr}=0.092$, respectively. As a result, the trough between the two band boundaries at $\mathrm{Sr}=0.03$ and $\mathrm{Sr}=0.06$, which failed to materialise in Figure 5, was clearly identifiable. While the tallest peaks in Band 1 in the measurements fell somewhat higher on the spectrum than the peak of the reference curve, this can be attributed to overall poorer performance at low frequencies.

These positive findings allowed us to realise the key advantage SPOD offers over conventional methods, i.e., the production of spatial modes that directly magnify the utility of spectral analysis. To elaborate on the brief introduction given in Section 2.2, the modes were spatial representations of the characteristic eddies, which the decomposition algorithm identified as energy-optimal basis vectors of the original sequence. Each mode oscillated periodically at a specific Strouhal number; the PSD supplied information regarding which frequencies and, consequently, which modes contained the highest energy content. We can, thus, reconstruct a reduced-order model of the fluctuating component of the fluid system by superimposing the leading mode at several peaks in the power spectrum.

The decomposition algorithm produced the spatial modes at a random phase ${\phi }_0$. This instantaneous slice frequently failed to convey the topology of the mode; it was, therefore, preferable to inspect the modes at different phase shifts $\mathsf{\Delta }\phi$. As we considered only the absolute magnitude of the velocity, the period of the oscillation was reduced from $2\pi$ to $\pi$, which we discretised into steps of $\pi /3$ for the visual presentation. Such a discrete progression of each mode through their $\pi$-periodic cycle is presented in Figure 9, Figure 10, Figure 11 and Figure 12.

The two dynamically active bands detected, through spectral analysis of each feature, a geometrically distinct characteristic eddy, which repeatedly emerged at significant peaks. Data in the Y00 and Z15 planes can be combined for a more complete phenomenological description of the three-dimensional fluid flow.

Representative spatial modes for the first band are presented in Figure 9 and Figure 10. The Y00 section was characterised by two parallel lobes, which roughly followed the contour of the annular jet with a clearly defined trough between them. While the bottom lobe, most prominently visible at $\mathsf{\Delta }\phi =2\pi /3$, remained predominantly straight, the top lobe fanned out and formed an identifiable stagnation area at the wall between $x/D=1.8$ and $2.2$. Most importantly, the phase evolution of the structure was clearly alternating: the top and bottom lobes swelled and receded periodically, one giving way to the other.

The Z15 section revealed a complementary image. For most of its cycle, the eigenmode evolved within the jet’s outer layer, i.e., the tall elliptical structure between $y/D=-0.2$ and $0.2$ that surrounded the core and featured the greatest mean flow speed gradient in the domain. The mode’s cycle began with an elongated lobe stretching along the shear layer to the right of the core; see Figure 10a. As the phase advanced, the lobe migrated upward and to the left, where it assumed a more compact, rounded shape at $\mathsf{\Delta }\phi =\pi /3$. By $\mathsf{\Delta }\phi =2\pi /3$, the lobe moved to the shear layer to the left of the core.

Having gathered information about the spatial characteristics of the oscillations in this frequency band in two perpendicular planes, we may attempt to reconstruct the three-dimensional structure of the phenomenon at hand. Both spatial modes in Y00 and Z15 exhibited decidedly directional motion surrounding the jet; while the analysis in Y00 consistently found elongated lobes moving perpendicular to the jet, emerging in the shear layer and dissipating in it again after traversing the core, Z15 indicated that the dominant motion in the cross-section plane was counterclockwise rotation throughout the shear layer. These factors aligned with the suggested presence of a counterclockwise precessing component of the jet, induced by the transition from an annular structure at the inflow gap. Cafiero et al. [29], in their investigation of a three-dimensional turbulent jet, registered precession at Strouhal numbers between $0.01$ and $0.025$; this aligned with the range detected in our experiment.

The modes associated with the second band differed significantly from the patterns discussed above. In Figure 11, the domain is split into two wide, staggered bands along the core of the jet. The bands progressively expanded, achieving their maximum amplitude at different $\mathsf{\Delta }\phi$; while the top band was most prominent at $\mathsf{\Delta }\phi =0$, the bottom band dominated at $\mathsf{\Delta }\phi =2\pi /3$. The corresponding structure in Z15 was found largely in the shear layer, though it expanded into the core of the jet further than the lobes seen in Figure 10. We can successfully match the lobe position in Z15, where the Y00 lobes intersected the perpendicular plane, e.g., Figure 11a and Figure 12a.

The most-notable difference from the previous three peaks was the directionality of the time evolution. While the spatial patterns of the first band exhibited predominantly rotary motion, modes in Band 2 exhibited cyclic growth parallel to it. This observation helped substantiate the hypothesis that the higher frequency peaks belonged to a distinct flow phenomenon from the previously discussed patterns.

7. Concluding Remarks

The focus of this work was an exploration of the utility of spectral proper orthogonal decomposition (SPOD) as an instrument for primary analysis of complex flow fields and the detection of organised motion therein. As an example of a largely unknown system, with practical engineering applications in mind, an abstracted model of a radial slice of a turbine centre frame (TCF) was selected to investigate its internal flow structures.

As a first measure of gauging the system’s dynamics, a hot wire anemometer was deployed in the hub cavity. Spectral analysis of the 1D signal at various positions confirmed the existence of organised motion outside of the range typically associated with turbulence mechanisms in open jets, justifying further investigation with more sophisticated equipment. Beyond that, the dominant unsteady structures were found to be stable throughout the entire range of Reynolds numbers at which the experiment was conducted.

The base procedure was found to perform poorly in scenarios where the sample rate was high while the total length of the recorded sequence was short; this excess in time resolution led to unnecessarily wide frequency ranges in spectral analysis, i.e., power spectra that were dominated by broadband turbulence, while the regions of interest were represented by only a few sample points. To remedy this shortcoming, a modification of the batch SPOD algorithm was developed. Instead of conducting Welch averaging with blocks split from the original sequence, the blocks were assembled by subsampling, i.e., assembling every nth frame into a series of new sequences with a lower effective sampling frequency $f_{\mathrm{s}}^{*}=f_{\mathrm{s}}/n$. This procedure achieved maximum possible spectral resolution at the expense of a lower Nyquist frequency, which in most applications is a favourable trade-off.

Armed with S2POD and high-speed particle image velocimetry (PIV) data in two perpendicular planes, we successfully identified two distinct flow phenomena:

• A precessing structure in the shear layer surrounding the hub cavity jet, characterised by counterclockwise motion at several peaks between $\mathrm{Sr}=0.01$ and $0.03$.
• A modulation of the jet’s flow speed at several peaks between $\mathrm{Sr}=0.07$ and $0.1$.

This analysis suggested that the inclusion of subsampling successfully broadens the applicability of SPOD. The dilemma between sampling frequency and recording time is inherent to all PIV experiments, even though it can be amortised through more sophisticated equipment in some cases. S2POD offers a way out without the need for hardware upgrades. Similarly, existing experimental or numerical data may not have ideal parameters for decomposition if they were originally recorded with different post-processing methods in mind. $\mathsf{\Delta }t$ tends to be very small, especially in numerical solutions, making it time-consuming and expensive to compute for extended periods of physical time. Gaining insight into the dynamics of such approaches through S2POD allows it to be tuned to resemble real systems beyond matching time-averaged conditions.