Polynomial element velocimetry (PEV): a technique for continuous in-plane velocity and velocity gradient measurements for low Reynolds number flows

Chaminda R Samarage; Josie Carberry; Gregory J Sheard; Andreas Fouras

doi:10.1088/0957-0233/23/10/105304

1. Introduction

The particle image velocimetry (PIV) technique, which has been under development for over two decades [1, 2], has been gaining popularity as a standard flow measurement technique, especially for medical applications [3–6]. In such applications, due to the difficult imaging conditions [3, 7], small flow structures are difficult to resolve, making accurate flow gradients difficult to determine [7]. Fluid mechanics properties derived from velocity gradients, such as the shear stress and vorticity, are important quantities that are required to make a physical interpretation of the flow.

Both the accuracy and spatial resolution of the flow measurements are important in determining an accurate representation of velocity gradients within these flows. There is a trade-off between the accuracy of PIV in measuring the velocity and its associated gradients and the spatial resolution of the PIV velocity measurements [8, 9]. The general PIV technique offers measurements that are discretized into sampling windows and there is a limit to which the size of a given sampling window can be reduced, i.e. physically a unit pixel. Increasing the sampling window size reduces the random error of the measurement (due to the increased amount of information that contribute to the cross-correlation map) while reducing the spatial resolution of the measurements [10]. The accuracy also depends on factors such as the particle seeding density [10] and image noise [11]. Decreasing the sampling distance between these measurements/sampling windows increases the spatial resolution of the measurements, but leads to random errors in the estimated flow gradients [8, 9]. Any underlying flow gradients within the sampling window are also captured within the cross-correlation map. However, with the general PIV technique, only the location of the maximum signal from the cross-correlation map is utilized. This maximum signal represents the modal velocity, and all other information stored in the cross-correlation map is discarded.

Existing techniques for improving the accuracy in the estimated velocity gradient include recursive local correlation [12], single pixel methods [13–15], image deformation methods [16–18], least-squares pattern matching methods [19–21] and more recently correlation-based least-squares methods [22–25]. At least three methods in the literature [19, 20, 32] are also capable of simultaneously and directly measuring the gradients in the flow in addition to the velocity. For completeness, a selection of the aforementioned methods is reviewed here.

1.1. Recursive local correlation

Hart [12, 26] proposed a technique to improve spatial resolution using a recursive correlation technique with a correlation-based error correction method. This method [26] involves multiplying adjacent correlation maps during PIV interrogation to remove spurious vectors. This error correction allows for discrete measurements at higher spatial resolutions using the recursive correlation technique proposed by Hart [12]. Due to the product operator in the correlation-based correction process, all of the information within the cross-correlation map is utilized to obtain the velocity measurement. However, the quality of the information within the cross-correlation map is greatly deteriorated in cases with low tracer seeding and/or high levels of image noise.

1.2. Single-pixel methods

Single-pixel resolution methods pioneered by Westerweel et al [13] use 1 × 1 pixel interrogation windows to obtain time-averaged single-pixel resolution velocity measurements from an image sequence. This method is capable of resolving small flow structures at the cost of utilizing a large number of image pairs, typically over 1000. Scharnowski et al [15] improved the two-point ensemble correlation [14] to determine Reynolds stresses at single-pixel resolution. This method determines the stresses based on the shape of the cross-correlation peak. A Gaussian fit-function is applied (as opposed to the optimization routine described in this paper) to determine the shape of the correlation peak. Using the peak shape characteristics, the probability distribution function (PDF) of the velocity is computed and from this, the Reynolds stresses. Although this technique is capable of resolving small-scale flow structures, single-pixel ensemble is applied to 20 000+ image pairs, and the technique is sensitive to the number of image pairs utilized in the ensemble correlation average.

1.3. Image deformation methods

In recent times, image deformation has become a part of the standard PIV interrogation process. The fundamentals of this technique date back to Huang et al [27] who proposed a particle image distortion technique intended to improve the reliability of the PIV technique in flows with high velocity gradients. Image deformation, as it is now known [16, 17], addresses some of the weaknesses of PIV relating to spatial resolution and peak locking [28]. This technique involves conducting standard PIV on the raw images and fitting the resulting velocity data to create a predictor field. This predictor field is used to distort the images and interrogation is conducted again at a finer resolution. The process is repeated iteratively until a convergence criterion is met. The method is sensitive to the interpolation schemes used to create the predictor fields and the image distortion. At least two studies have investigated optimizing these routines [29, 30]. Scarano [18] introduced a correction method that involves optimizing four correlation maps separately to directly determine the second-order derivatives for the two displacement components in x and in y. The derivatives are used to correct the displacement obtained from the cross-correlation analysis. This method demonstrated improved accuracy in determining the velocity gradient of an experimental wall jet flow with CCD noise and under-resolved particle images.

1.4. Least-squares image/pattern matching methods

Least-squares matching is a similar concept to image deformation and involves an optimization routine to match image segments by geometric deformation. This optimization can be run globally, using a routine to ensure global continuity of the solution such as that implemented in the correlation image velocimetry (CIV) technique [19, 31], or locally. Upon convergence, measures of the local velocity and the deformation rate of the image are simultaneously determined. Ruan et al [20] proposed a method (referred to as direct measurement of vorticity) to directly determine the vorticity of particle images using a local optimization routine. Kitzhofer et al [32] applied a similar local least-squares matching routine to 3D cuboids of volumetric data.

1.5. Advanced correlation-based methods

Recently, a set of PIV techniques using a least-squares correlation map matching approach have obtained 3D velocity measurements using the depth information inherent in micro-PIV flows [22, 24] and holographic flows [25]. These methods are established on a concept that the convolution of the probability density function of the velocity and the auto-correlation of the particle image closely resembles the shape of the cross-correlation peak [15, 22, 33]. To decode the depth information, holographic correlation velocimetry [25] uses the particle image diffraction pattern, while volumetric correlation velocimetry [22] uses the out-of-focus effects on the particle image. While these methods are capable of offering high-resolution measurements in the out-of-plane axis, the in-plane velocity measurements are still discretized into the interrogation sub-regions in which the total cross-correlation map was calculated.

1.6. Summary

The above-mentioned methods are optimized for their respective applications. The general PIV technique offers a single measurement at each interrogation region, which is obtained from the highest signal in the cross-correlation map representing the modal velocity within the underlying region of the flow. The accuracy in estimating velocity gradients is affected by the uncertainty in measuring the velocity and the spatial resolution of these measurements. Single-pixel methods offer measurements with a very high spatial resolution at the cost of a large number of image pairs. Hence, single-pixel methods trade temporal resolution for spatial resolution. Least-squares matching methods directly determine a continuous measurement for the velocity and its associated gradients by matching reconstructed particle volumes or image segments. While the accuracy of least-squares matching methods in estimating the gradient is not affected the same way as with general PIV, it has yet to be thoroughly studied.

A background on the factors that affect the accuracy of estimating flow gradients using standard PIV is presented in section 2. The technique that is the subject of this manuscript is described in section 3. In section 4, the proposed method is applied to synthetic computer-generated images for flow in a channel and flow past a square cylinder. The simulations are validated by applying the proposed method to two laboratory-based flows in a channel and the wake behind a circular cylinder at a low Reynolds number [34].

2. Accuracy of flow gradient measurements

PIV, in its standard form, involves dividing a pair of consecutively acquired images into several discrete domains called sampling windows. Upon subsequent PIV interrogation, each sampling window yields a two-component velocity measurement that closely represents the velocity of the tracer particles in the underlying flow. The accuracy in estimating the velocity gradient is dependent on factors such as the imaging resolution and the tracer particle seeding density. However, this accuracy is also dependent on the computational method used to determine the flow gradients. The velocity gradients can be estimated through finite-difference techniques [35], polynomial fitting [8] on a discrete grid of velocity samples or direct estimation using a least-squares image/pattern matching approach [19]. The accuracy of the estimated velocity gradient is dependent on the spatial resolution and the accuracy of the PIV velocity measurements. The sampling distance between the velocity measurements determines the spatial resolution of PIV measurements. The uncertainty in the PIV velocity measurement is dependent on factors such as image noise and the particle seeding density. The strength of the cross-correlation signal is mainly dependent on the image seeding density and it has been shown that at least seven particle images are optimal [10]. However, in most cases it is difficult to conduct the experiment multiple times to optimize seeding densities. In such cases, it is advisable to use the sampling window size as a means to effectively improve the accuracy of the PIV velocity measurements.

A Monte Carlo simulation utilizing computer generated synthetic images was conducted to study the accuracy of PIV in estimating flow gradient parameters such as the vorticity, as a function of the sampling window size and sampling distance. Synthetic images of 2048 × 2048 pixel² were generated for a Lamb–Oseen vortex flow [36], where the tangential velocity, u_t, and out-of-plane vorticity, ω_z, are given by

$\begin{equation} \mu _t r = \frac{\Gamma }{{2\pi r}}\left[ {1 - \exp \left( {\frac{{ - r^2 }}{{L^2 }}} \right)} \right], \end{equation} \tag{ 1 }$

$\begin{equation} \omega _z r = \frac{{ - \Gamma }}{{4\pi \nu t}}\exp \left( {\frac{{ - r^2 }}{{L^2 }}} \right), \end{equation} \tag{ 2 }$

where r² = X² + Y² and (X, Y) are global in-plane Cartesian coordinates, Γ is the circulation, ν is the kinematic viscosity and t is time. The core radius, L, of the Lamb–Oseen vortex is defined as

$\begin{equation} L = \sqrt {4\nu t} , \end{equation} \tag{ 3 }$

which describes the long-timescale viscous decay of the vortex.

The numerical analysis has been conducted for a circulation, Γ, of 1000 px² s⁻¹ and a core radius, L, of 100 pixels. All images are generated with 256 gray levels. The PIV vorticity error in the measurement is the difference between the estimated vorticity from PIV velocity measurements and the analytical solution from equation (2). The sampling window size, W, is varied from 8 pixels to 256 pixels. Varying the sampling distance, Δ, varies the window-overlap ratio and is varied from 8 pixels to 256 pixels. Figure 1 shows a contour plot of the RMS of the PIV error in estimating the in-plane vorticity, as a function of the sampling window size, W, and the sampling distance, Δ, normalized to the vortex diameter, 2L. The vorticity field is calculated from the in-plane PIV velocity measurements using the least-squares polynomial fit method described in Fouras and Soria [8].

The RMS error is optimized over a limited range of sampling window size and sampling distance. Under the given flow conditions, a PIV interrogation scheme utilizing 64 × 64 pixel² sampling windows (W/2L = 0.32) with a sampling distance of 16 pixels (Δ/2L = 0.08 or 75% window overlap) optimizes the estimated PIV vorticity. There is almost no variation in the error with sampling distance below Δ/2L = 0.16, and the error is primarily a function of the sampling window size, W. As the sampling distance is increased above Δ/2L = 0.32, the spatial resolution of velocity measurements used to estimate the vorticity is poor. There is a significant loss in accuracy of the estimated vorticity in this region of the parameter space.

As evident in the contour map shown in figure 1, the accuracy of the estimated flow gradient is dependent on the spatial resolution of the PIV data points as well as the sampling window size used for the PIV interrogation scheme. The effects of the sampling window size and sampling distance on the accuracy of estimated flow gradient can be explained with the concepts of random errors and bias errors. It is well known that the total error in estimating the velocity gradient is the sum of random errors and bias errors. Figure 2(a) illustrates the random error, bias error and the total error in estimating the velocity gradient as a function of the sampling distance between adjacent velocity measurements. This line plot is indicative of a single slice in the contour map in figure 1. Consider the optimal center of the graph where the total error is minimized. As the sampling distance is decreased (i.e. at higher spatial resolutions), the flow variations are fully captured by the numerous measurements and bias error in the gradient estimate is minimized. In this case, the velocity gradient error is primarily a function of the random error in the numerous velocity measurements contributing to the gradient estimation [8]. As the sampling distance is increased (i.e. at low spatial resolutions), the uncertainty in the measured velocity or random error transmitted in to the velocity gradient estimation method is reduced. The flow variations are not represented well by the low resolution and introduce bias errors to the gradient estimate [8]. The random error and the bias error cannot be minimized simultaneously resulting in an optimal spatial resolution as shown in the graph. In addition, random errors and bias errors introduced in to the velocity gradient estimate are also a function of the sampling window size as evident in figure 1.

Figure 2(b) shows a conceptual representation of the velocity gradient error map (depicted by the gray contours) presented in figure 1. Increasing the sampling window size improves the signal strength of the cross-correlation map and the shape of the cross-correlation peak is distorted due to underlying flow variations [37]. While the flow variations are well represented within the cross-correlation map, subsequent PIV interrogation, which utilizes the signal peak, yields the modal velocity over the sampling domain and bias errors are introduced into the velocity gradient calculation. Decreasing the sampling window size produces a cross-correlation map with poor signal strength and random noise due to the inadequate number of particle image pairs that contribute to the cross-correlation [10]. Typically at least seven particle image pairs are required for the cross-correlation map to have adequate signal strength for ideal measurements [10]. PIV interrogation on a cross-correlation map with noise artifacts leads to a higher probability of the interrogation routine selecting a noise peak instead of the signal peak. This leads to random velocity measurements or outliers [1]. While bias errors are minimized, there is a higher transmission of random error in to the gradient estimation process. The level of random error that is transmitted to the gradient estimation method can be reduced by vector validation routines [38, 39] to remove outliers in the velocity measurements. These underlying processes discussed here are the key factors that contribute to the elliptical shape of the contours that is evident in figure 1, and conceptually presented in figure 2(b).

We propose a method capable of velocity measurements with high spatial resolution utilizing a single pair of images. This method is capable of directly determining the velocity gradient while minimizing the aforementioned random errors and bias errors simultaneously.

3. The proposed technique

A novel method to directly determine a piece-wise grid of two-dimensional polynomials defining the entire flow is presented. The method, which the authors term polynomial element velocimetry (PEV), solves for a tensor-product of polynomials by modeling an estimate for the cross-correlation map and performing a nonlinear minimization with the measured cross-correlation map obtained from standard PIV.

Graphical representations of the steps involved in the PEV method are illustrated in figures 3–5. The overall flow field is discretized into a grid of regular regions called elements (I), as shown in figure 3. These elements are analogous to elements in spectral element methods used in computational fluid dynamics (CFD) [40–43]. Standard PIV with smaller sampling windows is conducted over the larger PEV elements (II). PIV measurements local to the region of flow inside the element are extracted (III) to determine a first approximation for the underlying flow. Although not a necessary condition for convergence, this first approximation reduces the time to convergence for the solver. Figure 4 shows the PEV process implemented over a single element in the flow field. The measured cross-correlation and auto-correlation maps are computed using two image snapshots for each PEV element (IV). Two-dimensional polynomials are fitted to the extracted PIV measurements to describe the local horizontal and vertical components of the flow within the element (V).

**Figure 3.** Graphical representation of the element discretization process. (I) The entire flow field is discretized into a grid of regular regions called elements. (II) Standard PIV with smaller sampling windows is conducted over the entire flow field and the PIV data local to the element is extracted to be utilized toward a first guess for the solver (III). *Note: Image data, PIV data and cross-correlation data shown in the figure are for illustrative purposes and are not a representation of experimental data.*
Download figure:
Standard image

**Figure 5.** Graphical representation of the process used to determine an estimate for the cross-correlation map (as depicted by process (VI) in figure 4). (VI.I) The estimated PDF, P, is determined from flow polynomials and the image intensities within the element. (VI.II) The PDF is convolved with the measured auto-correlation map, A, to determine an estimate for the cross-correlation map, C. *Note: Image data, PIV data and cross-correlation data shown in the figure are for illustrative purposes and are not a representation of experimental data.*
Download figure:
Standard image

The flow polynomials are used to estimate for the cross-correlation map (VI) within the element as shown in figure 4 and further expanded in figure 5. The PEV method applied to flows in this study uses a two-dimensional bi-cubic polynomial for each component of the local displacement. The flow polynomials for the horizontal component (Δx(x, y)) and the vertical component (Δy(x, y)) of the local displacement are given by

$\begin{eqnarray} &&\Delta x(x,y) = \sum\limits_i^3 {\sum\limits_i^3 {a_{ij} x^i y^j } } \;{\rm and}\;\Delta y(x,y) = \sum\limits_1^3 {\sum\limits_1^3 {b_{ij} } x^i y^j } ,\nonumber\\ \end{eqnarray} \tag{ 4 }$

where 0 ⩽ i + j ⩽ 3, and x and y are spatial coordinates in the cross-correlation map that range from 0 to 1 over the element. a_ij and b_ij are the coefficients of the polynomial that describes the horizontal component and vertical component of the local velocity, respectively.

The polynomials determine the ability of the technique to represent the flow and the order can be varied depending on flow conditions or gradients. For complex flows or special cases exhibiting spatial symmetries, Lagrangian polynomials or Fourier polynomials could instead be utilized. For example, it may be possible to represent discontinuities found in supersonic/hypersonic flows through the use of step functions. The polynomials for the displacement for Δx and Δy (equation (3)) can be divided by the time-interval between two successive snapshots to determine polynomials that define a map for the horizontal component (u) and the vertical component (v) of the local velocity. These displacement polynomials (equation (3)) can also be directly differentiated to obtain a continuous measure for the velocity gradients or vorticity, which is given by

$\begin{equation} \omega z \cong \frac{1}{{\Delta t}}\left[ {\frac{\partial }{{\partial x}}(\Delta y) - \frac{\partial }{{\partial y}}(\Delta x)} \right], \end{equation} \tag{ 5 }$

or directly integrated to determine quantities such as the volumetric flow rate.

It has been demonstrated in Fouras et al [33], and more recently in Scharnowski et al [15] that the cross-correlation, C, is the convolution of the auto-correlation map, A, of the image and the PDF, P for the displacement. In PEV, the PDF of the local flow velocity is constructed using the flow polynomials and convolved with the image intensity map, I, within the element as depicted by process (VI.I) in figure 5. This modification to the original concept allows for particle images to contribute more to the PDF reducing the effect of regions in the flow with no particle images. The cross-correlation is estimated by PEV using

$\begin{equation} C(x,y) = \left[ {I(x,y) \otimes P(x,y)} \right] \otimes A(x,y). \end{equation} \tag{ 6 }$

The estimated cross-correlation, derived from the initial fit to PIV data and the method above (VI.II), typically differs from the measured cross-correlation map. A nonlinear minimization routine is used to find the local flow polynomials that minimize the error between the estimated and measured cross-correlations within the element, shown by process (VII) in figure 4. In this study, the Levenberg–Marquardt routine [44, 45] with a Tikhonov regularization method [46] is used to compute the nonlinear optimizations and is repeated until local convergence criteria are satisfied (VIII). This combination of a nonlinear solver with regularization has been previously applied for least-squares routines [22–24]. The solver is used to minimize the cross-correlations within all the elements that define the flow. Since the flow should be continuous through the elements, a regularization method [46] is implemented during the nonlinear minimization process of each element. This method relies on adding weightings to minimize discontinuities in both the velocity components and their spatial derivatives at element interfaces. If necessary, with complex flows, this method allows for specific boundary conditions to be set at element interfaces. For a single element, there is some ambiguity in the resulting flow polynomials since the PDF can have the same shape for different flow fields inside the element (specifically, the direction of any rotating component is unknown). However, in a practical circumstance the flow field spans over a number of elements. As a result of the communication between element boundaries, this ambiguity is greatly reduced.

As explained above, PEV optimizes for the measured cross-correlation map over the sampling region (element) to determine a continuous map for the velocity and its associated gradients. With image deformation methods, the image is iteratively deformed using subsequent cross-correlation analysis and PIV data. These methods offer high spatial resolution measures while improving the accuracy of PIV over flows with high gradients. With higher order optical flow methods such as the CIV method [19] an optimization routine is performed directly on the image to determine the velocity and the velocity gradient without the need for added tracer particles.

To validate the proposed method, PEV is applied to synthetic flows as well as laboratory flows. PEV is applied to a single pair of images to obtain instantaneous measurements for the velocity and the vorticity and compared to instantaneous PIV measurements obtained for the same image pair. Similar to PIV, time-resolved PEV measurements are available by PEV analysis on subsequent image pairs. Table 1 summarizes the flow cases that are used to validate the measurement accuracy of PEV in comparison to that of PIV.

Table 1. Details of flow cases that are used in this report to validate PEV against PIV.

Case	Nature	Description	Flow conditions
A	Synthetic	Channel flow	u = 0 pixels, v_max = 10 pixels, Parabolic flow profile
B	Synthetic	Flow past a square cylinder	Re = 30, D = 512 pixels
C	Experimental	Channel flow	Re = 0.2, Q = 5 ml min⁻¹
D	Experimental	Flow behind a circular cylinder	Re = 30

4. Validation using synthetic flows

4.1. Synthetic image generation

To validate its accuracy, the proposed PEV method is applied to two synthetic flows. As shown in table 1, case A represents the fully-developed two-dimensional flow between closely-spaced parallel plates and case B is the flow around square cylinder at a Reynolds number of 30. With case A, the image pair generated has a resolution of 64 × 1024 pixels with a noise ratio of 0.05%. White noise is added to the image by generating a random number that has a Gaussian distribution with a standard deviation that is 0.05% of the maximum pixel intensity. The tracer particles used were assumed to have a Gaussian shape with a particle diameter of 4 pixels. The particle image diameter is defined as the full-width-half-max of the Gaussian function used to generate the particle image. The image is seeded with tracer particles at a seeding density of 0.0625 particles per pixel.

For the cylinder flow (case B), a synthetic image pair of resolution 1920 × 1920 pixels was generated with 0.05% image noise and a particle seeding density of 0.02 particles per pixel. In the image sequence, the square cylinder has a side-length measuring 512 pixels. The particle images have a diameter of 4 pixels. The flow around a square cylinder at a Reynolds number of 30 was simulated using in-house CFD software and the particles were displaced based on the resulting CFD solution. Figure 6 shows the synthetic image pair used for this flow case. The top and bottom halves of the image show a single snapshot and the overlapped snapshots, respectively. The reader is referred to the papers by Sheard et al [41–43], for further information on the CFD software and the meshes used to model the cylinder with square cross-section used in this study.

**Figure 6.** The synthetic image pair utilized for case B (table 1). The top half of the image shows a single snapshot of the image pair, while the bottom half of the image shows the two snapshots overlapped to help visualize the displacement of the particle images. Tracer particle images, assumed of Gaussian shape, are displaced based on a CFD solution for the flow past a square cylinder at Re = 30.
Download figure:
Standard image

4.2. Channel flow

PEV is applied to a synthetic flow in a channel and its performance compared to measurements obtained with standard PIV on the same flow. The particle images are displaced vertically with a parabolic profile with a maximum displacement of 10 pixels in a channel that has a radius, R = 32 pixels. Further details of the image generation process can be found under section 4.1. For this case the particles have no horizontal displacement (Δx = 0). The vertical displacement (Δy) and the vorticity (ω_z) are given by

$\begin{equation} \Delta y = 10\left( {1 - \frac{{(x - r)^2 }}{{R^2 }}} \right), \end{equation} \tag{ 7 }$

$\begin{equation} \omega _z = \left( {\frac{{ - 20}}{{R^2 }}} \right)(x - r), \end{equation} \tag{ 8 }$

where R is the radius of the channel, and r is the distance from the center of the channel.

Figure 7(a) shows a single snapshot of the synthetic image pair utilized for this case. Figures 7(b) and (c) show the in-plane vorticity contours and displacement vectors calculated using standard PIV (7b) and PEV (7c). The standard PIV interrogation involves performing a multi-grid recursive scheme from 64 × 64 pixel sampling windows down to 32 × 32 pixel sampling windows. A window of overlap ratio of 75% (sampling distance of 8 pixels) is used. As discussed earlier, the error in estimating velocity gradients, and in turn, estimating the vorticity is dependent on the window size and sampling distance. Therefore, the settings used here were determined to be the most optimal PIV interrogation settings based on the study discussed in section 2. Image deformation is implemented with bilinear velocity interpolation and bi-cubic image interpolation as carried out by Fouras et al [47]. The image deformation process is iterated for two passes. The in-plane vorticity is estimated using the least-squares method described by Fouras and Soria [8]. The flow is discretized to 1 × 3 elements, i.e. a sampling window size of 64 × 64 pixels. PEV is applied using a bi-cubic set of flow polynomials. The in-plane vorticity within an element is estimated directly using the derivatives of the flow polynomials.

Comparisons between the measured displacement and the estimated vorticity using the two methods, are conducted at y = 192 pixels, illustrated in figures 7(d) and (e). The PIV interrogation process yields five displacement measurements across the channel width. For clarity, only every second displacement measurement in the streamwise direction is shown in figure 7(b). For ease of comparison with figure 7(b), the continuous PEV field shown in figure 7(c) is sampled at the same locations as the PIV measurements. Random noise is present in the measurements obtained from PIV (figures 7(b), (d)). These measurements are instantaneous and extracted at a single slice across the channel at y = 192 pixels.

With standard PIV, the closest velocity measurement to the channel wall that can be obtained is within half of the sampling window size of the wall. However, for improved measurements in the wall-normal direction and for fair comparison, the PIV displacement measurements (⧫) are linearly extrapolated to the channel boundaries based on the no-slip boundary condition. The extrapolated measurements are depicted with open symbols (◊). Within the PEV solver, zero Dirichlet boundary conditions are imposed at the edges of the channel to restrict the solver from overfitting the underlying velocity profile. This is implemented using the regularization procedure explained in section 3. PEV analysis yields measurements (blue line) that are continuous across the channel width in comparison to the discrete PIV measurements. However, the PEV method overestimates the velocity near the channel wall and underestimates the velocity at the center of the channel. While this is expected of PIV analysis with large sampling windows, no such bias exists with PEV. The RMS error between the measured displacement using PIV and the theoretical flow profile is 0.532 pixels, while the RMS error between the displacements measured using PEV and the theoretical flow profile is 0.475 pixels. This accounts for an error reduction of 11% in measuring the displacement using PEV over PIV. On estimating the in-plane vorticity, the RMS error obtained with PIV is 0.115 s⁻¹, while the RMS error with PEV is 0.026 s⁻¹ resulting in a reduction of 77% in the RMS error with PEV over PIV. As evident in figure 7(c), for this first flow case, PEV produces a measurable improvement over PIV for the velocity measurement; however, PEV produces a significant improvement over PIV in estimating the vorticity, especially near the channel walls.

Obtaining unbiased, highly resolved flow measurements using PIV is difficult since the channel is only 64 pixels wide in the wall-normal direction. While PIV has been applied to flows similar to this [7], PIV is not ideal in similar harsh conditions where there are insufficient pixels to obtain a sufficient number of measurements to resolve structures in the flow. PEV, however, is well suited for this application as it optimizes for the shape of the cross-correlation map to yield a map for the local velocity and its associated gradients.

Figure 8 shows the comparison between the cross-correlation maps obtained for the region in the channel, denoted by B in figure 7. The 2D contour map of the measured cross-correlation map is compared against the cross-correlation map derived using the estimated flow polynomials from PIV data, and the final cross-correlation map from the PEV solver after convergence.

Distortion of the cross-correlation peak due to flow gradients in the underlying flow is evident on the 2D contour map in figure 8(a). To reduce the convergence time of the PEV solver, a series of polynomial functions are fitted to the PIV data to determine an initial starting point for the solver. Figure 8(b) shows the 2D contour map of the cross-correlation reconstructed using these estimated polynomials. While this step reduces correlation noise in the cross-correlation map, much of the peak definition, i.e. information on the underlying flow, is lost. Figure 8(c) shows the contour map of the cross-correlation map after the PEV solver has converged. The estimated cross-correlation map has reduced levels of noise and closely resembles the measured cross-correlation map. In contrast to PIV where a single window yields a single measurement, upon convergence, PEV yields a detailed map of the local velocity.

4.3. Flow past a square cylinder

The proposed PEV method is applied to a more complex flow, namely the flow past a square cylinder at a Reynolds number of 30. The images are generated using the methods described in section 4.1. Figures 9(a) and (b) show contours for the in-plane vorticity and vectors for the displacement as measured by standard PIV with image deformation and PEV. For the purpose of clarity, only every 16th vector in x and every 4th vector in y are shown in figure 9(a). With the plot in figure 9(b), the PEV measurements are sampled at the same locations as the PIV measurements to be consistent with the plot in figure 9(a). The PIV interrogation scheme utilized for this flow case is the same as that for case A since it is the set of parameters that optimize the accuracy of the estimated velocity gradient as described in sections 2 and 4.2. The PEV analysis is conducted over 15 × 15 elements (i.e. 128 × 128 pixel sampling windows). Comparison line plots for the displacement and vorticity at x/D = 0 are shown in figures 9(c) and (d), respectively.

**Figure 9.** Contour plots of the estimated in-plane vorticity using standard PIV (a) and PEV (b). The vectors show the displacement and have been skipped for clarity. (c, d) Comparisons of the measured displacement (c) and the in-plane vorticity (d) profiles for PIV ( $\protect\epsfbox{images/mst428424fun01.eps}$ ), PEV (blue line) and the theoretical flow (red line) are shown for slice A–A at *x/D* = 0. PEV yields displacement measurements that are 45% more accurate over PIV, and vorticity measurements that are 15.5% more accurate over PIV.
Download figure:
Standard image

$\protect\epsfbox{images/mst428424fun01.eps} $ — **Figure 9.** Contour plots of the estimated in-plane vorticity using standard PIV (a) and PEV (b). The vectors show the displacement and have been skipped for clarity. (c, d) Comparisons of the measured displacement (c) and the in-plane vorticity (d) profiles for PIV ( $\protect\epsfbox{images/mst428424fun01.eps}$ ), PEV (blue line) and the theoretical flow (red line) are shown for slice A–A at *x/D* = 0. PEV yields displacement measurements that are 45% more accurate over PIV, and vorticity measurements that are 15.5% more accurate over PIV.
Download figure:
Standard image

PIV and PEV both estimate the flow displacement well in the flow regions further from the cylinder walls. However, there are some slight discrepancies with the PIV measurements nearer the wall. The RMS error in measuring the horizontal displacement with PIV is 0.305 pixels while the proposed PEV method reduces this measurement accuracy by 45% to 0.169 pixels. As evident in figure 9(d), the estimates for vorticity using PIV yields measurements that are noisy. The RMS error between the estimated vorticity using PIV and the theoretical vorticity from the CFD package is 0.009 73 s⁻¹. PEV yields an estimate for the vorticity field with an RMS error of 0.008 22 s⁻¹ resulting in a 15.5% reduction in the RMS error. With both methods, there exist discontinuities in the velocity gradient measurements that are visible on the 2D contour plots in figures 9(a) and (b). In the case of PIV, reducing the spatial resolution of the velocity measurements would greatly reduce the random noise in the estimated vorticity measurements. However, as discussed in section 2, the accuracy of the vorticity estimate is reduced due to bias errors. With PEV, the elements can be overlapped in a similar method to window overlapping in PIV, to reduce the discontinuities in the estimated velocity gradient (refer figure 14). In the instance that time-averaged results are required, random fluctuations in the measurements can be reduced with temporal averaging processes [48].

Figure 10 shows the effect of smoothing the measurements obtained from PIV at x/D = −0.25 where the vorticity is stronger than at x/D = 0. A moving-average filter is applied to the PIV velocity measurements, denoted by $\protect\epsfbox{images/mst428424fun01.eps}$ , using a 5 × 5 kernel to greatly reduce the random noise fluctuations in the PIV velocity measurements. The smoothed PIV measurements are denoted by Δ in figure 10(a). For the PIV and the smoothed PIV velocity measurements, the in-plane vorticity is calculated using the method proposed by Fouras et al [8] and is presented in figure 10(b). While smoothing reduces random errors in the velocity measurements, bias errors are still present in the estimated PIV velocity gradient. This is evident in the vorticity measurements near the cylinder wall in figure 10(b). The PIV data could be further smoothed with multiple passes but introduces additional bias errors to the PIV measurements. PEV estimates the in-plane vorticity with reduced random errors, which is inferred by the smooth, noise free measurements as evident on the line plot. The vorticity estimates from PEV also have minimal bias errors in contrast to PIV and its smoothed counterpart.

5. Application to laboratory flows

5.1. Channel flow

To be consistent with the synthetic experiments, PEV is applied to two similar laboratory flows: a simple flow through a circular channel and the flow behind a circular cylinder. The channel flow is generated using a Perspex model of a circular channel with a 9.25 mm inner diameter. To match the refractive indices between water and the Perspex material, glycerine is mixed with the water. The flow is seeded with a glycerol solution (75% glycerol, 25% water) with 35 micron hollow glass spheres and illumination is achieved through a continuous Nd:YAG laser (Melles Griot). The flow is imaged with a PCO2000 (PCO AG) with a 105 mm lens (Nikon) at 14.7 frames per second (full frame) with an exposure time of 20 ms. The solution is pumped using a Harvard Apparatus syringe pump generating a flow rate of 5 ml min⁻¹. The flow in the channel has a Reynolds number of approximately 0.2.

Figure 11(a) shows a single image of the channel after average image subtraction to remove reflections. The top half of the image has relatively good image signal to noise ratio in comparison to the bottom half of the image, making it easier to compare the performance of both PIV and PEV methods in harsh imaging conditions. A standard PIV interrogation scheme uses an iterative scheme with 64 × 64 pixels sampling windows to 32 × 32 pixels sampling windows with 75% window overlap (i.e. a sampling distance of 8 pixels). A 2-pass image deformation process involving a bilinear velocity interpolation and a bi-cubic image interpolation is applied during interrogation. PEV is conducted with 1 × 3 elements (each 192 × 192 pixels) with a 2D bi-cubic polynomial describing the flow in each element.

Figures 11(b) and (c) show contours for the in-plane vorticity and displacement vectors, as estimated by PIV and PEV, respectively. For the purpose of clarity, only every fourth vector in y is shown in figure 11(b). For comparison, only the PEV measurements sampled at the same locations as the PIV measurements are shown in figure 11(c). Comparisons between the measurements using PIV and PEV are made at y = 96 pixels (figures 11(d) and (e)) and at y = 480 pixels (figures 11(f) and (g)). The two locations were chosen to be studied due to the large variation in the image signal to noise ratios between the two halves of the image.

As in section 4.2, the number of PIV displacement measurements in the wall-normal direction is increased by extrapolating the actual measurements, depicted here by filled symbols (⧫), to 0 at the channel wall. The extrapolated measurements are depicted by open symbols (◊). The displacement measured by PIV is noisy at y = 96 pixels (figure 11(f)) with poor image signal-to-noise ratio and poor seeding reducing the signal-to-noise ratio of the cross-correlation map. This random noise in the displacement measurements contributes to erroneous vorticity measurements (figure 11(g)), especially near the channel wall. However with PEV (blue line), as the element size is 192 × 192 pixels, the measured cross-correlation map captures more information within the image leading to a cross-correlation with an adequate signal-to-noise ratio as shown by 2D contour map in figure 12(a). The 2D contour map of the final cross-correlation map from PEV evaluation is shown in figure 12(b). PEV yields a more accurate measurement for the displacement and estimates the vorticity with dramatically reduced measurement noise, which we can infer to be a greatly reduced random error transmission. The expected distortion of the peak is captured well and has reduced levels of noise in contrast to the measured cross-correlation map thereby contributing to PEV's low-pass filter like feature.

**Figure 12.** (a) 2D contour map of the measured cross-correlation map at the base of the channel flow in figure 11(a). Although the image quality in this region is poor, the cross-correlation map captures underlying flow features, which is evident by the peak's distortion. The 2D contour map of the estimated cross-correlation map following PEV evaluation is shown in (b).
Download figure:
Standard image

At higher image signal-to-noise ratios and higher seeding densities, as is the case at y = 480 pixels (figure 11(d)), PIV yields improved displacement measurements. As evident in figure 11(d), the particle seeding density limits the amount of oversampling that can be done before PIV yields adjacent displacement measurements that are similar. PIV yields estimates for the vorticity that are dramatically improved over that in figure 11(f). Both PIV and PEV yield similar estimates for the vorticity at the channel walls (figure 11(e)).

5.2. Flow behind a circular cylinder

PEV is applied to a steady wake behind a circular cylinder. The circular cylinder was placed in a shallow water table and simultaneous measurement of the velocity field and the surface topography were made as described in Fouras et al [34]. The wake is steady with a Reynolds number of 30. Due to the depth of the water table, the flow is not directly comparable to the flow behind a fully submerged cylinder. Further details of the experimental procedure can be found in Fouras et al [34].

Figures 13(a) and (b) show the contours for the vorticity as estimated by PIV and PEV. The vectors, showing displacement, have been skipped for the purpose of clarity. For the PIV measurements in figure 13(a), only every 21st vector in x and every third vector in y is shown, while only every PEV measurement sampled at the same locations as the PIV measurements is shown for the PEV measurements in figure 13(b). PIV interrogation is carried out using the same scheme used on the synthetic flow past a square cylinder, as described in detail in section 4.3. PEV analysis is carried out over 16 × 16 elements (i.e. 128 × 128 pixel sampling windows) spanning the 2048 × 2048 pixels image. For comparative purposes, PIV was conducted using equivalent sampling windows to PEV (i.e. 128 × 128 pixel sampling windows) with a sampling distance of 8 pixels.

**Figure 13.** Contours plots showing the in-plane vorticity obtained with standard PIV using the optimal interrogation settings (a); PEV (b); and PIV using equivalent sampling window sizes as PEV (c). Comparisons between the measurements from PIV ( $\protect\epsfbox{images/mst428424fun01.eps}$ ) and PEV (blue line) are made at two slices: A–A at x/D = 1 (d, e) and B–B at x/D = 5 (f, g).
Download figure:
Standard image

Measurements of the flow displacement and the estimates of in-plane vorticity using standard PIV and the proposed PEV method are compared at two locations in the wake of the cylinder. Section A–A (figures 13(c) and (d)) is the location just behind the circular cylinder where there are high velocity gradients, and section B–B (figures 13(e) and (f)) is a location downstream in the wake with weaker velocity gradients. These regions were chosen to allow comparison of the performance of PIV and PEV with different levels of velocity gradients (due to the different velocity gradients that both methods need to compensate for in order to determine accurate measurements). The results shown in figure 13(c) demonstrate that both PIV and PEV are capable of capturing the sharp flow gradients in the cylinder wake. However, the vorticity estimates from PIV are noisy due to higher random error transmission from the displacement measurements as evident in figures 13(d) and (f). The random error transmission could be greatly reduced by using larger sampling windows for PIV interrogation. However, there would be insufficient resolution to fully capture structures in the flow. While reducing the sampling distance increases the number of velocity measures, the estimated vorticity field is far from ideal (figure 13(c)). While the measurements are smoother, i.e. reduced random error in the measurements, the vorticity measures have bias errors as evident in figures 13(e) and (g).

PEV yields a continuous estimate of the vorticity at both locations downstream of the wake with slight discontinuities at the element interfaces. These discontinuities arise from the regularization procedure that ensures global continuity across the entire solution. The elements in PEV can be overlapped in a similar fashion to that done as standard with PIV. Figure 14 shows how the initial PEV estimate from figure 13(f) can be improved greatly by conducting PEV with 50% element overlap.

**Figure 14.** Line plot of the estimated vorticity, at x/D = 5 (see figure 13), for different levels of element overlap. The blue line depicts the standard PEV solution, while the green and red lines depict the PEV solution curves with element overlapping of 50% and 75%, respectively.
Download figure:
Standard image

6. Conclusions

A technique to determine continuous velocity and velocity gradient measurements with low noise is presented. This process, PEV, optimizes for the measured cross-correlation map obtained from image data to determine a series of polynomial functions that describe the underlying flow.

With the general PIV technique, the velocity measurement is obtained from the signal peak in the cross-correlation map that represents the modal velocity. Estimated flow gradients using PIV velocity measurements are affected by random errors and bias errors. The error in estimating the flow gradient is dependent on the spatial resolution of flow measurements considered for calculation and the accuracy of the flow measurement itself. A Monte Carlo simulation was conducted in silico to study the effect of the PIV sampling distance (resolution of measurements) and the sampling window size (accuracy of the velocity measurement) on the accuracy in estimating flow gradients. The errors are coupled and therefore cannot be minimized simultaneously resulting in a single optimal combination.

The PEV method has been applied to a single instantaneous measurement of four sets of flows; two sets of laboratory flows and two sets of computer generated synthetic flows. Time-resolved PEV measurements are available by applying the proposed method to subsequent image pairs in the image sequence. The proposed method has been shown to yield continuous measurements with low levels of measurement noise for the velocity and vorticity in comparison to PIV. In the cases presented in this study, the PEV method was shown to reduce the velocity measurement error by up to 45% and the vorticity estimation error by up to 77%.

The random errors that arise in the estimated flow gradient measure due to increased amounts of oversampling of PIV data are dramatically reduced with the proposed PEV method due to the use of continuous polynomial functions to represent the flow. Bias errors in the estimated flow gradients are reduced as PEV optimizes for the entire cross-correlation map, taking into account the shape or distortion of the peak due to underlying flow gradients.

Acknowledgments

AF and CRS acknowledge the insightful discussions with Stephen Dubsky. Support from the Australian Research Council under grants DP0987643 and DP110101498 are gratefully acknowledged. GJS is supported by a Faculty of Engineering Seed Funding Grant, the NCI National Computational Merit Allocation Scheme, the VLSCI Resource Allocation Scheme, and Australian Research Council Discovery Grant DP120100153. CRS is a recipient of an Australian Postgraduate Award.

Polynomial element velocimetry (PEV): a technique for continuous in-plane velocity and velocity gradient measurements for low Reynolds number flows

Article metrics

Submit

Permissions

Author e-mails

Author affiliations

Dates

Abstract