A continuously scanning spatiotemporal averaging method for obtaining volumetric mean flow measurements with stereoscopic PIV

In this study, a variation of the stacked stereoscopic PIV technique is proposed to perform fully volumetric (3-dimensions, 3-components) measurements of average flow fields within a single experiment through the usage of an automated traversing system that continuously scans the SPIV light sheet over a linear path. The simultaneous measurement of the traverse location and the laser Q-switch pulse enables the automated assignment of instantaneous PIV fields to known physical coordinates, enabling spatiotemporal averaging in post-processing to obtain volumetric measurements of a flow field. This method provides a trade-off between spatial resolution of the volume measurements and statistical convergence of the spatiotemporal averages, enabling volumetric measurements under challenging experimental conditions where only stereoscopic PIV is viable. A comparison with the more traditional temporal averaging method and planar PIV is presented to demonstrate the capabilities and limitations of this technique in realistic, challenging experimental setups. It is found that the spatiotemporal averaging convergence behavior differs slightly from the traditional temporal averaging for the wake of a bluff body model, however relative errors lower than two standard deviations can still be attained. Thus, this technique presents a viable alternative for rapid 3D reconstruction of averaged flow fields that can provide invaluable insight of various flow topologies.


Introduction
Experimental volumetric flow field data can be extremely valuable to understand the flow physics of complex aerodynamic problems and to provide databases to validate computational models. Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) remain the most competitive techniques for acquiring quantitative measurements of flow fields, and several variations of PIV and PTV have been developed over the decades to increase the amount of information about the flow gathered in a single experiment. The leading techniques on volumetric flow measurement are tomographic PIV (Scarano 2013) and 3D PTV (Maas et al. 1993), and have been highly influential since the first introduction of tomographic PIV by Elsinga et al. (2006), enabling the volumetric visualization of very complex flow fields under a variety of conditions. Implementing tomographic PIV is very challenging and requires a large amount of expertise, limiting its general application. Although a complex web of trade-offs must be considered when performing tomographic PIV, the key physical constraint that prevents the probing of large volumes with tomographic PIV (or any light-volume particle based technique) is attaining a sufficiently large particle signal in the imaging cameras to overcome sensor noise, which is required to achieve a high-quality volume reconstruction. This is particularly the case in high-speed air flows, where the seed particle size cannot be increased to achieve a larger Mie scattering signal, as the excessive inertia that larger particles contain prevents them from following the flow in regions with large velocity gradients, such as shocks and vortices. Assuming a cuboid-shaped volume with the thickness t being its shortest side, a camera corrected by a Scheimpflug adaptor requires an f-number N f which is proportional to t (Scarano 2013) to ensure there is sufficient depth of field to maintain all particles within the volume in focus. Since the light collected by a stopped lens with an f-number N f is proportional to N −2 f , 1 3 56 Page 2 of 14 it imposes a light collection constraint that is proportional to t −2 . Furthermore, the light sheet thickness needs to be increased to illuminate the entire volume, meaning the irradiance produced by the laser in the illuminated volume is also proportional to t −1 , therefore reducing the particle brightness by the same amount. The two constraints combined culminate in a particle signal that is proportional to t −3 , which is a significant limitation. Such constraints would need to be compensated by either (1) more sensitive cameras, (2) laser pulses with larger energy, or (3) larger particles, which greatly limit the application of tomographic PIV to a subset of the experiments that can be performed with more forgiving PIV techniques such as 2D-2C (two-dimensional, two-component, or planar) or 2D-3C (two-dimensional, three-component, or stereoscopic) PIV. A simple numerical example with experimental dimensions typically encountered in small-scale high speed wind tunnel tests illustrates this limitation: A light sheet used in stereoscopic PIV in aerodynamic applications is typically t ∼ 2 mm in thickness. Assuming a volume of t = 80 mm needs to be visualized to understand the flow topology of interest, the thickness t would need to be increased by a factor of 40, which would reduce the particle intensity seen by a tomographic PIV setup by a factor of 40 3 =64,000. Thus, a trivial setup for stereoscopic PIV may become incredibly challenging, if not impossible, for a tomographic PIV system of same capabilities.
A great advantage of tomographic PIV and 3D-PTV is the fact the particles in the volume are illuminated and recorded simultaneously. This enables reconstruction of three-dimensional velocity fields for each instantaneous snapshot taken. In some applications, knowing the instantaneous velocity vectors may not be strictly necessary, and measuring the first and second order statistics in the volume may be sufficient to elucidate information about the flow field. By relaxing this constraint, previous researchers (Ostermann et al. 2016;Ostermann et al. 2017;Zigunov et al. 2022a;Meyer et al. 2002;Nakagawa et al. 2016;Cardano et al. 2008;Sellappan et al. 2018;Kumar et al. 2022) were able to reconstruct volumetric flow fields for important flow problems through the combination of multiple stereoscopic PIV (SPIV) planes in a desired volume. The approach utilizes a simple SPIV setup attached to a traversing mechanism, which can move the light sheet and cameras to take different slices of the flow field. At each slice, the traversing mechanism takes a number of images necessary to acquire a converged average to a desired uncertainty threshold. The traversing mechanism makes multiple stops at different cross-plane (z) coordinates, enabling the acquisition of average 2D-3C velocity fields that are later stitched to form a volumetric mean velocity field. The movement of cameras, laser and laser optics are simultaneous, either through the mechanical constraint of the acquisition assembly or through the synchronization of multiple traverse mechanisms. This constrained movement eliminates the necessity for re-focusing the cameras and permits the usage of a single calibration for all planes acquired, enabling the acquisition of multiple z planes. It is important to note, however, that the calibration can only be assumed to be constant if the test section walls are planar and the movement of the traverse mechanism is parallel to the walls the cameras are looking through, especially in the case of experiments in water channels or with other fluids with high index of refraction. A large number of images ( ∼100-500) for each plane is generally utilized to obtain a converged average, which is currently the most constraining limitation in the resolution of this type of technique. The experimental time required to complete this type of volumetric measurement is typically very long, as a non-negligible amount of time is required to acquire and download the images to the acquisition computer, typically of the order of several minutes per plane.
Other noteworthy efforts to perform volumetric reconstructions of flows include the usage of a rotating/ scanning mirror (Brücker 1995;Kitzhofer et al. 2009;Hess et al. 2011;Partridge et al. 2019) or a rotating prism (Hoyer et al. 2005) to deflect the laser light sheet such that it scans a volume; and then utilizing a sufficiently fast camera to acquire the PIV double-frame images. Each image in the sequence then corresponds to a different position of the laser sheet, allowing for a combination of the PIV vector fields into a volume. However, this type of scanning light sheet technique still suffers from the depth of field ( t 2 ) constraint previously discussed, as the camera would have to keep the particles in focus throughout the scan. This means particle illumination intensity would likely still be a problem for large volumes and small particles. Efforts to overcome this depth of field limitation, such as speaker-actuated lenses to scan the field of view (Chaves and Brucker 2010) have been considered in the past, however an integrated approach with synchronized actuation on the laser optics and camera lenses has not yet been attempted to the knowledge of the authors.
Recently, Rousseau and Ancey (2020) demonstrated that a planar PIV system mounted on a scanning traverse can be used to reconstruct the average 3D-2C flow over a bed of randomly packed glass spheres. They obtained a reconstruction of the 3D mean flow volume using spatiotemporal averaging of the planar PIV planes obtained, as the traverse is continuously scanning the volume instead of stopping at each plane to acquire a converged average.
In this work, further details and expansion of the continuously scanning PIV technique are considered. By using a stereoscopic PIV system for the acquisition of the continuously scanning plane, the reconstruction of the average 3D-3C volumetric flow field is possible by spatiotemporally averaging the SPIV planes recorded. Furthermore, a measurement of the location of the traverse for every plane is also proposed to ensure all instantaneous SPIV planes are correctly labeled, greatly improving the uncertainty in the location of the relative position of the instantaneous planes, as in practical experiments the acceleration of the traverse can distort the resulting volumetric flow field. The continuous scanning of the light sheet enables the rapid reconstruction of the average 3D flow field volume in a very short experimental campaign, enabling a trade-off between statistical convergence and experimental run time, which is very useful when experimental run time is limited such as in blow-down wind tunnels. Further exploration into the uncertainty and bias errors introduced when performing such a scan is presented to further understand the limitations of this technique. Three example flow fields of interest to the experimental fluid dynamics community are presented, demonstrating the capacity of this technique to provide novel insight on the structure of complex flow fields.

Hardware concept and timing
The general concept of the scanning stereoscopic PIV (S-SPIV) system utilized in the experiments described herein is presented in Fig. 1a. A traditional stereoscopic PIV setup including laser and optical components to generate a light sheet is mounted on a one-dimensional, linear traversing mechanism. In the specific setups utilized in this work, a Parker 506121-SLH traverse was used. A structure is also mounted on the traverse to allow the two SPIV cameras (LaVision sCMOS) to be mounted where required. The cameras and laser are then rigidly linked such that when the traverse moves, all optical components move. If the mechanical setup is such that the structure will not be sufficiently rigid, the cameras can be mounted in independent traverses, which would then require a means to tightly synchronize their motion. This possibility is not explored in this study.
The fact the traverse mechanism and cameras are mechanically linked together means that the laser sheet will remain focused through the motion of the traverse mechanism. This key feature enables a single focusing step and a single calibration (and posterior self-calibration) procedure that represents all z positions of the traverse, permitting the acquisition of a very large number of z planes.
In this experimental concept, the laser/camera system is traversed continuously throughout the experiment, after experimental conditions have been established. The laser sheet scans a predefined linear path spanning the coordinates z start and z end over a predefined time T. Assuming a laser firing rate r (double-frames per second), the distance between consecutive PIV planes is: This is illustrated in Fig. 1c. Evidently, each PIV plane is an instantaneous snapshot of the velocity field at a slightly different z coordinate. In order to compute a mean velocity field at an arbitrary center z coordinate, one needs to average a number n planes of fields. Under this concept, a spatiotemporal averaging method is defined in the z spatial direction, over the distance Δz: A trade-off between convergence of the statistical quantities ( n planes ) and effective spatial resolution (Δz) is then established. For a fixed data set, one can use more n planes for computing the averages and improve the statistical convergence of the average field calculated, however the result will be a "blurred" version of the true flow field, over a larger Δz .
(1) Δz = z end − z start rT (2) Δz = n planes Δz However, it is important to consider that SPIV is already a measurement that suffers from spatial averaging of the flow field through the interrogation of a finite-thickness laser sheet over a finite time. Thus, flow features smaller than the laser sheet thickness already are unable to be fully resolved (Raffel et al. 2018). With a sufficiently large experimental time T and firing rate r, Δz may be of order t or smaller for an acceptable n planes such that a converged mean velocity field can be attained within an uncertainty bound. In other situations where run times are very short, a partially-converged average flow field may be acceptable when it provides information about the flow that would not be available otherwise.
This concept is applied in three challenging experimental problems at different facilities at the Florida Center for Advanced Aeropropulsion (FCAAP), located in FAMU-FSU College of Engineering: After aligning the laser sheet and calibrating the cameras in this traversing setup, it is possible to position the light sheet at any z coordinate and acquire SPIV measurements at the corresponding plane. To know the precise position of the traverse system as it is in motion, a linear potentiometer is installed in the traversing mechanism. A data acquisition (DAQ) device then simultaneously acquires the Q-switch pulse train from the first laser pulse (Q1) and the linear potentiometer voltage, as shown in Fig. 1b. To ensure the DAQ device, which sampled the signals at 5 kS/s, does not miss the short Q-switch pulse (3 s long in our experiment), a Stanford Research DG535 delay generator (DG) is set to trigger from the Q1 pulse coming from the programmable timing unit (PTU), lengthening the pulse to 1 ms. The specific length of the pulses is not crucial, as long as it does not exceed the laser firing period (1/15 s in our experiment) and it is longer than a few DAQ sampling periods such that a rising edge can be reliably detected in post-processing.
As displayed in Fig. 1b, the simultaneously acquired signals from the DG and the linear potentiometer allow for accurate assignment of each PIV plane to a specific z location through the detection of the rising edges from the DG signal. Thus, in the event that the PIV system misses a frame due to computer processing limitations or camera download bandwidth issues, the location of the cameras and laser sheet can still be found accurately throughout the experiment.

Dynamic masks
In typical PIV applications, it is customary for the user to apply a manually-specified polygonal mask around the regions of interest in order to eliminate regions of the image that do not contain particles or contain excessive laser reflections. As the images acquired with this technique are taken at different z coordinates, the model position as seen by the cameras will be slightly different on every frame, which means manually masking each image becomes exceedingly burdensome. To alleviate this burden, a guided user interface (GUI) app was developed in MATLAB to generate a dynamically moving mask for every frame in the image, interpolating between way-points manually selected by the user at selected keyframes. The dynamic masking GUI, displayed in Fig. 2, is available on Github (Zigunov 2022), which also features an executable file for machines that do not have MATLAB installed.
As the traversing mechanism motion is linear throughout the scanning process, only a handful of keyframes are required to fully define a set of masks that will successfully delineate the PIV region of interest. The user loads the Fig. 2 GUI application utilized to generate the dynamic masks that exclude the moving model from the region of interest prior to PIV processing multi-frame files recorded from the cameras into the app, and an initially rectangular 4-vertex polygon mask is applied by default to the first and last frame. The user can then drag, add and delete vertices to the polygonal mask interactively, which is done by clicking on the polygonal mask vertices shown in Fig. 2. A linear interpolation is performed between vertices of two adjacent keyframes. Multiple keyframes can be added by simply selecting a frame and interactively dragging the vertices of the polygon, enabling the generation of complex mask motions while allowing the user to ensure intermediate masks are also adequate. Different polygons can be defined in case of multi-frame images, where different camera views are available, by changing the "camera/ frame selector" slider in Fig. 2. The masks are then saved as binary black-and-white images that can be utilized in the PIV pre-processing step to enable vectors only in the region of interest specified.

Volumetric velocity field computation through spatiotemporal averaging
After the acquisition of the data and processing of the masked camera images into SPIV vector fields, a stack of instantaneous vector fields with a total number of planes N planes = rT is available at a variety of z coordinates between z start and z end . The z coordinates are obtained from the measurements of the linear potentiometer connected to the traverse mechanism, which can readily be found for each instantaneous SPIV field by finding the rising edge of the laser Q-switch signal as shown in Fig. 1b through thresholding and differentiation. To obtain the mean velocity field ⃗ U(z) at a given plane, a spatial averaging operation is performed using the instantaneous vector fields ⃗ u(z � ): where z 1 = z − Δz∕2 , z 2 = z + Δz∕2 and n is the number of planes where the condition |z − z � | ≤ 0.5Δz (i.e., vectors within a z tolerance of ±0.5Δz ) is satisfied. In the ideal case with no frames dropped, n = n planes as defined in Eq. 2, but if frames are dropped then n < n planes . Only the instantaneous vector fields ⃗ u that satisfy this condition are used in the averaging process. This has the effect of spatially averaging, or "blurring" the vector field in the out-of-plane direction. The averaging period for ⃗ U(z) is t avg = n planes ∕r , and it must also be significantly larger than the slowest time scale of interest in the flow to obtain a converged result.
Note previous researchers (Ostermann et al. 2016;Ostermann et al. 2017;Zigunov et al. 2022a;Nakagawa et al. 2016;Cardano et al. 2008;Sellappan et al. 2018) obtained volumetric (3D-3C) velocity fields in past works by computing ⃗ U(z) for all snapshots taken at a fixed z location. This is in contrast with the proposed technique, where ⃗ U(z) is obtained for a fraction of the snapshots taken in the experiment, which is contained between two streamwise coordinates (|z − z � | < 0.5Δz) . The relaxation of the constraint that z must be constant for the averaging process to take place enables the trade-off between spatial resolution and statistical convergence of ⃗ U(z) for a given experiment, giving this averaging process a spatiotemporal aspect. More importantly, all instantaneous snapshots ⃗ u(z � ) can be captured within a single experiment, which presents a great advantage from the practical standpoint. More specifically, the requirement to label and organize the data sets for each z plane can be performed automatically through the simultaneously acquired linear potentiometer and Q-switch signals, and no pauses between z plane data acquisitions is required to save/ organize/label each data set corresponding to each z plane. Evidently, the time necessary to traverse a distance of Δz still must be significantly longer than the largest time-scale of interest in the flow, to ensure convergence of low-frequency behavior is achieved. This is, however, not a significant concern for experiments in high-speed aerodynamics.
As is typical with instantaneous PIV fields, some vectors within the field may be missing or contain spurious values due to locally poor particle density. To reduce the effect of vector outliers, the multivariate outlier detection (MVOD) method demonstrated by Griffin et al. (2010) for PIV fields is performed on a 5 ×5×n subset of the vectors. The averaging operation is then performed using only the vectors labeled as non-outliers. If faster processing is necessary, simpler outlier detection schemes may also be used.
Finally, it is worth noting that for moderately high speed PIV systems ( r > 200 Hz) it may be possible to achieve a very large number of averaging planes n even for Δz ≪ t , eliminating the spatial resolution trade-off previously mentioned. An experiment at such higher pulse rates was not performed in this work due to unavailability of materials, being left as a future work.

Mean 3D velocity fields
In this section, three example velocity field reconstructions using the technique described herein will be presented as preliminary results to demonstrate the capabilities of the proposed technique. As the focus of this manuscript is on the spatiotemporal SPIV technique, the discussion on the flow physics will be limited. Further exploration of the flow physics for the problems shown will be presented in future publications.

Compressible wake of the slanted afterbody
The first 3D-3C mean velocity field reconstructed utilizing this technique was the wake of the D = 60 mm, = 45 • slanted afterbody under moderate compressibility conditions ( M ∞ = 0.3 , Re D = 830, 000 ), which is displayed in Fig. 3 as a volume rendering utilizing the volumetric data visualization tool Paraview (Ahrens et al. 2005). The model consists of an elliptical forebody, a cylindrical body and a slanted afterbody. The slanted afterbody model with sharp edges has been previously studied by the authors under incompressible conditions (Zigunov et al. 2022a), and two distinct wake regimes were identified when the slant angle is 45 • as a function of Reynolds number. The addition of a variableradius fillet to this model seeks to represent more closely the geometrical features of a typical cargo aircraft afterbody, where the separation point is not defined by a sharp edge but by an adverse pressure gradient on a curved surface. A comparison is made in Fig. 3 between the wake of the sharp-edged slanted cylinder model (a) and the rounded edge model (b). The rounded edge model includes a variableradius fillet, which is defined with a radius of 0.31 diameters at the top edge and 0.24 diameters at the side edges of the model, while maintaining the trailing edge sharp.
Through highlighting the vorticity magnitude contours as a volume rendering, it is possible to fully observe the topology of the strong vortex pair formed in this wake. The vortex pair is also evidenced by the path taken by the streamlines seeded at key points through the flow field. It can be noted that the development of the vortex core is significantly affected by the presence of the rounded edge in Fig. 3b, resulting in smaller vortices that are significantly closer to one another when compared to the sharp-edged case. It is also possible to note that both sharp and rounded cases have a connection between the two vortex cores through the presence of an attached separation bubble at the centerline separation point.
In order to build this volumetric vector field, the traversing system described in Sect. 2.1 was installed and the processing technique proposed in this manuscript was applied. As the PSWT facility is a blow-down wind tunnel, only 60 s of run time could be afforded for an experimental test at M ∞ = 0.3 . The traversing mechanism started scanning at z start = 90 mm, at the furthest downstream location where the vortex pair is most developed. A traversal towards z end = 0 mm was then performed. A total of 695 planes were captured between z start and z end , averaging an effective laser firing rate of r = 11.6 Hz. This effective rate was lower than the intended firing rate (15 Hz) because frames were dropped as the camera link speed was insufficient to keep up with the data rate. This issue, however, does not affect the reconstruction as the Q-switch rising edges are measured by the DAQ system. Using Δz = 3 mm, an average of n planes = 23 spaced by Δz = 0.13 mm was employed for the reconstruction of the vector fields shown in Fig. 3. Although the number of planes used for averaging may not be sufficient for a fully converged 3D velocity field, it is sufficient to elucidate a very significant portion of the details about the flow topology studied. The averaging period for this case was t avg = 1.53 s, which is much larger than the time scales of the vortex wandering in this flow ( St D ∼ 0.6 ; f ∼ 1 kHz, details in Zigunov et al. (2022b)).
It is important to note that these 3D flow field reconstructions under challenging compressible conditions on a short run-time facility would not have been possible with the current tomographic PIV equipment available in our laboratory. Not only the volume measured is large (80×150×90 mm), but the seeding density control in the facility is difficult, greatly increasing the challenges related to utilizing a complex technique such as tomographic PIV, where seeding . Volume rendering shows normalized vorticity magnitude density variations can prevent the particle volume reconstruction process. By using the technique proposed herein, an invaluable 3D-3C mapping of the volumetric flow field was attainable, significantly improving the understanding of the flow physics of the problem studied that would otherwise not be available.

Generic ship wake
In this experiment, the wake of a generic destroyer ship (Government of Canada 2021) is examined with the spatiotemporal averaging technique herein described. More specifically, the wake of the ship superstructure above the helicopter landing spot under a headwind and no pitch/roll/ yaw are examined for a free stream velocity of W ∞ = 20 m/s, corresponding to a Reynolds number based on the total ship length of Re L = 570, 000.
As the superstructure wake constitutes a bluff-body wake with high unsteadiness (Tinney and Ukeiley 2009;Palm 2022), more SPIV planes were taken over the course of the experiment to attain a better convergence of the spatiotemporal average flow field. A total of 5,000 SPIV planes were taken at r = 7.5 Hz, distributed over a 150 mm streamwise distance, starting at z start = 15.7 mm and z end = 165.7 mm, where z = 0 is defined as the bottom corner of the landing deck. The landing deck is 85.7 mm long in the model examined, meaning approximately half of the SPIV measurements are within the landing deck and the other half is over the simulated sea surface, which in this experiment is a simple flat plate. With this setup, Δz = 0.03 mm, meaning for Δz = 3 mm a number of n planes = 100 was possible for the computation of the spatiotemporal averages. Thus, the averaging period for this case was t avg = 13.3 s, whereas the expected shedding frequency based on the ship height is 0.16 < St < 0.25 , or 46 < f < 72 Hz (Shukla et al. 2019).
The volumetric reconstruction of the ship wake is presented in Fig. 4a for a headwind of W ∞ = 20 m/s. A volumetric rendering of the regions with low velocity |⃗ u| < 5 m/s indicates the presence of separation bubbles behind the ship superstructure and between the ship deck and the sea surface. A representation of a generic helicopter model hovering at the dimensions provided in the Government of Canada (2021) report is also included, indicating the helicopter lands in the recirculation zone behind the superstructure for the conditions described herein, which is more evident in the center plane cross-section of Fig. 4b. Figure 4b shows contours of velocity magnitude overlaid by line integral convolutions (LIC) to highlight the patterns formed by the streamlines. Note the laser sheet used for SPIV is placed in the spanwise-vertical x − y direction, meaning the cross-section displayed in Fig. 4b utilizes all planes in the reconstruction. Note the helicopter model is added to the 3D plot after the reconstruction of the ship wake vector field is completed to display the position where the helicopter lands on the deck. There is no physical helicopter model in the experiments. Hence, the downwash created by the helicopter rotor It also can be noted in Fig. 4a that, despite our best alignment efforts, a small yaw angle is still present in the model, which can be noted by the asymmetry in the shape of the lower separation bubble between the ship and the sea floor. The yaw is lower than 0.5 • , as measured against the wall of the wind tunnel. Although an undesired feature of the results observed, the capturing of this small misalignment demonstrates the ability of the spatiotemporal reconstruction technique to accurately reproduce the volumetric velocity field studied.

Double-fin SBLI
Another difficult flow problem tackled using the spatiotemporal averaging technique herein described is the supersonic 3D shock-boundary layer interaction (SBLI) generated about two sharp fins. This configuration represents supersonic and hypersonic engine inlets with side-wall compression (Adler and Gaitonde 2019), presenting significant flow threedimensionality and interactions between crossing shock waves that are very difficult to understand with traditional diagnostics. Three dimensional velocity measurements for a single-fin SBLI were performed recently by Mears et al. (2021) using tomographic PIV. Due to the limitations of illumination intensity mentioned in Sect. 1, Mears et al. (2021) needed to perform stitching of two limited volumes with 10 mm thickness each, within a problem where flow features have sizes of the order of 30-60 mm, thus providing only a partial view of the flow field.
The double-fin configuration examined herein is described in detail by Seckin et al. (2022). The configuration where the wedge angle is 10 • is further examined now using the spatiotemporal averaging reconstruction technique. For reference, the fin height is h fin = 35.6 mm, the distance between fins is L = 63.5 mm (represented in Fig. 5a), and the total fin length in the streamwise direction is 102.6 mm. The free stream Mach number is M ∞ = 2 and the measured free stream velocity is W ∞ = 530 m/s. The incoming boundary layer thickness is = 3.5 mm.
In this experiment, the laser sheet is expanded to capture spanwise-vertical ( x − y ) planes, being traversed from z start = −5.8 mm to z end = 94.2 mm, for a total traversing span of 100 mm. In this coordinate system, z = 0 corresponds to the fin tips, increasing in the downstream direction. This traversal occurred over the course of 60 s, which at the rate of r = 15 Hz yields a total of 900 planes with Δz = 0.11 mm. An averaging distance of Δz = 2 mm was used to obtain the average vector fields, resulting in n planes = 18 . In order to attain a larger number of planes for improved convergence, ensemble averaging was utilized. Three experiments under the same conditions were performed and all planes were combined to produce an ensemble average data set with an effective n planes = 18 ⋅ 3 = 54 , which was considered adequate to converge most features of the flow. Eighteen planes result in an averaging period of t avg = 1.2 s. The dynamics of this double-fin interaction are broadband in nature, but low-frequency unsteadiness was identified at Strouhal numbers as low as St = 0.002 , corresponding to a physical frequency of f = 285 Hz, meaning at least 340 cycles of the lowfrequency SBLI unsteadiness are included in the averages.
A reconstruction of the 3D mean flow topology of this double-fin problem at M ∞ = 2 is presented in Fig. 5a along with a summary of the key features observed. Figure 5a shows a 3D view highlighting regions in the volume containing high positive (red) and negative (blue) divergence. The negative divergence regions highlight the oblique shock waves (1) and (2) produced in this flow, which present an intricate three-dimensional structure, starting with a flat wave near the fin tips and bowing in the streamwise direction due to a vertical pressure gradient induced by the open boundary at the top of the two fins. This bowing of the shocks is also highlighted as a velocity discontinuity in Fig. 5c, which presents contours of streamwise velocity at the center plane. Also at the center plane, at the bottom where the two oblique shocks merge, a Mach stem can be observed (indicated as (4)), showing the foot of the lambda shock that merges as a triangular surface at the corner of the two main oblique shocks produced by the fins. At the portion where the fins become straight again, two expansion fans form (indicated as (3) as positive, red divergence regions), which also have a complex 3D topology highlighted in Fig. 5a. At certain heights in the centerline, a slip-line shear layer (indicated as (5) in Fig. 5b) is observed as the larger strength of the oblique shock at the center plane decelerates the flow significantly more than the two oblique shocks off-center. Finally, it is also possible to observe the weak Mach waves coming from small seams between the test section and the nozzle block located upstream of the test section, as indicated in Fig. 5b, (6). The detection of these weak features is a hallmark of the ability of the spatiotemporal S-SPIV technique to resolve small velocity gradients. Such complexity in this flow topology is observed in this level of detail for the first time in experimental work, and are only possible with the equipment available in our laboratory through the usage of the spatiotemporal averaging technique herein described. As previously discussed, a similar flow was analyzed in the same facility by Mears et al. (2021) with tomographic PIV using the same equipment, obtaining only partial observations of the flow volume.

Comparison with planar PIV
To assess the ability of the S-SPIV technique to provide reliable flow field measurements, a direct comparison with the well-established planar PIV technique is herein presented. Two data sets are taken at the same conditions ( M ∞ = 0.3 , Re D = 830, 000 ) for the rounded edge slanted cylinder model, which was studied in the PSWT facility. As discussed in Sect. 3.1.1, the PSWT facility is a blow-down wind tunnel, meaning the run time available was limited to 60 s for the conditions examined, resulting in n planes = 23 for Δz = 3 mm. For Planar PIV, 1000 planes were captured to obtain a converged average.
A center-plane slice of the streamwise velocity field W at x = 0 is presented for the S-SPIV reconstruction experiment in Fig. 6b and compared against the planar PIV average field at matching conditions, shown in Fig. 6a. Within the window where both experiments have data, it is evident that all mean flow features were correctly captured by the S-SPIV shown in Fig. 6b, despite the relatively small number of planes used for averaging. More specifically, the low-velocity region at the separation bubble, the centerline velocity reduction as the free stream turns as well as the velocity defect of the trailing edge shear layer are all captured by both planar PIV and S-SPIV datasets. For further understanding of the underlying data used to obtain the mean field of Fig. 6b, the original instantaneous snapshots captured in the experiment are also displayed in Fig. 6c. Note that the previously described features are also distinguishable within the instantaneous snapshots.
A limitation of the setup utilized for the S-SPIV experiment is evident as a significant cropping above the surface were removed through the dynamic masking presented in Sect. 2.2, however they prevented the correlation of vectors near the surface, preventing the observation of the separation bubble recirculation zone that is easily observed in the planar PIV experiment of Fig. 6a. A potential improvement to the experimental setup would be to traverse the acquisition system at a 45 • angle, allowing for the positioning of the cameras from a more favorable viewpoint; although that may not be viable in all experimental setups. A more quantitative comparison between the planar PIV and S-SPIV fields is presented as a line plot in Fig. 7. The streamwise velocity is plotted as a function of z for the points along the dashed magenta line shown in Fig. 6 (a,b). Comparing the two lines, a good agreement is observed, however a difference between the two curves can be noted. As expected, the S-SPIV curve is more noisy, owing to the partially-converged averages from 23 planes. The underlying instantaneous snapshot data is also shown as green dots in Fig. 7, for comparison. Note the disagreement between the blue and black curves in Fig. 7 is greatest close to the minimum streamwise velocity value, where the standard deviation of the instantaneous snapshots (green dots) is the largest. This is a typical behavior of the mean convergence of a random variable, and is further explored in Sect. 3.3.
Furthermore, a part of the bias error between the two curves may also be attributed to the uncertainty in the positioning of the laser plane with respect to the model, since the two experiments require different laser/camera alignments. In the planar PIV case, the outline of the model is visible in the raw images, and therefore the vector fields can be located with respect to the model with an uncertainty of a few pixels (i.e., O(0.1 mm)). Conversely, the SSPIV volume is located with respect to the model by measuring the location of the calibration plate with respect to the physical model at z start = 90 mm. This introduces uncertainties related to the positioning of the calibration plate with respect to the model, which are of the order of O(2 mm) for this experiment in all axes. This spatial offset of the volume may be partially responsible for the bias error observed in Fig. 7.

Spatiotemporal average convergence
To further assess the convergence of the vector quantities calculated using the spatiotemporal averaging technique proposed in this study, a comparison is made with the purely temporal averaging method typically utilized in PIV experiments for the ship model. Five hundred SPIV snapshots were taken at a constant spanwise plane at z∕h = 1.8 ( z = 40 mm) downstream of the landing deck edge, as indicated in Fig. 4b by a blue dashed line. This plane is chosen because it includes regions of low unsteadiness (close to the free stream) as well as high unsteadiness (inside the separation bubble). The selected plane also displays a larger gradient in the traverse scanning direction z, allowing the exploration of the limits of the spatiotemporal averaging method. The purely temporal average obtained through this method is compared with the average obtained through the spatiotemporal method proposed in this study at the same z location for different averaging windows Δz = n planes Δz , all centered around the plane z∕h = 1.8.
Defining the normalized error (n s ) of the average of the variable , computed using a partial sample size n s : where N = 500 is the maximum number of samples acquired for the temporal averaging experiment, which is regarded as the most accurate estimate of the true mean. n s < N is the number of samples used for convergence assessment, and ̄(m) is the average of the variable computed using m samples. The variable can be a velocity component or its derivative. The normalized error (n s ) is normalized by the local standard deviation (N) , as the standard deviation of the mean estimator is ∕ √ N . (Rao 2002; Bendat and Piersol 1971) Figure 8a shows the spatiotemporal average field for the ship experiment at z∕h = 1.8 ± 0.065 , corresponding to Δz = 3 mm and n planes = 100 . The flow field obtained through the spatiotemporal averaging process greatly approximates the (4) (n s ) =̄( n s ) −̄(N) (N) Fig. 8 Comparison of the temporal (a) and spatiotemporal (b) averages obtained for the zero-yaw ship wake test at z∕h = 1.8 . In (c), the 95 th percentile of the normalized error for the velocity components is presented. For the spatiotemporal analysis, Δz = 0.03 mm and Δz max = 6 mm. In (d), the normalized error for the velocity derivatives is presented only for the spatiotemporal case ▸ reference flow field shown in Fig. 8b, obtained through a temporal averaging with n s = N = 500 samples. The asymmetry in the wake of the model, owing to the small yaw angle of the ship model ( ∼ 0.5 • ) is also correctly captured by the spatiotemporal averaging method. In Fig. 8c, the normalized error (n s ) computed using Eq. 4 is presented for the three velocity components in the field. As the field data includes a large number of vectors, the 95 th percentile of the normalized errors is presented and compared with the behavior expected from a Gaussian distribution at two standard deviations: First analyzing the behavior of the purely temporal averages, shown as solid lines for u, v and w, it is possible to note the normalized error trends as a function of ∼ 1∕ √ n s , as expected, always staying below the dashed black line of two standard deviations. The spatiotemporal averages, shown as dotted lines with dot symbols in Fig. 8c, present a similar behavior for Δz < 1 mm. Further increases in the number of samples continue reducing the magnitude of the error, however deviating increasingly further from the ∼ 1∕ √ n s trend. Thus, the spatiotemporal averaging technique reaches a point of diminishing returns at about Δz ≈ 3 mm for the parameters chosen for this experiment.
Similarly, the convergence behavior of the spatiotemporal averages for the in-plane average derivatives is presented Fig. 8d. The behavior of the convergence resembles the convergence of the average velocity components of Fig. 8c, however reaching a plateau about Δz ∼ 1 mm. Furthermore, the error in the derivative quantities exceeds the 1.96∕ √ n s line for the derivatives of the vertical velocity dv/dx and dv/dy, as well as for the out-of-plane velocity dw/dx and dw/dy, indicating a bias error is introduced due to the spatial averaging in the z direction.
The appearance of a bias error due to spatial averaging is easily explained for the continuous case by performing the averaging of the function (z) over a window spanning {z 0 − Δz∕2;z 0 + Δz∕2}: of the convergence behavior can be explored with future developments by sampling at a higher frame rate, increasing the experimental time, changing the size of the scanned volume, or performing ensemble averaging over multiple experimental runs.