Time-Volume Estimation of Velocity Fields From Non-Synchronous Planar Measurements Using Linear Stochastic Estimation

The work presented in this paper combines multiple non-synchronous planar measurements to reconstruct an estimate of a synchronous, instantaneous flow field of the whole measurement set. Temporal information is retained through the linear stochastic estimation (LSE) technique. The technique is described, applied and validated with a simplified combustor and FSN geometry flow for which 3component, 3-dimensional (3C3D) flow information is available. Using 3C3D data set, multiple virtual ‘planes’ may be extracted to emulate single planar PIV measurements and produce the correlations required for LSE. In 1 Corresponding author ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 2 this example, multiple parallel planes are synchronised with a single perpendicular plane that intersects each of them. As the underlying data set is known it therefore can be directly compared to the estimated velocity field for validation purposes. The work shows that when the input time-resolved planar velocity measurements are first POD (proper orthogonal decomposition) filtered, high correlation between the estimations and the validation velocity volumes are possible. This results in estimated full volume velocity distributions which are available at the same time instance as the input field – i.e. a time resolved velocity estimation at the frequency of the single input plane. While 3C3D information is used in the presented work, this is necessary only for validation; in true application planar technique would be used. The study concludes that provided the number of sensors used for input LSE exceeds the number of POD modes used for pre-filtering, it is possible to achieve correlation greater than 99%.


INTRODUCTION
The understanding of downstream velocity behavior is crucial to the design and development of many flow devices such as fuel injector nozzles. Great effort both in numerical modelling and experimental measurements is involved, with the state of the art in both areas becoming increasingly complex and expensive as more detailed information is sought.
The work presented in this paper will focus on the characterisation of the downstream flow of a radially-fed single stream air swirler which provides a generic swirling flow case that is applicable to fuel swirl nozzles (FSN) and combustor flows in gas turbine engines. A common approach in this application is to make use of techniques such as large eddy simulation (LES), for example, work by Dunham et. al. [2]. However, mainly due to the significant computational cost, there is an increasing use of statistical techniques. For example, Treleaven et. al. [3,4] uses proper orthogonal ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 3 decomposition (POD) to identify and study acoustic instabilities in a lean-burn FSN. This use of POD shows that not only can these techniques reduce processing time -but they may also increase the output of useful information from a data set; be it experimental or computational.
In the experimental measurement of flow fields, planar particle image velocimetry (PIV) has been commonplace for many years and it is possible to obtain high yield, high quality -both in terms of spatial and temporal resolution -data in relatively Extensions of the PIV technique such as holographic PIV and tomographic PIV allow 3C3D volume flow measurements to be carried out. Holographic PIV is generally considered the first true volumetric PIV technique [5][6][7] but its use is typically limited by its complex set-up [8] and fundamental imaging issues [9]. Spencer et. al. [1]  sufficiently illuminating a volume (rather than a sheet), positioned around a highly optically-accessible test rig. Lastly, the authors suggest that tomographic PIV requires greater smoothing of the data and evaluate an error in comparison to the SPIV measurements. Despite these difficulties, the technique is increasingly popular for those applications that afford good optical access proving the appetite for 3C3D velocity information and additional flow characteristics that may be obtained. Raffell et. al. [10] describes an overview of current 3D PIV techniques.
As computational / processing power continues to improve, particle tracking techniques are becoming more accessible to provide 3C3D velocity measurements. A growing body of literature makes use of a 3D Lagrangian particle tracking velocimetry (PTV). Examples such as shake-the-box (STB) [11,12] promises high quality 4D-PTV (spatial and time resolved) over a wide range of applications, which helps reduce some of the optical access requirements. However, as with techniques already discussed, there remains a significant hardware investment and increase in optical access requirement which may just not be possible in some applications such as gas turbine combustors.
This trade-off between representative geometry and optical access requirements is well established in the experimental community and has led to increased effort to make use of alternative applications of planar PIV measurements. Scanning PIV [13] can provide a quasi 3C3D measurement using a light sheet which is stepped through the depth of the volume in time-steps which are sufficiently smaller than the PIV interframe ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 5 time (i.e. volume sample frequency). More recently this type of system has been used to generate a moving reference plane allowing the study of vortex dynamics [14].
An alternative approach has been to make use of statistical techniques and multi-plane (asynchronous) 2(or3)C2D measurements to extract information or even reconstruct a 'measurement' volume. Volpe et. al. [15] demonstrate how a series of ten intersecting planes arranged five vertically and five horizontally may be used to study the wake topology behind a square-back Ahmed body (applicable to the automotive external aerodynamics). To understand the bi-stable wake behind the body, the temporal average of each stable state is considered by calculating two conditional averages. These are informed by a pressure signal that indicates which side of the body the wake is currently on, allowing a left-right distinction to be made. Further application of this technique is presented by Perry et. al. [16] where the authors compare the information from four 2D planes in a similar arrangement to tomographic PIV measurements at the same conditions. In that work, the authors also carry out POD analysis on the planar data to reveal bi-modal behavior; which is further revealed by investigating the tomographic PIV measurements.
In the swirling flows applicable to FSN and gas turbine combustors, the time average of planar data would not sufficiently describe the behavior of the highly three- demonstrate the revealing of thermoacoustic oscillations that would not be possible with ensemble averaging of the planes. However, unsteady features are still not fully revealed using this approach since a phase-averaged cycle will not resolve cycle-to-cycle variations.
In conclusion, it is highly desirable to obtain time resolved volumetric velocity measurements (3C3D). Whilst techniques do exist and are fast producing high-yield, high-quality measurements, the application of these is still restricted to sufficiently optically accessible test rigs and require significant hardware investment over more traditional PIV methods. There exist a handful of statistical techniques that allow a quasi-3C3D reconstruction, some including phase information. But these do not provide a true temporal resolution.
The work presented in this paper introduces the statistical technique, linear stochastic estimation (LSE), to application of volumetric velocity reconstruction using time resolved planar PIV (3C2D) as an input. It is an improvement over phase averaging as it does not require an explicit phase indicator on which to conditionally average the data. In this way LSE offers a way to improve on understanding the cyclic variation of rotating structures rather than recreating a single cycle average and importantly it is not limited to rotating quasi-periodic flows.

LINEAR STOCHASTIC ESTIMATION
The technique of stochastic estimation, introduced by Adrian [19][20][21], proposes that the velocity at one location can be stochastically estimated from conditionally correlated data at an unconditional source at the same instance in time, i.e. two-point ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 7 statistics. In the field of fluid mechanics, this presents a promising assistance to the experimental acquisition of velocity data and has been steadily developed in the time since its introduction by Adrian. The effect of higher order, quadratic stochastic estimation, QSE was soon after investigated, but found to give nearly identical estimations to the linear case, LSE [22,23].
Often, the sources for the estimation are from microphone (for example [24]) or pressure sensor arrays [25][26][27][28] as these measurements typically are less intrusive and can provide increased temporal resolution for a given experimental effort. However, disadvantages to the use of these typically wall-based sensors are described by Arnault et. al. [29] who show the limitations in the estimation of smaller turbulent structures which may not correlate well with these measurements. The study does however progress to show that it is possible to optimize the sensor locations to improve correlation using an algorithm developed by Muradore et. al. [30].
There are examples in literature that demonstrate the use of velocity measurements as the input sources for LSE. Work by Kerhervè et. al. [31] uses a hotwire rake with high temporal resolution of 30kHz correlated with high spatial resolution, low temporal resolution PIV measurements, generating a high spatial and temporal resolution estimation driven by the hot-wire.
The LSE technique is described mathematically by Equations 1 and 2 where u is velocity, p represents the sensors (or sources), where there are NR sensors and a is the correlation matrix. This terminology will be used through the presented work. This approach may be further extended into the 3rd spatial dimension to generate an estimation of 3C3D velocity. Figure 1 shows the velocity component contour on the YZ master plane with several example XY slave planes superimposed ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 9 (faded). Sensors are highlighted for the upper most XY plane showing that they should lie on the intersection between each slave plane and the master plane.
The steps involved in estimation of 3C3D velocity (in the presented orientation) is then as follows: 1.
Capture multiple slave XY planes (3C2D) at a range of (Z=const) values of interest. These should be long time averages to ensure convergence of the 2-point statistics calculated in that plane.

2.
Taking each slave plane independently, apply Equation 1. First take the multiple ( ) values of from along a (X=const) line, where this X location remains consistent between planes. Knowing ( ) in this plane then allows the respective correlation matrix, to be determined as the unknown. The calculation is carried out for each velocity component independently.

3.
For those velocity conditions to then be estimated using LSE, capture a single master YZ plane at the (X=const) location from the previous step. Now in this data, each (X=const) line may be combined with the corresponding a matrix to generate a velocity estimate in the intersecting XY plane at that value of Z.

Applying the technique in this manner means that all estimated velocity fields in
XY planes will now be statistically synchronized to the same time instance as the master YZ plane. Any estimations of these planes may therefore be resolved at the frequency of the captured master YZ plane, even if this is greater than the frequency of the captured slave XY planes as was demonstrated by Kerhervè et. al. [31]. In fact, it is preferable to calculate the correlation matrix over a long time period to ensure statistical convergence.
In the application described in this work the orientation of the frames was arbitrarily chosen, but the technique could equally as well be applied with a different choice of slave / master orientation.

PROPER ORTHOGONAL DECOMPOSITION
It is quite common for the LSE technique to be used in combination with another It is widely accepted that the spatial early modes -when ordered by energy as is the norm -contains the coherent motions and the latter modes contain the turbulence.
There exists a great number of techniques for choosing the 'cut-off' between these two sets. Many studies use a somewhat arbitrary value of 90% energy for example in the work by Graftieaux et. al. [34]. Butcher et. al. [35] shows how cross-correlation of the spatial modes may be used to identify the cut-off. In that work, the cut-off is used to ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 11 determine which POD modes are representative of coherent motions which vary cyclically. The selection of the number of POD used for the reconstruction will have a significant impact on the 'smoothing' effect of the filtering.
Approaches to the combination of POD and LSE fall into two categoriesprimarily based on whether LSE is applied to a POD filtered data set, or if the estimation is of the POD modes themselves. Podvin et. al. [36] demonstrates how coarse measurements of a single velocity component may be used in this manner to recreate POD modes and therefore reproduce the large-scale motions in the flow field on a finer mesh grid.

POD-LSE coupling
It is shown that when the LSE technique is applied to POD filtered data, there exists a coupling effect between the number of POD modes used for the initial filtering / reconstruction and the number of sensors used for the estimation of velocity fields [37].
They found that provided the number of independent sensors was sufficiently greater than the number of POD modes used for the filtering of the underlying velocity fields, a small error (<5% error vector magnitude) may be achieved. Figure 2 summarizes the relationship found between the number of sensors required for a desired maximum permissible error (up to a maximum of 100 sensors). In that work, the location of the multiple 'master points', , were randomly selected rather than along a straight line.
This mitigated the potential detrimental impact of selecting sensors in very close proximity because closely aligned locations of master points may not provide additional/independent statistical information. Therefore, it would follow that the work presented here may not follow this exact trend, given that the points selected are next to each other, as was required along the lines of intersection between the planes. Hence it is possible more sensors would be required than Figure 2 indicates.  controlled via a throttle valve restriction and is set to provide a nozzle exit Reynolds number of 7x10 4 , where Reynolds number is defined using the swirler exit diameter.

ACQUISITION OF TEST CASE DATA
Using the geometrical exit area, a bulk mean velocity of 1.70 ms -1 is observed.
Illumination was provided by a Litron LDY303HE laser with each of the two cavities set to 1kHz. This was directed through LaVision volume optics to provide a collimated illuminated volume of cross-section 70mm x 25mm. An aperture was applied to ensure sharp edges to the illuminated volume. Images were acquired using four  taken through the domain in the YZ and XY orientations respectively -each taken at mid-length. In general, these planes will be used for discussion and illustration throughout the remainder of the paper, as it is difficult to efficiently and effectively represent 3D data in print, particularly comparisons between 2 sets of 3D data.
However, all calculations will be carried out on the full volume.

DEFINITION OF PLANES
Using the presented test case data, it is necessary to define which data will be used as slave fields and which will be used as the master plane (to drive the subsequent estimations). Firstly, of the 1023 PIV fields (time-steps) captured, let 1000 be used for the generation of the correlation matrices and contain slave planes. The remaining 23 fields shall not be used in the generation of the correlation matrices -and only the master plane from this data set shall be extracted, containing all of the intersection points (sensors) velocity. In addition, the remaining volume of this set shall later be used to allow comparison to the generated estimations.

Slave planes
Given the orientation depicted in Figure 1, the XY planes shall provide slave planes. The intersection line will be defined at (x = 0) of this plane; although it should be To calculate the correlation matrix for each slave plane, the plane velocity matrix (68x64x1000) is collapsed to a (4352x1000) matrix. This may then be used together with the sensor velocities for that plane (a matrix of 64x1000) in Equation 1 to generate the correlation matrix for that plane.
Each slave fields' correlation matrix effectively represents the correlation between each spatial location in the plane and every sensor. The dimensions of each correlation matrix are therefore equal to the total number of spatial locations in the plane by the number of sensors. Figure 8 presents an example of a spatial correlation matrix for the Ux velocity component. There are a few noteworthy features of this plot.
Firstly, one may expect that the correlation should be similar to the spatial correlation about that point, with a peak at the same location. However, there exist several peaks of this (normalized) distribution, none of which are near. This is due to the velocity field correlation being an ensemble of contributions from several sensors (64 in this case).
However, if there were only one sensor used, then the distribution would mimic a spatial correlation.
For the presented example application there are a total of 21 slave planes (i.e. equal to the Z dimension), each with three correlation matrices; one for each velocity component. However, in practice the reconstructions could be achieved with fewer (or greater) slave planes, depending on the permissible reliance on interpolation between estimated planes.

Master plane
Following the same orientation, the master plane shall be defined as a midlength (in the X-direction) YZ plane. Each (X=const) line (examples presented in Figure 9) is taken from the master plane set and used with the appropriate correlation matrix and Equation 2 to generate the volume velocity field.

ASSESSMENT OF GENERATED ESTIMATIONS AND DISCUSSION
The process is carried out as described using the MATLAB program for each of the 23 validation master planes (time-steps) using 60 POD mode pre-filtering. Firstly, the ensemble average of the reconstruction of the estimated fields is presented in Figure 10 and can be compared to the ensemble average for the PIV data presented earlier ( Figure   5).
It can be seen in Figure 10 that there are slight differences to the PIV data posted in Figure 5. In fact, the validation shows agreement between the estimation and the validation data set. Differences arise as Figure 10 is the ensemble of only the estimations from the 23 master planes extracted from the validation set; whereas Figure 5 shows the ensemble of the 1000 slave velocity fields. Further exploration of the quality of the estimation is given in the remainder of this paper.
The purpose of carrying out the LSE technique -rather than the earlier discussed ensemble average reconstructions -is to provide time-volume velocity fields. And as ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 18 explained earlier, using fewer POD modes in the master field reconstruction than number of LSE sensors leads to less than 1% error in the estimation. The following discussion is regarding 60 POD mode reconstruction; the reader is reminded that there are a total of 64 LSE sensors at each intersection of planes, and therefore this condition is satisfied. By selecting an extract of an estimated XY plane from this data set; a direct comparison can be made to the equivalent field extracted from the PIV validation set (at the same time instance). This is presented in Figure 11a & b respectively.
Whilst the comparison is presented in this way (Figure 11), the two fields have very small disparity between them, << 1%. This confirms the finding in literature [37] regarding POD pre-filtering.
As part of the investigation, all three velocity components, across all 23 validation fields were compared to their validation measurement counterparts. A similar finding was observed in every condition. Therefore, the study concludes that the application of this technique may effectively estimate time-volume velocities from a single master plane measurement set.
Further, an investigation is carried out to establish the level of error introduced to the estimation if the number of POD modes is increased beyond the number of LSE sensors. For this purpose, five conditions were investigated, based on the number of POD modes used for pre-filtering: 70, 100, 130, 160 and no POD pre-filtering.
Considering the first mentioned condition, 70 modes, they may be directly compared by overlaying the data-set -this was not previously possible due to there being no discernable discrepancy between the two. Figure 12 shows how now that the ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 19 number of sensors is greater; albeit by only around 10%, there is some discrepancy between the estimation and the validation data. However, note the scale change in Figure 12 required before this difference is visible -showing there is still a strong similarity between the two velocity fields. This would be expected to become more significant as less pre-filtering (i.e. more POD modes) are included in the master planes.
To efficiently assess the accuracy of the estimations, a quantitative approach is taken. Two velocity fields may be directly compared by a correlation between the two.
Each of the pre-filtered cases are estimated and the correlation between estimations and validation time-volume calculated. Correlations for these conditions are shown in Figure 13.
The correlations in Figure 13 show that reasonable estimations may be made in the cases where the number of POD modes used for pre-filtering exceeds the number of LSE sensors, but this is dependent on what the acceptable level of error is deemed to be. Also included is the condition with no POD filtering (represented by the dashed line).
Even in the worst-case scenario, where there is no POD pre-filtering of the data, greater than 95% correlation exists between the estimated fields and the validation fields. It is ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 20 suspected that while the no-filtering cases did not utilize POD, it was necessary in their capture for gaussian filtering of the data due to noise on the tomographic measurements which may have impacted this result -leading to an unexpectedly high correlation in the no POD filtering condition. A side effect of this initial filtering (applied to all data in previous works [1]) will have been to remove some of the smaller structures that would otherwise have been represented by the higher order modes.
Therefore, effectively less information is required to estimate the velocity fields than may otherwise be the case.

METHOD CRITIQUE AND FUTURE WORK
This paper provides a validation of the LSE technique using a 3D-3C tomographic PIV data set where direct measurements are available to test the accuracy of the estimation results. This has been done carefully to ensure a priori knowledge available in this data set has not been used to assist the methodology. The origins, quality and uncertainty in this data has been previously described in [1]. Based on the divergence of the velocity field which should be zero, in this incompressible flow, it was estimated the instantaneous velocity uncertainty from the tomographic PIV was better than ±10% to 95% confidence.
The tomographic set up required for this validation does require a high degree of optical access. However, the technique is aimed at enabling 2D planar 2or3C PIV data to be volume reconstructed on an instantaneous basis. Without such an approach only time or phase averaged statistics can be meaningfully reconstructed into a volumetric picture. The proposed LSE method is aimed for use with a standard mono or stereo PIV ASME Journal of Engineering for Gas Turbines and Power GTP-19-1350 Butcher 21 system. Thus, the optical access to employ the technique would ordinarily be commensurate with a standard mono or stereo PIV systems with sufficient flexibility to traverse the measurement plane through the volume of interest.
Alignment of the measurement planes is an issue that has not been directly addressed in this validation. For the approach to work, the line of intersection in the master and slave planes will need to be identified with some degree of accuracy, otherwise the sensors from the master and slave plane will not be at identical locations.
Additionally, there may be some differences in velocity uncertainty from the two planes through alternate calibration of the PIV planes. There are solutions to these issues and is the subject of current work.
Where there is not good PIV vector yield (less than ~95% first choice vectors) then additional steps are required to ensure the LSE methodology provides mathematically robust estimates. Again, this is the subject of future work, but it would be expected that uncertainty of velocity from the stereo PIV would be reduced to typically around half of that of this tomographic PIV data, and this should help to improve the LSE results.
It is seen in this validation exercise that there can be very good correlation in the estimated and actual instantaneous velocity fields. This is unlikely to be as good if large scale coherent structures do not exist in the turbulence field. LSE can work with small scale turbulence, but sensor locations would need to be more uniformly distributed around the slave plane so that their typical spacing was not too far from the integral length scales that exist in the flow. Where sensors come from a line of intersection between two planes the user would not have the luxury of choosing the required sensor distribution and this application of LSE would be less valuable. However, there would be little to gain of instantaneous volume reconstruction as the flow field becomes more and more self-similar throughout, as in the case of near-homogeneous turbulence.
As a point in the estimated flow field becomes more remote from the sensor locations located on the plane of intersection then the magnitude of the 2-point correlations becomes smaller and the quality of the reconstruction would be expected to become poorer. This is demonstrated in Figure 14.
The profiles in Figure 14 are along two secant lines through a of the near circular swirl cone (z/Ds=0.02). They run parallel to the diametral line on which the sensor points are taken from the master plane. In Figure 14 the profile is close to the diametral line. In Figure b the profile is at about a quarter diameter further away. The two peaks in the turbulent profile associated with the swirl cone and the inner shear layer helical vortex structures are thus thickened on the far field chord line as it cuts the swirl cone obliquely. In this flow where the turbulence field has a strong coherence to it, it can be seen the LSE estimation captures the flow features, and the higher order statistics associated with them, well. The quality of the second order statistic degrades as the LSE estimated flow field becomes more remote from the sensor locations. In the example above the tke is up to 40% lower using LSE directly. However, using POD filtered data prior to applying the LSE reconstruction as advocated in this paper, then in the example presented, the turbulence statistics are captured with as good fidelity as allowed by the POD reconstruction, where that POD reconstruction is generated by as many spatial modes as there are sensor locations available in the master plane.

CONCLUSION
A method is developed and presented allowing a time-volume velocity estimation to be generated based on a set of sensors of planar measurements. To demonstrate the methodology, it has been applied to a simplified gas turbine combustor and FSN flow for which existing time-volume velocity measurements were available to allow for effective validation. Although the full volume information was available, virtual planes have been extracted and used to emulate planar measurement inputs. The advantage of this approach is a priori knowledge of the flow field allows a direct comparison between the estimations and the measured velocities.
During the first step of the technique, akin to a calibration procedure, the spatial resolution of the final estimations is determined by the number of slave fields to be captured, and subsequent correlations generated against the master plane intersection.
The temporal resolution of the estimations is then determined by the frequency of the master plane measurements under test conditions. This highlights a potential further use of the technique; if high-speed planar information is available, then the full volume estimation will also be at this frequency.
The link between POD filtering and LSE reconstructions suggested in literature has been investigated in this application. The findings generally agree strongly, but there is some evidence that other PIV filtering techniques have smoothed the data and removed some of the smaller structures of the flow.
Finally, the application of this technique does not necessarily require new data sets, it may be retrospectively applied to any data where intersecting non-synchronised measurements have been carried out. Indeed, these may generally be velocity fields from PIV data but the mathematical approach would work well with other quantities such as pressure or scalar fields from LIF in both non-reacting or reacting flow fields.