OpenPIV-Matlab -- An open-source software for particle image velocimetry; test case: birds' aerodynamics

We present an open-source MATLAB package, entitled OpenPIV-Matlab, for analyzing particle image velocimetry (PIV) data. We extend the PIV analysis with additional tools for post-processing the PIV results including the estimation of aero/hydrodynamic forces from the PIV data of a wake behind an immersed (bluff or streamlined) body. The paper presents a detailed description of the packages, covering the three main parts: generating two-dimensional two component velocity fields from pairs of images (OpenPIV-Matlab), spatial and temporal flow analysis based on the velocity fields (Spatial and Temporal Analysis Toolbox), and wake flow analysis along with the force estimates (getWAKE Toolbox). A complete analysis with a variety of post-processing capabilities is demonstrated using time-resolved PIV wake data of a freely flying European starling (\emph{Sturnus vulgaris}) in a wind tunnel.


Motivation and significance
Particle Image Velocimetry (PIV) is a state-of-the-art optical flow measurement technique. Using PIV, two-and three-dimensional (2D/3D) velocity fields are obtained with high spatial and temporal resolution, in a non-intrusive manner. PIV is widely used in research and industrial applications. For a detailed description of PIV principles and methodology, please refer to the book by Raffel et al. [1].
OpenPIV is an international scientific community that develops free and open-source software that performs PIV analysis to obtain the velocity fields from images and a variety of post-processing tools to elucidate the physics of the investigated flow. OpenPIV houses several software packages, schematically shown in Fig. 1. Originally, OpenPIV was developed to analyse PIV images to yield a 2D velocity field using Matlab platform [2]. Later on, the spatial and temporal analysis toolbox was developed for the characterization of the measured flow field as well as provide turbulence characteristics and some features of spectral analysis based on PIV data. Two additional toolboxes were later added: Pressure and Proper Orthogonal Decomposition (POD) [2,3]. The pressure toolbox uses Poisson formulation of the flow equation to estimate the pressure field, assuming the boundary conditions of the field is known; using an iterative technique [2]. The POD toolbox was developed to extract energetic modes from velocity or vorticity fields to study coherent motions in turbulent flows [3,4]. Another OpenPIV toolbox extended the post-processing capabilities with calculation of the density fields using the background-oriented schlieren technique [5]. Most recently, we added a toolbox to estimate aero/hydrodynamic loads based on PIV data acquired in a body wake, entitled getWAKE [6].
In parallel to Matlab toolboxes, OpenPIV team developed also C++ and Python branches (see Fig. 1). Detailed description of the C++ version can be found in Taylor et al. [7]. The Python branch has additional features, including adaptive meshing, correlation of 3D-PIV images, real-time processing, among others. A paper of the OpenPIV-Python version is out of scope of this work and will be released separately. The set of tools described hereinafter, provides an open-source "full-stack" package for analyzing and post-processing PIV data measured in the body wake. The significance is in the variety of tools, together providing a complete solution supporting new PIV users and experts that study aero/hydrodynamics based on wakes of dragged or self-propelled bodies.
We demonstrate OpenPIV-Matlab capabilities using the PIV near wake flow data measured behind a freely flying European starling (Sturnus vulgaris) in a wind tunnel [8] and present how we can shed light on the role of unsteady aerodynamics in natural flyers locomotion. The wake signature behind birds for estimating their aerodynamic performance is a very active field of research due to rapid technological advancement of high-speed lasers, cameras and data transfer solutions [9,10,11,12].
The OpenPIV-Matlab stands-out due to the advanced analysis toolboxes with its capabilities beyond the other open-source PIV software, such as PIVlab [13], MATPIV, among others.

Software Architecture
OpenPIV-Matlab uses the MATLAB (Mathworks Inc.) language to provide fast advanced PIV processing and post-processing tools and straightforward development process.
• opevpivgui.m is a GUI comprised of several subroutines that allow to import PIV images, pre-process them, analyse using fast Fourier transformbased cross-correlation algorithm, filter and interpolate the flow field, and export the velocity vector maps as ASCII files.
• spatialbox.m is a GUI that allows the user to load a series of velocity maps created by openpivgui.m and calculate various flow characteristics (mean and turbulent, velocity derivatives, energy terms, auto-correlation, etc.), plot them as contours, vector fields and spatial distribution profiles and storing the processed data in Matlab MAT format.
• wake.m is a GUI which loads the MAT file exported by spatialbox.m. For the case of time-resolved wake data, it enables reconstruction of a complete wake signature behind the body using a cross-correlation algorithm that overlaps consecutive velocity maps. Based on the reconstructed wake data, wake.m also estimates aerodynamic body forces such as profile drag and cumulative circulatory lift.
In the following we describe the functionality of each toolbox.

Software Functionalities
An example of a flowchart of analysis of the body wake PIV data is depicted in Fig. 2.

OpenPIV
The PIV raw images acquired at the wake downstream of a body are imported into Matlab GUI using File->Load. The user defines the main PIV parameters, including: magnification or scale (i.e.: pixels/meter), time interval between laser pulses, (∆t), image pre-processing function (e.g. contrast enhancement or inversion of shadowgraphy images), ROI (region-of-interest) for the PIV analysis, interrogation window size and spacing/overlap size (pixels), and filters: signal-to-noise ratio (S/N) type and threshold and an outlier velocity threshold.
The OpenPIV utilizes a cross-correlation algorithm that yields a displacement vector map by correlating two consecutive PIV images. The output data in pixels of displacement is converted to physical units (meters/second) using the ∆t and magnification (scale). The main openpivgui.m loop consists of the following steps: • Cross-correlation algorithm is applied to sub-image square or rectangular interrogation windows. OpenPIV utilizes an FFT-based cross-correlation algorithm to process pairs of PIV images to yield the velocity field maps [1], therefore sizes are typically of 2 n ×2 n size (32×32, 16×64, etc.). The spacing/overlap value controls the spatial resolution of the grid x, y at which we estimate horizontal and vertical velocity components (u, v). Larger interrogation window size reduces resolution but less affected by the background noise.
• The ratio of the maximum cross-correlation peak to the average crosscorrelation value, or the ratio of the peak to the second highest peak are used as a measure of the signal-to-noise ratio (S/N) [14,15,16]. Selection of S/N type is based on the data quality (i.e., strong contrast images). The choice of S/N threshold value that marks erroneous vector is obtained through a trial and error procedure, based on manual assessment of the vector field.
• After the cross-correlation, additional filters are applied for validation and removal of outliers marked by the S/N ratio or removed based on statistics of the flow field. A so-called global filter removes vectors with length that are larger than the mean of the flow field plus N times its standard deviation. The outlier filter parameter in the GUI indicates the value of N .
• Local filter is performed on small neighborhoods of vectors using the 3 × 3 kernel, removing vectors that are more than 3 times local standard deviation distant from the local mean of the 8 nearest neighbor vectors.
• It is desired to complete the analysis with less than 5-10% of erroneous vectors. After removing the outliers, the missing values are filled using iterative interpolation, based on the valid neighborhood vectors.
The velocity field result is stored in three ASCII files: the raw vector results (dataName_noflt.txt), filtered results (dataName_flt.txt) and interpolated data (dataName.vec). Output files have headers that define the list of variables, the units (e.g. pixels/dt or m/s) and the size of the field in terms of rows and columns, followed by the 5 columns of data: x, y, u, v, and S/N.

Spatial and Temporal Analysis Toolbox
The Spatial and Temporal Analysis Toolbox GUI (spatialbox.m) loads series of dataName.vec files into a 3D flow velocity array where the 3rd dimension is the number of the flow field in the ensemble (or time for the time-resolved case). The data is automatically decomposed (using the Reynolds decomposition) into mean and turbulent fluctuations and provides both qualitative and quantitative visualization tools of a large variety of flow properties. Qualitatively, the toolbox offers several options to show the calculated flow properties in the form of colored contour maps, colored contour lines and vector representation. Quantitatively, it offers the user to plot the properties in a 2D profile format, by selecting regions of interests, or cross-sections. The flow properties include the spatial velocity derivatives and the properties such as vorticity, rate of strain, along with the turbulent parameters such as turbulent intensity, Reynolds stress, turbulent kinetic energy, production, dissipation and enstrophy. Furthermore, the toolbox includes an additional feature allowing to perform some basic spectral analysis using auto-correlation functions applied to the velocity field. The flow characteristics computed in this toolbox for the velocity maps are then exported as a binary Matlab MAT file.
For more details, the reader is referred to the following illustrative example and documentation listed in Table 1.

getWAKE Toolbox
The getWAKE toolbox is designed for wake data analysis, and specifically for the case where several consecutive flow velocity maps contain time evolution of the flow in the wake or motion of the same vortices as they shed behind the body. The toolbox allows to define multiple experimental setup related parameters, required for the reconstruction (i.e. combination) of the wake signature behind a body, followed by the forces estimate procedure. The parameters include the PIV parameters, body parameters (dimensions and weight) and an incoming flow parameters (free stream velocity, fluid density and dynamic viscosity) in SI units and the motion type (stationary or cyclic flapping motion, etc.). For the flapping case, for instance, the parameters relate to three phases of downstroke, transition and upstroke. In some cases, a sub-region from the velocity map is selected mainly to remove noise at the edges. For the wake reconstruction, velocity maps sequence parameters and the cross-correlation flow parameter (different then the one used in the OpenPIV analysis, applied to the instantaneous or fluctuating velocity fields) are defined. Further details regarding the wake reconstruction scheme are available in Appendix A.
The reconstructed wake signature is presented using the velocity fluctuations vector fields (u , v ) and colored with the normalized spanwise vorticity field ω z c/U ∞ , where ω z = ∂v/∂x − ∂u/∂y, c is the characteristic length scale of the body and U ∞ is the freestream velocity. Multiple visualization options are available, including contour threshold, vorticity or swirl strength [17] contours, with or without Gaussian smoothing, etc. The reconstructed wake can be shown and exported as a static image or as an animated movie.
The wake data enables the estimation of the profile drag and cumulative circulatory lift coefficients. The lift is calculated using either i) Panda and Zaman method [18] or ii) a direct summation of the circulation values. Further details regarding the forces estimation are given in Appendix B.
Forces can be presented as a function of the normalized streamwise wake distance (x/c, where c is the body size, or chord length of the wing) or time (t). If the wake is generated due to a flapping wing cyclic motion, one can correlate the wingbeat phases (downstroke, upstroke) with the velocity maps in the wake by setting the flapping wingbeats kinematics inputs. Thus, the force coefficients can be presented for a specific wingbeat or a sequence of wingbeats that correspond to the reconstructed wake signature, where the different wingbeat phases are highlighted.

Illustrative Example
We exemplify the analysis of PIV wake data measurements obtained in the wake of a freely flying European starling (Sturnus vulgaris) in a wind tunnel; see Kirchhefer et al. [8] and Ben-Gida et al. [19] for details. The experimental setup utilized for PIV measurements behind the bird is depicted in Fig. 3. Experiments were conducted in a climatic closed-loop wind tunnel at the Advanced Facility for Avian Research (AFAR), at Western University (Canada). The wind tunnel test section with cross-sectional area of 1.2m 2 is specifically designed for simulating the flight conditions experienced by birds during long distance migratory travel. The flight conditions reported in this work correspond to atmospheric static pressure, a temperature of 15 • C, and relative humidity of 80%.
A European starling (Sturnus vulgaris) had been trained to fly in flapping flight mode in the wind tunnel. At the time the experiments were performed the bird had a mass of 76g, where its wings had an average chord of c = 6cm, a maximum wingspan of b = 38cm and an aspect ratio (wingspan squared divided by the wings lifting area) of AR = 6.4. A typical cruising flight speed of U ∞ = 13.5m/sec was chosen for the experiments [8].
Flow measurements were taken using a long-duration time-resolved PIV system, developed by Taylor et al. [7], consisting of a 80W double-head diodepumped Q-switched Nd:YLF laser at a wavelength of 527nm and two CMOS cameras (Photron FASTCAM-1024PCI) with a spatial resolution of 1MP at a sampling rate of 1kHz (see Fig. 3). Olive oil aerosol particles with an average size of 1µm [20] were introduced into the wind tunnel using a Laskin nozzle from the downstream end of the test section. One camera was used for measuring the PIV wake data (in the streamwise-normal plane) whilst a second one was used for recording the bird kinematics simultaneously with the PIV, where the PIV wake measurements were taken 2 chord lengths behind the starling. The fields of view of the PIV the bird kinematics cameras were 2c × 2c and 9c × 9c, respectively. The PIV system sampled pairs of wake images at a rate of 500Hz (2msec intervals), which allows sufficient resolution for temporally resolving the wake of the bird that flap its wings at a frequency of 15Hz.
The distance between the bird wing and the location of the laser sheet was determined by placing a set of IR detectors that trigger the laser once the bird passed the laser location about 3-4 chord length upstream (see Fig. 3).

OpenPIV analysis
The velocity fields were computed from the PIV wake data behind the bird with 32 × 32pixel 2 interrogation windows and 50% overlap (16 × 16pixel 2 ), thus yielding a spatial resolution of 32 vectors per average chord of the bird (c), equal to 1.8 vectors per millimeter. Here, the second type of the S/N type was used (with a threshold value of 1) and an outlier filter value of 100 was chosen. An example of velocity vector maps computed by the OpenPIV is depicted in Fig. 4. The raw velocity map (non-filtered) is depicted in Fig. 4a, where the outlier vectors are colored in red. The filtered velocity map (after applying global and local filters) is depicted in Fig. 4b. The final velocity map is shown in Fig. 4c, where the interpolated vectors are colored in green.  The final velocity vector maps in the wake, which were computed using the OpenPIV (opevpivgui.m) and exported as .vec files, are then imported into the Spatial and Temporal Analysis Toolbox (spatialbox.m) for the computation of the required flow characteristics. Fig. 6 depicts an example of the instantaneous spanwise vorticity fields computed using the Spatial and Temporal Analysis Toolbox, from the velocity vector maps depicted in Fig. 5, in the near wake behind the starling during flight.   The flow characteristics and the velocity vector maps are then exported from the Spatial and Temporal Analysis Toolbox as a binary MAT file, and imported into the getWAKE Toolbox GUI (wake.m). The characteristics of the experimental setup (PIV system, bird and freestream flow) are set in the getWAKE GUI. Using these inputs, the wake signature behind the starling is reconstructed for a sequence of velocity maps, which corresponded to flapping wingbeat phases: downstroke and upstroke, as depicted in Fig. 7. The bird appears to fly from right to left; thus, the downstream wake essentially occurred earlier (corresponding to the downstroke phase), while the upstream wake occurred later (corresponding to the upstroke phase). Here, we utilized the velocity fluctuations (u , v ) for cross-correlating the velocity maps. The time variation of the profile drag and cumulative circulatory lift forces exerted on the bird is depicted in Fig. 8. We demonstrate the forces calculations over a single and multiple wingbeat phases. The time variation of the forces for a single wingbeat are given as a function of the non-dimensional time, t/T ; where T is the time period of the flapping wingbeat cycle. The white shaded regions in Fig. 8 correspond to the downstroke phases of each flapping wingbeat, whereas the gray shaded region corresponds to the upstroke phases. Herein, C d , total drag coefficient is the summation of the C d0 and C d1 components, which are the steady and unsteady components of the drag coefficient, respectively. The cumulative circulatory lift coefficient was computed based on Panda and Zaman method [18], see the procedure in Appendix B).

Impact
The main contribution of this work is to provide an open source code for PIV data analysis and post-analysis. OpenPIV-Matlab is a complete 2D PIV data analysis software package that can handle PIV and its extensions (i.e., time-resolved, stereo) images and perform a detailed flow analysis. To the best of our knowledge, OpenPIV-Matlab is the most complete set of tools for the flow analysis based on data acquired using optical measurement techniques in fluid dynamics, providing thousands of lines of code as an open source. The toolbox is used in a variety of fields associated with fluid dynamics such as such as mechanical, aerospace, civil and environmental engineering, biology, chemistry, earth sciences and so forth. The package presented herein offers the only known open-source software for advanced post-processing of time-resolved PIV wake data behind immersed bodies, enabling visual description of the wake evolution over time and space estimation of the aero/hydrodynamic forces.

Conclusions
We present a Matlab (Mathworks Inc.) software package for analysis and post-analysis of PIV data. OpenPIV is comprised of three toolboxes that can analyse PIV images to yield a 2D2C (two-dimensional, two-component) velocity vector fields using cross-correlation technique, post-analysis the velocity fields calculating various flow properties including turbulence and spectral analysis and estimation of the aero/hydrodynamic forces exerted over immersed bodies in fluids. As test-case study of the toolbox applicability to PIV time-resolved data, we present flow analysis of the near wake region behind a freely flying bird in a closed-loop wind tunnel. We demonstrate the functionality and usefulness of the toolboxes as applied to a data collected in a wake region that is unsteady and turbulent. The main contribution is a free, flexible, extendable and validated platform for both basic and advanced PIV related analysis and post-processing.

Conflict of Interest
We wish to confirm that there are no conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Appendix A. Wake Reconstruction
The wake reconstruction (see an example wake measured behind a flapping starling in Fig. 7) is based on time-resolved PIV images taken from a stationary camera yielding Eulerian observation of the flow field behind an immersed body (either stationary or undergoing flapping wing motion). The wake reconstruction method described herein was originally developed by the first author and later on certain aspects of the methods appeared in [21]. Thus, there are similarities in terms of its core principles, however, the utilization and output of the reconstruction scheme detailed below are different.
For the wake reconstruction, we assume the body's position did not change much relative to the measurement plane. Based on Taylor's frozen turbulence hypothesis [22], we assume that the turbulent flow remains relatively unchanged as it passes through the measurement plane. This hypothesis implies that there is no significant variation of a spatial velocity distribution over the timescale required for observation, supported by previous studies, e.g. Zaman and Hussain [23].
As a minimal requirement, sufficient temporal resolution is needed to track the flow patterns as they propagate from one velocity map to another (see Fig. A.1). For example, if the wake is measured behind a bird in flapping flight with a flapping frequency of 10Hz, the PIV velocity vector maps must be sampled at a rate of 100Hz or higher (≥ 200Hz PIV raw images sampling rate); thus, enabling a given flow structure to be tracked by at least three consecutive velocity maps.
The wake composite image is generated by offsetting each consecutive PIV velocity map with a calculated instantaneous convection velocity and then overlap the images, while keeping the mid region of each instantaneous PIV velocity map (to avoid overlapping 'noisy' data from the velocity map edges). The instantaneous convection velocity, which determines the offset if each PIV velocity map, is calculated based on a cross-correlation algorithm that examine the match of the velocity vector fields (u, v) or the fluctuating velocity vector fields (u , v ) of two consecutive velocity maps. The cross-correlation coefficient, which determines the match of two consecutive PIV velocity maps, is calculated as follows (with respect to a given flow property s): If C = 1, the two velocity maps are identical (perfectly matched), whereas if the two images are different, then C < 1. (I, J) is the size of the overlapping area and σ s is the standard deviation of the flow property s. Here, s is to be replaced with (u, v) for computing the cross-correlation coefficients with respect to the velocity vector components (C u , C v ), or with (u , v ) for computing the cross-correlation coefficients with respect to the velocity fluctuations vector components (C u , C v ). The spatial shift (X, Y ) of any instantaneous PIV velocity map is computed based on overlapping that achieves the maximum correlation If the cross-correlation failed during the wake reconstruction process (e.g., when the PIV data is low in quality or if the computed spatial shift is larger than the velocity map boundaries), the spatial shift of each instantaneous PIV map is determined from the freestream advection velocity U ∞ and the time difference ∆t; i.e., X = U ∞ ∆t.

Appendix B. Forces Estimation
The variation of the drag force coefficient C d can be estimated from the PIV wake data as shown in Ben-Gida et al. [19]. The analysis is based on the flow momentum equations that with appropriate assumptions leads to the following formulation of the drag coefficient for the two-dimensional wake case: where (x, y) is the Cartesian coordinates system used for the streamwise-normal plane in the wake; the x-axis indicates the downstream direction and the y-axis is the normal direction. Conceptually, the steady drag coefficient component C d0 is proportional to the velocity deficit at the wake, whereas the unsteady drag coefficient component C d1 is associated with the unsteady flow motion. While the steady drag term can be obtained from the near wake velocity field, the unsteady drag term requires information regarding the entire control volume (surface in the case of 2D PIV where we assume volume per unit length) surrounding the body over time. We assume most of the unsteady disturbances generated by the unsteady motion are obtained from the velocity field at the near wake where both unsteady contribution and viscous effects have not dissipated yet. Therefore, we approximate the full surface integral of the unsteady term to include only the velocity field obtained from the PIV experiments in the body near wake. Here, U ∞ is the mean undisturbed streamwise velocity, and h and l are the vertical and horizontal extent of the computed velocity field in the wake, respectively. For the steady drag coefficient, we average the various profiles along the streamwise extent of each PIV map, thus yielding a single C d0 value to represent each PIV map. It is noteworthy that drag coefficient is computed for a sequence of velocity maps, and therefore it does not require a priori information of the reconstructed wake. The variation of the cumulative circulatory lift coefficient ∆C lcirc can be estimated from the PIV velocity vector fields based on Stalnov et al. [24], Ben-Gida et al. [25] and Nafi et al. [26]. Any fluid motion around a body is accompanied by the shedding of vortices into the wake. By analyzing these vortical patterns in the near wake, one can estimate the unsteady lift exerted on the body. Herein, the estimation of increment in the time-dependent lift throughout the wake is evaluated from the PIV velocity fields by utilizing Wu's viscous flow approach [27], which was later expressed by Panda and Zaman [18]. This method can only be applied for near wake flow measurements behind bodies, where it is assumed that the wake has not deformed yet and interactions between the vortices (which results from the tip and root regions) shed into the wake are not significant. Assuming two-dimensional, incompressible flow and neglecting added mass effects, the time-dependent circulatory lift force exerted on a body can be evaluated from the near wake flow field, as follows [18]: x-moment of the vorticity field In the above equation, the first integral from left is the first x-moment of the vorticity field (A), with ω z (t) as the instantaneous spanwise vorticity field. The second integral from left is the contribution from the viscous term (diffusion), where ν is the kinematic viscosity and ρ is the fluid density. Applying Taylor's hypothesis (as introduced above for the wake reconstruction), dx = U ∞ dt, one can transform the spatial derivative in the left-side integral of Eq. (B.2) into a temporal one. Moreover, by interchanging the leftside integral in Eq. (B.2) with the time derivative (using Leibniz integral rule), one can re-write Eq. (B.2) as follows: Therefore, the change in the lift in time δτ can be expressed accordingly: where the corresponding change in the circulation δΓ is given by: Since at the beginning of the unsteady motion the lift is unknown, we shall refer to the estimated lift component as an increment in the circulatory lift [24,25] that is generated from the beginning of the motion. Based on Eq. (B.3), the cumulative circulatory lift at time t, ∆L circ (t), is computed accordingly: with the vorticity flux term ζ(t) being expressed as follows: Here, u c is the advection velocity at which the characteristics of the wake collectively travel downstream. The cumulative circulatory lift coefficient at time t is therefore expressed as: ∆C lcirc (t) = 2 cU ∞ˆt 0 ζ(t)dt = 2Γ(t) cU ∞ (B.8) As detailed in Section 2.2.3, the getWAKE Toolbox suggests two options for computing the cumulative circulatory lift coefficient in the wake. The first option is based on Panda and Zaman method [18], where the vorticity flux ζ(t) defined by Eq. (B.7) is computed for each individual PIV velocity map obtained in the wake behind the body (as function of time), after applying a threshold on the vorticity contours. For each velocity map (or time t), the advection velocity is defined as the mean streamwise velocity component along the x-direction, u c ≈ u x (y), and the spanwise vorticity field is estimated as the local mean spanwise vorticity along the x-direction, ω z (x, y, t) ≈ ω z x (y, t). Moreover, for each velocity map, the second order derivatives of u are computed using a least squares differentiation scheme and then approximated with their local mean value along the x-direction; ∂ 2 u(x, y, t)/∂x 2 ≈ ∂ 2 u/∂x 2 x (y, t) and ∂ 2 u(x, y, t)/∂y 2 ≈ ∂ 2 u/∂y 2 x (y, t). All the above vectors, representing each instantaneous velocity map, are then integrated over the y-direction to yield the vorticity flux ζ(t), as presented in Eq. (B.7), and the time variation of the cumulative circulatory lift coefficient in the wake (see Eq. (B.8)). Using the option described above, the computation of the cumulative circulatory lift coefficient for a sequence of velocity maps does not require a priori information of the reconstructed wake. Nevertheless, if needed, in the getWAKE Toolbox, the user can choose to compute the vorticity flux ζ(t) directly from the reconstructed wake signature and not from individual velocity maps; with or without a threshold applied on the vorticity contours. In doing so, the advection velocity is approximated as the freestream velocity, u c ≈ U ∞ and the spanwise vorticity field array is taken directly from the reconstructed wake, ω z (x, y), where the x-direction and time t are interchangeable by applying dt = U ∞ dx. The second order derivatives of u are computed using a least squares differentiation scheme directly from the reconstructed wake; ∂ 2 u(x, y)/∂x 2 and ∂ 2 u(x, y)/∂y 2 . All the above arrays, representing the full reconstructed wake map, are then integrated over the y-direction (at each x location) to yield the vorticity flux ζ(t), as presented in Eq. (B.7).
In the second option, the cumulative circulatory lift coefficient in Eq. (B.8) is computed based on a direct summation of the circulation values throughout the reconstructed wake signature. Thus, ζ(t) is computed as the integral over the spanwise vorticity field of each individual PIV velocity map, accordingly: