Time-Supersampling 3D-PIV Measurements by Vortex-in-Cell Simulation

Particle image velocimetry (PIV) is currently the dominant technique for measurement of the instantaneous velocity field in unsteady and turbulent flows relevant for aircraft aerodynamics. The study of such flows is essential for realisation of improvements to propulsion systems and the external configuration of aircraft and to allow for reduction of their fuel consumption, emissions and noise levels. By the advent of solid-state diode-pumped lasers and fast CMOS imagers, time-resolved PIV presently allows for study of unsteady and turbulent flows at low to moderate flow velocities (V ∼ 10 m/s). However, at high measurement rates, laser pulse energy available is strongly limited and consequently, the feasible repetition rates are insufficient for time-resolved measurement of unsteady flows at speeds relevant for aircraft aerodynamics (V ∼ 100 m/s). In order to relax the measurement rate requirements currently limiting the study of unsteady and turbulent flows, the present thesis work investigates a novel approach combining techniques from computational fluid dynamics with PIV velocity measurements. A numerical approximation to the Navier-Stokes equations is proposed to increase the temporal resolution of correlated PIV measurements, done at a measurement rate below what is historically dictated by the Nyquist criterion for reconstruction of the velocity fluctuations. The principle of the time-supersampling method is that the spatial information available by the measurements can be leveraged to increase the temporal-resolution. Feasible measurement rates are especially critical for tomographic PIV measurements and therefore the present work focusses on application to the three-dimensional datasets obtained from tomographic PIV experiments. The solution of the governing equations is based on the Vortex-in-Cell (VIC) method (Christiansen, 1973) under the hypothesis of incompressible flow. The unsteady numerical simulation of the dynamic evolution of the measured flow is applied within the 3D measurement domain and time-integration is performed between pairs of consecutive measurements. Initial conditions are taken from the first measurement field and time-resolved boundary conditions are approximated between the two fields by an advection based model. Temporal continuity of the velocity field is obtained by imposing a weighted-average of forward and backward time-integrations. The accuracy of the proposed time-supersampling method is studied with two experimental datasets obtained from time-resolved tomographic PIV measurements: a turbulent wake, and a circular jet. The results are compared to linear interpolation, advection-based supersampling, and measurement data at high sampling rate. The proposed method improves upon the M.Sc. Thesis Jan F.G. Schneiders

The principle of the time-supersampling method is that the spatial information available by the measurements can be leveraged to increase the temporal-resolution. Feasible measurement rates are especially critical for tomographic PIV measurements and therefore the present work focusses on application to the three-dimensional datasets obtained from tomographic PIV experiments. The solution of the governing equations is based on the Vortex-in-Cell (VIC) method (Christiansen, 1973) under the hypothesis of incompressible flow. The unsteady numerical simulation of the dynamic evolution of the measured flow is applied within the 3D measurement domain and time-integration is performed between pairs of consecutive measurements. Initial conditions are taken from the first measurement field and time-resolved boundary conditions are approximated between the two fields by an advection based model. Temporal continuity of the velocity field is obtained by imposing a weighted-average of forward and backward time-integrations.
The accuracy of the proposed time-supersampling method is studied with two experimental datasets obtained from time-resolved tomographic PIV measurements: a turbulent wake, and a circular jet. The results are compared to linear interpolation, advection-based supersampling, and measurement data at high sampling rate. The proposed method improves upon the Robustness and sensitivity to outliers and precision errors is studied. Blocks of outliers present in the data are smoothed by the proposed technique, but remain in the reconstructed velocity fields and for accurate time-supersampling the measurement data must first be subjected to outlier detection and removal. On the other hand, a stable reconstruction is demonstrated when the measurement data is corrupted by both spatially uncorrelated and correlated precision errors up to levels at which the experiment would normally be discarded. At high precision error levels, the evaluation of the weighted average of a forwards and backwards integration allows for a reduction in precision error. For assessment of the reconstruction quality in real-world situations and in the absence of reference measurements, an a-posteriori error estimate is proposed.
The present thesis demonstrates the sampling rate requirements for tomographic PIV can be strongly reduced when the data is time super-sampled using the proposed Vortex-in-Cell method, thereby extending the range of application of time-resolved 3D-PIV by tomographic PIV or similar techniques. On the one hand the technique allows for larger measurement domains and on the other hand it opens doors to time-resolved measurement of high speed flows (V ∼ 100 m/s) by tomographic PIV. Potential wider applications of the proposed time-supersampling technique in the fields of noise reduction and evaluation of acceleration and pressure from an instantaneous PIV measurement are illustrated using data from a tomographic PIV experiment on a jet and further development and validation of the latter is recommended as a subject for future work. of the tiny tracer particles a high-power light source is required and for air flows typically lasers are both required and used. These high power requirements limit the temporal resolution that can be achieved, as the laser pulse energy drops when increasing the measurement rate. For tomographic PIV, repetition rates up to 5 kHz have been reported in a recent topical review by Scarano (2013), which is largely insufficient to temporally resolve aerodynamic problems, where flow velocities may approach 100 m/s or more .
Intuitively one easily agrees with the statement that when measurements are performed with a sampling frequency below the frequency dictated by the Nyquist criterion, accurate reconstruction of the velocity fluctuations between measurements is only possible in some special cases. However, is this really true? Consider a tennis player throwing up a tennis ball for his service at time T 0 . A photo of this situation at time T 1 > T 0 is schematically given in Figure 1.1a. Say the height h and, analogously to a PIV measurement, the velocity v of the tennis ball at T 1 are measured. On inspection of the situation, one can directly imagine what happens for t > T 1 ; the height of the tennis ball will increase, until the ball is fully decelerated due to the presence of gravity, after which depending on the skills of the tennis player it will be served at its highest point. In this case, from only one single velocity measurement, almost the complete path of the tennis ball is determined.
The location of the tennis ball can be calculated using either a linear approximation (Figure 1.1b) or a non-linear model (Figure 1.1c), where only the latter allows for estimation of the highest point of the tennis ball. These approximations are possible because a model for the movement of the tennis ball is available through Newton's second law of motion. Consider now the case of a real tomographic PIV velocity measurement; it is an instantaneous snapshot of the flow, but analogously to the case of the tennis ball, the flow governing equations are known, which may allow for calculation of the velocity fluctuations in a short time interval around the measurement. Ideally, based on simulation of the flow governing equations and leveraging the available spatial information, the temporal velocity fluctuations in between two consecutive measurements can be reconstructed.
For air flows relevant for aircraft aerodynamics, the flow governing equations are the Navier-Stokes equations (Anderson, 2006). Based on the principle of frozen turbulence, an approxim-ation of these has been applied to PIV measurements by Scarano and Moore (2011) in order to increase the temporal resolution of the measurements, or as it was appropriately called by Scarano and Moore (2011), to time-supersample the measurements. This technique, solving the linear advection equation, is analogous to the linear approach sketched in Figure 1.1b and was demonstrated to accurately reconstruct velocity fluctuations between correlated measurements in flows where the assumption of frozen turbulence holds. However, in flows with strong curvature of streamlines (e.g. in shear layers) the reconstruction quality degrades rapidly.
To improve upon the advection model, in the present thesis a novel time-supersampling approach is proposed, which is by taking into account non-linear flow physics analogous to the sketch in Figure 1.1c. It combines both computational and experimental fluid dynamics, not to validate one or the other, but to yield flow fields with high temporal resolution, that could not be readily calculated or measured using either computational or experimental fluid dynamics alone. At the time of writing this thesis, PIV is considered the most appropriate and dominant technique for velocity measurements with high spatial resolution (Westerweel et al, 2013) and therefore the present work focusses on application to time-resolved PIV measurements. The investigation is however essentially applicable to any measurement technique that allows for measurement of velocity fields with high spatial resolution.
This thesis is divided in three main parts. The first part considers a concise literature review, with first in chapter 2 an overview current state-of-the-art PIV techniques and literature on application of CFD to PIV measurement data. From this review, vortex methods are selected as CFD technique for application to measurement data and therefore these methods are reviewed in more detail in chapter 3. In the second part of this thesis, the time-supersampling method is proposed, developed and validated. First, in chapter 4 the method is discussed in detail, after which in chapter 5 its numerical implementation is presented. The proposed technique is validated with tomographic PIV experiments in chapter 6 and additionally the robustness and sensitivity of the method to measurement noise is analysed in chapter 7. In the final part of this thesis, the scope is broadened and alternative applications of the proposed technique are discussed in chapter 8. Finally, chapter 9 summarises the conclusions and recommendations of the thesis work.
In the introduction, the requirement for a technique to reconstruct velocity fluctuations in between time-resolved PIV measurements was highlighted. In this chapter this is elaborated in more detail by means of a concise literature review. First in section 2.1 state-of-the-art PIV techniques are reviewed, after which section 2.2 summarises literature on post-processing procedures for PIV using techniques from the field of computational fluid dynamics. The results from this literature survey are synthesised in section 2.3.

State-of-the-Art Particle Image Velocimetry
Particle image velocimetry (PIV) is currently the dominant method for measurement of velocity fields (Westerweel et al, 2013;Thurow et al, 2013). As discussed in the introduction, it allows for measurement of velocity through correlation of particle image pairs. For a discussion on the working principles of the method, the reader is referred to recent review papers by Westerweel et al (2013) and Scarano (2013), or for a more in-depth discussion to the book by Raffel et al (2007). For the purpose of the present thesis, in this section the state-of-the-art techniques are summarised, with particular focus on PIV techniques allowing for time-resolved measurements.
In contrast to point-wise measurements, PIV captures the velocity field in a finite measurement domain at many points, typically 10 3 − 10 5 , with sufficient spatial resolution to permit calculation of instantaneous vorticity and rate of strain (Westerweel et al, 2013). Where historically PIV measurements could only be made in a plane, following the introduction of tomographic PIV by Elsinga et al (2006), measurements of all three velocity components in a three dimensional measurement volume have become feasible. Notably, measurement volumes of up to 50 cm 3 in air flows and up to 200 cm 3 in water tanks have been reported at limited measurement rates on the order of 1 Hz (c.f. Table 2.1). Furthermore, tomographic Cylinder wake 10 10 (3.7 × 3.6 × 0.8) 400 1 (Elsinga et al, 2006) Turbulent boundary layer 10 13 (4.2 × 3.0 × 1.0) 400 1 (Elsinga et al, 2007) Backward facing step 7 52 (7.0 × 9.3 × 0.8) 400 (double-pass) 1 (Schröder et al, 2011) Turbulent boundary layer 3 32 (6.0 × 6.0 × 0.9) 250 (double-pass) 1 (Atkinson et al, 2011) TR Cylinder wake (in water) 0.2 160 (8.0 × 10 × 2.0) 200 (double-pass) 7  TR Trailing edge 14 18 (4.7 × 4.7 × 0.8) 10 (multi-pass) 2,700 (Ghaemi and Scarano, 2011) TR Turbulent boundary layer 7 19 (3.4 × 3.0 × 1.9) 21 (multi-pass) 5,000 (Schröder et al, 2009) TR Turbulent boundary layer 10 3 (2.0 × 4.1 × 0.4) at 1 kHz: 25 10,000 (Pröbsting et al, 2013) (multi-pass) PIV experiments at sufficient spatial resolution for evaluation of spatial derivatives of velocity in such measurement domains have been reported in literature. For example, Worth and Nickels (2011) measured by tomographic PIV fine-scale structures in homogeneous isotropic turbulence and evaluated the velocity gradient tensor for calculation of the fields of enstropy and dissipation rate, and Schröder et al (2013) evaluated vorticity in the flow following a backward facing step. It should be remarked in turbulent flows at high Reynolds numbers, spatial resolution becomes more critical (Tokgoz et al, 2012). Nevertheless, tomographic PIV is to date the only established experimental technique for measurement of velocity fields in such flows in a relatively large measurement volume (Westerweel et al, 2013).
The advent of diode-pumped solid state lasers and fast CMOS imagers has allowed for timeresolved PIV measurements. Laser power available for each pulse however drops significantly when increasing repetition rate. For uncorrelated measurements at 1 Hz, pulse energies of 400 mJ are reported, which is in stark contrast to pulse energies as low as 10-20 mJ for time-resolved measurements (Table 2.1). To obtain better illumination of the measurement domain, dual-pass or even multi-pass systems are used, which essentially redirect laser light towards the measurement region through an arrangement of mirrors (Scarano, 2013). However, even using such techniques pulse energy becomes increasingly critical at higher repetition rates. Additionally, CMOS cameras can presently record 1024 × 1024-pixel images at frame rates of 5 kHz (Westerweel et al, 2013), but to allow for higher frame rates the image format needs to be reduced. It should be remarked the limitation on laser pulse energy is especially critical for time-resolved tomographic PIV measurements, where a relatively large three-dimensional volume needs to be illuminated. For this reason, the present thesis work focusses in particular on methods applicable to 3D datasets obtained from tomographic PIV.
The feasible repetition rates quoted above are mostly sufficient for water flows and low speed air flows, but largely insufficient for aerodynamic problems with flow velocities of 100 m/s or more. To overcome the measurement rate limitations, dual-PIV systems have been proposed by Jakobsen et al (1997) and Liu and Katz (2006). These systems allow for calculation of the material derivative of velocity at a single time instant without the requirement for timeresolved system. However, even in the most ideal case these systems only provide information at uncorrelated time instances, as the limitations on repetition rates are still present.
When the measurements are temporally under-resolved, i.e. at a measurement rate below what is dictated by the Nyquist criterion, reconstruction of the velocity fluctuations using pointwise interpolation techniques between consecutive measurements is not possible. Under the assumption of flow quasi-periodicity, Druault et al (2005) proposed an interpolation technique based on a proper orthogonal decomposition for time interpolation of PIV measurements. This approach is however as expected only valid for the large-scale structures of the flow (Bouhoubeiny and Druault, 2009). Similarly, synthetic temporal derivatives can be obtained by phase averaging measurements (van Oudheusden et al, 2005). This does not require a high repetition rate, but does require phase synchronisation and most importantly strong periodicity in the measured flow field.
From the discussion above follows that no state-of-the-art tomographic PIV measurement technique allows for time-resolved velocity measurements at repetition rates sufficient for high-speed flows. Recall now the tennis ball example from the introduction. Where the aforementioned techniques fail to reconstruct velocity fluctuations from temporally under resolved measurements in turbulent flows, a physics based interpolation technique is expected to improve upon them. This is investigated in more detail in the next section, by means of a literature review on current methods applying CFD to velocity measurements.

Combining Computational Fluid Dynamics and PIV
Direct application of techniques used in computational fluid dynamics (CFD) to measurements is not trivial, due to for example the sometimes strong stability requirements of CFD codes and relatively noisy measurements (detailed in chapter 7). Additionally, as will be shown in the experimental validation, the physics based interpolation technique proposed in this thesis relies heavily on the availability of three dimensional measurements of all three velocity components. The latter only recently became feasible after introduction of tomographic PIV by Elsinga et al (2006). Possibly because of this, only limited literature is available on the topic of increasing temporal resolution by simulation of the Navier-Stokes equations.
In the following sections a summary of the relevant state-of-the-art literature available is presented. First, three existing approaches for increased temporal resolution based on the combination of numerical simulations and PIV measurements are given in chronological order in sections 2.2.1 to 2.2.3. Afterwards, the scope of the investigation is broadened to discuss three techniques not aimed at increasing temporal resolution of PIV measurements, but which use CFD techniques respectively for pressure calculation (section 2.2.4), filling spatial gaps (section 2.2.5) and in the field of computer vision (section 2.2.6). The main conclusions of this literature survey are summarised in section 2.3. Hayase and Hayashi (1997) proposed to improve numerical simulations of the Navier-Stokes equations by adding an artificial body force f to the momentum equation,

Measurement-Integrated Simulation
where the body force is proportional to the difference between the numerical simulation and a velocity measurement multiplied by a gain factor K. This general framework is called Measurement-Integrated (MI) simulation and allows for improvement of numerical simulations through measurements even in a small part of the computational domain. However, as discussed in Yamagata et al (2008), no satisfactory improvement of the temporal resolution is obtained when the velocity measurements are temporally limited. A similar method was proposed by Suzuki et al (2009a,b), who present a hybrid algorithm where velocity is computed as the weighted average of a Large Eddy Simulation (LES) and a particle tracking velocimetry (PTV) measurement. Similar to the MI simulations by Hayase and Hayashi (1997), this technique requires simulation of the full flow field and uses measurements in part of the domain to improve the numerical simulation.
A different approach, which can also be considered a measurement-integrated simulation is proposed by Ma et al (2003). In their paper it is proposed to extract eigenmodes from PIV measurements and use those to improve a numerical simulations. Again, the basis of such a method is the numerical simulation and the measurements are solely used for improvement of the simulation.

Linear Advection Model
In contrary to the MI simulations, which employ measurements to improve numerical simulations, Scarano and Moore (2011) proposed a technique aimed at improvement of measurements by reconstruction of velocity fluctuations between correlated measurements. A procedure the authors appropriately named 'time-supersampling.' The technique is based on the assumption of frozen turbulence and calculates velocity at time t in between two PIV velocity measurements u m at times T 1 and T 2 by solving with ∆T = T 2 − T 1 and where the locations x 1 and x 2 are found by solving the linear advection equations 3) The convective velocity u c is estimated following an iterative procedure, with an initial estimate based on the local velocity. In contrast to the MI simulation techniques, which require simulation of the full flow field, the advection model takes as input only the velocity measurements and computations are performed on the same computational domain and grid as the measurements. Also, the linear advection model ensures continuity of the solution by taking a weighted average a forward and backward projection between two measurements, which is inexpensive compared to more advanced variational data assimilation techniques employed by the methods discussed in the next section. Consequently, and also because a relatively simple linear model is employed, the method is relatively inexpensive and readily applicable to three-dimensional datasets from tomographic PIV, which has been demonstrated in Scarano and Moore (2011) and Ghaemi and Scarano (2011).
In flows where the assumption of frozen turbulence holds to a large extent, such as for example in a turbulent boundary layer where small turbulent fluctuations are transported at a relatively high convection speed, the method is able to accurately reconstruct the temporal velocity fluctuations between correlated measurements. Significant reductions of the required measurement frequencies for accurate reconstruction have been demonstrated by Scarano and Moore (2011) for such flows. The model however deteriorates rapidly in flows with strong curvature of stream lines and acceleration, such as for example in a shear layer.
In the next section more advanced iterative adjoint-based methods are discussed. Taking into account non-linear flow dynamics, these promise to improve upon the linear advection model in non-linear flow regimes.

Iterative Adjoint-Based Navier-Stokes
Improvement upon the linear advection model is expected when the full Navier-Stokes equations are simulated on the measurement domain. Lemke and Sesterhenn (2013) consider a such an approach, where the squared difference between a numerical simulation of the Navier-Stokes equations and the measurement is minimised by means of an iterative adjoint-based framework. A similar adjoint based data assimilation technique was proposed by Gronskis et al (2013), where measurements are used to estimate inflow and initial conditions for a direct numerical simulation.
The approach proposed by Gronskis et al (2013) focusses mostly on usage of measurements to improve direct numerical simulations. This is illustrated by Figure 2.1, which Gronskis et al (2013) show in experimental validation of their method. The PIV measurements used in the assimilation procedure (green diamonds) already have a sufficiently high temporal resolution to capture the vortex shedding event, considering it is captured by approximately five measurements. It is argued however that the result of the numerical simulation and assimilation offers a higher spatial resolution than the PIV measurements. Note this is not the topic of the present thesis, where focus lies on increasing temporal resolution of measurements done at repetition rates below those dictated by the Nyquist criterion. The paper by Gronskis et al (2013) however has a very relevant notion on the outflow boundary condition, which is both unsteady and unknown in between subsequent measurements. To allow for the numerical simulation to advance to the next time step, Gronskis et al (2013) estimate the unsteady outflow boundary conditions using a convection equation, where u c is an estimate for mean convection velocity of the main structures in the flow. The problem of unknown boundary conditions in between measurements also applies to the time-supersampling approach proposed in the present thesis and is discussed extensively in chapter 4.
The adjoint-based techniques discussed in this section can be seen as more advanced extensions of the MI simulations by Hayase and Hayashi (1997) discussed earlier. The procedures are essentially very generally applicable, however, due to the iterative nature of minimisation of the cost function, are also significantly more expensive than the techniques discussed before. Remarkably, Lemke and Sesterhenn (2013) specifically note "that the technique demands high computational costs." Gronskis et al (2013) make a similar remark and note that, wile theoretically possible, extension of the method to three-dimensional datasets is "computationally much more intensive". Consequently, to date application of the methods to 3D datasets has not been demonstrated.
It should be remarked Gronskis et al (2013) note an interesting future application of their study may be to perform the computation of complex external flows, by considering boundary and initial conditions far downstream of the body, in contrast to the MI simulations discussed before, which considered the full flow domain. Indeed, measurements can possibly be employed to only simulate the region of interest, circumventing the requirement for simulation of the full flow domain.

Momentum Equation for Pressure Field Calculation
While not strictly a time-supersampling technique, CFD techniques have been used for calculation of pressure fields from PIV measurements. As will be discussed in detail in section 8.3, the momentum equation allows for calculation of the pressure gradient, when the spatial and temporal derivatives of velocity are known (van Oudheusden, 2013). For steady flows or calculation of the mean pressure field this approach is frequently used, as in those cases the temporal derivatives are not required and spatial derivatives are provided by PIV (see for example Regert et al 2011).
More advanced techniques for pressure calculation using the momentum equation and its divergence have been proposed by for example Hosokawa et al (2003), Regert et al (2011), de Kat andvan Oudheusden (2012) and de Kat and Ganapathisubramani (2013). These however all rely on measurements with high temporal resolution or the assumption of steady flow. This shows however that a general time-supersampling technique also has potential for application in the field of pressure calculation, as it will be able to calculate temporal derivatives from the reconstructed velocity fluctuations. This is elaborated in more detail in chapter 8.

Filling Gaps in PIV Data with Navier-Stokes Simulation
Another approach that does not focus on time-supersampling is the technique proposed by Sciacchitano et al (2012) to fill spatial gaps in PIV measurements (e.g. shadow regions) using a finite volume incompressible Navier-Stokes solver. The solver is applied in the gap region, with only a small buffer layer around it. This gives rise a boundary condition problem, because the time step required for forward time integration is smaller than the time step between two subsequent measurements. To make the problem well-posed, Sciacchitano et al (2012) use the advection model (section 2.2.2) to calculate boundary conditions at integration times between subsequent measurements.
When the flow measurements are sufficiently well resolved in time, or when the advection model provides a good reconstruction of the velocity fluctuations to provide boundary conditions at times between measurements, the approach proposed by Sciacchitano et al (2012) allows for a significant improvement over standard interpolation techniques in spatial gaps. This is illustrated in Figure 2.2, which shows the reconstruction of an artificial spatial gap in the PIV measurement of a circular jet by both linear interpolation (middle) and the numerical simulation technique (right). Clearly the numerical simulation shows significant improvement over the linear interpolation. However, the finite volume Navier-Stokes solver used by Sciacchitano et al (2012) suffers from significant numerical dissipation. For example, the the negative peak in the horizontal velocity component (Figure 2.2) is reduced by 40%. An alternative and more conservative discretisation of the Navier-Stokes equations based on the vorticity formulation of the momentum equation is used in the field of computer vision, as discussed in the next section.

Alternative Approaches in the Field of Computer Vision
In the field of computer vision, physics based methods have been proposed for fluid motion estimation. The developments have been surveyed in a recent review article by Heitz et al (2010). Notably, Corpetti et al (2002) improved motion estimation using the continuity equation and Cuzol and Memin (2005) have employed a vortex particle discretisation of the momentum equation. Cuzol and Memin (2005) estimate by a very limited set of vortex particles the vorticity and velocity fields as illustrated in Figure 2.3, which shows a measurement of a vortex at the tip of an airplane wing (left figure) and the corresponding velocity field (middle figure) estimated from a limited number of vortex particles (right figure). This reduced-order model of the flow using a limited number of particles allows for use of an advanced data assimilation technique, such as the Ensemble Kalman Filter (Evensen, 2003;Papadakis et al, 2010), as simulation is relatively inexpensive compared to full models. Related to this is also the work of Ross et al (2010), who propose a method for extraction of vortices from velocity measurements in order to define a reduced-order model of the flow field by for example a limited set of point vortices.
Motion tracking using these reduced-order models has a relevant application for flows where a small number of large scale structures dominate the flow and when there is limited interest in smaller structures. In general this is certainly not the case for the flows considered with PIV and it is desired that the result of time-supersampling has the same spatial resolution as the original PIV measurement. It is however relevant to note the choice for the vortex particle discretisation of the Navier-Stokes equations employed in the references given in this section. As will be discussed extensively in the next chapter, such a Lagrangian discretisation with vortex particles leads to a very stable and conservative simulation. This makes vortex methods an interesting alternative to classical grid-based Navier-Stokes solvers, which often suffer from significant artificial diffusion and strong stability conditions, requiring small time steps and mesh spacings (see for example Sciacchitano et al, 2012).

State-of-the-Art Time-Supersampling
Following the literature review presented in the previous section, it can be concluded that only the linear advection model by Scarano and Moore (2011) has proven, in flows were assumptions made by the model hold, to be able to reconstruct velocity fluctuations between measurements done at a rate below that dictated by the Nyquist criterion. Importantly, application to three-dimensional datasets has also been demonstrated. The model however deteriorates rapidly in flows where the assumption of frozen turbulence does not hold.
Alternative and more advanced approaches discussed have not been been applied to threedimensional datasets and most importantly, at the time of writing no substantial increase of the temporal resolution of the original measurements is demonstrated by any of the alternative methods. Therefore, while limited to flows dominated by a relatively large advection velocity and limited streamline curvature, the linear advection model is considered the current stateof-the-art technique for time-supersampling of PIV measurements. This highlights one of the main motivations for the present work, in which a time-supersampling technique which allows for accurate reconstruction for flows in general is developed.
Recall the discussion in the introduction, where it was proposed that an improved timesupersampling technique can possibly be obtained by considering a more complete description of the flow: the Navier-Stokes equations. Following the literature review in the previous sections, a vortex particle discretisation of the momentum equation in vorticity-velocity formulation is considered promising for application to measurement data, considering it is currently successfully applied in the field of computer vision to measurement data. As noted in for example Cottet and Koumoutsakos (2000) and references therein, vortex methods have excellent conservation properties and because of their Lagrangian and mesh-free character are not limited by conventional stability criteria, which is in stark contrast to the finite volume Navier-Stokes solver employed by Sciacchitano et al (2012) for reconstruction of the velocity field in spatial gaps in PIV datasets. Additionally, good accuracy of vortex methods in under-resolved cases is found by Cottet et al (2002), who compared spectral and vortex methods in three-dimensional flows. This is especially relevant considering the present work does not aim for improvement of the spatial resolution of PIV measurements, but only the temporal resolution.
It should be noted that vortex methods are essentially only applicable to incompressible flows (Eldredge et al, 2002). However, this assumption is accepted, as there is already a realm of applications within the field of incompressible flows. Indeed, most of the references in the present chapter consider incompressible flows only. Also, in the field of computer vision vortex methods have presently been applied to two-dimensional datasets and as mentioned before, no extension to three-dimensional datasets is demonstrated. Therefore, in the next chapter literature on vortex methods is reviewed to determine suitability of different families of vortex methods for application to three dimensional measurement data and to select the most suitable implementation for the present thesis work.
In the present chapter the advantages and disadvantages of different vortex methods are discussed, considering the specific application to PIV measurement data. First, in section 3.1 a fundamental introduction into vortex methods is given. From this introduction, it follows a distinction can be made between mesh-free and hybrid vortex methods, which is discussed in detail in section 3.2. Subsequently, extending the discussion to three-dimensional flows allows for a further distinction between vortex particle and filament methods, as detailed in section 3.3. Finally, in section 3.4, treatment of solid interfaces in the measurement domain is discussed.

Introduction to Vortex Methods
Vortex methods can be found in literature dating back to the early 1930s, with most notably the hand-calculations by Rosenhead (1931), who was the first to compute the evolution of a vortex sheet by discretising it into elemental vortices. For a complete historical overview of the development of vortex methods following the work of Rosenhead, the reader is referred to the book on vortex methods by Cottet and Koumoutsakos (2000). The 'modern' vortex methods discussed in this chapter originate from the late 1970s and first review papers including discussions on numerical simulation of three-dimensional flows using vortex methods were written by Leonard (1980), Leonard (1985) and Anderson and Greengard (1985).
Vortex methods are based on discretisation of the vorticity transport equation in a Lagrangian form. The vorticity transport equation is a vorticity-velocity formulation of the incompressible Navier-Stokes equations, Figure 3.1: Illustration of surface S i containing the black particle at x i .

Overview of Vortex Methods
x j x j + ∆t ⋅ u(x j ,t) with Equation (3.1) follows directly when the curl operator is applied to the incompressible momentum equation. For completeness, its derivation is given in Appendix A.
Vortex methods are very suitable for calculation of viscous flows, as discussed extensively in for example Cottet and Koumoutsakos (2000) and Barba et al (2005). The viscous terms are however neglected in the present work to allow for reverse time integration. The desire for the latter will become apparent in chapter 4. For the purpose of the present discussion, it is important to note that neglecting viscosity does not limit the range of application of the time-supersampling technique to inviscid flows only. The measured flow field itself may and in most relevant cases will have been affected by viscous processes. It is only for the relatively short integration time in between two subsequent measurements that the viscous terms are considered negligible in the measurement domain. This same assumption is typically also also made in the field of pressure evaluation from PIV measurements (e.g. van Oudheusden et al 2007 and Ghaemi et al 2012), as will be discussed in more detail in section 8.3.
In order to illustrate the Lagrangian character of vortex methods, in the remainder of this section a two-dimensional flow in the absence of body forces is considered. In the case of two-dimensional flow, the vorticity vector reduces to ω = [0, 0, ω z ] T and therefore the vorticity transport equation under the aforementioned assumptions reduces to Two-dimensional vortex methods differ from grid-based Navier-Stokes solvers by the discretising vorticity not on a grid, but by a sum over N vortex particles, where δ is the two-dimensional Dirac delta function and α h i the (constant) particle strength. The latter is an estimate of the circulation around a particle located at x i . For example, when the particles are initially regularly distributed at distance h in each direction (Figure 3.1), with S i = h 2 . To satisfy equation (3.4), the movement of each vortex particle is governed by an ordinary differential equation (Leonard, 1980), This implies vortex particles are advected with the local flow velocity as illustrated in Figure 3.2. Velocity and vorticity are related through a Poisson equation, which follows from application of the curl operator to equation (3.3), The solution of equation (3.8) allows, with appropriate boundary conditions, for calculation of the divergence free flow field. It should be remarked it is not trivial that the solution of equation (3.8) yields a divergence free velocity field and for a proof thereof the reader is referred to Daube (1992) and Clercx (1997). Starting from an initial vorticity field ω(x, 0) = ω 0 and vortex particle distribution, the divergence free flow field can be calculated by equation (3.8), which is used to update the vortex particle locations and accordingly vorticity, after which the integration procedure is repeated.
A distinction between different classes of vortex methods can be made based on the method used for solution of equation (3.8). In case the flow field is at rest at infinity, this equation can be solved using an infinite domain Green's function as discussed in Leonard (1980). Alternatively, Christiansen (1973) proposed to solve the Poisson equation on a grid, which results in a hybrid method where the grid values of velocity are interpolated to the Lagrangian vortex particles. These different solution methods correspond respectively to 'mesh-free' and 'hybrid' vortex methods. In the next section the differences between these methods are explained in more detail.

Mesh Free and Hybrid Methods
Assuming the flow is at rest at infinity, a truly mesh-free and infinite domain vortex method can be obtained when calculating velocity using the Biot-Savart integral, which for a 2D flow and in the absence of interior boundaries reads (Leonard, 1980) u where e z = [0, 0, 1] T . It should be noted a similar integral can be derived for 3D flows (see for example Leonard, 1985). Substitution of equation (3.5) into (3.10) allows for calculation of the velocity at a location x by a summation of the velocities induced by each vortex particle, (3.11) To illustrate this equation, the velocity field induced by a single vortex particle and by two vortex particles is illustrated in Figure 3.3, which shows the expected spinning motion of the flow around an imaginary axis through the vortex particles. This formulation is especially useful for infinite domain problems, with a relatively small number of vortex particles. For example, in He and Zhao (2009) the mesh-free formulation is used to model rotor wake dynamics. Indeed, particles only need to be initialised in parts of the domain occupied by vorticity.
To allow for efficient solution of equation (3.11), instead of looping over all vortex particles, the Fast Multipole Method is employed in modern mesh-free vortex methods (see for example Cheng et al, 1999 andCocle et al, 2008). This technique essentially groups vortex particles that lie close together and treats them as a single source. It should be noted that vortex particles tend to gather in regions with strong velocity gradients and therefore need to be frequently reinitialised on regular locations to preserve accuracy for long-time simulations (Cottet and Koumoutsakos, 2000).
An alternative vortex method is the Vortex-in-Cell (VIC) method proposed by Christiansen (1973). This method solves equation (3.8) on a computational grid using a Poisson solver and subsequently interpolates the velocity values from the grid to the Lagrangian vortex particles. Note the latter interpolation step does not lead to increased computational cost compared to the mesh-free methods, considering in modern vortex methods particles are typically regridded to the grid node locations at each integration time step (Cottet et al, 2002). Also, despite the development of fast summation techniques for mesh-free vortex methods, in bounded domains the solution of equation (3.8) on a grid using a Poisson solver based on the Fast Fourier Transform is in many practical situations at least one order of magnitude faster (Cottet and Poncet, 2003). Simulations by Cottet and Poncet (2003) show that only when vortex particles occupy approximately less than 10% of the computational box required by a grid-based Poisson solver, mesh-free vortex methods become more effective than hybrid vortex methods. Similar conclusions have been reported in other recent literature by for example Brown (2000), Chatelain et al (2008) and Ossmani and Poncet (2010).
Considering the relatively small measurement domains of tomographic PIV, the regions occupied by vorticity generally span large part of the measurement domain, if not the full domain. The measurement domain often does not reach up to the the free stream flow (see for example Ghaemi and Scarano, 2011), violating one of the main assumptions in derivation of equation (3.11). Therefore, hybrid vortex methods which solve Poisson equation (3.8) on a computational grid are considered most appropriate for application to PIV measurement data.
Before moving to the next section, for completeness it should be remarked Christiansen (1973) originally related velocity to vorticity through the vector potential ψ, This has in 2D the advantage that only one Poisson equation needs to be solved for the scalar stream function, where equation (3.8) needs to be solved for both velocity components. In 3D, both formulations however require solving three Poisson equations. Because formulation of velocity boundary conditions on the PIV measurement domain boundary is more natural with equation (3.8), where the velocities can be specified directly, this formulation will be used.
In the next section the extension of vortex methods to 3D flows is discussed. In 3D, the distinction between mesh-free and hybrid methods and the discussions in this section are still valid, but also an additional distinction between vortex particle and vortex filament methods can be made.

Extension to Three Dimensions
Short after the introduction of two-dimensional vortex methods, the methods were extended to three-dimensions (see for example Leonard, 1985). Where the vortex stretching term is zero in the two-dimensional case, it is non-zero in three-dimensional flows and can be identified at the right hand side of the vorticity transport equation, again in the absence of body forces, (3.14) Analogously to the two-dimensional situation, vorticity can be discretised by a sum over vortex particles (Leonard, 1985), which in three-dimensional flows carry the particle strength vector α h i (t). Analogously to the 2D situation (Figure 3.1), when the particles are initially located on a regular grid, the initial particle strength vector is calculated from an initial condition on vorticity ω 0 (see for example Giovannini and Gagnon, 2006), Because the vortex stretching term can produce local vorticity intensification and reorientation, the particle strength vector is time dependent and to satisfy the vorticity transport equation is governed by an ordinary differential equation, It should be remarked various alternative formulations of equation (3.17) can be found in literature. For example, the observation that leads the the so-called the mixed scheme (Cottet and Koumoutsakos, 2000), which allows for computational savings because of symmetry of the matrix (∇ + ∇ T )u, and the transpose scheme used by for example Winckelmans and Leonard (1993), which has the advantage of conserving total vorticity. Note the original formulation does not guarantee conservation of total vorticity as after discretisation, vorticity is not guaranteed to be divergence free (Winckelmans and Leonard, 1988). A second conservative scheme used by for example Cottet (1996), Walther and Koumoutsakos (2001) and Cocle et al (2008) is, Recent numerical results discussed in Cottet and Koumoutsakos (2000) show the transpose and mixed schemes actually yield less accurate results than the classical scheme (3.17), which is expected to be caused by the fact that the transpose scheme may amplify initial disturbances of the divergence free constraint (Speck, 2011). Furthermore, as discussed in Cottet and Koumoutsakos (2000), it is expected that when particles are regridded frequently, the other two schemes yield similar accuracy. For the present application to measurement data, the classical scheme (3.17) is chosen over the formulation of (3.21), as the latter requires evaluation of spatial derivatives of both velocity and vorticity, whereas the classical formulation only requires first order spatial derivatives of velocity, which is preferred to limit measurement noise amplification, considering the present application to measurement data.
Up this point, the methods discussed discretise vorticity by a sum over Lagrangian vortex particles. In three dimensions, an alternative discretisation can be made using a finite number of vortex lines, which can be either closed or extending to infinity (Leonard, 1985). The strength carried by each vortex filament is a scalar, which remains constant over time, similar to the 2D scalar particle strengths and therefore such vortex filament methods can be seen as a more straightforward extension of two-dimensional vortex methods. It was however already noted by Chorin (1980) that there is no obvious method to generate filaments at a solid boundary in a consistent manner. This is also discussed by Cottet and Koumoutsakos (2000), who note the Lagrangian character of the filaments is purely an inviscid property and that there is no clear-cut way to simulate diffusion with a filament method. In contrast, diffusion can readily be implemented in vortex particle methods, using for example techniques described in for example Cottet and Koumoutsakos (2000) and Barba et al (2005). Additionally, note that the vortex filaments are infinite or closed lines. As mentioned earlier, tomographic PIV measurements typically consider relatively small measurement volumes, requiring special treatment for vortex filaments crossing the domain boundary. Because of this, and because the possibility to include effect viscosity in the future is desirable, the vortex particle discretisation is selected for the present application.

Treatment of Solid Interfaces
Obtaining reliable PIV measurements near solid walls is not trivial due to the presence of reflections and the possible appearance of ghost particles in the case of tomographic PIV (Stanislas et al, 2003;Theunissen et al, 2008). Time-supersampling of PIV measurements including solid wall boundaries, either inside the measurement domain or along one of the domain boundaries, has to deal with increased measurement errors near the solid walls and possibly even with a spatial gap between the first reliable velocity vector and the interface. The latter is essentially a problem of gappy PIV (Sciacchitano et al, 2012) and a simulation may be required to fill the gap region near the interface. Because the present thesis work can be considered as the first attempt to improve upon the linear advection model by simulation of the Navier-Stokes equations, the problem of solid wall boundaries is considered outside the scope of the present thesis. It is however desired that the method can be extended to such problems in the future and therefore in this section possibilities for extension of the VIC method to include solid interfaces in the computational domain are briefly summarised.
When the computational domain contains a solid interface, the assumption of inviscid flow needs to be revisited, considering viscosity is essential for correct implementation of the no-slip condition. As mentioned earlier, vortex methods are especially suitable for modelling of viscous flows and the assumption of inviscid flow is made mainly to allow for reverse time integration. As will be discussed in the next chapter, the backwards time integration is weighted with the result of a forward time integration to obtain a continuous solution in time. Because the vorticity transport equation is not reversible when viscosity is included, an alternative method for obtaining continuity of the solution is required when viscous flows are simulated.
For inclusion of an interface in the computational domain, a body fitted grid can be used with appropriate modification to the particle-grid interpolation functions near the boundary (see for example Poncet, 2002). An interesting alternative is posed by the advent of immersed boundary methods (Mittal and Iaccarino, 2005). In general, this is a family of methods for simulation of viscous flows on grids that do not conform to the shape of the immersed and solid interfaces. More specifically, the idea is incorporated in vortex methods by enforcing the boundary conditions on interfaces through the addition of body forces (Cottet and Poncet, 2003;Morgenthal and Walther, 2007). The immersed boundary techniques for vortex methods are especially interesting for the present application, as they do not require generation of body fitted meshes and allow for relatively easy incorporation of complex boundary conditions. Their accuracy can not always compete with body-fitted methods (Cottet and Poncet, 2003), but in view of the present application to PIV data, which often is of lower quality near interfaces, this may be considered acceptable. Therefore, in case the proposed super-sampling technique is to be extended for application to measurement data including solid interfaces, it is recommended that the immersed boundary technique is investigated.
From the discussions in the present chapter it is concluded that a hybrid vortex particle method, the Vortex-in-Cell (VIC) method, is the most appropriate vortex method for application to PIV measurements for the purpose of time-supersampling. In next chapter, the time-supersampling procedure using VIC is outlined.

Time-Supersampling by VIC Simulation
Presented at the 10 th International Symposium on PIV (Schneiders et al, 2013) Journal paper submitted to Experiments in Fluids currently under review In this chapter the outline of the time-supersampling technique using Vortex-in-Cell (VIC) simulation, selected in the previous chapter, is proposed. First, in section 4.1 the procedure itself is outlined, after which in section 4.2 estimation of the domain boundary conditions in between measurements is discussed. Finally, in section 4.3 a method to a-posteriori estimate the reconstruction accuracy is proposed.

Time-Supersampling Procedure
Examples of tomographic PIV measurements with sufficient spatial resolution for evaluation of the velocity gradient tensor were given in chapter 2. Considering such measurements are possible, a PIV measurement can be used directly as initial condition for VIC simulation on a computational domain equal to the measurement domain, as the initial vorticity field can be calculated from the three-dimensional velocity measurement. This is considered a much more direct and computationally efficient approach to time supersampling than simulation of the complete flow field as done by Measurement-Integrated simulations (section 2.2.1), but requires estimation of boundary conditions in between measurements as will be discussed in section 4.2.
With an initial velocity measurement at time T 1 as initial condition for the VIC simulation, the measured velocity field can be integrated forward in time until the next measurement at time T 2 is reached, as illustrated in Figure 4.1 (top figures). Subsequently integrating from time T 2 to a third measurement at T 3 , using the measurement at time T 2 as initial condition leads to a discontinuity at time T 2 . To obtain continuity of the solution, a simple approach is to apply a correction, linearly weighted over time, based on the difference between the solutions backward integration forward integration backward integration forward integration at each measurement time. On the other hand, significantly more complex alternatives are available from the field of data assimilation. In chapter 2 it was however already discussed that such techniques come at high computational costs. For example, the Ensemble Kalman Smoother (Evensen and van Leeuwen, 2000) is a powerful technique allowing for incorporation of estimates of both measurement errors and errors in the numerical simulation, but requires a large number of simulations for propagation of the uncertainties on each velocity vector through the non-linear equations.
An interesting and rather inexpensive alternative is suggested by Scarano and Moore (2011), who consider the weighted average of both a forward and backward projection of the advection equation between two subsequent measurements (section 2.2.2). Analogously, when considering a VIC simulation, the weighted average of a forward and backwards time integration between two subsequent measurements can be used to obtain continuity of the solution. This only requires two simulations between consecutive measurements, but is considered a significant improvement over a single forward integration, as will also be shown in the experimental validation (chapter 6). Recall from the previous chapter that backward time integration is possible thanks to the assumption of inviscid flow. Therefore, the weighted average of a forward time integration u h f and a backward time integration u h b between two measurements at times T 1 and T 2 is calculated to obtain continuity, for T 1 ≤ t ≤ T 2 and with Note that the resulting velocity field u h w is divergence free, because the divergence operator is linear and both u h f and u h b are divergence free. It should be remarked that in the remainder of this thesis work the subscripts are dropped and unless stated otherwise u h refers to the weighted average of both a forward and backward integration.
The time-supersampling procedure proposed above is limited by a maximum integration time, i.e. a maximum time separation ∆T between two subsequent measurements. A theoretical maximum integration time can be derived from the simplified case of scalar convection, as the time it takes for the information contained in the measured flow field to be convected out of the measurement domain, where L is the length of the measurement domain in direction of the dominant convective velocity V c . The inverse directly yields the requirement on the minimum acquisition frequency, Whether or not the theoretical maximum integration time can be reached depends largely on the quality of the boundary conditions, which need to be estimated in between consecutive measurements as discussed in the next section.
It should be remarked that when time-integration cannot be done up to the next measurement because the time separation is too large (e.g. uncorrelated instantaneous measurements), the method can still be applied to obtain information about temporal velocity fluctuations in a short time τ around an instantaneous measurement, as also illustrated in Figure 4.1 (bottom figures). As will be discussed chapter 8, when validated this may allow for calculation of pressure fields from instantaneous tomographic PIV measurements.

Domain Boundary Condition Treatment
The value of velocity is known at the boundary of the domain only at the times measurements are available. Therefore, the governing equations are ill-posed unless an estimate of the intermediate boundary conditions is formulated. To investigate this problem in more detail, consider forward time integration starting from a measurement at time T 1 . As discussed in the previous chapter, from the initial measurement, vorticity ω h (x, T 1 + ∆t) in the measurement domain Ω can be calculated without the need for boundary conditions. It is for solution of Poisson equation (3.8) that velocity at time T 1 + ∆t needs to be prescribed on the domain boundary ∂Ω.

Because the Poisson equation is linear, it can be split into four separate equations. With
the solution of equation (3.8) at T 1 + ∆t equals the sum of the solutions of The sum of u 1 and u 2 describes the initial flow field at time T 1 , and u 3 and u 4 the change in velocity in the measurement domain over one integration time step, due to respectively the change in vorticity in the domain and the change in velocity on the domain boundary. Clearly, equations (4.6) and (4.7) can be solved directly as both velocity and vorticity at T 1 are known. Furthermore, Poisson equation (4.8) can be solved as ω h (x g , T 1 + ∆t) is calculated using the VIC procedure. On the other hand, equation (4.9) accounting for changes in velocity on the domain boundary cannot be readily solved, because u h (x, T 1 + ∆t) is not measured.
When the velocity at a boundary is known, such as in case of uniform flow (i.e. free stream condition), it can be applied directly. When the velocity is varying slowly along the boundary a point-wise linear interpolation between consecutive measurements can be considered an acceptable approximation of the velocity boundary condition at an intermediate time in between two measurements. If unsteady conditions involve the boundaries of the domain, a possibly better approximation can be made using the advection model by Scarano and Moore (2011), outlined in section 2.2.2. The advection model was used previously for estimation of unsteady boundary conditions by Sciacchitano et al (2012) for simulation of the Navier-Stokes equations to fill spatial gaps in PIV measurements (section 2.2.5) and a similar method was used by Gronskis et al (2013) for estimation of unsteady outflow boundary conditions (section 2.2.6). Analogously, the advection model is also used in the present project for estimation of unsteady boundary conditions.
To illustrate an important limitation of the approach outlined above, consider an inviscid and incompressible flow, accelerated and decelerated by a piston moving right and left in an infinite tube (Figure 4.2). In this case both a linear interpolation and the advection model are unable to reconstruct the changes in velocity on the domain boundary (dashed in Figure 4.2), as they depend entirely on the motion of the piston between the measurements at T 1 (left figure) and T 2 (right figure). Also, in this situation, ω = 0 in the entire measurement domain and hence by equation (4.8) velocity u 3 = 0. On the other hand, u 4 = 0, but unknown, which implies that in this rather extreme case the velocity fluctuations in the measurement domain cannot be reconstructed. The piston case is a rather extreme case, however, a similar situation can occur for example in case of experiments on a model with moving boundaries, such as the robot dragonfly DelFly (de Clercq et al, 2009). The problem is an information theoretical one, which can only be correctly solved when the model is extended with, in the piston case, information about movement of the piston, or in general, with information about the external flow dynamics. The proposed time-supersampling technique simulates the interior flow dynamics and exterior effects at a frequency higher than what can be reconstructed from the measurements are not included in the reconstructed velocity fields. It should be remarked this can be considered 'the best possible' when only the measurements are taken into account. Indeed, the state-of-the-art advection model (section 2.2.2) suffers from the same limitation.
Considering there is a degree of uncertainty in assessment of the quality of the estimated boundary conditions, a method to assess the reconstruction quality in the computational domain is desired. In the next section such an estimate is proposed, based on the difference between the forward and backward time integrations.

Reconstruction Error Estimation
The proposed vortex method for time-supersampling is a physics-based interpolation technique. The accuracy of the reconstructed velocity field is dependent upon a number of factors: the measurement accuracy of the initial conditions and of the estimated boundary conditions (see the previous section); the quality of the flow model (e.g. the validity of neglecting viscosity); and the numerical error in the model discretization.
The accuracy of the reconstruction is estimated using the Euclidean norm of the difference between the solutions produced by forward and backward time marching. Such an a-posteriori error estimate has already been proposed by Scarano and Moore (2011), however, an improved estimate is proposed by averaging differences at all reconstructed time steps within one measurement interval, weighted according to proximity to a measurement plane: where u h f and u h b are respectively the solutions of forward and backward time integration and where the weights w n are calculated from the tent function, where t n is an integration time between two consecutive measurements at T 1 and T 2 and ∆T = T 2 − T 1 .
In the present work this error estimate is evaluated and compared to the actual reconstruction error calculated using reference measurements, as part of the experimental validation of the proposed time-supersampling technique in chapter 6. But first, in the next chapter the numerical implementation of the proposed technique is presented.
Presented at the 10 th International Symposium on PIV (Schneiders et al, 2013) Journal paper submitted to Experiments in Fluids currently under review The numerical implementation of the proposed time-supersampling technique using VIC is presented in this chapter. The time-integration procedure is outlined in section 5.1, after which the criteria to set the integration time step are given in section 5.2. The numerical implementation of the VIC method is verified by simulation of an analytical vortex ring in section 5.3.

Time-Integration Procedure
The VIC method is explained in detail in literature, see for example Cottet and Koumoutsakos (2000). In this section the numerical treatment will be explained, focusing on the changes to the method for application to measurement data.
Starting from a measurement u m (x g , T 1 ) at time T 1 on a grid with N nodes at locations . . , x g N }, the time integration procedure consists in summary of the steps given in Figure 5.1. In the sections below, the steps in this figure are discussed in detail.
Step 1: Initial Condition The initial condition for the VIC procedure is the vorticity ω h (x g , T 1 ) at the initial measurement time T 1 , determined from the initial PIV measurement u m (x g , T 1 ). The proposed technique does not try to increase the spatial resolution of the measurements and the computational grid used is chosen equal to the measurement grid with nodes at locations x g . As discussed by for example Foucaut and Stanislas (2002) and Raffel et al (2007), three Step 1.
Calculate ω h (x g , t + ∆t) on grid u h (x g , t + ∆t) different techniques allow for calculation of vorticity from velocity u m (x g , T 1 ) measured by PIV. The first one is based on finite difference operators applied to neighbouring grid points, the second one is based on the analytical derivation of an analytic function fitted through the velocity field (see for example Abrahamson and Lonnes, 1995) and the third one is based on calculation of the circulation. Alternatively, following the method proposed in Ruan et al (2001), the vorticity vectors may be calculated directly from digital particle images measured by PIV. This has the advantage over the aforementioned methods that use adjacent velocity vectors for gradient evaluation, that uncertainty in the vorticity values originates directly from the particle image noise, instead of from uncertainty in the velocity vectors. Such a technique is however not well established and in the present study, the vorticity field will be calculated from velocity data to keep focus on the proposed time-supersampling technique.
From the study by Foucaut and Stanislas (2002) follows that good accuracy can be obtained using a second-order central difference scheme, in comparison to the alternative techniques for calculation of vorticity from velocity vectors. This is linked to the fact that the central difference scheme presents a cut-off frequency equal to that of the PIV data. Furthermore, higher order schemes are not desired as their respective noise amplification factor (Foucaut and Stanislas, 2002) increases with increasing order. Therefore, to set the initial condition, the vorticity vectors ω h (x g , T 1 ) are calculated using second order central differences from u m (x g , T 1 ).
Step 2a: Divergence Free Velocity Field The divergence free velocity field u h (x g , T 1 ) is calculated by solving Poisson equation (3.8) using the Fast Fourier Transform (FFT). This is in VIC literature a standard and efficient technique for solution of the Poisson equation on Carthesian grids (Cottet and Koumoutsakos, 2000) and the technique is well documented in for example (Sbalzarini et al, 2006). At the  initial time T 1 boundary conditions are taken from the measured velocity field at time T 1 . At subsequent times, boundary conditions are approximated as explained in section 4.2.
It should be noted that the calculation of the divergence free velocity field is intrinsic to the VIC procedure and that it is not calculated for noise reduction or data pre-conditioning. Several studies have proposed to impose the divergence free condition for reduction of measurement noise (see for example Schiavazzi et al 2013 andde Silva et al 2013). There is however at the time of writing no established approach for minimisation of the divergence error in PIV measurements. The effect of obtaining the divergence-free reconstruction by the VIC procedure on the measurement noise is considered outside the scope of the present study and will not be further discussed here. On the other hand, the effect of measurement errors on the time-integration procedure is assessed in chapter 7.
Step 2b: Particle Initialisation The measurement domain is divided into N equal volumes V p i centred around each grid node, as illustrated in Figure 5.2 for the 2D case with a Carthesian grid. In case of a uniform Cartesian 3D measurement grid, V p i = h 3 , where h is the distance between grid nodes.
A vortex particle is initialised in each volume V p i at the location of the grid node, As mentioned by for example Kudela and Kosior (2011), clustering of particles tends to occur in regions with high velocity gradients. To prevent that particles gather together, the vortex particles are reinitialised on the grid nodes at each integration time step. As discussed in chapter 3, this also means that grid values do not have to be interpolated to the particles, as at each integration time step their locations initially coincide with the grid nodes.
Step 3: Particle Advection and Strength Update The vortex particle location x p i (T 1 + ∆t) at the next time step is calculated by solving equation (3.7) using first-order Euler integration, where the integration time step ∆t is set as discussed in section 5.2. After particle advection, the particles no longer coincide with the grid nodes, as was illustrated in Figure 3.2.
Particle strength is updated to account for vortex stretching. The updated particle strength, α h i (T 1 + ∆t), is calculated by solving equation (3.17) using first-order Euler integration, where the spatial velocity derivatives are calculated using a second-order central-difference scheme.
Step 4: Interpolation from Particles to Grid Vorticity ω h (x g , T 1 + ∆t) is calculated on the grid nodes by , withφ an interpolation kernel. Originally Christiansen (1973) used the Cloud-in-Cell interpolation technique (Birdsall and Fuss, 1969), where vorticity is assigned to the grid points surrounding the particle through a piecewise bilinear interpolation kernel. It was however noted already by Christiansen (1973) that better results could be obtained with a smooth interpolation and by taking into account more grid points. Similar conclusions were found by Ebiana and Bartholomew (1996), who assessed the performance of various interpolation kernels. Significantly improved results have been reported in the past decade using the M 4 interpolation kernel introduced by Monaghan (1985), which has been used for VIC simulations by, amongst others, Walther and Koumoutsakos (2001), Ould-Salihi et al (2000), Cottet et al (2002), Giovannini and Gagnon (2006), Morgenthal and Walther (2007) and Kosior and Kudela (2013). This third order interpolation kernel, given in Figure 5.3 for completeness, conserves total circulation and both linear and angular impulse and will also be used for the present VIC implementation.  (2001); for example, the black particles assign their vorticity to different grid nodes.
The interpolation kernel (5.5) is compact, chosen such that each particle influences 64 surrounding grid points (a 4×4×4 point kernel). Because one particle is initialised at each grid point, multiple particles will influence the vorticity at each grid point. Solving (5.4) directly by looping over all particles results in an inefficient loop. As suggested by Walther and Koumoutsakos (2001), the procedure can be vectorized when particles are ordered such that consecutive particles are at least 5 grid points apart in each spatial direction ( Figure 5.4). In this order, particles stored consecutively in memory assign their strength to different subsets of the mesh and the assignment procedure can be vectorized.
To correct for the fact that vortex particles near the domain boundaries redistribute their vorticity to grid points outside of the computational domain, Cottet and Koumoutsakos (2000) suggest an iterative method. To allow for a more computationally efficient solution of the equations, in the present study a vorticity boundary condition is overlayed over grid points on the domain boundary. The vorticity on the boundary is calculated from the velocity boundary conditions, which are approximated as explained in section 4.2.
Step 5: Updated Velocity Field Calculation The velocity u h (x g , T 1 + ∆t) is subsequently calculated from ω h (x g , T 1 + ∆t) by solving equation (3.8) as in step 2 above. The integration procedure is repeated to calculate the velocity field at time T 1 + 2∆t from the velocity field at T 1 + ∆t and onwards.
The procedure for backwards time-integration between measurements at T 1 and T 2 is similar to the procedure outlined in this section, however, it starts from the second measurement at T 2 and integrates backward in time with time step −∆t.

Integration Time Step
Inviscid vortex methods are not limited by classical CFL type stability conditions. As noted by for example Cottet and Poncet (2003) and Chatelain et al (2008), Lagrangian particle advection is linearly unconditionally stable and for non-linear stability it is sufficient that particles do not collide. This imposes a limit on the time step related to the rate of strain of the flow, This condition does not involve the mesh size and is often less restrictive than a classical CFL condition.
To allow for vectorisation and to simplify implementation of the assignment procedure, where vorticity is calculated on the grid from advected vortex particles (step 4), it is assumed that each particle does not move more than the mesh cell size h in each spatial direction over the integration time ∆t. This leads to a restriction on the maximum time step, which is identical to the standard CFL condition, A third pragmatic constraint on the maximum time step follows simply from the desired sampling frequency: In view of the stability conditions discussed in this section, it should be remarked that reducing the mesh spacing h does not result in growth of instabilities when the conditions above are satisfied. However, when the spatial resolution is insufficient for evaluation of the spatial derivatives, this generally leads to underestimation of the initial vorticity and consequently to underestimation of the amplitude of the reconstructed velocity fluctuations. This poses a requirement on the spatial resolution of the tomographic PIV measurement, which can be challenging in the turbulent flow regime at high Reynolds numbers (Westerweel et al, 2013).

Numerical Verification of the VIC Code
The proposed time-supersampling technique is developed for application to experimental datasets. However, with appropriate initial and boundary conditions, also pure numerical simulations are possible with the VIC implementation presented in this chapter. In this section, the implementation of the presented VIC method, in a Matlab program written by the author, will be verified by numerical simulation of a translating vortex ring. This test case was already considered for vortex methods in Leonard (1985) and recently VIC simulation results have been reported on the translational velocity of a vortex ring by Kudela (2012, 2013), which will be employed as independent reference. The analytical thin-cored vortex ring is described in section 5.3.1 and the results of the numerical simulations are given in section 5.3.2. Subsequently, in section 5.3.3 the results of a convergence study are given, in order to verify the method is indeed first order in time and second order in space.
It should be remarked that this section is solely presented as verification of the correct implementation of the VIC method presented in this chapter. Complete validation of the full time-supersampling procedure is done using experimental datasets in chapter 6.

Analytical Vortex Ring
The simulation of a thin-cored inviscid vortex ring in a quiescent fluid has been used recently as a verification experiment for parallel and GPU based VIC methods by Kudela (2012, 2013). In the present study, this test case has been chosen, as it allows for verification of the full 3D inviscid code and comparison of the results to both an analytical approximation and another recent VIC implementation. As discussed in for example Green (1995), when the initial vorticity distribution of a thin-cored vortex ring is chosen to be uniform with ω φ = K, the translational velocity of the ring can be approximated by Kelvin's formula (Kelvin, 1867), where Γ 0 is the initial circulation of the vortex ring, R 0 the ring radius and 0 the core radius, with 0 /R 0 1 (see also Figure 5.5). The constant K is set such that the initial circulation of the vortex ring equals,  Kosior and Kudela, 2012 Results from the present VIC code Figure 5.6: Velocity of the vortex ring at t = 6 as a function of initial circulation, in comparison to results from Kudela (2012, 2013) and Kelvin's formula.
The vortex ring is modelled on the computational grid by setting the initial vorticity to (5.11) and to zero elsewhere in the domain. Note that in derivation of equation (5.9), it was assumed that the vortex ring translates steadily and without change of form (Green, 1995). It is however expected that during the VIC simulation, the initially uniform vorticity distribution in the vortex core will change over time, especially considering the discontinuous initial condition and because the vortex ring is an unstable structure (Green, 1995;Kosior and Kudela, 2013). Equation (5.9) therefore serves only as an indication of the translational velocity of the ring and should not be taken as an exact reference.

Simulated Vortex Ring Velocity
The parameters for the VIC simulation are set equal to those used by Kudela (2012, 2013). The computational domain is a 2π × 2π × 2π box with the origin in its centre at x = [0, 0, 0] T . Furthermore, the ring radius R 0 = 1.5 and core radius 0 = 0.3. The initial circulation is varied, Γ 0 ∈ [0.25, 1.5], and the simulation is run up to T = 6 with ∆t = 0.01. However, instead of periodic boundary conditions, no through-flow boundary conditions are used, as the code was developed for use with Dirichlet boundary conditions only. Also, because of limited computational resources, a grid with 59 nodes, instead of 129, in each direction is used. Equation (5.9) shows the translational velocity of the vortex ring increases with increasing initial circulation. This is confirmed by the simulation results given in Figure 5.6. However, it can be seen that the results from the VIC codes are lower than the result from equation (5.9), especially for higher values of circulation. The same result was found by Kudela (2012, 2013), who explained the higher translational velocities predicted by Kelvin's formula by the fact that the formula was derived for constant and uniform vorticity in the core, hence having a discontinuity in vorticity at the surface of the ring, whereas during the simulation the vorticity distribution smooths and and changes over time.
The results of the present simulation are within 5% of the results obtained by Kudela (2012, 2013). This is considered acceptable given the differences in especially the boundary conditions and mesh size. Furthermore, kinetic energy dropped over the duration of the simulation less than 1%, compared to 2% in Kosior and Kudela (2012).
The results in this section show good correspondence with the reference results, indicating the VIC code has been implemented correctly. Before turning to a more quantitative validation of the time-supersampling procedure using experimental datasets, in the following section spatial and temporal convergence of the VIC code is studied.

Convergence
In this section both temporal and spatial convergence of the method is investigated. Results with different time steps ∆t and grid node spacings h are compared to a reference solution u r generated using respectively a small step and a fine mesh. The discretization error is calculated as the root mean squared error over the domain at the final time T , All simulations are done using an initial circulation of Γ 0 = 0.5 and with the same computational domain and vortex ring parameters as in the previous section. For the study of temporal convergence the time step was varied between ∆t = 2 · 10 −1 and 10 −3 and results were compared to a reference solution computed with ∆t = 10 −4 . For the reference solution, the vortex ring was simulated up to T = 0.2 on a computational grid with 49 nodes in each direction. Time stepping is done using first order Euler integration and therefore first order convergence is expected. This is confirmed by the results plotted in Figure 5.7 (left figure), where the slope of the error plot equals 1 on the log-log scale.
Spatial convergence is studied by varying the number of grid nodes in each spatial direction between 19 and 99. The simulation is run over 100 time steps with ∆t = 0.1 (T = 1). The reference solution is calculated using 109 nodes in each direction, with fourth order central differences for all spatial derivatives. The M 4 kernel used for interpolation is of the third order (see for example Giovannini and Gagnon, 2006;Morgenthal and Walther, 2007), however, the proposed VIC method uses second order central differences in the interior domain and therefore second order convergence is expected. This is directly confirmed by Figure 5.7 (right figure), where the error in the interior domain is plotted by the solid dots and the slope of a line through the points on the log-log plot is 2.
In order to verify third order convergence of the interpolation kernel, the simulation was also done using fourth order central differences. For those simulations, the interpolation kernel is limiting the order of spatial convergence and therefore third order convergence is expected. This is confirmed directly by Figure 5.8, where the slope of the line through the points is 3.
This concludes the numerical verification of the implementation of the VIC code presented in this chapter. It should be stressed that this is only a preliminary verification and full experimental verification of the time-supersampling technique is presented in the next chapter.
Preliminary results were presented at the 10 th Int. Symp. on PIV (Schneiders et al, 2013) Journal paper submitted to Experiments in Fluids currently under review The proposed time-supersampling technique has been implemented as a Matlab program and in the present chapter the technique is validated against two datasets from existing tomographic PIV experiments, following the assessment method explained in section 6.1. The first validation experiment is that of a turbulent wake behind a NACA airfoil (section 6.2) and the second is a circular transitional jet (section 6.3). After the method has been validated for application to tomographic PIV measurements, applicability of the method to two-dimensional datasets is assessed in section 6.4.

Assessment Method
In section 4.3, an a-posteriori error estimate was proposed, based on comparison of the forward and backward integrations. However, for validation of the time-supersampling method, a comparison of the reconstructed velocity u h with the true velocity is desired. Because time-supersampling is essentially reconstruction of velocity fluctuations in between PIV measurements, by definition these are not measured by the PIV system. To allow for validation despite this, Scarano and Moore (2011) suggested for validation of their linear advection model, a validation based on super-sampling a sub-sampled dataset of PIV measurements.
Scarano and Moore (2011) consider reference tomographic PIV measurements u m as the 'ground truth.' The PIV measurement rate f m is subsequently decreased artificially to f s by sub-sampling the original data set (Figure 6.1) according to a Sub-Sampling Factor (SSF), sub-sampled (f = fs < fm) super-sampling by linear interpolation Figure 6.1: Illustration of reference measurements sub-sampled with SSF=6; the dashed line represents a linear interpolation through the sub-sampled measurements.
Subsequently, time-supersampling is applied to the sub-sampled velocity fields and the resulting velocity fields are compared to the original time series. This is first done qualitatively by visual comparison with the measured velocity fields in the spatial and temporal domains. For a quantitative comparison, the metric used for the reconstruction uncertainty in the reconstructed velocity is the Euclidean norm of the difference in reconstructed velocity with the reference measurements and normalised by a reference velocity, where t n are the times at which reference measurements u m are available. Additionally the method is compared with two other reconstruction techniques: the first is point-wise linear interpolation; the second makes use of the advection equation . The latter is considered the current state-of-the-art time-supersampling technique (see chapter 2) and is therefore selected as benchmark technique.

NACA Airfoil Trailing Edge Flow
In this first case a fully developed turbulent boundary layer leaves the trailing edge of a NACA 0012 airfoil. The setup of the tomographic PIV experiment is illustrated in Figure 6.2. Measurements are performed on a 40 cm chord airfoil at a free-stream velocity of 14 m/s. The resulting Reynolds number based on chord length is Re C = 386,000 and based on the boundary layer momentum thickness is Re θ = 1,300. Full description of the experiment is given in Ghaemi and Scarano (2011) and salient features of the dataset are given in Table 6.1.
To fulfil the light demand of the tomographic PIV experiment, the laser beam is introduced in a multi-pass light amplification system to increase illumination intensity in the measurement domain (Ghaemi and Scarano, 2011). The measurement volume extends for approximately 5 times the boundary layer thickness along the streamwise and spanwise directions. The volume thickness of 8 mm is slightly less than that of the boundary layer (δ 99 = 10 mm) before the  trailing edge. The tomographic recording system is composed of four Photron Fastcam CMOS imagers. Despite this advanced setup, the tomographic PIV system could not be used at a measurement rate beyond 2,700 image pairs/second, which is at the limit to capture the temporal velocity fluctuations. At this measurement rate, a fluid parcel at the outer edge of the boundary layer travels approximately 5 mm (0.5δ) between subsequent measurements and will appear approximately 10 times in the measurement domain before exiting it downstream. Accordingly, the maximum SSF calculated by equation (4.3) from the minimum frequency of the sub-sampled dataset equals approximately SSF max = 10.
The flow in the near wake of the airfoil features the typical pattern of fully developed turbulent boundary layers with elongated low and high speed streaks (Figure 6.3). Most vortical structures, cane vortices and hairpin packets, concentrate across the low-speed regions, where ejections are more pronounced (Tomkins and Adrian, 2003). In the near wake region behind  the airfoil, Taylor's hypothesis of 'frozen turbulence' holds reasonably well and supersampling by advection has been demonstrated to be an accurate approach for increasing temporal resolution . Therefore, this experiment has been chosen to verify whether the VIC method is able to reconstruct the time-resolved velocity fields in conditions where the advection-based model can be considered as a reference technique.

Divergence Free Velocity Field
The first two steps in the time-supersampling procedure are respectively calculation of vorticity ω h (x, T 1 ) and the divergence free velocity field u h (x, T 1 ) from a PIV measurement u m (x, T 1 ).   It should be recalled this is done as an essential part of the VIC time-intergration procedure and not for pre-conditioning of the data. Because of measurement and model errors, generally u h (x, T 1 ) = u m (x, T 1 ). To assess quality of the divergence free reconstruction, contour plots for all three velocity components in the spanwise plane at time t = T 1 = 0 have been plotted in Figure 6.4, with in the top row the reference measurements and in the bottom row the divergence free reconstruction. From visual inspection of this figure it can directly be concluded that the velocity fields are in good correspondence. Note boundary conditions for solution of Poisson equation (3.8) at time t = 0 are taken directly from the initial measurement. The latter is not divergence free and therefore it is expected that in a small region near the domain boundaries, the divergence of velocity in the reconstructed field is non-zero. This is readily confirmed by Figure 6.5, which shows the velocity divergence in the spanwise plane for both the measured and reconstructed velocity field. As expected, divergence of velocity in the interior domain is practically zero in the reconstructed field.
For a more quantitative assessment of the reconstruction quality, Figure 6.6 shows the reconstruction error, calculated using equation 6.2 from the divergence free reconstruction of 100 reference measurements (right figure). The reconstruction error is also shown for each separate velocity component and all the reconstruction errors have been normalised with respect to the free stream velocity (respectively the first three figures). On the boundary, velocity is taken from the measurements and therefore the error goes to zero on the boundaries of the domain. Furthermore, typically in tomographic PIV the velocity component in direction of the laser sheet thickness is least reliable (de Silva et al, 2013). This is confirmed by Figure 6.6, which shows that the errors in the streamwise and spanwise velocity components are similar (≈ 0.8%) and that slightly larger errors are present in the w-component (≈ 1%).
The error calculated from all three velocity components is approximately 1.4%. As this is the difference between the measured velocity field and the divergence free velocity field, and because the low-speed wake flow is expected to be divergence free, this is considered to be an estimate for the measurement error. It should be remarked this estimate is only slightly in excess of the measurement error of 1.3% predicted in Scarano and Moore (2011)  measurements conducted on the same experiment. The slightly larger measurement error is expected, considering the errors in tomographic PIV are expected to be higher.
In chapter 2 it was noted the spatial resolution of PIV measurements is generally sufficient to allow for calculation of the spatial derivatives of velocity. When spatial resolution is limited, vorticity peaks are expected to be underestimated, leading to underestimation of the amplitude of the velocity fluctuations. This is confirmed by Figure 6.7, which shows the divergence free reconstruction when the velocity field is spatially sub-sampled such that the mesh spacing becomes 3h (red / grey lines). With this severely limited spatial resolution, the second order finite difference scheme is unable to correctly evaluate the velocity spatial derivatives and the resulting divergence free velocity field is strongly smoothed. On the other hand, Figure 6.7 shows, considering the estimated measurement error, at the original mesh spacing sufficiently accurate reconstruction (black lines) is possible.
Based on the discussion in this section, steps 1 and 2a of the VIC procedure discussed in the previous chapter are validated for this first test case, as indeed a divergence free velocity field is calculated with an acceptable reconstruction error compared to the reference measurements. The actual time-supersampling procedure is validated in section 6.2.3, but first in the next section the integration time step is determined.

Integration Time Step
The minimum integration time steps set by the different criteria discussed in section 5.2 have been calculated from the measurement data and are listed in Table 6.2. The limit on the maximum vortex particle displacement is most critical, requiring the minimum integration time step to be ∆t min = 0.095 ms, which is approximately (20f m ) −1 . In comparison, Sciacchitano et al (2012) quote a time step requirement of (2, 500f m ) −1 for simulation of the shadow regions in a flow over the same airfoil, but with a lower chord based Reynolds number of 100,000.
To assess the reconstruction quality at this minimum integration time step, the tomographic PIV dataset is sub-sampled with a factor 8 (SSF = 8) and subsequently one single forward integration between two consecutive measurements from the sub-sampled dataset is done. The resulting velocity time history is shown in Figure 6.8 for a single point in the centre of the measurement domain. In this figure, also a linear interpolation between the subsampled measurements is shown. As discussed in section 6.1, for these simulations, only the measurements at times t · f m = 0 and t · f m = 8 have been used and clearly the linear interpolation is unable to reconstruct any of the relevant velocity fluctuations. On the contrary, it can be seen in the time history plots that already using only a forward time integration and with the minimum integration time step, rather good correspondence with the reference measurements (solid circles) in between the sub-sampled measurements is obtained.
Full assessment of the reconstruction quality is given in the next section. In the present section, first the effect of reducing the integration time step is studied. Simulation results with a reduced integration time step are also plotted in Figure 6.8. Visually the results do not change substantially when reducing the time step and it is expected that the reconstruction error, calculated with respect to the reference measurements, for each of the forward integrations is significantly larger than the differences between the solutions with the different integration time steps. Indeed, the reconstruction error for all single forward integrations at SSF = 8 is approximately 3%, relative to the free stream velocity, whereas the difference between the solutions with ∆t/∆t min = 1/2 and 1/4 is approximately 0.3%. Because visually the results  are converged with ∆t/∆t min = 1/3 and because the reconstruction errors are an order of magnitude larger than the truncation errors made with time integration, no further reduction of the integration time step is considered necessary. Furthermore, it should be remarked that to limit storage requirements, the results are stored with ∆t = 0.05 ms, corresponding to a super-sampling frequency of 20 kHz.

Results of Time-Supersampling
The time history of the streamwise velocity component at a point in the centre of the measurement domain is shown in Figure 6.9 for the cases of SSF = 4 and SSF = 8. The temporal fluctuations of the streamwise velocity component have a typical time scale slightly less than 1 ms (Ghaemi and Scarano, 2011). In this respect the current measurement rate of 2.7 kHz is just sufficient to describe these fluctuations and it is therefore expected that even at low SSF the reconstruction by linear interpolation will fail to reconstruct temporal fluctuations. This is readily verified by inspection of the time series in Figure 6.9, where it is shown that when the measurements are sub-sampled from 2.7 kHz to 675 Hz (SSF = 4), the linear interpolation fails to reconstruct the temporal velocity fluctuations. Most fluctuations of time scales smaller than 2 ms are clipped. In contrast, both the advection and VIC method appear to be adequate, as they return a temporal evolution of the velocity fluctuations with a good qualitative agreement with the measured reference data series. It should be noted that the reconstruction by advection passes exactly through the sub-sampled data. In contrast, the reconstruction obtained with the VIC method does not, since a divergence free velocity field is produced from the actual datasets. For reference, the divergence free reconstruction has been calculated at each time instant and is also plotted in Figure 6.9 (blue crosses). As expected, the VIC reconstruction a passes through these exactly at each time instant corresponding to a sub-sampled measurement.
When the sampling frequency is reduced further to 337.5 Hz (SSF = 8), the reconstruction quality of both the advection and VIC models degrades, which is confirmed by the summary of reconstruction errors given in Table 6.3. In calculating these relative errors the free stream velocity is used as reference velocity and a dataset of 100 snapshots is considered. For comparison to the reconstruction error, the last column shows the a-posteriori error estimate ε VIC,a−post. , which is based on the backward and forward integration only and calculated using equation (4.12).
To illustrate the calculation of the a-posteriori error estimate, both forward and backwards time-integrations are plotted in Figure 6.10, where the time histories of the streamwise and The results show that as expected, for higher SSF, the difference between the forward and backward integrations (shaded areas in Figure 6.10) becomes larger. Accordingly the a-posteriori error estimate is larger for higher SSF, as also confirmed by the results in Table 6.3. Comparison of the error estimate to the direct error shows the values are in good agreement, which indicates that the a-posteriori error estimate from equation (4.12) may be taken as an estimator for the reconstruction quality. For low SSF, the error estimate is smaller than the direct error, which is for a large part due to the presence of measurement noise.
Even for a potentially exact reconstruction by supersampling, the measurement noise gives a lower bound on the error. In the present experiment, the measurement uncertainty has been estimated to be approximately 1.4% of the free-stream velocity based on the analysis of the velocity divergence (section 6.2.1). It is therefore expected that the reconstruction error is overestimated by a similar amount. Comparing the reconstruction errors in Table 6.3 to the estimated measurement error shows the additional error introduced by super-sampling using a linear interpolation is significantly larger than the measurement error. The situation is improved substantially by using both the advection and vortex models. For super-sampling using the linear advection model it remains in the same order up to approximately SSF = 4 and the additional error from super-sampling using the vortex model remains smaller than the measurement error for SSF < 6.
As expected from the previous discussion, the spatial distribution of the reconstruction error of VIC and linear advection (Figure 6.11) are comparable for SSF = 4. With a coarser time resolution (SSF = 8) the direct error estimate calculated for the VIC technique is slightly lower than that observed for the advection model. The situation near the inflow and outflow boundaries shows an increase of the error and a more equivalent result between the two methods, as the boundary conditions used for the VIC reconstruction are equal to the solution of the advection equation on the boundaries. Inside of the measurement domain the vortex method improves upon the advection model, indicating that perfect frozen turbulence conditions cannot be assumed. Especially closer to the trailing edge of the airfoil, the VIC model improves upon the advection model. This is explained by the velocity deficit of approximately 30% in the inflow region of the measurement domain, which invalidates the assumptions made by the advection model.
In summary, the results by Scarano and Moore (2011) have been confirmed and indeed the advection model yields accurate results for this experiment where the assumption of frozen turbulence holds to a large extent. The VIC model has been shown to yield approximately equivalent results up to SSF = 4 and slightly better results for higher SSF. In the next section a case where the assumption of frozen turbulence does not hold is considered, to verify if also in such a case the VIC method is able to improve upon the advection model.

Transitional Circular Jet
The second case is that of a transitional circular jet ( Figure 6.12), which flows from a nozzle with exit diameter D = 10 mm installed at the bottom wall of a water tank. The tomographic PIV measurements are performed with a jet exit velocity of 0.45 m/s, yielding a Reynolds number Re D = 5, 000 based on the jet diameter D. The relevant experimental parameters are listed in Table 6.4 and a full discussion of the experimental parameter setup is available  in Violato and Scarano (2011).
For the jet flow, less accurate results were reported by Scarano and Moore (2011) using the advection model and therefore it is chosen as a second test case for the VIC method. The free shear layer is dominated by the Kelvin-Helmholtz instability and strong vortices are shed at approximately 30 Hz (Strouhal number of 0.67). The hypothesis of frozen turbulence is not valid in this region, which as discussed in Scarano and Moore (2011) invalidates the assumptions made by the advection model.
The flow field is visualised in Figure 6.12 (left figure), which shows isosurfaces for the different velocity components and vortical structures as identified by the Q-Criterion (Hunt et al, 1988). The tomographic measurements are done with a repetition rate of 1,000 Hz, whereby the vortex shedding (at 30 Hz) is captured with very high temporal resolution. Based on the jet exit velocity, the theoretical maximum SSF is estimated by equation (4.3) to be approximately SSF max = 90. However, as will be shown, the practical maximum SSF equals 60. This value still represents quite an extreme as it corresponds to a sampling frequency of only 16.7 Hz, which is about one quarter of the frequency required by the Nyquist criterion for reconstruction of the temporal velocity fluctuations.
Following the same outline as in the previous section, first the calculation of the divergence free flow field is verified (section 6.3.1). Subsequently, the integration time step is chosen (section 6.3.2) and the results of time-supersampling are presented and analysed (section 6.3.3).

Divergence Free Velocity Field
The divergence free reconstruction of an instantaneous velocity field is visualised by the Q-Criterion (Hunt et al, 1988) at time t = 0 in Figure 6.12 (right figure), next to the original reference measurement (left figure). Visually the results are nearly identical, indicating sufficient spatial resolution of the reference measurements. To allow for a more quantitative  comparison, the average reconstruction error relative to the jet exit velocity in the axial plane has been calculated with equation (6.2) and plotted in Figure 6.13. The errors are largest in the shear layer, where they are in the order of 1-2.5% of the jet exit velocity, corresponding to 0.1-0.25 voxels per measurement. Because the difference between the divergence free reconstruction and the measurements is largest in the shear layer, the divergence in velocity for the measurements is also expected to be largest in the shear layer. This is readily confirmed by Figure 6.14, which shows the average divergence of velocity in the axial plane.
To illustrate the relatively larger error levels at x/D = ±0.5, Figure 6.15 shows the divergence free reconstructions of the radial (left figure) and axial (right figure) velocity components along the line x/D = −0.5 in the axial plane. This figure shows that when mesh spacing is reduced to 4h, the velocity field is excessively smoothed (red / grey lines), but that with the original mesh spacing h only minor differences can be observed (black lines). Especially in the inflow region, y/D < 1.5, the divergence free reconstruction is a more smooth representation of the measured flow field. Also, the amplitude of the axial velocity fluctuations is slightly reduced by the divergence free reconstruction, with a reduction in the axial velocity peak at y/D ≈ 1.5 of approximately 4%. For the purpose of the present project, this error is accepted, considering good correspondence between the reference measurement and the divergence free reconstruction in large part of the domain (see also Figure 6.12). Improvement of this first step in the time-supersampling procedure is however envisaged to lead to improved reconstruction quality and could possibly be achieved by pre-processing the measurements using techniques for minimisation of the divergence error proposed in Schiavazzi et al (2013) and de Silva et al (2013).

Integration Time Step
The minimum integration time steps set by the different criteria discussed in section 5.2 are listed in Table 6.5. As for the previous experimental test case, the limit on particle displacement is most critical, setting the minimum time step to ∆t min = 0.75 ms. Analogously to the investigation for the previous test case, the effect of decreasing the time step further is investigated. Figure 6.16 shows the result of forward VIC simulations in a point in the shear layer with SSF = 30. Further reduction of the integration time step below the minimum time step shows minor improvement and the results are visually converged with ∆t/∆t min = 1/4. As argued in the previous section, further reduction of the integration time step is not considered necessary as the reconstruction error is significantly larger than the temporal discretisation   error. The integration time step is therefore set to ∆t/∆t min = 1/4 for all further simulations presented in this section. To reduce storage requirements, the simulation results are stored with a time separation of 1 ms, making the sampling frequency of the super-sampled dataset equal to the measurement frequency.

Results of Time-Supersampling
The time histories of the reconstructed radial and axial velocity components are plotted in Figure 6.17 at a point within the shear layer (x/D = -0.5, y/D = 2, z/D = 0), for SSF = 30 (top figure) and SSF = 45 (middle figure). This location has been selected as it corresponds to a maximum in reconstruction error of the divergence free field (Figure 6.13) and because it is sufficiently far away from the boundary to exclude boundary effects, as discussed later in this section. The bottom figure additionally shows the velocity time history at the same x and z coordinates, but closer to the outflow region (y/D = 3.5) with SSF = 45. Considering that the period of the temporal fluctuations equals approximately 30 ms (corresponding to vortex shedding) when the data is sub-sampled with SSF = 30, the resulting frequency becomes half of that dictated by the Nyquist criterion. Not surprisingly, the linear interpolation is unable to reconstruct the velocity fluctuations as illustrated in Figure 6.17. Furthermore, the advection model completely fails to reconstruct the fluctuations at this low sampling frequency, as the assumption of linearized trajectories according to frozen turbulence does not hold in this flow case.
In contrast, the VIC method is able to accurately reconstruct the velocity fluctuations, even though the sampling frequency is reduced to the vortex shedding frequency. Because the divergence free step slightly damped the the axial velocity peak around y/D = 1.5 and considering the jet exit velocity of approximately 0.5 m/s, a minor reduction in amplitude of the axial velocity peak is expected at y/D = 2, approximately 10 ms after the first measurement. This is confirmed by the results in Figure 6.17, which shows a reduction in the axial velocity peak of approximately 4%. This corresponds to the figure found earlier in the divergence free reconstruction and indicates numerical dissipation in the time-integration procedure is limited, especially considering Sciacchitano et al (2012) reported amplitude reductions of 40% with a finite volume Navier-Stokes solver to fill gaps in PIV measurement. Further reduction of sampling frequency to 22 Hz (SSF = 45) shows some further reduction in amplitude of the temporal velocity fluctuations for the VIC method. However, no phase distortion is noticeable.  To analyse the spatial patterns produced by the different reconstructions, the case of SSF = 30 is taken and the temporal evolution of the jet is depicted in three time instants (Figure 6.18) by visualisation of the Q-Criterion (green isosurfaces). The top plots show the reference measurements. The bottom three respectively show a reconstruction using linear interpolation, the advection model and the VIC method. At t = 0 ms the linear interpolation and advection model are equal to the reference measurements, whereas the VIC reconstruction shows minor difference as it returns the divergence free flow field. Already 5 ms later both the linear interpolation and advection model cannot reconstruct the fluctuations as expected from the previous discussions. The VIC method however captures the correct phase and location of the vortex rings. Further at t = 15 ms the reconstruction is at the middle time between two sub-sampled measurements. At this time instant the error is expected to be largest, because both forward and backward integrations are used. The advection model breaks apart the vortex rings and is clearly inadequate to reconstruct this velocity field with strong flow shear and rotation. In contrast, the VIC method appears to reconstruct the velocity field rather accurately without any visible distortion.
The reconstruction error, normalised using the jet exit velocity, is given in Table 6.6 at the position in the free shear layer corresponding to the first two plots in Figure 6.17. From the analysis on the divergence free velocity field, the largest errors are expected to be found in the shear layer. This is confirmed by Figure 6.19, which shows the spatial distribution of the reconstruction error for SSF = 30. From both the  maintains reconstruction errors within 3-4% up to SSF = 40 (f s = 25 Hz), confirming its suitability for these type of flows. Only for very low SSF ≤ 5, the advection model gives smaller reconstruction errors than the VIC model. This is expected as the advection model is forced to the reference measurements and instead of to the divergence free reconstruction. Because the advection model is employed to provide the unsteady outflow boundary conditions, the reconstruction errors of the VIC model are larger near the outflow boundary. Considering the solution is the weighted average of both a forward and backward integration, an estimate of the region affected by the incorrect boundary conditions can be made by where L e is the distance from the boundary affected by the boundary condition. This equation is derived by assumption of convection of turbulent fluctuations on the boundary with a constant convective velocity V c . Setting V c to the jet exit velocity yields for SSF = 30 that L e /D = 0.7. A line drawn at this distance from the outflow boundary divides the domain in a region expected to be unaffected and and a region affected by the outflow boundary condition, as shown in Figure 6.20. In the same figure the results for other SSF have also been plotted. As can be seen, equation (6.3) is a rather accurate estimate of the regions eroded by the incorrect boundary conditions. Notably, because equation (6.3) can be solved before actually performing an experiment, it can be used to obtain an estimate for the required minimum measurement domain.
The three-dimensional velocity fields were reconstructed from the PIV particle images using the fluid trajectory correlation (FTC) algorithm as introduced by Lynch and Scarano (2013), because the technique yields time accurate measurements with significant error reduction. The time-supersampling procedure benefits from the relatively low measurement errors, as illustrated by Figure 6.21, which shows only small differences between the forward and backward integrations at the mid-times in between measurements. The a-posteriori error estimate calculated from the differences between the forward and backward time-integrations is listed in Table 6.6, The estimate underestimates the direct reconstruction error by approximately one percent point, which is for a large part due to the presence of measurement errors. The spatial distribution of the a-posteriori error estimate plotted in Figure 6.22 resembles the actual direct reconstruction error. Notably, the higher reconstruction errors near the outflow boundary are accurately estimated, indicating the error estimate can be used for a-posteriori estimation of the eroded region near the boundaries, in addition to the estimate provided by equation (6.3).
As discussed earlier, the maximum value of SSF that can be theoretically applied to this experiment is approximately 90. The practical maximum value of SSF however cannot be increased beyond 60 due to the uncertainty in the boundary conditions, which in this case cannot be estimated accurately by the advection model. Increasing SSF from 50 to 60 already increases the reconstruction error substantially (Table 6.6). However, this rather extreme SSF of 60 corresponds to a measurement rate reduced from 1,000 Hz to only 16.7 Hz, which is approximately half the frequency at which vortices are shed and one quarter of the required sampling frequency according to the Nyquist criterion.
Up to this point in this chapter, it has been demonstrated that the proposed method yields slightly more accurate results than the advection model in flow cases where frozen turbulence can be assumed and that the proposed method is able to accurately reconstruct velocity fluctuations in shear layers with strong curvature of streamlines. In such flows the advection model deteriorates quickly and significant improvement can be obtained by the proposed VIC method. Notably, it has been shown that a stable and accurate reconstruction is possible when reducing the sampling frequency to the vortex shedding frequency (SSF = 30).
The results in the present section already demonstrate a significant improvement of the proposed model with respect to the state-of-the-art advection model. To complete the validation, a sensitivity and robustness analysis is presented in the next chapter. First however, in the next section applicability of the method to two-dimensional datasets is discussed.

Two-Dimensional Measurements
Planar PIV remains one of the most popular measurement techniques for fluid flows (Westerweel et al, 2013) and even though laser power requirements are less critical, application of the timesupersampling technique to planar PIV measurements is considered relevant. In this section the VIC method is applied to two-dimensional slices extracted from the two tomographic PIV experiments discussed in the previous sections in order to assess applicability of the method to two-dimensional datasets. This approach is chosen over validation against planar PIV experiments, as it allows for direct comparison of the results to the results presented above, which considered the full three-dimensional datasets. Before the time-supersampling results for planar data are given, first the required changes to the governing equations are discussed in section 6.4.1. Subsequently the actual results are presented in section 6.4.2.

Implications for the Governing Equations
Time-supersampling is essentially only relevant for unsteady flows. Turbulence in such flows is a strongly three-dimensional phenomenon and as discussed extensively in Boffetta and Ecke (2012), two-dimensional turbulence is strictly speaking never realised in nature. Planar PIV measurements of flows are therefore in most practical applications never a full description of the flow field, as they lack the out-of-plane components. In this section the 2D VIC model is derived from the full 3D equations to show how this can lead to reduced reconstruction quality.
By definition, in truly two dimensional flows the 'out-of-plane' velocity component w is zero, where it should be remarked that This illustrates that the divergence free flow field can only be calculated from planar data, by solution of this Poisson equation and under the assumption of two-dimensional flow, if both | ∂ωy ∂z | | ∂ωz ∂y | and | ∂ωx ∂z | | ∂ωz ∂x | and if the out-of-plane velocity gradients of u and v are small. It should be noted hence not only the flow itself needs to be sufficiently two-dimensional, but also the measurement plane needs to be aligned such that the aforementioned inequalities hold.
Modelling the flow with a two-dimensional version of the VIC method has not only implications for calculation of the divergence free flow field, but also for the time-integration steps. The inviscid vorticity transport equation for the scalar vorticity ω z becomes because the vortex stretching term (ω · ∇) u = 0 in 2D, implying particle strength remains constant over an integration time step and particles are only advected with the local flow velocity. The out-of-plane velocity component w = 0 and therefore particles are advected only in-plane, . (6.9) The discussion in the present section exposes that, obviously, it is not trivial or even possible to model a three-dimensional flow using a two-dimensional model. With the current approach the out-of-plane velocity and the out-of-plane derivatives are neglected entirely. Approaches to improve upon this, possibly using optimisation procedures treating the out-of-plane components as unknowns, are considered outside the scope of the present thesis. It is therefore not expected that the proposed technique allows for accurate time-supersampling of 2D measurement data, considering most practical applications are in unsteady and turbulent flows which are inherently three-dimensional. To assess this, in the next section the reduction in reconstruction quality is evaluated when considering only a two-dimensional slice from the tomographic PIV experiments discussed earlier in this chapter.

Time-Supersampling in 2D
First the case of the turbulent wake described in section 6.2 is considered. The spanwise plane at y = 0 is extracted from the tomographic PIV dataset and serves as 2D dataset. From the full 3D dataset, all three components of the vorticity vector have been calculated at time t = 0 and are given in the histograms shown in Figure 6.23. The histograms show directly  that all three components of the vorticity vector are approximately equally distributed in the chosen plane, indicating that the turbulent wake flow is, as expected, 3D and that none of the components may be neglected. Indeed, the divergence free reconstruction calculated using only ω z with the two-dimensional Poisson equation given in the previous section the lacks essential velocity peaks as shown in Figure 6.24 (right figure) in comparison to the reference measurement (left figure). For completeness, the divergence free field calculated using the full 3D dataset is also given (middle figure), which shows good correspondence as already discussed in the previous sections.
Calculating the difference between the divergence free reconstruction and the measurements in the planar case using 1% when the plane is reconstructed using the full three-dimensional dataset. Note the latter is lower than the reconstruction error reported in section 6.2.1, as errors in the out-of-plane component are not included. The reconstruction quality when time-supersampling is therefore expected to be reduced by at least a factor three. This reduction in reconstruction quality is clearly visible in Figure 6.25, which shows the super-sampling results with SSF = 4 using both the 2D and 3D VIC model in the centre of the measurement domain. As discussed earlier in this chapter, the 3D VIC code allows for accurate reconstruction of the temporal velocity fluctuations at this SSF. On the contrary, the 2D code is able to reconstructions some of the velocity fluctuations, but fails to capture all fluctuations accurately. The linear advection model has also been implemented in two-dimensions and the timesupersampling results are also plotted in Figure 6.25. Interestingly, this model is able to yield rather accurate results for this flow case and the model does not seem to suffer much from absence of the full 3D velocity data. This is explained by the fact that the flow and plane considered is one where the assumption of frozen turbulence holds very well, the dominant   When time-supersampling slices extracted from the axial plane in the transitional jet experiment, the advection model is expected to fail also in 2D. On the contrary, a more reasonable VIC reconstruction is expected because in the axial plane of the jet, the in-plane vorticity components are significantly smaller than the out-of-plane component (see Figure 6.26). This is readily confirmed by Figure 6.27, which shows the time-supersampling results at a point in the shear layer with SSF = 30. Note however that the axial velocity fluctuations are damped substantially, which leads to a significantly higher reconstruction error in the 2D case ( V IC,2D = 7.6%) than when the full 3D dataset is considered ( V IC,3D = 2.3%).
The results in this section show that the proposed technique can be applied to two-dimensional datasets, but that quality of the results depends strongly on whether or not the flow case considered is sufficiently two-dimensional. The reconstruction quality of the two-dimensional model using planar data deteriorates quickly in three-dimensional flows and application of the proposed technique to planar data is not recommended, considering most relevant flow cases are inherently three-dimensional.
In this chapter the proposed time-supersampling technique has been validated using experimental tomographic PIV datasets. When applied to three-dimensional datasets, the method has shown an improvement over the state-of-the-art advection model in all flow cases. Most notably, in a shear layer with strong streamline curvature the proposed technique is able to accurately reconstruct velocity fluctuations, even when the dataset was sub-sampled to below the frequency dictated by the Nyquist criterion. No stability issues were found, despite real PIV datasets including measurement noise were used. However, to assess the effect of increased measurement errors, in the next chapter the robustness of the method to measurement errors is assessed.
The proposed time-supersampling technique was validated in the previous chapter using real-world tomographic PIV experiments. These, being measurements, included measurement errors and it was already noted the measurement errors form a minimum to the attainable reconstruction quality. In this chapter the effect of measurement noise is studied in more detail to assess robustness of the proposed technique. First, in section 7.1 the assessment method is outlined. Subsequently, sensitivity to two main types of measurement errors are studied: firstly outliers (section 7.2), and secondly precision errors (section 7.3).

Assessment Method
In the recent review paper by Westerweel et al (2013) it is noted that "the study of ghost particles (Elsinga et al, 2011) is currently the only available theoretical framework for understanding measurement errors in tomographic PIV." Because very limited literature is available on errors in tomographic PIV measurements and a specific error investigation is considered outside the scope of this thesis, this chapter focusses on the two main types of measurement errors expected to occur: outliers (section 7.2) and precision errors (section 7.3).
In the previous chapter, reference measurements u m were considered as ground truth, but in reality included measurement errors m . Hence with u t the true velocity field. To assess the effect of measurement noise, in this chapter artificial errors a are added to the transitional jet tomographic PIV measurements, sumption that m a , i.e. an artificial error is added, which is substantially larger than the measurement error. The jet case is selected, because it was shown to have a relatively small measurement error, as it could be processed with the FTC technique (section 6.3). Time-supersampling of sub-sampled 'noisy' datasets is subsequently done with the proposed technique and the results are compared to the original reference measurements in order to assess the effect of the artificial error level on the reconstruction quality.

Sensitivity to Outliers
Velocity vector fields obtained by PIV can contain outliers, as discussed in for example Raffel et al (2007). These outliers are spurious velocity vectors, deviating in both amplitude and direction from nearby velocity vectors (Westerweel, 1994) and can occur in blocks when interrogation window overlap is used when processing the PIV images. They occur when there is no dominant correlation peak and consequently the magnitude of a component of a spurious vector can range anywhere between ± 1 2 W s , with W s the interrogation volume size. Outlier detection and replacement is listed as one of the first post-processing steps for PIV data in Raffel et al (2007), and can be done using for instance the techniques discussed in Westerweel (1994) and Westerweel and Scarano (2005).
As discussed in Raffel et al (2007), in challenging experimental situations spurious vectors in 5% of the domain can be expected. In this section it will be investigated whether the VIC technique is able to handle unprocessed velocity fields, still containing such a significant number of outliers. Figure 7.1 shows part of the jet measurement with 5% spurious vectors added to the data, both as single vectors and in 3 × 3 × 3 blocks. The artificial spurious vectors are generated by selecting a velocity between ± 1 2 W s for each component of the velocity vector with uniform probability.
The outliers cause, with their random direction and large magnitude, strong discontinuities in the velocity field. From the results in the previous chapter it is expected that the velocity field is smoothed in the second step of the VIC procedure, where the divergence free reconstruction of the initial measurement is calculated from the measurement data. This is confirmed by the velocity vector fields in Figure 7.1, which shows the VIC reconstructions of the velocity fields at the initial time step. In the case of single (1 × 1 × 1) outliers, the velocity vectors are largely corrected in the high speed flow region, where the difference between the outlier magnitude and the flow velocity is small. Figure 7.1 shows however the outliers have a noticeable effect in the region outside of the jet (x/D > 0.75). Furthermore, when the outliers appear in blocks of 3 × 3 × 3 vectors, the VIC code is unable to reconstruct the reference velocity field starting from the initial time step. From this analysis, it can be concluded that the method in its current form should not be used for outlier replacement, as it is unable to remove the blocks of outliers. Because the VIC method is chosen for its good conservation properties, the vortices induced by large blocks of outliers are expected to remain in the solution throughout the time integration procedure. This is confirmed by Figure 7.2, which shows the result is severely corrupted by the outliers already at t = 0 (middle figures) and continues to degrade when integrating forward in time (right figures).
The analysis in this section shows that the proposed time-supersampling method is unable to correctly replace especially blocks of outliers. Therefore, to obtain a correct velocity field reconstruction, the PIV measurement data must first be subject to outlier detection and replacement, following techniques proposed by for example Westerweel (1994) and Westerweel and Scarano (2005). It should be remarked that outlier detection and replacement is, while not trivial, at the moment a standard first step in post-processing PIV data (see for example Raffel et al, 2007). Further investigation in use of the VIC procedure and especially the divergence free velocity field calculation step for correct replacement of outliers is left as a topic for further research.

Sensitivity to Precision Errors
The robustness of the VIC model in cases where the measured flow field contains a large precision error is investigated in the present section. In contrary to the outliers discussed in the previous section, which were spurious vectors that could have any amplitude allowed by the interrogation volume size, precision errors are relatively small random errors affecting the magnitude of each velocity component. For modelling purposes, in this section, it is assumed to be normally distributed, and temporally uncorrelated. The assessment is performed by super-sampling velocity fields with SSF = 30 (f s = 33 Hz) and therefore the latter assumption is justified, as the time between subsequent measurements considered is sufficiently large.
The analysis in this section considers both spatially uncorrelated and correlated precision errors. Similar as in the case of outliers, spatially correlated precision errors can occur when interrogation volume overlap is used in processing of the PIV particle images. Scarano (2013) shows in an overview of recent tomographic PIV experiments interrogation volume overlap of 75% is used most frequently. With such overlap, the precision errors of seven neighbouring vectors in one direction are expected to be correlated, with five neighbouring vectors relatively strongly correlated by sharing at least 50% of the interrogation volume used by the midvector. For the present investigation, spatially correlated precision errors are simulated by adding uncorrelated noise to a coarse grid and interpolating this to the finer measurement grid. The resulting autocorrelations in one spatial direction are given in Figure 7.3, showing as expected the wide correlation peak in the case of spatially correlated noise (right figure).
In Carr et al (2009) it is discussed that spatial correlation of precision errors essentially leads to a more smooth and continuous error distribution. This is confirmed by Figure 7.4 which shows the radial velocity along the line x/D = −0.33, z/D = 0 in case uncorrelated noise (left figure) and correlated noise (right figure) is added, in both cases such that ε m = 10% of the jet exit velocity. This error ε m is defined as the root mean squared difference between the reference measurements and the measurements including artificial noise, similar to the reconstruction error defined by equation 6.2. The VIC divergence free reconstruction is also plotted in Figure 7.4 for both cases of correlated and uncorrelated errors. Remarkably, the divergence free reconstruction handles the case of uncorrelated noise very well and significantly reduces the measurement error by 75% to approximately ε = 2.5%. On the other hand, in case of spatially correlated noise the divergence free reconstruction rather accurately reconstructs the velocity field with the artificial noise. Consequently, the error compared to the reference measurements remains rather high, ε ≈ 6.5%, but is still 35% lower than the error in the original measurements with artificial noise.
From the discussion above it is expected that the proposed method is able to handle uncorrelated  precision errors rather well. This is confirmed by the contour plots in Figure 7.5, which shows the time-supersampling result for SSF = 30 at t = 15 ms (right figure) when starting from an initial condition with uncorrelated artificial precision noise (left figure, ε m = 10%). Despite the simulation started from a measurement corrupted by a significant amount of noise, the VIC procedure is stable and able to relatively accurately reconstruct the reference flow field at t = 15 ms (middle figure). The time history of radial and axial velocity at a point in the shear layer plotted in Figure 7.6 confirms these results and shows relatively accurate reconstruction is still possible when the measurements with a significant noise level are sub-sampled to below the Nyquist frequency. The acceptable amount of precision noise is not unlimited however. When the artificial noise level is increased such that ε m = 26% of the jet exit velocity, the simulation is no longer stable as shown in the time history in Figure 7.7. However, with measurement error levels of a quarter of the jet exit velocity, the measurement would normally be discarded and not even considered for time-supersampling.
Expectations for accurate and stable reconstruction are lower for the case of spatially correlated noise following the results in Figure 7.4. The accuracy of the divergence free reconstruction, obviously, suffers from the smooth correlated noise, which is hard to detect even visually. The flow field with correlated noise (ε m = 10%) used as initial condition for supersampling with SSF = 30 is given in Figure 7.8 (left figure). Remarkably, even at this rather extreme noise level and despite significant distortions to the correct flow field (Figure 7.4), the simulation is stable and yields very reasonable results at the mid-time in between two measurements (right figure). Especially at the mid-time between two consecutive measurements, reconstruction is especially good as the reconstruction benefits from evaluation of the weighted average from the backward and forward time integration.
To investigate this in more detail, the reconstruction errors of both the individual forward and backward integrations and the weighted average have been evaluated at each time instant in between two measurements and plotted Figure 7.9. This figure shows for all artificial noise levels and for both the forward and backward integrations the expected behaviour of increasing reconstruction error over integration time. However remarkably, for high artificial error levels (ε m > 6%) the reconstruction error of the individual forward and backward integrations remains approximately constant over a short time, indicating the proposed time-supersampling technique can allow for noise reduction as will be explored further in chapter 8. Furthermore, it was already noted that the weighted average yields significant improvement in reconstruction quality, which is readily confirmed by the results in Figure 7.9. For relatively low noise levels (ε m < 6%) the error is highest at the mid-time in between two sub-sampled measurements, but for high noise levels (ε m > 6%) the error reaches a minimum in between two measurements. The VIC reconstruction with such high noise levels benefits strongly from the weighted average of the two independent forward and backward integrations, as the measurement errors are Similar to the case of spatially uncorrelated noise, there is a maximum noise level that still allows for stable simulation. For spatially correlated noise this is found at approximately ε m = 15% of the jet exit velocity. At such a noise level the measurement would normally be discarded however. Indeed, Figure 7.4 illustrates ε m = 10% is already an extreme error level. Therefore, the simulation proves to be stable also for all realistic levels of correlated precision errors.
The results in this section show that the proposed method is able to produce a stable solution, even when the precision error level is increased beyond realistic levels. Indeed, the Lagrangian vortex particle discretisation used in the time-integration procedure proves to be very stable with excellent abilities to handle measurement noise. Additionally, for high artificial noise levels the procedure reduces the error level and benefits greatly from the weighted average of the forward and backward integration. It should be remarked that from the analysis in the previous section it follows that the method is not able to correct outliers effectively. Time integration is still possible, but the outliers remain in the solution, degrading reconstruction quality. Therefore, outlier detection and replacement procedures should be applied to the measurements before time-supersampling as discussed in the previous section.
This concludes verification and validation of the proposed time-supersampling technique. In the previous chapter it was demonstrated that in all three-dimensional test cases the method yields better results than the, now previous, state-of-the-art advection model. Furthermore, in this chapter, robustness to handle significant and even unrealistic amounts of precision errors is demonstrated. In the next chapter, the discussion focusses on alternative applications of the proposed technique. Possibilities for reduction of measurement noise were already noted briefly in this chapter and will be discussed in more detail in the next chapter. Furthermore also application to instantaneous measurements is discussed, which may allow for calculation of temporal derivatives and pressure from instantaneous velocity measurements.
All discussions up to this point focussed on time-supersampling, i.e. reconstruction of temporal velocity fluctuations in between consecutive measurements. It was however already hinted at in the previous chapters that other there are also alternative applications for the proposed technique. This chapter takes a step back from the focus on time-supersampling to allow for a discussion of these alternative applications. First, section 8.1 focusses on measurement noise reduction by VIC simulation. Second, application to instantaneous and uncorrelated measurements is discussed in section 8.2. Third, in section 8.3 possibilities to calculate unsteady pressure fields from both correlated and uncorrelated measurements using VIC simulation are discussed.

Measurement Noise Reduction
In the previous chapter it was concluded that whereas the method in its current form cannot effectively replace outliers, it has very good capabilities to handle precision errors. Notably, the reconstruction error remains approximately constant over a short time when integrating between measurements (recall Figure 7.9). Following up on this, noise reduction by timesupersampling is envisaged through a sliding technique as illustrated in Figure 8.1. The top figure illustrates that when the data is sub-sampled and subsequently super-sampled with SSF = 2, starting from both the first and second measurement three approximations of the velocity at each time instant can be calculated: the divergence free field at T n and the results from respectively a forward and backward integration between T n−1 and T n+1 . Under the assumption that the measurement noise is temporally uncorrelated, averaging the three solutions is expected to lead to a noise reduction of T 10 T 11 T 12 T 13 SSF = 3: T 10 T 11 T 12 T 13 where ε m is the error in the divergence free reconstruction calculated in the first two steps of the VIC procedure (chapter 5). In general this equation becomes The errors ε are here defined as the root mean squared error following equation 6.2. To illustrate the proposed approach for noise reduction, Gaussian noise is artificially added to the transitional jet experiment as described in the previous chapter and subsequently the sliding time-supersampling technique is applied.
The resulting radial velocity field when Gaussian noise with a rather extreme standard deviation of 15% (ε m = 26%) of the jet exit velocity is added to each velocity component is shown in Figure 8.2 (middle figure). Especially in the region y/D < 2, with this noise level it is very difficult to even visually identify the shed vortices. The aforementioned noise reduction procedure has been applied with SSF = 8 to yield the noise reduced velocity field shown in  To quantify these results, the reconstruction errors are calculated by equation 6.2. When artificial noise with σ = 15% of the jet exit velocity is added to each velocity component, the artificial measurement error becomes ε m = 26%. Calculation of the divergence free flow fields already significantly reduces this error to ε m = 7%. Application of the time-sliding VIC procedure outlined above is subsequently expected by equation (8.2) to reduce the error, in case of SSF = 8, by a factor √ 15 ≈ 4 to 1.8%. As given in Table 8.1, the actual error of 1.9% is slightly in excess of this, indicating the reconstruction error is not exactly constant when performing the simulations with SSF = 8. This is expected, as in the previous chapter it was noted that the reconstruction error generally increases over time when integrating over longer It should be remarked that the above analysis by no means provides a complete validation and that the present chapter only proposes alternative areas where it is envisaged that the proposed VIC technique for time-supersampling can be applied. Also, it should be remarked the procedure proposed in this section works on the level of velocity fields obtained from cross-correlation analysis of the PIV particle images. At a lower level considering the particle images, Lynch and Scarano (2013) and Jeon et al (2013) have demonstrated a significant extension of the dynamic velocity range (Adrian, 1997) based on fluid trajectory correlation and its integration with ensemble averaged correlation. Investigation of possible integration with recent time-resolved PIV algorithms for processing particle images is recommended as a further research topic, in case noise reduction with VIC is pursued.
This section took a small sidestep and presented a noise reduction technique for time-resolved datasets with relatively high temporal resolution, by calculation of independent velocity fields at single time instants, allowing for error reduction by averaging the independent fields. Recall however the time-supersampling technique proposed in this thesis was developed for application to correlated measurements with limited temporal resolution. In the next section focus is back on measurements at low repetition rates. However, instead of considering correlated measurements, the possibility for application to instantaneous and uncorrelated measurements is discussed.

Material Derivative from Instantaneous Measurements
Intuitively, one expects time resolved measurements to be required for calculation of temporal derivatives. In the present thesis it has however been shown that velocity fluctuations could be reconstructed even if measurement data was sampled below the frequency dictated by the Nyquist criterion. In this section it is proposed to go further, by application of the proposed method to instantaneous and uncorrelated measurements. As already suggested in section 4.1 and illustrated in Figure 4.1, this may allow for reconstruction of the velocity fluctuations in a small time interval τ around the instantaneous measurement and consequently for calculation of the temporal velocity derivative. When successfully validated, this may pose a much simpler alternative to dual and multi-pulse PIV systems proposed by for example Liu and Katz (2006) and Lynch and Scarano (2013).
The procedure is essentially the same as the one for time-supersampling described in chapter 5. However, instead of integrating in between two measurements, only one integration time step is simulated with an instantaneous measurement as initial condition. To make the problem well-posed, boundary conditions are determined by a single forward time projection of the advection based model by Scarano (2013), analogously to the full time-supersampling procedure (section 4.2).
It has been recently shown that for time-resolved PIV measurements, evaluation of the material derivative from a Lagrangian viewpoint gives important advantages in terms of reduced truncation errors (Lynch and Scarano, 2013). However, in the present case the integration time step ∆t can be chosen arbitrarily small. In these conditions the Eulerian evaluation of the temporal derivative does not suffer from significant truncation errors. Therefore, with the updated velocity field, after one integration time step from a measurement at T m , the acceleration is calculated using a single-sided finite difference scheme, which allows for calculation of the material derivative by For the present discussion, the code has been applied to an instantaneous tomographic PIV measurement of the transitional jet discussed in section 6.3. Two types of boundary conditions have been applied, first constant boundary conditions (i.e. no acceleration on the domain boundary) and second advection type boundary conditions. Figure 8.4 shows the resulting temporal derivatives in the plane x/D = 0. For comparison, the temporal derivative has been evaluated from the full time-resolved dataset (left figures). As expected from the experimental validation, in the interior domain the VIC simulation can accurately calculate the temporal derivatives already when constant boundary conditions are assumed (middle figures). Furthermore, some improvement near the boundary is achieved by the advection boundary conditions (right figures). This is confirmed by the plots in Figure 8.5, which The error in velocity acceleration by the VIC simulation near the outflow boundary follows from the information theoretical problem that the velocity boundary conditions are only known at the measurement time instant. To illustrate this, recall the discussion on domain boundary conditions in section 4.2. Laplace equation (4.9) needs to be solved for u 4 , which accounts for changes in the velocity in the interior domain due to the change in velocity on the domain boundary over one integration time step ∆t. Consider now the Taylor's expansion where the temporal derivative of velocity follows from the momentum equation, Substitution of equation (8.5), with (8.6), into Laplace equation (4.9) allows for splitting of u 4 into u 4 = u 4,a + u 4,b , where u 4,a and u 4,b are solutions of From the PIV measurements, Laplace equation (8.7) can directly be solved. However, the pressure gradient on the domain boundary is not measured and equation 8.8 cannot be solved, leading to an error in the reconstructed velocity field. Similar to the case of time-supersampling, this procedure for evaluation of the material derivative from an instantaneous measurement therefore simulates only the interior flow dynamics and exterior effects are not taken into account, as they are not measured. Hence the procedure essentially is valid only when exterior effects are negligible in the interior domain. It should be remarked this is an information theoretical limitation that is applicable also when more advanced optimisation techniques are employed and it can essentially only be solved by providing the model with information about the external flow dynamics.
The rather remarkable result that from a single tomographic PIV snapshot (two particle images), the acceleration of velocity may possibly be approximated is especially relevant for pressure calculation from PIV velocity measurements. This is discussed in more detail in the next section.

Pressure from PIV using VIC
In only a decade, techniques that determine the fluid flow pressure from PIV velocity measurements have emerged (Baur and Konteger, 1999) and come to a degree of maturity that justifies their application in practical problems. These developments have been surveyed in a recent review article by van Oudheusden (2013). Obtaining detailed information on the pressure field in flows is of increasing interest, as surface pressure fluctuations play a major role in aeroelastic and aeroacoustic phenomena (e.g. Koschatzky et al, 2011). For example, for assessment of broadband noise emissions, determining the spectrum of pressure fluctuations in boundary layers or at the trailing edge of an airfoil is recognised as an experimental challenge (Pröbsting et al, 2013;Morris, 2011 (Liu and Katz, 2006;van Oudheusden, 2013). Ghaemi et al (2012) have directly evaluated those terms in a turbulent boundary layer, showing that they typically are two orders of magnitude smaller than the other terms.
Considering the differential form of the material derivative (equation 8.4) exposes the requirement for time resolved measurements through appearance of the temporal velocity derivative. It should be remarked that under the assumption of incompressible flow, the velocity field is divergence free and the acceleration term in equation (8.4) drops out of the Poisson equation. Still, especially considering the relatively small tomographic PIV measurement domains, the temporal derivative is typically required for determination of Neumann boundary conditions for pressure by equation (8.10). Traditionally, calculation of the temporal derivative of velocity requires time-resolved measurements with high temporal resolution or dual-PIV systems. The time-supersampling procedure proposed in the present thesis reduces this requirement by allowing for reconstruction of velocity fluctuations between correlated measurements. In addition, in the previous section, application of the proposed technique to instantaneous measurements was discussed and the possibility to calculate the temporal derivative of velocity from an instantaneous measurement was illustrated. It is envisaged that this may allow for calculation of the instantaneous and unsteady pressure field from an instantaneous measurement.
To demonstrate this application to instantaneous measurements, the unsteady pressure field is calculated from the single tomographic PIV snapshot from which in the previous section the temporal derivative of velocity was calculated (Figure 8.4). can be accurately calculated in the interior domain, pressure is also evaluated accurately. Also the more quantitative plot of pressure along the line x/D = 0, z/D = 0.33 in Figure 8.7 shows good correspondence of the pressure field calculated from an instantaneous measurement using VIC to the reference pressure calculated from the time-resolved dataset.
The results in this chapter show promising alternative applications for the time-supersampling approach proposed in this thesis work. It was noted earlier already that the preliminary results presented here in this chapter do not provide a full validation, such as presented in the previous chapters on time-supersampling. The discussions however highlight possibilities for further research. It is recommended that the latter focusses on both further development of the proposed techniques and on validation by independent reference measurements (e.g. pressure transducers). This and other recommendations following from the present work are detailed in the next chapter, together with the main conclusions from the present thesis work.
To reduce challenging measurement rate requirements for time-resolved tomographic PIV experiments, in the present thesis work the Vortex-in-Cell (VIC) technique for reconstruction of velocity fluctuations between correlated measurements is proposed, developed and subsequently validated against datasets obtained from existing tomographic PIV experiments. In the present chapter, the conclusions from this thesis work are first summarised in section 9.1. Afterwards, the recommendations for further work regarding improvement of the method and alternative applications are summarised in section 9.2.

Conclusions
Feasible measurement rates for time-resolved velocity measurements by tomographic PIV are strongly limited by available laser pulse energy. Following a literature review into stateof-the-art PIV techniques, it is concluded that the linear advection model by Scarano and Moore (2011) is currently the only validated technique for time-supersampling of 3D datasets obtained from tomographic PIV experiments. By time-supersampling, velocity fluctuations between measurements are reconstructed, reducing measurement rate requirements. However, the advection model deteriorates rapidly in flows with strong streamline curvature (e.g. shear layers). To improve upon this, an interpolation technique based on simulation of the Navier-Stokes equations is investigated, which leverages the available spatial resolution to increase the limited temporal resolution.
The Vortex-in-Cell (VIC) technique is proposed for solution of the incompressible and inviscid Navier-Stokes equations between consecutive measurements. The numerical solution is evaluated on a domain that corresponds to the 3D measurement volume and time-integration is performed between pairs of consecutive measurements to achieve a temporal resolution not available from the original measurements. The measured data serves as initial conditions for the VIC method simulating dynamical evolution of the measured fluid flow and continuity of the solution is obtained by the weighted average of a forward and backward time-integration. In this way the spatial resolution available by the measurements is projected into time.
The proposed VIC method is compared to linear interpolation and the advection model. Two experimental test cases are selected to evaluate the applicability of the method to real measurement conditions; the turbulent wake behind a NACA 0012 airfoil and the case of a transitional jet. In the case of the turbulent wake, where Taylor's hypothesis of frozen turbulence holds to a large extent, both the advection model and VIC method are able to accurately reconstruct temporal velocity fluctuations. Furthermore, for very low sampling frequencies, 6-10 times smaller than the original sampling frequency of the measurement, the VIC model improves slightly upon the advection model. Secondly, in the case of a transitional jet, where the free shear layer is dominated by strong vortices that roll-up under the effect of the Kelvin-Helmholtz instability, the advection model essentially fails to represent the strongly non-linear physics caused by the vortices in the shear layer. In contrast, the VIC method accurately reconstructs the velocity time series, even when the sampling rate is reduced to a fraction of the Nyquist frequency.
An a-posteriori error estimate is proposed to estimate the reconstruction quality. Good correspondence of the error estimate with the direct reconstruction error is found for all test cases. Also, the error estimate is found to give an estimate of the eroded region near unsteady domain boundaries, which follows from uncertainty in the boundary conditions between consecutive measurements.
Robustness and sensitivity to two main types of measurement noise is studied; outliers and precision errors. It is found that the reconstruction quality suffers from data corrupted by outliers with a strong deviation in both direction and magnitude from the true flow, especially when the spurious vectors occur in blocks. Consequently, measurement data should first be subjected to outlier detection and replacement before application of the proposed technique. On the other hand, excellent robustness to precisions errors was found. Both in the cases of spatially uncorrelated and correlated errors the simulation was accurate, relative to the noise level introduced in the measurements, and stable up to extreme and unrealistic error levels.
The study demonstrates that the sampling rate requirements for tomographic PIV can be strongly reduced when the data is time super-sampled using the Vortex-in-Cell method. This has two main advantages. First, reduced repetition rate requirements allow for higher laser pulse energy per measurement, making larger measurement domains feasible. Second, it is anticipated that in high speed flows (e.g. V ∼ 100 m/s), time-resolved information can still be obtained by time-supersampling measurements at a limited measurement rate with the proposed VIC technique, thereby extending the range of applicability of time-resolved 3D-PIV by tomographic PIV or similar techniques.
The thesis work additionally illustrates the potential for a wider application of the proposed VIC method. A noise reduction technique is envisaged by sliding time-supersampling, which calculates multiple approximations of the velocity at a single time instant from a time-resolved dataset. On the other hand, also application to instantaneous and uncorrelated measurements is envisaged, which can allow for evaluation of the temporal derivative of velocity and the unsteady pressure field from a single snapshot. These alternative applications have been illustrated by test cases based on the transitional jet experiment and their validation is a possible extension of the present thesis work, as detailed in the next section.

Recommendations for Future Work
Throughout the present thesis report, different possibilities for future projects are noted. These recommendations for future work can be split into two main categories: (i) improvement of the time-supersampling method and (ii) development and validation of alternative applications for the proposed technique.
First recommendations for improvement of the proposed technique are discussed. Here, the most important recommendation for a future project is extension of the method to handle flow cases including solid walls in or along the boundary of the measurement domain. It should be remarked that because of reflections, PIV measurements are typically of lower quality in the near wall region and appropriate adjustments to the procedure to account for this may be required. Also, viscosity is required for correct implementation of the no-slip boundary condition and extension of the method to include viscous effects may be necessary. When extension of the method to incorporate solid wall boundary conditions is pursued, immersed type boundary conditions are considered very promising in combination with the VIC method, considering they allow for solution of the equations on a Cartesian grid and can be used together with efficient Poisson solvers based on the Fast Fourier Transform.
A second recommended procedural improvement to the method is related to improvement of the initial condition used for the simulation. As noted in the validation, quality of the initial condition determines to a large extent the reconstruction quality. Presently the initial divergence-free flow field is calculated by the VIC procedure, but it is envisaged that preprocessing of the flow field to obtain a better estimate of the divergence free flow field leads to increased reconstruction quality. Here especially the works of Schiavazzi et al (2013) and de Silva et al (2013) for estimation of the solenoidal velocity field from PIV measurements are relevant.
The second category of recommendations for further work considers development and validation of alternative applications. As discussed extensively in the previous chapter, it is envisaged that the method can be used for noise reduction of time-resolved datasets and for calculation of the velocity material derivative and pressure from single tomographic PIV snapshots. The latter is considered most relevant, noting alternative and promising techniques functioning on the level of the PIV particle images are being developed for time-resolved datasets. Calculation of unsteady pressure fields from instantaneous and uncorrelated velocity measurements, when validated by preferably independent pressure transducers, can allow for radical simplification of pressure evaluation in comparison to alternative techniques such as time-resolved and multi-pulse PIV.
Vortex methods solve the vorticity transport equation by means of a particle discretisation and the equation forms therefore the basis for the proposed method. For completeness, its derivation from the full Navier-Stokes equation is given in this appendix equations to allow for better understanding of the made assumptions.
Under the assumption of incompressible flow the energy equation decouples from the mass and momentum conservation equations. Additionally, under the assumption of constant viscosity the mass and momentum conservation equations reduce to 0 = ∇ · u, 0 = ρ ∂u ∂t + ρ (u · ∇) u + ∇p − µ∆u + ρf .
Dividing the momentum equation by the density ρ and application of the curl operator yields The right hand side of this equation is treated term by term to rewrite it into the vorticity transport equation. Starting with the first term, because the curl and temporal derivative commute In order to rewrite the second term, note the vector relation (it should be remarked all the following vector relations are taken from Chorin andMarsden 1990 andSpeck 2011) and therefore (with ω = ∇ × u) Note that because ∇ × ∇φ = 0, (A.5) for any scalar field φ, which, because ∇ · u = 0 and ∇ · ω = ∇ · ∇ × u = 0, reduces to