Dense velocity reconstruction from tomographic PTV with material derivatives

Abstract

A method is proposed to reconstruct the instantaneous velocity field from time-resolved volumetric particle tracking velocimetry (PTV, e.g., 3D-PTV, tomographic PTV and Shake-the-Box), employing both the instantaneous velocity and the velocity material derivative of the sparse tracer particles. The constraint to the measured temporal derivative of the PTV particle tracks improves the consistency of the reconstructed velocity field. The method is christened as pouring time into space, as it leverages temporal information to increase the spatial resolution of volumetric PTV measurements. This approach becomes relevant in cases where the spatial resolution is limited by the seeding concentration. The method solves an optimization problem to find the vorticity and velocity fields that minimize a cost function, which includes next to instantaneous velocity, also the velocity material derivative. The velocity and its material derivative are related through the vorticity transport equation, and the cost function is minimized using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm. The procedure is assessed numerically with a simulated PTV experiment in a turbulent boundary layer from a direct numerical simulation (DNS). The experimental validation considers a tomographic particle image velocimetry (PIV) experiment in a similar turbulent boundary layer and the additional case of a jet flow. The proposed technique (‘vortex-in-cell plus’, VIC+) is compared to tomographic PIV analysis (3D iterative cross-correlation), PTV interpolation methods (linear and adaptive Gaussian windowing) and to vortex-in-cell (VIC) interpolation without the material derivative. A visible increase in resolved details in the turbulent structures is obtained with the VIC+ approach, both in numerical simulations and experiments. This results in a more accurate determination of the turbulent stresses distribution in turbulent boundary layer investigations. Data from a jet experiment, where the vortex topology is retrieved with a small number of tracers indicate the potential utilization of VIC+ in low-concentration experiments as for instance occurring in large-scale volumetric PTV measurements.

1 Introduction

The spatial resolution of tomographic particle image velocimetry (PIV) and particle tracking velocimetry (PTV) measurements is directly related to the seeding concentration of the tracer particles. For tomographic PIV, interrogation volumes typically include approximately 5–20 particles. For PTV analysis, the scattered velocity measurements are obtained at the particle locations only. Maximum seeding concentration and thereby spatial resolution is typically dictated by the maximum particle image density, N _p , in particles per pixel (ppp) that can be dealt with by particle triangulation or tomographic reconstruction in view of the ghost particles phenomenon (Elsinga et al. 2006; Lynch and Scarano 2015). With a four camera system, the maximum particle image density, N _p , for tomographic reconstruction is approximately 0.05 particles per pixel (ppp, Elsinga et al. 2006). When three cameras are employed, the particle image density should not exceed 0.025 ppp. The particle image density scales linearly with volume thickness, hence this reduction corresponds to a measurement volume thickness reduction of 50 % in comparison with a four camera setup. Opacity of the medium comes as an additional constraint especially for measurements in water, where excessive seeding concentration leads to diffusion of the laser light and loss of image contrast (Michaelis et al. 2010; Scarano 2013). Also, inhomogenous seeding could result from non-ideal tracing behavior and centrifugal effects in vortex cores (Raffel et al. 2004). Finally, in the case of large-scale PIV experiments in wind-tunnels that make use of helium-filled soap bubbles (HFSB, Bosbach et al. 2009; Scarano et al. 2015), a poor concentration of tracers is reported due to the limited rate of production (Caridi et al. 2015).

Particle tracking techniques could potentially offer a higher spatial resolution than tomographic PIV by avoiding the spatial filtering effect inherent to the cross-correlation analysis used for tomographic PIV. Sitou and Riethmuller (2001) discussed further the possibility of super-resolution PIV, originally introduced by Keane et al. (1995), by refining the PIV result using a PTV algorithm. For volumetric measurements, 3D-PTV poses, however, a strong limitation on the seeding concentration to allow reliable particle triangulation, with a maximum particle image density reported of 0.005 ppp for a three-camera system (Maas et al. 1993). This limit value is raised when using tomographic reconstruction instead of triangulation, which as noted above can be done up to an image source density of approximately 0.025 ppp using a three-camera system (Elsinga et al. 2006). In the case of time-resolved tomographic PTV, particle trajectories can in addition be leveraged to increase accuracy of the velocity measurements, using procedures inspired for example by the time-resolved PTV algorithm by Malik et al. (1993). Also, using the advanced iterative particle reconstruction (IPR, Wieneke 2013) triangulation method, the very recently introduced Shake-the-Box technique (Schanz et al. 2016) has demonstrated accurate Lagrangian particle tracking at seeding concentrations on the order of those used for tomographic PIV. The velocity measured with particle tracking techniques is returned at scattered locations corresponding to the instantaneous particle positions. Although some data post-processing operations can be performed on scattered data (Neeteson and Rival 2015), it is common to transport the information onto a Cartesian uniform grid before performing more general data post-processing and visualization operations. A straightforward approach is to apply tri-linear interpolation (Fig. 1, left) or to average the velocity vectors in an interrogation window analogously to PIV (Fig. 1, middle).

Fig. 1

Velocity reconstruction (blue vectors) from PTV measurements (red vectors). Linear interpolation between two PTV velocity measurements (left), PIV interrogation window (shaded blue) approach (middle) and VIC+ interpolation between two PTV particle trajectories (right, the orange vector is D u /Dt). The gray lines indicate the reconstruction grid

Tri-linear interpolation does guarantee continuity of the velocity field, but it yields a discontinuous (piecewise constant) distribution of the velocity gradient field, which prevents an accurate estimation of the vorticity field. Several interpolation methods have appeared in the literature that deal with the problem of reconstructing scattered PTV data onto a Cartesian grid with the constraint of smoothness. The adaptive Gaussian windowing technique (AGW, introduced by Agüí and Jiménez 1987) has been widely used and recently radial basis function regression (RBF, among others Casa and Krueger 2013) has been proposed. The choice of a length scale is necessary for both AGW and RBF. The length scale should be larger than the inter-particle distance such to suppress measurement noise, at the cost of low-pass filtering the result. Zhong et al. (1991) and Vedula and Adrian (2005) proposed to impose the velocity field to be divergence free (viz., incompressibility constraint). Furthermore, very recently Gesemann et al. (2016) presented the ‘FlowFit’ method. This method makes a B-spline reconstruction of the velocity field that minimizes a cost function, which can include, similar to the method proposed in this manuscript, a weighted sum of differences between measured particle velocity and acceleration, and other differences such as nonzero divergences of velocity. The weights are set based on the flow properties (e.g., incompressibility) and measurement quality (e.g., measurement uncertainty and seeding concentration). It should be remarked that techniques invoking the incompressibility constraint are applicable to volumetric measurements only, as from planar PIV no information on the out-of-plane components of the velocity gradient is available and the in-plane velocity field is rarely divergence free in applications involving turbulent flow fields.

For PIV measurements where instead of spatial resolution, temporal resolution is insufficient (i.e., sampling frequency below the Nyquist limit), Scarano and Moore (2012) have proposed a technique that leverages information available by instantaneous velocity measurements in the spatial domain to increase resolution in the temporal domain. Using a linearized advection equation, the technique reconstructs the velocity temporal evolution in between temporally under-sampled measurements. This procedure was christened ‘pouring space into time.’ It was demonstrated to reconstruct past the Nyquist limit in flows where the assumption of frozen turbulence holds. The concept was generalized by Schneiders et al. (2014) to three-dimensional incompressible flows by extending the linearized advection equation to the vorticity transport equation using the vortex-in-cell (VIC) technique. Analogously, information in the temporal domain has been employed recently for increased accuracy of measurements in the spatial domain. Fluid trajectory correlation (FTC, Lynch and Scarano 2013) and ensemble-averaged cross-correlation (FTEE, Jeon et al. 2014) have demonstrated improved correlation for time-resolved tomographic PIV. Tomographic reconstructions using multiplicative algebraic reconstruction (MART) were improved by leveraging the temporal information with motion tracking enhanced MART (MTE, Novara et al. 2010). The recently proposed sequential MTE method (SMTE, Lynch and Scarano 2015) improved the original MTE method by using a time-marching instead of a time-sliding approach. In the field of discrete particle tracking methods, Shake-the-Box (Schanz et al. 2016) has demonstrated accurate Lagrangian particle tracking at high particle image densities (similar to the limits of SMTE). These methods use the temporal history and expected continuity of the flow as a way to leverage the flow physics. Variational methods, proposed recently in the field of optical flow, apply in addition a regularization based on the Navier–Stokes equations, as discussed in the review paper by Heitz et al. (2010).

These techniques have not yet dealt specifically with increasing spatial resolution of the instantaneous velocity measurements, in the cases where the spatial resolution is limited by tracer particle seeding concentration, as discussed in the first paragraph. As information in the spatial domain has proven to allow for an increase in temporal resolution, the reversed principle of pouring time into space potentially offers an increase in spatial resolution. In the present study, it is investigated whether next to the instantaneous velocity measured by time-resolved tomographic PTV, also the velocity material derivative can be used to increase spatial resolution of the reconstructed velocity field (Fig. 1, right). The velocity material derivative is obtained by tomographic PTV from the second temporal derivative of a particle trajectory. The spatial and temporal information is linked through the flow governing equations, i.e., Navier–Stokes. The relevant applications mentioned above consider flow cases at speeds well below Ma = 0.3, and thus the incompressible Navier–Stokes equations can be used. These, however, still include pressure, which is not measurable by PIV. The pressure term is avoided in the present study by using the vorticity transport equation and the vortex-in-cell framework, as applied earlier to tomographic PIV measurements by Schneiders et al. (2014), (2016a).

The proposed method is assessed using both a simulated experiment and two tomographic PIV experiments. In particular, recently published results of PIV measurements in a turbulent boundary layer give an indication of the issues related to spatial resolution, with the turbulent velocity fluctuations being underestimated of approximately 20 % (Pröbsting et al. 2013). This is a direct consequence of the spatial modulation associated with the cross-correlation analysis used for PIV when seeding concentration is limited. The numerical assessment in the present study considers the flow conditions of the above experiment in a simulated experiment, to study the dependency of the reconstruction of turbulent velocity fluctuations on the particle seeding concentration. For experimental assessment, the existing turbulent boundary layer experiment is used to assess to which extent the proposed technique is able to recover the turbulent velocity fluctuations in a real-world experiment. In addition, the case of a well-known circular jet experiment in water illustrates the potential of the method to reconstruct individual vortices at the lower limit of seeding concentration.

2 The VIC+ method

The input to the VIC+ method is the result of a volumetric Lagrangian particle tracking velocimetry procedure, as discussed in Sect. 2.1. After particle tracking, VIC+ is applied as discussed in Sect. 2.2. The method requires closure parameters as discussed in Sect. 2.3.

2.1 Lagrangian particle tracking

Consider a time-resolved volumetric velocimetry measurement. The proposed method applies to measurement results from techniques that follow a volumetric particle tracking approach over multiple time-instants, such as 3D-PTV, tomographic PTV and Shake-the-Box. For the present study, an existing tomographic PTV algorithm, which is reported in Schneiders et al. (2016b), is used. The tomographic PTV algorithm is based on the MART algorithm (Elsinga et al. 2006), which yields the particle intensity within a voxel discretized volume. Individual particles are identified using a peak-finding algorithm based on a kernel of 5 × 5 × 5 voxels. Sub-pixel position of the particle centroid is obtained using a 3-point Gaussian fit (Raffel et al. 2007). Particle trajectories are built using the multi-step algorithm by Malik et al. (1993), where in case multiple particle trajectories are found within a search window, a criterion is applied that selects the trajectory with minimum change in acceleration. The discrete positions of the particle along its trajectory are fitted with a polynomial function, which regularizes the particle position, velocity and acceleration (Novara and Scarano 2013). The analytical time-derivatives of the polynomial yield from each particle trajectory, γ, at time-instant t _i and location x _γ (t _i ) a measurement of velocity, u _m , and a measurement of the velocity material derivative, D u _m /Dt.

2.2 Dense velocity reconstruction using VIC+

The VIC+ method uses the velocity–vorticity formulation of the incompressible and inviscid Navier–Stokes equations as outlined in Sect. 2.2.1. The actual VIC+ optimization procedure is subsequently discussed in Sect. 2.2.2.

2.2.1 Velocity–vorticity formulation

Consider a computational grid with mesh spacing h, where h is set as discussed in Sect. 2.3.1. Vorticity, ω _h , is defined on the grid using radial basis functions (Sect. 2.3.1) and related to velocity, u _h , through a Poisson equation (Christiansen 1973; Schneiders et al. 2014),

$$\nabla^{2} \varvec{u}_{h} = - \nabla \times\varvec{\omega}_{h} ,$$

(1)

with boundary conditions on velocity. Second-order central differences are used for evaluation of all the spatial derivatives in the present work, except at the volume boundaries where first-order single-sided differences are used. The velocity material derivative, D u _h /Dt, is defined on the grid by,

$$\frac{{{\text{D}}\varvec{u}_{h} }}{{{\text{D}}t}} = \frac{{\partial \varvec{u}_{h} }}{\partial t} + \left( {\varvec{u}_{h} \cdot \nabla } \right)\varvec{u}_{h} .$$

(2)

The temporal velocity derivative in Eq. (2) is calculated from the temporal vorticity derivative by solution of a Poisson equation (Schneiders et al. 2016a, b),

$$\nabla^{2} \frac{{\partial \varvec{u}_{h} }}{\partial t} = - \nabla \times \frac{{\partial\varvec{\omega}_{h} }}{\partial t} ,$$

(3)

with boundary conditions on the temporal velocity derivative. The temporal vorticity derivative in turn is calculated from the vorticity transport equation,

$$\frac{{\partial\varvec{\omega}_{h} }}{\partial t} = \left( {\varvec{\omega}_{h} \cdot \nabla } \right)\varvec{u}_{h} - \left( {\varvec{u}_{h} \cdot \nabla } \right)\varvec{\omega}_{h} .$$

(4)

2.2.2 Optimization procedure

The velocity and velocity material derivative fields can be calculated at a single time-instant from Eqs. (1) and (2) if the radial basis coefficient weights for vorticity and boundary conditions on velocity, and its temporal derivative is known on the computational grid. Depending on the flow case, one or more boundary conditions could be known a priori from a no-slip condition or a free stream condition. The remaining vorticity values on the grid and boundary conditions are, however, unknown. Consider a vector, $\varvec{\xi}$, that collects the unknown radial basis function weights for vorticity and boundary conditions on velocity and its temporal derivative,

$$\varvec{\xi}= \left\{ {\varvec{\xi}_{\varvec{\omega}} ,\beta_{1}\varvec{\xi}_{{BC,\varvec{u}}} ,\beta_{2}\varvec{\xi}_{{BC,\partial \varvec{u}}} } \right\} ,$$

(5)

where coefficients β ₁ and β ₂ are scaling parameters (Sect. 2.3.2). The vector $\varvec{\xi}$ contains the optimization variables for the VIC+ optimization procedure.

The goal of the VIC+ optimization procedure is to find the optimization variables, $\varvec{\xi}$, that minimize a cost function, J, under the constraint that the motion equations in velocity–vorticity form (Sect. 2.2.1) are respected at a single measurement time-instant. The cost function is defined by Eq. (6),

$$J = J_{\varvec{u}} + \alpha^{2} J_{{D\varvec{u}}} ,$$

(6)

where α is a weighting coefficient (Sect. 2.3.3), J _u is given by Eq. (7) and J _{D u} is given by Eq. (8),

$$J_{\varvec{u}} = \sum\nolimits_{p} {\left\| {\varvec{u}_{h} \left( {\varvec{x}_{p} } \right) - \varvec{u}_{m} \left( {\varvec{x}_{p} } \right)} \right\|}^{2} ,$$

(7)

$$J_{{D\varvec{u}}} = \sum\nolimits_{p} {\left\| {\frac{{D\varvec{u}_{h} }}{Dt}\left( {\varvec{x}_{p} } \right) - \frac{{D\varvec{u}_{m} }}{Dt}\left( {\varvec{x}_{p} } \right)} \right\|}^{2} ,$$

(8)

where u _h and D u _h /Dt are calculated from Eqs. (1) and (2) and are evaluated at the particle locations, x _p , by linear interpolation from the computational grid. The cost function penalizes the difference between the PTV measurements and the velocity and material derivative at a single measurement time-instant calculated from the optimization variables. The optimization problem does not include time-integration of the vorticity transport equation. Time-integration is avoided because the cost function for the optimization (Eq. 6) only considers the instantaneous velocity and its material derivative. The optimization problem is solved iteratively using the limited-memory Broyden–Fletcher–Goldfarb–Shanno method (L-BFGS, Liu and Nocedal 1989), which is a commonly used method for large-scale optimization problems. The L-BFGS method is a quasi-Newton method technique that approximates the standard BFGS method using a limited amount of memory. It approximates the typically dense Hessian matrix to reduce the cost of storing and manipulating it. For details on the Hessian approximation and the step size determination, the reader is referred to Liu and Nocedal (1989). At each iteration, the method requires the gradient of the cost function with respect to the optimization variables, $\varvec{ }\partial J/\partial\varvec{\xi}$. This gradient could be approximated using finite differences, requiring N evaluations of the cost function, where N is the number of optimization variables. Following recent literature (Gronskis et al. 2013; Yang et al. 2015; Yegavian et al. 2015; Lemke and Sesterhenn 2016; among others), an alternative and more cost-efficient method to calculate the gradient is provided by the adjoint approach. This approach gives the exact gradient, but more importantly the computational cost of the adjoint approach is approximately equal to only one evaluation of the cost function. Therefore, within VIC+, the adjoint approach is used for evaluation of the gradient $\partial J/\partial\varvec{\xi}$. The steps to calculate this gradient using the adjoint approach are given in Appendix A. It should be remarked that the adjoint procedure within VIC+ is relatively simple in comparison with the adjoint procedure in the aforementioned literature, because many calculation steps within VIC+ are linear operations and because the procedure does not involve time-integration. Only the state at a single time-instant needs to be kept in memory for evaluation of the adjoint equations and therefore the method does not suffer from the typically large memory requirements for adjoint-based optimization techniques, which led to the proposal of for example a storage/recomputation strategy by Gronskis et al. (2013).

Appendix B contains a pseudo-code of the VIC+ method and Fig. 2 schematically summarizes the method. The PTV measurements (lines 5–7 in the pseudo-code; green box in Fig. 2) are used to obtain an initial estimate of the velocity field on the grid (line 10; gray box in Fig. 2), using for example the adaptive Gaussian windowing technique (AGW, Agüí and Jiménez 1987). This initial condition provides the initial estimate of the degrees of freedom vector, $\varvec{\xi}$ (lines 15–18). The optimization procedure is illustrated by the blue loop in Fig. 2 and spans lines 22–69 in the pseudo-code. It starts by calculation of velocity and its material derivative on the computational grid from the optimization variables (lines 24–29) and subsequently the interpolation of these quantities to the tracer particle locations (lines 32–33). The cost function (Eq. 6) is calculated at lines 36–38, after which the gradient of the cost function with respect to the degrees of freedom is calculated using the adjoint of the code (lines 41–53). The degrees of freedom are then updated using the step size and direction found by the L-BFGS algorithm (lines 59–61). When the system is converged (Sect. 2.3.2; lines 64–66 in the pseudo-code), the optimization loop is stopped and the output velocity field is calculated (lines 72–73).

Fig. 2

VIC+ reconstruction framework. From the PTV input (left, green box), an initial velocity estimate is made (gray box), which is input into the VIC+ iterative procedure (blue loop) to find the optimization variables that yield a velocity and velocity material derivative of minimum discrepancy with the PTV measurements (orange box)

The procedure is christened as ‘VIC+’. It is inspired by time-supersampling using VIC simulation (Schneiders et al. 2014), and the ‘plus’ symbol refers to the simultaneous use of the particle velocity and its time-derivative.

2.3 Problem closure

The VIC+ method introduces multiple parameters whose values need to be determined before performing the optimization. In recent literature on methods with similar weighting parameters, optimal setting is either left out for future research (e.g., Papadakis and Mémin 2008) or manually chosen on a case-by-case basis (Gronskis et al. 2013; Gregson et al. 2014; among others). In an optimization framework similar to VIC+, the ‘FlowFit’ method for interpolation of scattered data to an Eulerian grid (Gesemann et al. 2016) makes use of weighting parameters to penalize differences between the reconstruction and the velocity and material derivative measurements, and other differences such as nonzero curvature and divergence depending on the flow properties (e.g., incompressibility and measurement uncertainty). However, no direct relations are provided to obtain the values of these weights. To achieve a general treatment of experimental data and avoid case-dependent, or even user-dependent parameter selection by tuning, the present study includes a fixed approach to closure for the mentioned parameters. These parameters are discussed separately in the sections below.

2.3.1 Mesh spacing and radial basis functions

The VIC+ method uses finite differences for evaluation of derivatives on the computational grid, requiring a grid-spacing that makes truncation errors negligible. The grid node spacing, h, is set based on the particle concentration, C, using Eq. (9)

$$h = \frac{1}{4}C^{ - 1/3}$$

(9)

This equation implies that 64 grid nodes are employed per tracer particle. For tomographic PIV, the amount of grid nodes typically also exceeds the amount of tracer particles, as the window overlap factor is often set to 75 % to avoid unnecessary truncation when evaluating spatial derivatives (Scarano and Poelma 2009; Togkoz et al. 2012). Also, considering an optimization technique similar to VIC+, the recent study of Schanz et al. (2016) reports that in the ‘FlowFit’ interpolation method typically 5–20 B-spline cells are chosen per tracer particle. For the VIC+ method, at each grid node Gaussian radial basis functions, ϕ(r) are employed. The free parameter, σ, which defines the width of the Gaussian radial basis functions, is chosen slightly larger than the grid node spacing (σ = 1.1 h), and the resulting radial basis functions are defined by Eq. (10)

$$\phi \left( r \right) = e^{{ - \frac{{r^{2} }}{{2\sigma^{2} }}}} = e^{{ - \frac{{r^{2} }}{{2.4h^{2} }}}} .$$

(10)

Vorticity is calculated on the grid through Eq. (11),

$$\varvec{\omega}_{h} = \sum w_{i} \phi \left( {\left\| {\varvec{x} - \varvec{x}_{i} } \right\|} \right) ,$$

(11)

where w _i are the radial basis function weights. The boundary conditions of velocity and its temporal derivative are defined analogously. On average, 0.7 particles are present within a sphere with a radius of 2σ that is centered around a grid node. Because each particle provides measurements of both velocity and material derivative, each radial basis function includes on average 1.4 measurements (of three components) within such a sphere. In addition, the system is constrained by the motion equations in velocity–vorticity formulation (Sect. 2.2.1). As discussed in Sect. 2.2.2, the cost function of the nonlinear system of equation is iteratively reduced using a gradient-based quasi-Newton method (the L-BFGS method), where the optimization variables are the radial basis function weights for vorticity and the boundary conditions on velocity and the velocity temporal derivative.

2.3.2 Scaling of the optimization variables

The optimization variables ($\varvec{\xi}_{\omega }$, $\varvec{\xi}_{BC,u}$ and $\varvec{\xi}_{BC,\partial u}$) have different units, requiring weighting or non-dimensionalization of the optimization variables. A non-dimensional vector with the optimization variables $\varvec{\xi}^{\varvec{*}} ,$ can be defined by non-dimensionalization of the optimization variables with their standard deviation,

$$\varvec{\xi}^{*} = \left\{ {\frac{1}{{\sigma_{\omega } }}\varvec{\xi}_{\omega } ,\frac{1}{{\sigma_{BC,u} }}\varvec{\xi}_{BC,u} ,\frac{1}{{\sigma_{BC,\partial u} }}\varvec{\xi}_{BC,\partial u} } \right\} ,$$

(12)

where σ _ω , σ _BC,u and σ _BC,∂u are the standard deviation of $\varvec{\xi}_{\omega }$, $\varvec{\xi}_{BC,u}$ and $\varvec{\xi}_{BC,\partial u}$ as defined by Eqs. (13)–(15),

$$\sigma_{\omega } = \sqrt {\frac{1}{{n_{\omega } }}\sum \left( {\varvec{\xi}_{\varvec{\omega}}^{{}} - \bar{\varvec{\xi }}_{\varvec{\omega}} } \right)^{2} } ,$$

(13)

$$\sigma_{BC,u} = \sqrt {\frac{1}{{n_{BC,u} }}\sum \left( {\varvec{\xi}_{BC,u}^{{}} - \bar{\varvec{\xi }}_{BC,u} } \right)^{2} } ,$$

(14)

$$\sigma_{BC,\partial u} = \sqrt {\frac{1}{{n_{BC,\partial u} }}\sum \left( {\varvec{\xi}_{BC,\partial u}^{{}} - \bar{\varvec{\xi }}_{BC,\partial u} } \right)^{2} } ,$$

(15)

with n _ω the number of grids points on which vorticity is an optimization variable and n _BC,u and n _BC,∂u the number of boundary points on which velocity and its temporal derivative are optimization variables.

Multiplication of $\varvec{\xi}^{*}$ by σ _ω gives the vector of optimization variables, $\varvec{\xi}$, used for VIC+,

$$\varvec{\xi}= \sigma_{\omega }\varvec{\xi}^{*} = \left\{ {\varvec{\xi}_{\omega } ,\frac{{\sigma_{\omega } }}{{\sigma_{BC,u} }}\varvec{\xi}_{BC,u} ,\frac{{\sigma_{\omega } }}{{\sigma_{BC,\partial u} }}\varvec{\xi}_{BC,\partial u} } \right\} = \left\{ {\varvec{\xi}_{\omega } ,\beta_{1}\varvec{\xi}_{BC,u} ,\beta_{2}\varvec{\xi}_{BC,\partial u} } \right\} ,$$

(16)

where the weighting coefficients β ₁ and β ₂ are defined by Eqs. (17) and (18).

$$\beta_{1} = \frac{{\sigma_{\omega } }}{{\sigma_{BC,u} }} ,$$

(17)

$$\beta_{2} = \frac{{\sigma_{\omega } }}{{\sigma_{BC,\partial u} }}$$

(18)

2.3.3 Weighting coefficient α

A non-dimensional cost function, J ^*, can be defined by non-dimensionalization of J _u and J _Du by the variance of the PTV velocity and velocity material derivative measurements,

$$J^{*} = \frac{1}{{\sigma_{u}^{2} }}J_{u} + \frac{1}{{\sigma_{Du}^{2} }}J_{Du} ,$$

(19)

where σ _u and σ _Du are the standard deviation of the PTV velocity and velocity material derivative measurements,

$$\sigma_{u} = \sqrt {\frac{1}{{n_{PTV} }}\sum\nolimits_{p = 1}^{{n_{PTV} }} {\left\| {\varvec{u}_{p} - \bar{\varvec{u}}_{p} } \right\|}^{2} } ,$$

(20)

$$\sigma_{Du} = \sqrt {\frac{1}{{n_{PTV} }}\sum\nolimits_{p = 1}^{{n_{PTV} }} {\left\| {\frac{{D\varvec{u}_{p} }}{Dt} - \overline{{\frac{{D\varvec{u}_{p} }}{Dt}}} } \right\|}^{2} } ,$$

(21)

with n _PTV the number of tracer particles (i.e., the number of PTV measurements). Multiplication of the non-dimensional cost function J ^* by σ ² _u gives the cost function used in the present study,

$$J = \sigma_{u}^{2} J^{*} = J_{u} + \frac{{\sigma_{u}^{2} }}{{\sigma_{Du}^{2} }}J_{Du} = J_{u} + \alpha^{2} J_{Du} ,$$

(22)

where the weighting coefficient α is defined by,

$$\alpha = \sigma_{u} /\sigma_{Du}$$

(23)

2.3.4 Convergence criterion

A convergence criterion similar to, among others, Lemke and Sesterhenn (2016) is used, where an iterative procedure is stopped (converged) when the change in the cost function J becomes small. Here, the convergence level is selected at 10⁻⁴ % of the initial value of J.

3 Simulated turbulent boundary layer

3.1 Test case and data processing

The numerical assessment considers the time-resolved tomographic PIV experiment of a flat-plate turbulent boundary layer (u _∞  = 10 m/s, δ ₉₉  = 9.4 mm), simulated from a direct numerical simulation (DNS) by Bernardini and Pirozzoli (2011). The Reynolds number based on boundary layer thickness is Re _δ  = 8185. This numerical simulation was used previously for assessment of tomographic PIV processing techniques by Pröbsting et al. (2013) and Lynch and Scarano (2015). The salient boundary layer parameters are summarized in Table 1.

Table 1 Boundary layer parameters for the numerical assessment

The configuration of the simulated experiment is chosen similar to the existing real-world time-resolved PIV experiment by Pröbsting et al. (2013), whose data are also used in the present study for assessment. The experimental setup is schematically illustrated in Fig. 7. The simulated measurement volume has streamwise length L _x /δ ₉₉  = 2.1 and spanwise length L _z /δ ₉₉  = 4 corresponding to the real-world experiment, and the wall-normal length is chosen to capture the full turbulent boundary layer height (L _y /δ ₉₉  = 1.1). The procedure to generate the simulated experiment is chosen similar to the approach used by Novara and Scarano (2013). Tracer particles are randomly and uniformly distributed in space. Their motion is obtained by integrating the velocity over time using a fourth-order Runge–Kutta method. With the same voxel/mm ratio of 0.042 as for the existing experiment by Pröbsting et al. (2013), intensity volumes of 470 × 250 × 890 voxels are generated. A reconstructed particle diameter of two voxels is assumed, and the exact particle positions in the objects are purposely corrupted with normally distributed random noise with a standard deviation of 0.2 voxels to model a remaining registration error (Novara and Scarano 2013).

Objects are generated with a time-separation of 0.1 ms, corresponding to a measurement rate f _acq  = 10 kHz, which is equal to the existing real-world experiment by Pröbsting et al. (2013). The free stream particle displacement at this acquisition frequency is 20 vox. Similar to Novara and Scarano (2013), the approach followed here solely simulates the motion of particle tracers and tomographic imaging is not included in the simulation. Therefore, also the artifacts of tomographic reconstruction such as ghost particles (Elsinga et al. 2006) are not reproduced. This approach allows a study of the effect of tracer particle concentration on the potential velocity reconstruction quality, separately from that of the tomographic reconstruction quality. A range of seeding concentrations from C = 50 particles/δ ³₉₉ to C = 2000 particles/δ ³₉₉ is considered. Figure 3 shows extracts from the projected volumes from the intensity volumes at the seeding concentrations considered. For the lowest particle concentrations (C = 50 and 200 particles/δ ³₉₉ ), the particle image density is below N _p  = 0.005 ppp, which is a typical level for 3D-PTV. At C = 2000 particles/δ ³₉₉ , the particle image density equals N _p  = 0.045 ppp, corresponding to an optimal working condition for tomographic reconstruction (Elsinga et al. 2006) and advanced triangulation methods such as iterative particle reconstruction (IPR, Wieneke 2013).

Fig. 3

Projected volumes of the test cases at different particle concentration. The seeding density, C, for each strip is indicated in particles/δ ³₉₉ and the particle image density, N _p , in particles per pixel (ppp). The colormap is clipped at I = 2·10⁵ for clarity, and larger values of intensity are colored red

The objects are processed using the proposed technique, VIC+, and for comparison also using PIV (spatial cross-correlation), linear interpolation, adaptive Gaussian windowing (AGW, Agüí and Jiménez 1987) and VIC with no account of D u /Dt. The latter is deemed important to establish the specific contribution of the additional material derivative term in improving the velocity reconstruction. For tomographic PIV analysis, an iterative cross-correlation technique with multi-grid and window-deformation (Scarano and Riethmuller 2000) is used with four iterations and an interrogation volume size such that ten particles are included in the measurement box. A window overlap of 75 % is used. Following Pröbsting et al. (2013), the aspect ratio of the interrogation volumes is chosen as 2 × 1 × 2, with the smallest dimension along the wall-normal to account for the velocity gradient in the boundary layer. After each iteration, the vector fields are subject to universal outlier detection to remove spurious vectors (Westerweel and Scarano 2005). For the intermediate iterations, the vector fields were filtered using a second-order polynomial filter. No filter was applied to the final PIV velocity fields.

The other techniques require input from PTV analysis. To this extent, 7-snapshot tracks are built, through which a second-order polynomial is fitted to regularize the particle positions (Novara and Scarano 2013). The analytical time-derivatives of the particle tracks yield the velocity and velocity material derivative measurements at the time-centered position of each track. For both VIC without D u /Dt and VIC+, the no-slip condition is prescribed at the wall. On the other boundaries, the time-averaged velocity is prescribed, which can be obtained from cross-correlation analysis or more efficiently from time-averaging of the result of linear interpolation of the PTV results to a regular grid. To prevent this boundary condition from affecting the instantaneous velocity, the computational volume is padded by 30 % at all boundaries, except into the wall, similar to the padding procedure used in Schneiders et al. (2016a). The padded domain is not considered for further data analysis. For VIC+, the VIC without D u /Dt result is used as initial condition, for which in turn the AGW result is used as initial condition. Furthermore, the weighting coefficient α is calculated using Eq. (17), which gives a value of 4·10⁻⁴ for α. In the present study, computational time was not considered as a limiting factor. However, in future studies, a potential improvement of the initialization approach can be envisaged. Inspired by SMTE (Lynch and Scarano 2015), starting from the second time-step in a series the VIC+ method could potentially be initialized efficiently using the result from the previous time-instant advected forward in time over the time corresponding to the laser pulse separation.

3.2 Assessment of the results

First, the results are inspected by visual comparison of the high- and low-speed streaks and the vortical structures in the turbulent boundary layer. Isosurfaces of low- and high-speed streaks are plotted in blue (u′/u _∞ = −0.07) and yellow (u′/u _∞ = +0.07) for a single time-instant in Fig. 4. The vorticity magnitude isosurfaces at the same time-instant are plotted in Fig. 5. Figures 4a and 5a show the DNS result. The expected streamwise elongated low- and high-speed streaks in the turbulent boundary are visible from the blue and yellow isosurfaces. To aid the visual comparison, a sub-volume in the center of the volume is selected in Fig. 4. In Fig. 5, the pair of hairpin vortices on the low-speed streak in the volume center is extracted for visual comparison of the methods. Figures 4b and 5b show the organization of tracer particles in the sub-volumes at the seeding concentrations considered. From left to right, the seeding concentration ranges from 50 to 2000 particles/δ ³₉₉ . For clarity, the particles are colored by streamwise velocity. Figures 4c–g and 5c–g show the corresponding reconstructions using PIV (c), linear interpolation (d), AGW (e), VIC without D u /Dt (f) and VIC+ (g).

Fig. 4

Instantaneous reference flow field in the measurement volume (a), visualized by isosurfaces of u′/u _∞ = −0.07 (blue) and u′/u _∞ = + 0.07 (yellow). Instantaneous particle distributions at seeding a concentration of (from left to right) 50, 200, 500, 1000 and 2000 particles/δ ³₉₉ (b, particles colored by velocity magnitude). Flow field reconstruction using PIV (c), linear interpolation (d), AGW (e), VIC without D u /Dt (f), and VIC+ (g)

Fig. 5

Instantaneous reference flow field in the measurement volume (a), visualized by isosurfaces of u′/u _∞ = −0.07 (blue) and vorticity magnitude (|ω| = 2500 1/s, red). Instantaneous particle distributions at seeding a concentration of (from left to right) 50, 200, 500, 1000 and 2000 particles/δ ³₉₉ (b, particles colored by velocity magnitude). Flow field reconstruction using PIV (c), linear interpolation (d), AGW (e), VIC without D u /Dt (f), and VIC+ (g)

The PIV analysis (Figs. 4c, 5c) yields a spatially low-pass filtered version of the reference velocity field. The filtering is particularly pronounced at the lower seeding concentrations, whereas, for the concentrations above 1000 particles/δ ³₉₉ , the velocity field recovers the main high- and low-speed streaks. The pair of hairpin vortices (Fig. 5) is combined in one hairpin vortical structure at the two highest concentrations, indicating the effect of spatial filtering associated with the cross-correlation approach in this case. The linear interpolation result (Figs. 4d and 5d) does not suffer from the filtering effect of the cross-correlation approach. At the lowest concentration, it shows the low-speed streaks, which were not visible from the PIV analysis. At the highest seeding concentration, the velocity streaks reconstructed by the linear interpolation are visually in good correspondence to the DNS reference. The results in Fig. 5c indicate the presence of vortices starting from C = 500 particles/δ ³₉₉ . At C = 2000 particles/δ ³₉₉ , the pair of hairpin vortices just becomes visible. The reconstruction, however, exhibits a visible noise component with random fluctuations degrading the spatial coherence of the velocity and vorticity field. The AGW result (Figs. 4e and 5e) is essentially similar to linear interpolation, with some better clarity in the detection of the flow structures due to the increased spatial coherence introduced by the weighting function.

At seeding concentrations above 1000 particles/δ ³₉₉ , both the VIC without D u /Dt (Fig. 4f) and VIC+ (Fig. 4g) methods are not expected to show a visual improvement for the low- and high-speed streaks as the other PTV-based methods already capture the relevant low- and high-speed streaks in the volume. This is confirmed by the results in Fig. 4. Due to the choice of vorticity as basis for the VIC reconstructions, both VIC with and without D u /Dt are expected to reconstruct a smoother vorticity field than the reconstructions resulting from linear interpolation and AGW. Indeed, smoothness and coherence of the hairpin vortices are improved with the VIC methods in Fig. 5. For VIC without D u /Dt, the pair of hairpin vortices can be identified at concentrations above 1000 particles/δ ³₉₉ . When the D u /Dt term is used for the reconstruction with VIC+ (Fig. 5g), at half the seeding concentration (500 particles/δ ³₉₉ ) the hairpins can be identified.

In the introduction, it was noted that in recent volumetric particle velocimetry experiments a damping of the turbulent velocity fluctuations is found. To illustrate how the application of the velocity reconstruction techniques impacts estimation of the turbulent velocity fluctuations, velocity statistics are plotted in Fig. 6 for the range of seeding concentrations considered. This figure shows from left to right: u _mean , u′ _rms , v′ _rms , w′ _rms and (u′v′) _rms . From top to bottom, the seeding concentration increases from 50 particles/δ ³₉₉ to 2000 particles/δ ³₉₉ . The statistics were obtained using 100 snapshots (measurement duration of 10 ms), and the data ensemble along x and z directions were also used for averaging. As can be seen in Fig. 6, all techniques are able to capture already at the lowest seeding concentration the mean streamwise velocity profile above y/δ ₉₉ > 0.2. Closer to the wall, both PIV and AGW yield an overestimation of the time-averaged streamwise velocity due to the effect of the interrogation volumes and Gaussian windows not being able to capture the strong velocity gradient near the wall. As expected from the qualitative comparison in Fig. 4, the differences between PIV- and PTV-based techniques are further pronounced in the u′ _rms statistics. Differences between the PTV-based techniques remain small, however, for u′ _rms , as these fluctuations are associated with low- and high-speed streaks occurring as rather elongated structures, which correspond to relatively small wavenumbers. In contrast, the wall-normal velocity fluctuations corresponding to smaller flow structures such as hairpins and hairpin packets (e.g., Adrian et al. 2000) are expected to be poorly captured at low concentration of tracers, as their wavenumbers are filtered by the different reconstruction techniques. At the highest seeding concentration, the PIV result yields most significant damping at y/δ ₉₉ = 0.3 (~35 %). Linear interpolation reduces this to 20 %. The adaptive Gaussian windowing technique performs similar to VIC without the aid of D u /Dt at the highest seeding concentration case with 10 % damping. When the material derivative is used for the reconstruction with VIC+, damping of v′ _rms and w′ _rms is 4 %.

Fig. 6

Turbulence statistics calculated from the different velocity field reconstructions; from top to bottom, the tracer particle concentration, C, equals C = 50, 200, 500, 1000 and 2000 particles/δ ³₉₉

The VIC+ result remains within 5 % at y/δ ₉₉ = 0.3, also when the particle seeding concentration is reduced by an order of magnitude to 200 particles/δ ³₉₉ . At this seeding concentration, the differences between VIC+, VIC without the use of D u /Dt (10 % damping) and AGW (20 % damping) are more pronounced, as expected also from the study of the hairpin vortex reconstruction above. The cross term (u′v′) _rms requires most tracer particles for correct reconstruction, as indicated by the more significant damping at especially the low seeding concentration cases. The trend of the results remains, however, and, as also indicated by the isosurfaces of instantaneous u′ and vorticity in Figs. 4 and 5, at each seeding concentration case PIV yields significantly damped statistics and VIC+ the most accurate reconstruction. For the present study, both VIC+ and cross-correlation for tomographic PIV were implemented in Mathworks MATLAB 2015b. The computation time for VIC+ was found to be approximately four times longer than for cross-correlation.

4 Experimental assessment

The experimental assessment considers two tomographic velocimetry experiments. First, the findings of the numerical assessment are validated by an existing real-world experiment of a turbulent boundary layer (Sect. 4.1). Second, the method is applied to measurements of a circular jet flow (Sect. 4.2).

4.1 Turbulent boundary layer

To complement the numerical assessment on the turbulent boundary layer, the real-world tomographic PIV experiment of Pröbsting et al. (2013) in a similar turbulent boundary layer is considered here. The experiment considers a turbulent boundary layer over a flat plate at a free stream velocity of 10 m/s. The Reynolds number based on the local boundary layer thickness (δ ₉₉  = 9.4 mm) is Re _δ  = 6240. The measurements are performed with four LaVision HighSpeedStar CMOS cameras equipped with Nikon Micro-Nikkor 105 mm prime lenses and a Quantronix Darwin-Duo Nd:YLF laser. The cameras and laser are operated in single-frame mode at 10 kHz for a time-resolved recording. A multi-pass light amplification system is installed, following the indications of Ghaemi and Scarano (2010), to increase the illumination intensity. Salient details of the flow configuration and the experimental setup are given in Tables 2 and 3, respectively. Sketches of the setup are given in Fig. 7. For a more complete discussion of the experimental setup, the reader is referred to Pröbsting et al. (2013).

Table 2 Boundary layer parameters for the experimental assessment

Table 3 Tomographic PIV setup

Fig. 7

Schematic of the tomographic PIV experiment (schematics not to scale); top view (left) and back view (right); the measurement volume is indicated by the black dashed line in details A and B

The original tomographic PIV measurements by Pröbsting et al. (2013) were processed using a volume deformation iterative multi-grid technique with a final interrogation volume size of 32 × 16 × 32 vox at 75 % overlap. This resulted in a vector spacing of 0.16 mm along the wall-normal and 0.33 mm in the stream and spanwise dimensions. For increased correlation signal, a short-time sliding-average correlation technique with a kernel size k = 4 (Sciacchitano et al. 2012) was used. The no-slip condition is imposed at the wall to stabilize the interrogation and reduce the number of spurious vectors (Theunissen et al. 2008).

For the PTV-based techniques, the experiment is re-processed using tomographic PTV. A three-step particle tracking algorithm based on Malik et al. (1993) is applied to the tomographic reconstructions to build particle tracks over five snapshots. Particles are identified by peak-finding in a 5 × 5 × 5 voxel neighborhood and sub-pixel accuracy of particle location is obtained by fitting of a Gaussian through their intensity distribution. An adaptive search radius based on the estimated rms velocity fluctuations from PIV is used. Approximately 3200 tracks are found per measurement volume, corresponding to C = 2000 particles/δ ³₉₉ . To regularize the particle trajectories, following Novara et al. (2010), a second-order polynomial is fitted through the particle locations. Velocity and velocity material derivative are evaluated at the center time-instants of each polynomial fit by evaluating analytically the first and second temporal derivatives. No trajectories could be evaluated below y/δ ₉₉ = 0.06 due to reduced laser light intensity close to the wall, as visible in the z-intensity profile of the MART reconstruction (Fig. 8, left). This figure also shows the expected peak intensity value at the wall (y/δ ₉₉ = 0) due to the wall reflections. The normalized histogram of identified particle trajectories along z-depth is plotted in Fig. 8, right. The same settings as in the numerical assessment are used for all PTV-based reconstruction techniques.

Fig. 8

Profile of y-intensity of the tomographic reconstructions (left) and normalized histogram of identified particle trajectories along the wall-normal direction (right)

Analogous to the numerical assessment, first the results are compared qualitatively. Isosurfaces of instantaneous velocity magnitude and vorticity magnitude are plotted in Fig. 9, left. In the same figure, the top-views on the right show the instantaneous isosurfaces at two consecutive time-instants to allow for inspection of temporal coherence of the structures. From top to bottom, the result from PIV, linear interpolation, AGW, VIC without D u /Dt and VIC+ with D u /Dt reconstructions are plotted. Analogous to the results in the numerical assessment, all methods appear to capture the high- and low-speed streaks of streamwise velocity, which correspond to relatively small wavenumbers. For closer inspection of the structures corresponding to larger wavenumbers, a single vortical structure ‘A’ has been highlighted in Fig. 9. For the tomographic PIV result (Fig. 9a), the structure is marginally visible at times t = 0 ms and t = 0.4 ms, where it has convected downstream. As expected from the numerical assessment, the PTV-based techniques show an increase in vortical structures. The AGW result (Fig. 9c) shows most vorticity isosurfaces, but structure ‘A’ cannot be identified in the two time-instants plotted, indicating the result is affected by a significant amount of measurement noise. Both VIC without D u /Dt (Fig. 9d) and VIC+ (Fig. 9e) show a regularized result. In particular, the VIC + result shows the elongated structure ‘A’ at both time-instants t = 0 ms and t = 0.4 ms. This confirms in this real-world experiment the increased coherence of vortical structures with VIC+, that was also observed in the numerical assessment.

Fig. 9

Isosurfaces of velocity magnitude (|u| = 7 m/s, yellow) and vorticity magnitude (|ω| = 2500 1/s, blue). Perspective view at t = 0 ms (left) and top-view at t = 0 (middle) and t = 0.4 ms (right). From top to bottom: PIV (a), Linear Interpolation (b), AGW (c), VIC without D u /Dt (d) and VIC+ (e)

Parallel to the analysis in the numerical assessment, the turbulent velocity fluctuations are considered. The velocity statistics are plotted in Fig. 10. For the statistics, the measurement duration was 200 ms and 200 snapshots were used at a time-separation of 1 ms. The statistical evaluation of the turbulent flow properties is obtained considering the data ensemble along the homogeneous directions x and y. The near-wall region where no particle trajectories could be identified is masked gray in this plot. As expected from the numerical assessment, the mean velocity profile (Fig. 10, left) is captured by tomographic PIV in good correspondence (within 4 %) to the DNS simulation by Bernardini and Pirozzoli (2011). The statistics of u ^′_rms associated with to low- and high-speed streaks occurring as rather elongated structures and corresponding to relatively small wavenumbers are captured with damping up to 5 % in comparison with the DNS. On the other hand, statistics of v ^′_rms and w ^′_rms correspond to larger wavenumbers and are damped by approximately 20 % by tomographic PIV.

Fig. 10

Velocity statistics in the turbulent boundary as approximated using different processing techniques. No particle tracks were identified in the gray-shaded area

The time-averaged velocity outside the masked region is also for the PTV-based techniques in good correspondence with the DNS reference. Similarly, the streamwise rms velocity fluctuations of all methods are within approximately 5 % of the DNS value outside of the masked region. Inside of the masked region, no particle trajectories were found with tomographic PTV and the peak value of u ^′_rms in the masked region is not found by any of the PTV-based techniques. On the other hand, the PIV analysis retrieves the expected peak in u ^′_rms at y/δ ₉₉ = 0.03. The PIV result is based on the short-time average cross-correlation (sliding-average correlation, SAC, Sciacchitano et al. 2012) based on four frames, which adds the needed robustness to the analysis in this region of low light intensity. Instead, in the same region, although individual particles are occasionally detected, their systematic tracing is not possible due to the low peak intensity. As a result, no particle trajectory is identified in the near-wall region and consequently the PTV-based techniques do not capture the peak value of u ^′_rms in the masked region. The differences between the techniques, outside of the masked region, become more marked when the wall-normal and spanwise velocity fluctuations are considered. Linear interpolation shows slightly increased damping of v ^′_rms and w ^′_rms with respect to PIV. The AGW result closely follows the PIV result above y/δ ₉₉ = 0.15. Closer to the wall, AGW shows a significant increase in turbulent velocity fluctuations, which is ascribed to measurement noise, based on the discussion above of the isosurfaces (Fig. 9c). The isosurface visualization in Fig. 9 of the results using VIC with and without D u /Dt showed increased coherence of the vortical structures with respect to the other techniques, indicating a potential increase in rms velocity fluctuations is not due to increased noise. The statistics of v ^′_rms and w ^′_rms are with VIC without D u /Dt slightly higher than the tomographic PIV result. When D u /Dt is leveraged with VIC+, the rms fluctuations of wall-normal and span-wise velocity components are found within 5 % of the DNS reference, in comparison with 20 % damping by PIV. Inside of the masked region, both VIC results follow a decreasing trend of the velocity statistics. Because there are no particle trajectories available in this region, this trend is solely imposed by the no-slip condition used for VIC. Therefore, the relatively good comparison with DNS is artificial in this region. The results outside of the masked region validate the findings of the numerical assessment on a real-world tomographic PIV experiment and show that the velocity material derivative can be used using the VIC + technique to improve reconstruction quality.

It should be remarked that the present experiment was taken at a particle image density of N _p  = 0.07 ppp (Pröbsting et al. 2013) and thus at the limit of the capabilities of tomographic reconstruction (Scarano 2013). The measurement captured 45 % of the boundary layer thickness. Increasing volume thickness at the same seeding concentration would increase the particle image density beyond the limit for tomographic reconstruction (i.e., doubling volume thickness yields a particle image density of 0.14 ppp). This highlights the relevance of the VIC+ reconstruction technique for cases where the spatial resolution is limited by the tracer particle seeding concentration.

4.2 Jet flow

The second experimental case considers data from an available time-resolved tomographic PIV measurement of a circular jet in water by Violato and Scarano (2011). For full details of the experiment, the reader is referred to the latter manuscript and here only the salient details are repeated. The nozzle exit diameter is D = 10 mm, and the contraction ratio equals 56:1. The exit velocity of the jet in water is 0.5 m/s, corresponding to a Reynolds number based on the nozzle exit diameter of Re _D  = 5000. Neutrally buoyant polyamide particles of 56 μm diameter are used and dispersed homogeneously in the water tank. The particle seeding concentration is 0.7 particles/mm³, corresponding to 700 particles/D ³. Illumination is provided by a Quantronix Darwin-Duo Nd-YLF laser, and images are recorded using four Lavision HighSpeedStar 6 CMOS imagers. With this configuration, the particle image density equals 0.04 ppp at the jet axis. Measurements are acquired at 1000 Hz, corresponding to 35 samples per vortex shedding period. At this frequency, the particle displacement at the jet exit is approximately eight pixels along the jet axis. Further details of the experiment are summarized in Table 4 and available in Violato and Scarano (2011).

Table 4 Experimental parameters for the transitional jet following Violato and Scarano (2011)

Figure 11 shows the instantaneous organization of the four vortex rings in the measurement volume, visualized by isosurfaces of Q-criterion (Hunt et al. 1988), as obtained by tomographic PIV. The measurement is reprocessed using tomographic PTV following the procedure in Novara and Scarano (2013), who used a particle tracking algorithm to compute the velocity material derivative from this experiment for the purpose of evaluating the pressure field. Approximately 9000 particle trajectories are found per measurement volume, corresponding to the expected seeding concentration of 700 particles/D ³. Figure 11, middle shows the resulting 30-snapshot particle trajectories in the measurement volume, color-coded by velocity magnitude.

Fig. 11

Result from the tomographic PIV measurement of the jet in water (left, vortices detected by the Q-criterion, Q = 1100/s²). Reconstructed 30-snapshot particle tracks from tomographic PTV at the full seeding concentration of 700 particles/D ³ (middle). Sub-sampled particle tracks (right, 7 particles/D ³). The tracks are color-coded by velocity magnitude

In order to simulate a low relative seeding concentration, which can occur in several cases as discussed in the introduction, a small subset of the particle trajectories is randomly selected to examine the minimum case of approximately 7 particles/D ³. The resulting particle trajectories are shown in Fig. 11, right. With such a coarse distribution of information, the reconstruction of the instantaneous velocity and associated vortices (visualized by the Q-criterion) using scattered tri-linear interpolation (Fig. 12, left) is not possible.

Fig. 12

Reconstruction of vortices distribution (isosurfaces of Q-criterion, Q = 1100/s²) from sub-sampled tomo-PTV (95 particle tracks, 7 particles/D ³). Linear interpolation (left). VIC without D u /Dt (middle). VIC+ (right). The red dots indicate the instantaneous position of the sub-sampled particles

Following the above numerical and experimental assessment, it is expected that VIC without D u /Dt (Fig. 12, middle) improves upon the linear interpolation. As can be seen in Fig. 12, middle, this is indeed the case and now some parts of all four vortex rings can be recognized by the Q-Criterion visualization. The pattern of coherent and closed vortex rings is, however, not reproduced. Including in the minimization, the measured velocity material derivative (VIC+) yields in this case a significant addition to the reconstruction of the spatial velocity and vorticity distribution, as the vortex rings are retrieved more clearly as shown in Fig. 12, right. In particular, using VIC+, three vortex rings are reconstructed at a seeding concentration of 7 particles/D ³, which corresponds to a total of 95 tracer particles in the measurement volume.

5 Conclusions

The VIC+ method is proposed for reconstruction of instantaneous velocity from time-resolved volumetric PTV measurements, by leveraging the temporal information available by the measurements in form of the velocity material derivative. This is christened as pouring time into space. In both the numerical and experimental assessment considering turbulent boundary layer measurements, VIC+ is able to reconstruct turbulent velocity fluctuations at a fraction of the seeding concentration required for tomographic PIV. In the numerical assessment, coherent hairpin vortices could be identified at a concentration of 500 particles/δ ³₉₉ , whereas for tomographic PIV, a concentration beyond 2000 particles/δ ³₉₉ is required. When VIC is employed without the use of the velocity material derivative, a concentration beyond 1000 particles/δ ³₉₉ is needed. In the real-world experiment, turbulent velocity fluctuations that were damped by approximately 20 % by tomographic PIV could be restored to within 5 % by use of the VIC+ technique. The application to the tomographic experiment of a circular jet flow with a subset of tracer particle trajectories simulates a situation of very low seeding concentration (7 particles/D ³) where VIC+ reconstructs coherent vortex rings that could not be obtained by VIC without the use of the velocity material derivative.

Both the numerical and experimental assessment confirm the advantage of pouring time into space for reconstruction of instantaneous velocity from time-resolved volumetric flow measurements. The proposed technique offers its best performance when dealing with experiments where seeding concentration is limited, which is relevant to a number of cases, including large-scale experiments with HFSB. Moreover, the experiment in the turbulent boundary layer indicates that VIC+ reconstructs the instantaneous velocity field with less spatial averaging of the turbulent fluctuations than the cross-correlation approach used in tomographic PIV.

Appendices

Appendix A: Adjoint procedure

This appendix outlines the adjoint procedure for efficient and exact calculation of the gradient $\partial J/\partial\varvec{\xi}$ for the gradient-based optimization procedure. The equations are derived through the tangent-linear of the VIC+ code, as explained for general cases in Giering and Kaminski (1998). For a detailed discussion on how the adjoint code is derived, the reader is referred to the latter paper and here only the resulting equations are given. Following the naming convention of Giering and Kaminski (1998), adjoint variables are named corresponding to their counterparts in Sect. 2, but with the prefix δ ^*.

The procedure follows in reverse order the steps to calculate the cost function at each iteration (Sect. 2.2.2), and thus starts with the adjoint counterpart of the cost function,

$$\delta^{*} \varvec{u}_{h,p} = 2\left( {\varvec{u}_{h,p} - \varvec{u}_{m} } \right) ,$$

(24)

$$\delta^{*} \frac{{D\varvec{u}_{h,p} }}{Dt} = 2\alpha^{2} \left( {\frac{{D\varvec{u}_{h,p} }}{Dt} - \frac{{D\varvec{u}_{m} }}{Dt}} \right) ,$$

(25)

where the subscript ‘h,p’ indicates the variables on the computational grid, interpolated to the PTV measurement locations. The adjoint of linear interpolation is applied to move from these scattered locations back to the grid. Linear interpolation can be written as a linear system, $\varvec{u}_{h,p} = {\mathbf{A}}\varvec{u}_{h}$, where A is the (sparse) matrix containing the weights for the linear interpolation. The corresponding adjoint becomes (see also e.g., Claerbout 2014),

$$\delta^{*} \varvec{u}_{h} = {\mathbf{A}}^{\text{T}} \varvec{ }\delta^{*} \varvec{u}_{h,p} ,$$

(26)

$$\delta^{*} \frac{{D\varvec{u}_{h} }}{Dt} = {\mathbf{A}}^{\text{T}} \varvec{ }\delta^{*} \frac{{D\varvec{u}_{h,p} }}{Dt} ,$$

(27)

where the superscript ‘T’ means the transpose. Velocity was calculated from vorticity using a Poisson equation (Eq. 1). The Laplacian is self-adjoint, hence the existing Poisson solver can be used to solve,

$$\nabla^{2} \delta^{*} \varvec{R}_{\varvec{\omega}} = \delta^{*} \varvec{u}_{h} ,$$

(28)

where $\delta^{*} \varvec{R}_{\varvec{\omega}}$ is the adjoint variable corresponding to the right hand side of the Poisson equation, $\varvec{R}_{\varvec{\omega}} = - \nabla \times\varvec{\omega}_{h}$. The sub-script ‘h’ is dropped in the remainder of this appendix for clarity. Subsequently,

$$\delta^{*}\varvec{\omega}= - \tilde{\nabla } \times \delta^{*} \varvec{R}_{\varvec{\omega}} ,$$

(29)

where $\tilde{\nabla } = \left[ {\frac{{\tilde{\partial }}}{\partial x},\frac{{\tilde{\partial }}}{\partial y},\frac{{\tilde{\partial }}}{\partial z}} \right]$ is the adjoint of the gradient operator. Finite differences are used for the spatial derivatives, hence the corresponding adjoint operators are the transpose of the finite difference matrices. The radial basis function evaluation (Eq. 11) can also be written as a matrix multiplication, $\varvec{\omega}= {\varvec{\Phi}}\varvec{w}$, where Φ contains the value of the radial basis function (Eq. 10) at each grid point and w is the vector with the weights. The corresponding adjoint equation equals,

$$\delta^{*} \varvec{w}_{1} = {\varvec{\Phi}}^{\text{T}} \delta^{*}\varvec{\omega},$$

(30)

The velocity material derivative was calculated from Eq. (2). The adjoint counterpart of this step is,

$$\delta^{*} \frac{{\partial \varvec{u}}}{\partial t} = \delta^{*} \frac{{D\varvec{u}}}{Dt} ,$$

(31)

$$\delta^{*} \varvec{C} = \delta^{*} \frac{{D\varvec{u}}}{Dt} ,$$

(32)

where $\delta^{*} \varvec{C}$ is the adjoint variable corresponding to the convective part of the velocity material derivative. The adjoint variable $\delta^{*} \varvec{u}_{2}$, is calculated from $\delta^{*} \varvec{C}$ using,

$$\delta^{*} u_{i,2} = \frac{{\partial u_{j} }}{{\partial x_{i} }}\delta^{*} C_{j} + \frac{{\tilde{\partial }}}{{\partial x_{j} }}\left( {u_{j} \delta^{*} C_{i} } \right) ,$$

(33)

written in Einstein notation for clarity. The tilde again indicates the adjoint counterpart of the spatial derivative, which equals the transpose of the finite difference matrix. Subsequently, $\delta^{*} \varvec{w}_{2}$ is calculated from $\delta^{*} \varvec{u}_{2} ,$ analogously to $\delta^{*} \varvec{w}_{1}$ using Eqs. (28)–(30).

Analogous to Eqs. (28) and (29), the adjoint counterpart of the temporal vorticity derivative is obtained from $\delta^{*} \frac{{\partial \varvec{u}}}{\partial t}$,

$$\nabla^{2} \delta^{*} \varvec{R}_{{\partial\varvec{\omega}}} = \delta^{*} \frac{{\partial \varvec{u}}}{\partial t} ,$$

(34)

$$\delta^{*} \varvec{ }\frac{{\partial\varvec{\omega}}}{\partial t} = - \tilde{\nabla } \times \delta^{*} \varvec{R}_{{\partial\varvec{\omega}}} ,$$

(35)

The temporal derivative of vorticity is calculated from the vorticity transport equation (Eq. 4). The adjoint counterpart yields,

$$\delta^{*} \omega_{i,2} = \frac{{\partial u_{j} }}{{\partial x_{i} }}\delta^{*} \frac{{\partial \omega_{j} }}{\partial t} - \frac{{\tilde{\partial }}}{{\partial x_{j} }}\left( {u_{j} \delta^{*} \frac{{\partial \omega_{i} }}{\partial t}} \right) ,$$

(36)

$$\delta^{*} u_{i,3} = \frac{{\tilde{\partial }}}{{\partial x_{j} }}\left( {\omega_{j} \delta^{*} \frac{{\partial \omega_{i} }}{\partial t}} \right) - \frac{{\partial \omega_{j} }}{{\partial x_{i} }}\delta^{*} \frac{{\partial \omega_{j} }}{\partial t} ,$$

(37)

where again the Einstein notation is adopted for clarity. Subsequently, $\delta^{*} \varvec{w}_{3}$ is calculated from $\delta^{*}\varvec{\omega}_{2}$ using Eq. (30) and $\delta^{*} \varvec{w}_{4}$ is calculated from $\delta^{*} \varvec{u}_{3}$ using Eqs. (28)–(30). The gradient $\partial J/\partial\varvec{\xi}_{\varvec{\omega}}$ is then calculated by,

$$\frac{\partial J}{{\partial\varvec{\xi}_{\varvec{u}} }} = \delta^{*} \varvec{w}_{1} + \delta^{*} \varvec{w}_{2} + \delta^{*} \varvec{w}_{3} + \delta^{*} \varvec{w}_{4}$$

(38)

Analogously, the gradients $\partial J/\partial\varvec{\xi}_{{\varvec{BC},\varvec{u}}}$ and $\partial J/\partial\varvec{\xi}_{{\varvec{BC},\partial \varvec{u}}}$ for the boundary conditions are calculated.

Appendix B: VIC+ Pseudo-code