Robust Mode Analysis

In this paper, we introduce a model-free algorithm, robust mode analysis (RMA), to extract primary constituents in a fluid or reacting flow directly from high-frequency, high-resolution experimental data. It is expected to be particularly useful in studying strongly driven flows, where nonlinearities can induce chaotic and irregular dynamics. The lack of precise governing equations and the absence of symmetries or other simplifying constraints in realistic configurations preclude the derivation of analytical solutions for these systems; the presence of flow structures over a wide range of scales handicaps finding their numerical solutions. Thus, the need for direct analysis of experimental data is reinforced. RMA is predicated on the assumption that primary flow constituents are common in multiple, nominally identical realizations of an experiment. Their search relies on the identification of common dynamic modes in the experiments, the commonality established via proximity of the eigenvalues and eigenfunctions. Robust flow constituents are then constructed by combining common dynamic modes that flow at the same rate. We illustrate RMA using reacting flows behind a symmetric bluff body. Two robust constituents, whose signatures resemble symmetric and von Karman vortex shedding, are identified. It is shown how RMA can be implemented via extended dynamic mode decomposition in flow configurations interrogated with a small number of time-series. This approach may prove useful in analyzing changes in flow patterns in engines and propulsion systems equipped with sturdy arrays of pressure transducers or thermocouples. Finally, an analysis of high Reynolds number jet flows suggests that tests of statistical characterizations in turbulent flows may best be done using non-robust components of the flow.


Introduction
It is rarely known at sufficient detail that the operation of many technologically important processes relying on strongly driven fluid or reacting flows can be controlled or optimized through modeling. For example, due to chamber feedback and unsteady heat release [1,2], it is difficult to establish a priori how the fuel efficiency and power output of airplane engines or efficacy of industrial reactors can be enhanced while maintaining their integrity. Attempts to reduce the fuel/air ratio in an airplane engine can cause flame extinction [3]; strong driving of industrial reactors can initiate highly damaging, and spontaneously forming "hot spots" [4]. The ability to establish the nature of such flows and identify their changes can aid in design improvements to control or suppress undesirable facets of the flows. Such an approach is outlined here.
A major difficulty in control is incomplete knowledge of the underlying systems [1]. Although the Navier-Stokes equations modeling fluid flows are known [5][6][7][8], reactiondiffusion systems for combustion are typically only known approximately. Even in rare instances when the underlying dynamics are known, the absence of symmetries or other simplifying constraints in realistic configurations preclude the derivation of solutions. Difficulties can originate from nonlinearity of the underlying equations and the resulting sensitive dependence on initial conditions, as well as the preponderance of nonlinearly u(x, t) = ∑ a n (t)Φ n (x). (1) The most widely-used class of coherent structures is derived through proper orthogonal decomposition (POD) [15][16][17][18][19]. They are eigenfunctions of the matrix of pair-wise correlations of the field at different sites, arranged in decreasing order of their contribution to variations in the field. POD gives the most "efficient" decomposition in the sense that at any level of truncation, no other class of coherent structures captures a larger fraction of variations [16]. Typically, dynamics of a n (t)'s are coupled due to nonlinearity of the underlying system.
A second class of coherent structures is given by dynamic mode decomposition (DMD) [20][21][22][23][24]. A dynamic mode is a spatial structure that oscillates at a fixed (complex) frequency; that is, the coefficients in Equation (1) evolve as a n (t) ∼ exp[iΛ n t], for Λ n ∈ C. The spatio-temporal flow dynamics can be interpreted as the sum of oscillating dynamic modes. Under appropriate conditions, DMD is a form of Koopman mode analysis, originally constructed to explain the linearity of quantum mechanics even for classically nonlinear systems [22][23][24][25][26][27]. Unlike for POD modes, the evolution of a dynamic mode is independent of the others. In addition, dynamic modes are a natural class of coherent structures to expand flows containing, for example, periodic vortex shedding.
DMD offers means to partially address the question posed earlier on delineating flow dynamics from noise and other irregular facets. As in Fourier expansion of a noisy signal, DMD expands all facets of the signal u(x, t), including noise, in periodically evolving modes. The question then is which of these modes represent dynamical features of the flow. Robust mode analysis (RMA) identifies them as dynamic modes common to multiple, nominally identical, experimental realizations; the commonality is established through proximity of the dynamic modes and the associated eigenvalues. The dynamic mode decomposition and implementation of RMA are outlined in Section 2.
RMA differs from robust principal components analysis (RPCA), whose intent is to partition a spatio-temporal field into a low-rank constituent and one that consists of a collection of "sparse" modes [28,29]. RPCA can, for example, separate rapid changes in a traffic pattern from a slow varying background. In contrast, RMA does not require differences in time-scales or "energy content" to differentiate robust and non-robust features. We will show through examples that robust modes are not necessarily those with the highest energy.
Three applications of RMA are given in Sections 3-5. Section 3 analyzes PLIF data from an experiment in a reacting flow past a single symmetric bluff body. At sufficiently high Reynolds numbers, bluff-body flows exhibit symmetric vortex shedding, where pairs of vortices are shed periodically from either side simultaneously [3]. In lean mixtures, the flow also exhibits von Karman vortex shedding, where vortices are shed periodically from alternate sides of the bluff body [30]. RMA yields a collection of robust dynamic modes, which can be grouped by comparing downstream flow rates. The two classes so derived are found to be symmetric and von Karman vortex shedding. The results in Section 3 were adapted from the Ref. [31].
Section 4 re-analyzes reacting flows behind a single symmetric bluff body, now using time-series from a small number of transducers, rather than through spatio-temporal fields. Such analyses can prove important in following, for example, changes in flow patterns within airplane engines or in industrial scale reactors. Their failure is often caused by the (slow) progression of domains with micro-cracks or material degradation. Changes in flow patterns may precede catastrophic failure. In contrast to PIV or PLIF which require robust optical access, time-series measurements can be made continuously, since most engines and propulsion systems are equipped with sturdy arrays of pressure transducers and thermocouples. The analysis requires "enhancing" the spatial field using extended DMD [24,27], which in this application, is made using Chebyshev polynomials. RMA, implemented on the extended field, once again, yields one robust mode in rich mixtures and two in lean mixtures. Time-series of projections of the first to pairs of symmetrically-placed transducers can be used to infer that it is up-down symmetric (justifying its association with symmetric vortex shedding), while the second is not (von Karman shedding).
The final application presented in Section 5 involves turbulent jet flows which are quantified using statistical properties of the flow. The experimental data consist of highresolution, high-frequency (100 kHz), two-dimensional velocity fields in a plane of symmetry. The flow data can be used to extract distributions of single-point fluctuations (which are empirically found to be nearly Gaussian) and the variance-covariance matrix for velocity fluctuations at two sites. Moments of structure functions between two points are derived from this observation. The experimental flow fields are likely to be corrupted by perturbations, such as acoustic reflections from the laboratory walls that are not included in an analysis. We assume that perturbations extraneous to the experimental setup can be identified through RMA. The remaining flow fields are those whose statistics should be compared with predictions. In fact, we are able to show that structure functions of the non-robust velocity fields more closely follow the computed one than the original flow.
Section 6 summarizes some implications of our findings.

Dynamic Mode Decomposition
Dynamic mode decomposition (DMD) [20][21][22]24] provides a class of coherent structures. In contrast to POD, which categorizes coherent structures according to correlations, DMD assorts them through spatial structures that oscillate with a single (complex) frequency. Hence, it is a natural approach to extract periodic flow constituents, such as symmetric or von Karman vortex shedding. In addition, DMD has a theoretical basis in the Koopman operator theory [23][24][25], and the dynamics of the modes are decoupled from each other.
We provide a summary of DMD: suppose spatio-temporal dynamics of a field at time-step n is expressed as U n = u(x, t n ), u being a field recorded at site x and the entire "snapshot" is recorded at discrete time-points t n . The primary conjecture in DMD is that the operation to transform from U n to U n+1 is independent of the "initial" state U n , and hence that a sufficiently large number of snapshots can be used to find properties of the underlying dynamics. (Note that snapshots obtained during a transition between two states far-away in state space may not satisfy the requirement.) Under the assumption, when snapshots are recorded at sufficiently short uniform time-intervals (so linear expansion suffices), we could write where A is a representation of the flow dynamics within the duration ∆t = t n+1 − t n . When N is sufficiently large, singular value decomposition is used to find the spectrum of A.
The coherent structures (dynamic modes) are the eigen-vectors Φ n (x), which evolve as a n (t) = exp(Λ n · t), (the generally complex) Λ n being the associated eigenvalue. Thus, each dynamic mode evolves periodically with a period 2π/Imag(Λ n ). The field u(x, t) can be expanded in dynamic modes [20][21][22]24].

Robust Mode Analysis
Robust mode analysis (RMA) uses dynamic modes (1) to differentiate dynamical features of a flow from noise and other irregular facets, and (2) to identify flow constituents. Here, as in POD or Fourier analysis of a time-series, noise and other irregular (e.g., initialcondition differences amplified by sensitive dependence on initial conditions) features contribute to DMD as well. RMA is predicated on the expectation that while dynamical features of the flow will retain their form (i.e., dynamic modes and eigenvalues, although not necessarily the intensity) in different realizations of an experiment, irregular features will differ between experiments. Under this scenario, robust flow modes are those common to multiple realizations of an experiment.
To be specific, denote the dynamic modes and eigenvalues for the first experiment by {Φ n (x), Λ n }, and those of the next set of experiments by {Φ (g) n (x), Λ (g) n }. Dynamic modes from the different groups are assumed to be the same if they satisfy the following criteria: • Imaginary parts of the eigenvalues from multiple groups should be sufficiently close; specifically, where δ I is a cutoff. • The real parts of an eigenvalue identified from the first condition should not vary significantly among the subgroups; specifically, for a cutoff δ R . For work reported here, we limit consideration to δ R = δ I . • Dynamic modes from the different groupings that satisfy Equations (3) and (4) should j (x) are normalized eigenfunctions associated with proximate eigenvalues, then for a cutoff ∆. The minimization over the phase θ is implemented to account for the fact that eigenfunctions computed from different experimental realizations may differ in phase.
Robust mode analysis is conducted on the entire spatial field, which is typically a few thousand measurements. The cutoff values δ R , δ I , and ∆ depend on properties of an experimental system, including the level of noise, and need to be established heuristically. We limit consideration to δ R = δ I in order to restrict the parameter search. We illustrate the selection of the cutoff values using spatio-temporal flow fields in a turbulent jet flow of Reynold's number 21,000, described in Section 5. The data consist of axial and radial velocity fields on a lattice of size N = 67 × 51 at a rate of 100 kHz for a duration 0.08 s.
We note first that for unit vectors Φ i (x) and Φ Furthermore, the probability density function for the angle θ between two random vectors in N dimensions satisfies P(θ) ∼ sin N−2 θ. Thus, for ∆ sufficiently smaller than √ 2, the likelihood of the left side of the last inequality being satisfied with no underlying reason is minuscule. For fixed δ R , we search for robust modes beginning with a small ∆ (say ∼ 0.5), and ascertain increases in their number with increasing ∆. There is a proliferation of robust modes outside the Nyquist range (i.e., |ImagΛ| > π) when ∆ increases beyond ∼0.85. Furthermore, no new robust modes appear in (0.75, 0.85), motivating the choice of ∆ = 0.75. Next, we compute how the number of robust modes increases with δ R ; this is shown in Figure 1. These considerations suggest that there are 13 robust modes in the jet flow, three of which have real Λ's.
We studied the effect of noise on the choice of δ R , δ I , and ∆ by adding synthetic Gaussian "noise" to the experimental data. The magnitude of the synthetic noise was quantified using a measure of the inverse of the signal-to-noise ratio η, defined as the ratio of the standard deviation of the Gaussian noise to that of the field. As η increases beyond ∼0.3, the number of robust modes is found to decrease. The remaining modes can be associated with those for the noise-free case through proximity of the eigenvalues. In addition, some (although not all) of the robust modes lost on adding noise can be recovered by increasing δ R .
Note that cutoffs δ R , δ I , and ∆ do not depend on the (mean) magnitude of the experimental data fields. Since the eigenfunctions are normalized, the difference in phases is independent of the magnitude. In addition, since the eigenvalues of the matrix A are expressed as exp(Λ), the imaginary parts do not depend on the magnitude as well. Although the real parts do depend on the mean magnitude, differences between pairs of eigenvalues from different experiments are independent of the scale.
RMA is best used on data from a group of distinct experiments. However, we have used the following approach in cases where only a single experimental recording is available. The full set of snapshots is assigned as Group 1. The remaining groups are constructed using (1) the first half of snapshots, (2) the second half of snapshots, (3) even-numbered snapshots, and (4) odd-numbered snapshots. We have established robust modes in combustion flow behind a bluff body [31], in axisymmetric injector flow [32], and swirling combustion flows [33] and found flow constituents consistent with prior studies. The analysis also uncovered a novel broken-symmetry mode in injector flows [32]. Although DMD is highly sensitive to noise and subsets of snapshots used for analysis [29,34,35], we find that imaginary parts of eigenvalues of robust modes in different sub-groupings typically do not change by more than a fraction of a percent. We finally consider the construction of a robust flow constituent. As an example, a constituent may include a robust dynamic mode and its second harmonic, which itself may be robust. Our goal is to select robust modes that should be associated with the same flow constituent. For a given dynamic mode, we can use the evolution of the "phase" of the coefficient a n (t) of Equation (1) (i.e., tan −1 (Im a n (t)/Re a n (t))) to quantify its motion. However, we note that the phase of the second harmonic will advance at twice the rate of that of the fundamental. In order to correct this difference, we note that geometric features of the harmonic, such as the number of local peaks, will also be twice that of the primary (as illustrated through examples in the Ref. [31]). Thus, we can associate a phase for each robust mode via θ n (t) = 1 s n tan −1 Im a n (t) Re a n (t) , where s n is the number of local peaks in the dynamic mode. The "flow rate" of the robust mode is then defined as ω n = dθ n (t)/dt. With this definition, the fundamental and its harmonics will have the same flow rate. Once all robust modes with the same flow rate are identified, the associated flow constituent is reconstructed via Applications in the use of the flow rate ω n (t) to construct flow constituents are given in the next section [31].

The Experiment
Experiments were conducted within an optically accessible, atmospheric-pressure combustion test section that contains a bluff-body flame holder for flame stabilization. Air is delivered into a 152 mm × 127 mm rectangular test section at a constant rate of 0.32 kg/s. While the air rate is maintained constant, propane fuel is added and mixed upstream of the flame holder, which is a v-gutter with a width of 38.1 mm and an angle of 35 • . Additional facility details and flame-holder dimensions are provided in the Ref. [3].
The nature of instabilities in bluff-body flows depends on fractions of fuel and air in the mixture. It is quantified using the equivalence ratio φ defined as the fuel/air ratio of the mixture compared to that of its stoichiometric value (when all reactants are fully utilized). Fuel-rich mixtures are those with larger equivalence ratios, while lean mixtures are those with smaller equivalence ratios. In our configuration, flame blow-off occurs when φ < 0.55 and our experiments are performed for equivalence ratios between 0.6 and 1.1. Symmetric vortex shedding [36,37] is observed in the entire range of φ, while the asymmetric von Karman vortices are only observed in a lean mixture (φ 0.85) [3].
A schematic of the experimental setup is shown in Figure 2a [38], and two-dimensional images of hydroxyl (OH) behind the bluff-body v-gutter were acquired utilizing planar laser-induced fluorescence (PLIF), performed using a 10-kHz, diode-pumped, solid-state Nd:YAG laser and a tunable dye laser. The 532 nm output of the Edgewave laser was used to pump the dye laser for obtaining tunable laser output at ∼586 nm. This wavelength was then frequency-doubled at ∼283 nm to excite the Q 1 (9) rovibrational transition in the A 2 Σ+ → X 2 Π (1, 0) band of OH. The Q 1 (9) transition has a low Boltzmann-fraction sensitivity between temperatures of 1000 and 2400 K, minimizing the need for a temperature correction on Boltzmann distributions when extracting flame fronts. The PLIF signal of OH was collected employing a LaVision dual-stage, high-speed UV intensifier (IRO) coupled to a Photron SA-5 CMOS camera. The collected light was filtered using a Brightline Semrock filter with ∼90% transmission between 300 and 340 nm. Each recording contained 8000 frames separated by 0.1 ms. Only the last 4000 frames were used for the DMD analysis.
Below, we summarize results from RMA on a flow at equivalence ratio φ = 0.8, when both symmetric and von Kaman vortex shedding modes are active. These results were previously reported in the Ref. [31].  Figure 3a shows the eigenvalues of DMD modes and Figure 3b, those of the robust modes, which are identified using 2000-frame sections of the video. (Increasing the number of frames further introduces numerical errors in singular value decomposition, as can be inferred from emergence of a large number of eigenvalues with significantly different magnitudes from 1.) We compare dynamic modes from this full sequence with those from the first half, second half, odd-numbered snapshots, and even-numbered snapshots. The real parts of the robust modes vary between the subgroups and we show the mean and standard deviation for each robust mode. Note that the modes that remain past the decay of transients neither decay nor grow, and are thus expected to have zero real parts of their eigenvalues [22,23]; this is seen to be the case for the robust modes in Figure 3b. The mode numbers given in Figure 3b are ordered according to their energies. Robust modes 4 and 5 are symmetric about the x-axis, and modes 38 and 58 show a von Karman-like structure, as can be seen from Figure 4a,b, the real and imaginary parts of Φ 4 (x), and Figure 4c,d, the real and imaginary parts of Φ 38 (x).

Results from Robust Mode Analysis
The phase dynamics of these five modes are shown in Figure 5. The angular velocities ω 4 (t) and ω 5 (t) are approximately 0.073 and 0.078 radians/frame, respectively, suggesting that they belong to a single flow constituent. The angular velocities ω 38 and ω 58 (t) are approximately 0.055 and 0.051 radians/frame, suggesting that they belong to a second flow constituent. However, the difference in angular velocities between the two pairs of modes (and the symmetry of the first pair and asymmetry of the second pair) implies that the two constituents are distinct. The reproduction of the dynamics of Φ 38 shows that the constituent represents von Karman vortex shedding. The difference between the downstream flow rates of the symmetric and von Karman vortices means that the overall shedding is not periodic.
Finally, we present the reconstructions of the two flow constituents. The top row of Figure 6 contains several snapshots from the original flow. The next two rows show reconstructions of the symmetric and asymmetric constituents at the corresponding time-points. These can be interpreted as symmetric and von Karman vortex shedding. Here, evidence is found for the success of using DMD in separating the flow into distinct constituents. The symmetries of the constituents and the properties of their dynamics (for example, that the angular velocities are not identical) provide the information that is essential in the development of accurate low-dimensional models for the flow.  . Phase dynamics of robust modes for the reacting flow. "Angular velocities" ω 4 ≈ 0.073 radians/frame and ω 5 ≈ 0.078 radians/frame are nearly identical, suggesting that they belong to the same flow constituent. (ω 61 is slightly higher and is not used in the reconstruction.) Similarly, ω 38 ≈ 0.055 radian/frame and ω 58 ≈ 0.051 radians/frame are close, suggesting that they form a second constituent. Successive curves are shifted for clarity. This figure was adapted from the Ref. [31]. Mode 61 has a slightly higher angular velocity and appears less symmetric than modes 4 ad 5. For this work, we have chosen to not include it in the first flow constituent.

Bluff Body Flow Assessed by Discrete Devices: Extended DMD
Here, we once again implement RMA on combustion flows behind a symmetric bluff body, but using only pressure time-series from a small number of transducers. One of the key requirements in DMD is that the space of measurements are sufficiently large that eigenfunctions lie in their span. This is likely to fail when only measurements from a small number of transducers is all that is available. Consequently, a pre-processing step is implemented to extend [24,27,39] the "data-space"; we use Chebyshev polynomials of measurements at each transducer for the task. Robust modes computed from the analysis can be compared with those from Section 3 or derived from other approaches.
A broader goal belies these paradigmatic studies. Damage to engines or industrial reactors are often initiated by localized weakening of material (due to cyclic loading) through onset of microscopic fractures [40,41], which progress slowly until feedbackinduced acceleration precipitates catastrophic failure. New flow patterns may be expected to emerge within a system prior to failure, and if so, techniques to identify them could provide predictive signatures of failure.

The Experiment
Here, measurements of the flow in the experimental arrangement described in Section 3.1 are made using a set of pressure transducers located on at x = −17D, −8D, −4D, 0D, 1D, 5D, 7D, 8D, and at the exit 14D on the bottom surface and at 0D, 1D, and 5D on the top surface, where D = 38.1 mm is the width of the bluff body. Pressure measurements in each transducer were made for 20 s at a rate of 20 kHz. Figure 7 shows the arrangement schematically. The analysis below is for a combustion flow at 440 F of equivalence ratios φ = 0.95, when only symmetric vortex shedding is observed, and at φ = 0.85 when both symmetric and von Karman shedding are observed. The trailing edge of the symmetric bluff body is assigned to be at downstream locations being positive. Note that transducers 7, 8, and 9 are placed opposite to 4, 5, and 6, respectively. A detailed explanation of imaging optics and measurement techniques can be found in the Ref. [42]. Here, we confine our analysis to measurements from transducers 1-9.

Extended DMD and Robust Mode Analysis
Time-series from transducers 1-9 are "extended" synthetically using Chebyshev polynomials. The choice of Chebyshev polynomials is motivated by the fact that they provide finer temporal resolutions while remaining finite. Further, their use is computationally faster than the use of, for example, trigonometric functions. The number of polynomials used for the extension is made as follows: evaluate robust eigenvalues with increasing numbers of Chebyshev polynomials and search for convergence. For the bluff body data, the eigenvalues converge when the number is ∼8. We used 10 Chebyshev polynomials in our analysis. Such extensions, referred to as extended DMD, have been shown to successfully enhance the data space [24,27,39]; the spatial field is effectively enhanced to 90 sites. Robust mode analysis is implemented using five sub-intervals, each of 10,000 time steps (i.e., 500 ms), as distinct, nominally-identical experiments. Robust modes were identified using the criteria outlined before.
Only a single pair of robust modes Λ = −0.002 ± 0.66i and period 10.5 ms, is found in the first experiment, at the equivalence ratio φ = 0.95. Projections of the mean and the standard deviation of the robust mode to the space of the signals from the nine transducers are shown in Figure 8a. Observe that the values of the robust mode at transducers 4, 5 and 6 are nearly identical to those at transducers 7, 8, and 9 respectively. Since the pairs of transducers (4, 7), (5, 8), and (6, 9) are at the identical axial locations, this near-equality indicates the robust mode is up/down symmetric; that is, consistent with symmetric vortex shedding.
The flow at φ = 0.85 is significantly more disordered and contains six robust modes with periods 10.5, 6.3 and 3.0 ms. The last mode, whose projections to the space of signals from the transducers are shown in Figure 8b, does not exhibit the up/down symmetry, suggesting that the mode is not associated with symmetric vortex shedding.
To conclude, what we have shown in this section is that robust mode analysis of the time-series from a set of appropriately placed transducers can be used to extract specific details (in our case, the presence or absence of up/down symmetry) of flow constituents. As already pointed out, these determinations are made directly from data, and do not rely on theoretical considerations.  (4, 7), (5,8) and (6,9) seen in the mode of (a).

Non-Robust Constituents in Turbulent Jet Flow
Turbulent flows are characterized using statistical objects, such as time-delayed correlation functions, structure functions, and spectral energy densities [5][6][7][8]. The final example highlights that statistical features of turbulent fluid flows are closer to theoretical predictions when robust modes are eliminated from the flow fields.
We analyze an axisymmetric jet flow of Reynolds number ≈21,000 measured using particle image velocimetry (PIV). The PIV data permit us to compute velocity distributions at individual sites empirically and to establish the variance-covariance matrix between pairs of points. As shown below, we are thus able to compute structure functions of all orders between two sufficiently far-away points. The agreement between the flow data and an analytical expression is shown to improve when robust modes are eliminated from the flow.

The Experiment
The critical technology that facilitated the resolution of the turbulent flow in both spatial and temporal domains was the development of a long-duration (100-ms), highenergy burst-mode laser capable of producing pairs of pulses at a high frequency [43,44]. It permitted the extraction of two-dimensional velocity measurements through particle image velocimetry (PIV) at a rate of 100 kHz and a resolution of flow features (both spatially and temporally) down to the Taylor micro-scale. For example, at 100 kHz in the present study; ∼8500 snapshots of such velocity-fields were acquired. We report on a turbulent jet flow of Re ≈ 21,000 emanating from a cylindrical inlet of diameter D = 4.7 mm. As it moved downstream, the axisymmetric flow expanded (nearly) linearly. Velocity measurements were made on a grid of 51 × 67 of lattice size 0.268 mm placed symmetrically in an axial plane approximately 13D downstream of the inlet, as shown in Figure 9a. Figure 9b shows a snapshot of the flow velocity field. The PIV measurements were recorded once it was determined that the flow had stabilized and the seed density was sufficiently high. That statistical objects had converged was verified by computing them over different sub-intervals. The time average of the velocity field v(R; t) and its deviation from the mean are denotedv(R) and v (R; t), respectively, where R is measured from the center of the nearboundary. Mean velocity fields and their fluctuations about the mean are found to be neither homogeneous nor isotropic [45].

Robust Mode Analysis
Implementing RMA with δ R = δ I = 1.0 and ∆ = 0.75 gives 12 robust modes in addition to the mean flow, which together contain approximately 16% of the energy. Figure 10 shows the energy in the first 60 modes, illustrating that robust modes are not those with the highest energy. In particular, this observation implies that POD [15][16][17][18][19] or other energy-based modal decompositions cannot be used to identify robust modes. Flow patterns associated with one of the robust modes whose period is T = 1.28 ms is shown in Figure 11. The origin of the robust cyclic modes in a given configuration is difficult to gauge. They may originate from the flow dynamics, details of the experimental setup, such as imperfections in the inlet flow, or from extraneous factors, such as acoustic reflections from laboratory walls. Our goal is to test whether they need to be eliminated from the turbulent flow prior to computing statistical objects. Figure 11. Flow patterns in the pair of robust modes 6 and 7 whose eigenvalues are −0.013 ± 0.489i and whose period is T = 1.28 ms.

Velocity Fluctuations and Structure Functions
Theoretical studies on statistical properties of turbulent flows, derived on the basis of general considerations such as energy transmission between scales, rely on the isotropy and homogeneity of the flow [6, [46][47][48]. Unfortunately, turbulent jet flows are inhomogeneous and do not satisfy these scaling laws [45]. We derive a possible comparison that relies on the empirical observation, illustrated in Figure 12, that point-wise velocity fluctuations in the jet flow, both in axial and radial directions, are approximately Gaussian distributed. We focus on axial velocity fluctuations, v (R; t), which are distributed as Although standard deviations are location-dependent, the near-normality of pointwise fluctuations is found at all sites. The next step is to evaluate the joint distribution of v 1 ≡ v (R; t) and v 2 ≡ v (R ; t) of axial velocity fluctuations at points R and R . Although these velocity fluctuations are expected to be uncorrelated at large distances, this fails to hold for the short duration (∼6 ms) of the experiment. However, the variance-covariance matrix for fluctuations between R and R can be evaluated from the flow data. The joint distribution of v 1 and v 2 is found to be nearly a centered bivariate normal distribution whose form is evaluated from the time-series v 1 and v 2 .
Our goal is to compute structure functions [7,8,46] of all orders using the distribution of the velocity difference (v 1 (t) − v 2 (t)). The first step toward this goal is the diagonalization of the variance-covariance matrix. Denote the eigenvalues and eigenvectors of Σ by Λ 1,2 and u 1,2 . The new variables are distributed on uncorrelated normal distributions with variances Λ 1 and Λ 2 . Next, let the angle between the original basis and the (orthonormal) basis {u 1 , u 2 } be φ. The projections of a point (v 1 , v 2 ) to the new eigen-basis are (v 1 cos φ + v 2 sin φ, v 2 cos φ − v 1 sin φ). It can then shown that the difference with α 1 = cos φ − sin φ and α 2 = − cos φ − sin φ. Correlations between velocity fluctuations at neighboring sites are characterized using the family of structure functions where e ≡ e (R − R) is the unit vector in the direction (R − R). Increasing ζ emphasizes larger deviations from the mean. Structure functions are homogeneous on the symmetry axis, that is, they do not depend on the "origin" R [45], as confirmed empirically. Our computations of structure functions are implemented on a 60 ms (6000-frame) segment of the flow. Using Equation (10), it is found that Any deviation from this expression is due to the slight differences from normality of velocity fluctuations, see Figure 12. Figure 13a shows structure functions for two points separated by 10.7 mm (40 lattice points) along the symmetry axis. The computation is implemented for four such sets R and R . (For these points Σ 11 = 27.16 ± 1.54, Σ 22 = 24.34 ± 0.94, and Σ 12 = 0.03 ± 0.01 1.) The diamonds in Figure 13a are the averages and the error bars are the associated standard deviations of S ζ (R, R ). The data differ from the expression (12), shown by the dashed line.
The corresponding comparison for non-robust constituents of the jet flow is shown in Figure 13b. Now we find excellent agreement. (The dashed line here is computed from the non-robust flow.) This subtle test suggests the need to eliminate robust flow constituents prior to testing statistical theories of turbulent flows.

Discussion and Conclusions
We introduced a model-free methodology to extract, from experimental data, the primary constituents in a flow. The methodology is expected to be particularly useful in studying strongly driven fluid or combustion flows [1,2,[5][6][7][8], where nonlinearities of the underlying flow induce chaotic dynamics [9]. The lack of precise governing equations, absence of symmetries or other simplifying constraints in realistic configurations, and the presence of flow structures over a wide range of scales reinforce the importance of direct analysis of experimental data.
The proposed approach was predicated on the assumption that primary constituents in fluid or combustion flows are common in multiple, nominally identical realizations of an experiment [31][32][33]. Further, it is assumed that these flow constituents are combinations of dynamic modes that are common between the different experimental realizations, the commonality established via proximity of the eigenvalues and eigenfunctions. Flow constituents are constructed by summing over robust dynamic modes with the same phase dynamics. It needs to be emphasized that, while evaluations of dynamic modes are often fickle, depending sensitively on the precise collection of frames and sites selected for the analysis [29,34,35], robust modes do not suffer this fate.
Robust cyclic facets of a flow may consist of measurement errors, and factors (acoustic reflections from laboratory walls, etc.) extraneous to the system studied. In addition, small-scale features like eddies may be assigned a non-robust designation as well when their appearance, location, and dynamics depend on chaotic aspects of the global flow, that is, they are not cyclical.
We illustrated the use of robust mode analysis (RMA) in a lean combustion flow behind a v-shaped symmetric bluff body, where vortex shedding takes both symmetric and von Karman forms. RMA yielded five robust modes from which two flow constituents were constructed using their phase velocities. One constituent appears to describe symmetric vortex shedding, and the other von Karman shedding. RMA analysis in rich mixtures only yield the first constituent [31]. The approach outlined here was shown to be capable of extracting known flow constituents in shear coaxial jets with and without transverse acoustic driving [32] and in stable and unstable turbulent swirling reacting flows [33].
Difficulties in extracting spatio-temporal flow fields in working engines and reactors motivate us to inquire if RMA can be implemented using time-series from a limited number of point sources such as pressure transducers or temperature sensors, which are implanted in many such systems. We used time-series from nine transducers to study, once again, reacting flows behind a bluff body in lean and rich mixtures. The analysis required synthetically extending the data space [24,27,39,49], in our case, using Chebyshev polynomials. Once again, RMA yielded a single robust mode in rich mixtures and two in lean mixtures. Moreover, the values of one of the robust modes was (nearly) identical at symmetrically-placed transducers, suggesting that the mode represents symmetric vortex shedding.
The final issue we addressed pertains to turbulent flows, whose characterization requires statical objects [5][6][7][8]. We analyzed velocity field data in an axisymmetric turbulent jet flow of Reynolds number ≈21,000, which was shown to contain 13 robust modes. The question we raised is whether statistical characterizations apply to the full flow field or only to its non-robust constituents. We proposed a test that relied on an empirical observation that joint velocity fluctuations were distributed on centered bivariate normal distributions, from which the behavior of structure functions of all orders can be derived. The expression shows closer agreement with non-robust flow fields than with the full flow, suggesting that statistical comparisons of turbulent flows may be better tested on non-robust flow constituents.
The ability to extract robust flow constituents directly from experimental data and to identify their metamorphoses can have broad ranging applications, especially as a signature of impending catastrophic failure in working engines or industrial reactors. It would prove especially useful if the analysis can be implemented using data from easily embedded pressure transducers or thermocouples. Further studies are needed to test if observable changes in flow structures due to progressive material degradation/wear-and-tear precede catastrophic failure of a system.

Data Availability Statement:
The spatio-temporal flow velocity data that support the findings of this study, in Matlab format, have been deposited in "figshare.com", with the accession code "10.6084/m9.figshare. 12789926." Grid points x and y of the 51 × 67 lattice and the velocity fields u(x, y, t) and v(x, y, t) in the x and y directions are given for jet flows at Reynolds number 21000.