- Email: [email protected]
Adaptive PIV algorithm based on seeding density and velocity information
Flow Measurement and Instrumentation 51 (2016) 21–29
Contents lists available at ScienceDirect
Flow Measurement and Instrumentation journal homepage: www.elsevier.com/locate/flowmeasinst
Adaptive PIV algorithm based on seeding density and velocity information Kaikai Yu, Jinglei Xu n Jiangsu Province Key Laboratory of Aerospace Power System, Nanjing University of Aeronautics and Astronautics, Nanjing 21006, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 18 January 2016 Received in revised form 14 June 2016 Accepted 7 August 2016 Available online 8 August 2016
An adaptive particle image velocimetry (PIV) processing algorithm to increase local spatial resolution and obtain accurate results is presented. Adaptive sampling criterion is based on seeding density and the analysis of velocity information. The proposed methodology places more measurement points in the region where high seeding density exists or flow parameters vary drastically. This methodology is effective in PIV image processing, especially in non-optimal PIV experimental conditions. The working principle of the adaptive processing method in this paper is to generate spring force on the bias of the deviation between the ideal and actual distances of measurement points. Under artificial force, the sampling points are moved to desirable places to satisfy sampling density function. Sampling density function is important in processing and ingeniously conceived based on the combination of vorticity and velocity gradients. The viability of the method is elaborated through synthetic and experimental tests. The synthetic test verifies the effectiveness of the adaptive processing on the basis of flow information. The bias and random error can be reduced by 7.3% and 4.0%, respectively. The result of the experimental test shows that the sampling points and the interrogation window size can be arranged according to the spatial distribution of seeding density and the actual flow features. Compared with iterative window deformation method, spatial resolution is locally enhanced, and the robustness and reliability of the result are increased, particularly in poor seeding regions. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Particle image velocimetry Adaptive Sampling points Seeding density Velocity information Spring force
1. Introduction Particle image velocimetry (PIV) has become the universal measurement technique to obtain instantaneous velocity and its derivatives, study physical flow mechanism, and design the contour of vehicles. PIV images are analyzed by segmenting the images into small interrogation windows and using the correlation method based on fast Fourier transform to obtain the resulting velocity [1,2]. The most important parameters for the conventional PIV processing are the size of the interrogation windows, the window deformation strategy, and the overlap factor. The result is located in a discrete position, which depends on the PIV processing parameters set by the users (i.e., velocity information can be obtained only at a fixed instance with finite resolution). This strategy, which commonly uses a regular mesh grid in the processing method, is irrespective of the underlying signal [3]. Almost all the flow causes variations of local velocity and inhomogeneity of seeding density. Hence, the traditional PIV algorithm using the
n
Corresponding author. E-mail address: [email protected] (J. Xu).
http://dx.doi.org/10.1016/j.flowmeasinst.2016.08.004 0955-5986/& 2016 Elsevier Ltd. All rights reserved.
uniform interrogation window cannot satisfy the actual requirement completely. PIV has undergone considerable progress on the analysis algorithm over more than a decade, and many researchers have made efforts to improve the performance of the processing algorithm. One of the typical improvements, the iterative strategy, has been implemented in the PIV algorithm. Both the accuracy of the measurement and spatial resolution have been improved [4]. This method has been commonly adopted because of its significantly augmented measurement capabilities and high spatial resolution. Moreover, extensive research has been conducted on the iterative image deformation method to improve the performance of the algorithm. Astarita [5,6] studied the influence of interpolation schemes in the image deformation method and provided useful suggestions on selecting the interpolation scheme. Schrijer [7] researched the stability of the iterative process and concluded that a well-designed spatial filter would improve process stability and increase the reliability of the velocity result. The iterative method improves the quality of the PIV result, but the spatial resolution remains limited because the uniform sampling method ignored the variations of seeding density and velocity dataset. Hence, two approaches can be applied to obtain high resolution in the processing algorithm: Particle Tracking velocimetry (PTV) and its
22
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
developments, and the adaptive method for making PIV correlation on the basis of seeding density and velocity information. The disadvantages of the traditional PIV processing method are caused by the adopted uniform spatial interrogation windows, which resolve velocity. Therefore, the actual situation of particle seeding cannot be considered. PTV has the potential to obtain a much higher spatial resolution because the velocity is based on the seeding particle although it uses a similar method to acquire the experimental images as PIV. However, PTV has inherent disadvantages, such as low particle density; by contrast, particle density is much higher in PIV. A reliable heuristic method called super PIV has been proposed to solve this problem [8–10]. The method has been successfully applied in actual experiments, thereby demonstrating that the method can increase spatial resolution. However, the performance of this method has been limited by the predictor algorithm, particle identification process, and its localization algorithm. Further development along this direction was made by Westerweel [11] by applying a single-pixel ensemble correlation and has been extended further by Billy et al. [12] and Kähler et al. [13]. Recently, the method was applied to compressible flows with high Mach numbers [14]. This innovative method has also increased the spatial resolution significantly for many microscopic applications when the flow is mostly (quasi-) stationary. Another approach is the adaptive method. Rohaly et al. [15] proposed a hierarchical processing scheme, which increased the correlation area gradually until the signal-to-noise ratio satisfied the requirement. This scheme is the basis of the adaptive method that considers seeding density only. A universal method is to apply the 2D Gaussian function on the basis of the velocity information to the interrogation area to reshape the window [16,17]. Becker [18] followed this idea but designed the Gaussian function on the basis of an error model function. The adaptive method is further implemented by resampling the measurement points on the basis of seeding density and velocity information [19–21]. Theunissen et al. [19,20] relocated measurement points by using the 2D PDF method, whereas Yu et al. [21] used artificial force to adjust the locations of the measurement points. Adaptive PIV algorithm is a research topic that has been attracting considerable interest. The concept of adaptive processing method has been recently applied in the tomography PIV by Matteo Novara [22]. The reliability of the adaptive method was verified when more sampling points were dynamically redistributed to better resolve the region where more spatial resolution can be obtained or need to be processed. So an adaptive PIV processing algorithm to increase local spatial resolution and obtain accurate results is presented in this paper. The proposed method utilizes a different approach with high performance to implement the adaptive procedure. The concept of the implementation is based on the authors’ previous work [21]. But the seeding density and velocity information are combined to construct the sampling density function rather than based only on the seeding density in the previous work. The sampling points change gradually under the influence of the artificial spring force until they satisfy the requirement of the local sampling density. Using the proposed sampling method, more details of the flow can be resolved. The first section of this paper describes the detail and the procedure of the adaptive method. The sampling density of the velocity is also introduced in detail. The performance, reliability, and robustness of the proposed method are verified through processing synthetic and experimental images. The advantages of the proposed adaptive sampling method are also presented.
2. Adaptive algorithm 2.1. Construction of adaptive sampling criterion The adaptive algorithm described in this paper samples the measurement points and is designed as a self-contained subroutine. The method can be added to the existing PIV processing algorithm with minimal modification, which is convenient for users to utilize and make some improvements. The intrinsic processing algorithm can be implemented to obtain the final result once the positions of sampling points are obtained. High quality and resolution results can be obtained by concentrating sampling points on the positions where they are most needed at the expense of a reasonable computational cost. Hence, the two main components of the adaptive processing method are an optimal sampling criterion and a strategy to redistribute the points. From this perspective, the adaptive criteria should be designed carefully. The adaptive processing method is described in this paper. It can relocate the sampling points by adapting to the flow features and the particle density. Thus, the sampling criterion must be organized through a special method. The sampling criterion cannot be organized clearly given that researchers do not always know a priori what constitutes the complex flow structures of the problem. The computational errors of the cross-correlation first come to mind because cross-correlation is the principal method to process the PIV images. However, computational errors remain unclear and cannot instruct the sampling procedure directly. The computational errors of cross-correlation tend to increase in regions where seeding density is low or velocity changes rapidly. Thus, seeding density must be the main criterion to instruct the sampling procedure. Yu et al. [21] chose seeding density as the first candidate in their work to implement the sampling procedure, which demonstrated the reliability and robustness of the proposed adaptive method based on spring force model. The details are given in our previous work [21]. Velocity and its derivatives are considered in the construction of the sampling criterion in this paper. Velocity must be adapted after the first sampling procedure based on the seeding density because velocity information remains unclear before sampling the points. The adaptive sampling procedure based on the seeding density is the first step. The sampling procedure based on the velocity information is the second step. The adaptive mesh is used in computational fluid dynamics (CFD) to finely resolve flow feature. A choice of characteristic flow parameters can be detected in CFD, such as velocity, density, pressure, and temperature. Hence, the solution-based adaption method adopted in CFD is usually based on solution features, such as gradients, curvature, or specific features of the flow field (e.g., shock waves and expansion waves). The velocity and its derivatives should be the best candidates for the adaptation given that PIV algorithm obtains only the velocity results. Moreover, a combination of suitable parameters is necessary. Therefore, velocity gradients should be used to detect viscous regions and vorticity to detect shear layers. A new distribution of sampling points can be generated based on the solution obtained from the first step of seeding density-based sampling distribution. This distribution is more suitable for capturing flow features. Hence, the key point of the adaptation procedure is the construction of the sampling density specification function, which specifies the distances between sampling points through the experimental region. The distributions of sampling points can be generated by the influence of the spring force once the function has been established. Prior knowledge of the adaptive mesh in CFD presents that the combination of the flow feature should be considered. Otherwise, the sampling criterion may cause the sampling points to become excessively concentrated or sparse. Hence, an optional sampling
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
23
criterion can be constructed according to the choice of Theunissen [19] to use the standard deviation of the local variations of the velocity.
1 W2
σ k ( x , y) =
i = W /2
j = W /2
∑
∑
( uk ( x + i, y + j) − uk )2 , k = 1, 2
i =−W /2 j =−W /2
(1)
where x and y are the local coordinates and k are the index of the U and V component velocity. uk and uk are the local velocity and the mean velocity of the specified area W, respectively. A combination of the standard deviation is utilized in general.
σ=
1 2
2
∑ σk2
(2)
k=1
This sampling adaptive criterion has been constructed and applied in the experimental processing procedure by Theunissen [19]. Its robustness and reliability have also been verified. A different sampling criterion is adopted and constructed in this paper. Two different parameters which have been found to work well are related to the local gradients and the curl of the velocity.
λ1 (i, j ) = ‖ curl u‖
(3)
2
λ2 (i, j ) = ‖gradient u ‖
u =
∑ uk2
(4)
k=1
where u is the velocity vector and k are the index of the U and V component velocity. The reasons for this choice are the gradients of velocity can measure the local change of the flow and sense shock waves effectively and the curl criterion can use the local rotationality of the flow and performs well in locating shear layers. The standard deviations of the two parameters are computed, and the ratios of the local detected parameter to the relative standard deviation are used to characterize the refinement region.
ζk (i, j ) =
λ k ( i, j ) 1 XY
i=X ∑i = 1
j=Y ∑ j=1
(5)
where i and j are the coordinates and k are the index of the different parameters. X and Y are the matrix sizes of the parameters. The parameters are normalized by intermediate values to distinguish the relative value of the flow feature parameters. The function is defined as follows:
ξk (i, j ) = ζk (i, j )/median (ζk (i, j ))
(6)
After the relative value of the entire region is computed, additional locations where the sampling points should be placed can be determined. Given that the initial distribution of the sampling points has been determined, the function cannot be applied directly and must be transformed to another form. The above functions shows that refinement should be conducted when ξk is larger, which means the distances of the sampling points decrease as the value of ξk increases, and vice versa. Under this requirement, a specified function g should be defined as:
ρk (i, j ) = g (ξk ) k = 1, 2
(7)
The range of the sampling density is limited, which is from 0.5 to 1.5. When ξk → 0, g is 1.5. When ξk → ∞, g is 0.5. Moreover, g is equal 1 when ξk = 1. Hence, a typical function has been constructed.
⎧ − 0.5ξk 0.88 + 1.5 0 < ξk < 1 k = 1, 2 ρk = g (ξk ) = ⎨ ⎩ 0.4ξk−1.10 + 0.6 1 ≤ ξk ⎪
⎪
Fig. 1 displays the distribution of the sampling density function. This figure illustrates that the sampling ratio decreases as the value of ξk increases. Moreover, the max value of the sampling ratio is 1.5, and the minim one is 0.5. This setting can avoid excessive refinement or coarsening of the sampling points. From the foregoing section, two detected flow features can be selected as the adaptive criterion, which causes different distributions of the sampling points. A combined function is used in this paper to make better use of the two features.
(
ρfinal = min ρgradient , ρcurlz
)
(9)
The sampling density of the velocity information has been determined. The main procedure of distributing the sampling points and its applications in the processing synthetic and experimental images are introduced in the next subsections. 2.2. Adaptive sampling procedure
k = 1, 2
(λ k ( i, j ) − λ k )2
Fig. 1. Distribution of the sampling density as a function of ξk .
(8)
Measurement point sampling can be implemented once the sampling criterion is settled (including the seeding density and the criterion of the velocity). This subsection details the procedure of distributing the sampling points. The sampling criterion for the velocity and the seeding density can be obtained according to the method described in Section 2.1. The main procedure of the method is similar to that of Yu et al. [21]. The concept of artificial spring force is utilized. The initial distances of the sampling points are determined after initializing the distributions of the sampling points. The deviations of the initial and the desired distances of the points produce spring force. Spring force pushes the sampling points to the ideal place under the iterative procedure. The concept is first proposed by Persson and Strang [23] in the generation of CFD mesh with high quality. The mesh generated by the spring mode method is verified to be higher than that filtered by the Laplacian method. The detailed procedure is illustrated briefly as follows. 1. The number of the interrogation windows is set according to the number of identified particle images (N), the number of particles within the interrogation window (NI), and the mean overlap value (OF). Hence, the number of the interrogation windows can be defined as
Nw =
N (1 − OF )2NI
(10)
24
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
in Subsection 2.1. 5. The seeding density and the sampling density are attained. The velocity information causes the desired distance to change although the sampling points meet the requirement of the seeding density in the first sampling step. The sampling points move to the new iteration place after the iteration because of the spring force. The final positions and the number of the sampling points are then determined. Then, the final velocity is obtained by using the interrogation method. The procedure of the method can be summarized as follows: The seeding density is computed. The points are placed according to seeding density. The flow feature is computed, and the sampling density function is constructed. The points are placed according to the sampling density based on the position obtained from the first sampling step. The final velocity is computed. Fig. 3 schematically displays the procedure of the adaptive processing method. 2.3. Interrogation method Fig. 2. Typical 2D Delaunay triangulation and generation of the artificial spring force.
The sampling points can then be distributed randomly. The random algorithm in the paper adopts the method proposed by Kocis and Whiten [24]. 2. The sampling points are triangulated by Delaunay triangulation. Fig. 2 displays the typical 2D Delaunay triangulation. The neighboring sampling points of the specified point Pi can be recognized once the triangulation is generated. The desired distance between Pi and Pj, which is associated with the local seeding density, the deviation of the original, and the desired distances, can be obtained. The force is then generated.
⎧ if hj > L ij ⎪ k ( hj − L ij ) Fij = ⎨ ⎪ if hj ≤ L ij ⎩0
(11)
where Fij is the repulsive force to the node Pi from Pj, and hj is the ideal distance between Pi and Pj. The ideal distance is deduced from the seeding density in the first sampling step and from the combination of the distances from the first step and the sampling density in the second sampling step. Lij is the original distance between Pi and Pj. The resultant force of each point is the addition of Fij. Sd is seeding density (in particles per pixel 2)).
hk, j
⎧ NI ⎪ (1 − OF ) k=1 =⎨ Sd ⎪ k=2 ⎩ ρ ⋅L ij
(12)
3. After finishing the computation of the resultant force, the equilibrium positions where all the measurement points meet the requirement of seeding density and the sampling ratio are solved from the system F(p) ¼0. p is the location of sampling points. The details of the solution of the system are not introduced because they are beyond the scope of this work. Under the sophistic hypothesis, the new positions of the sampling points can be obtained by using the forward Euler method .
( )
pin + 1 = pin + Δt⋅F pin
(13)
where n is the iterative number. The sampling points move gradually to the equilibrium positions according to the equation. 4. After approximately 10–20 iterations, the sampling points in the first sampling step arrive at the equilibrium positions. On the basis of these points, the velocity can be obtained by using the cross-correlation method and the sampling density introduced
Fig. 3 schematically illustrates the procedure of the adaptive PIV processing method. The displacements of each sampling points are obtained by cross-correlation based on fast Fourier transform. An iterative procedure is adopted (i.e., iterative window deformation method [4]). Contrast-limited adaptive histogram equalization [25] is utilized before the core processing to enhance the contrast of the original PIV images and increase the reliability of the final result. A high-pass image enhancement method is applied in the computation of the seeding density to decrease the computational cost. For comparison, the traditional method utilized in the following two verification cases adopts the same way to obtain the displacements of the interrogation windows, including iterative window deformation method, the image preprocessing method and the verification method for the final dataset. And the final interrogation areas in iterative window deformation method are determined according to the actual conditions. The predictor needed by the construction of the sampling criterion is acquired by one-pass cross-correlation after the first sampling based on seeding density. Given that the sampling points are not uniformly distributed spatially, an interpolation scheme must be selected to obtain the gridded data. The velocity results are interpolated over the regular grid from the discrete sampling points using natural neighbor interpolation in this paper [26]. The sampling criterion is then computed and filtered by the method proposed by Garcia [27] to increase the stability of the algorithm. The sampling points are relocated using artificial spring force once the final sampling criterion is obtained. Finally, the velocity on the sampling point is computed by using the window deformation method. Moreover, the interrogation areas are decreased to the specified value iteratively. The final dataset is verified by the global standard derivation and the universal outlier detection method proposed by Westerweel [28]. For comparison with the traditional processing method, the scattered data can be interpolated to the uniform grid by natural neighbor interpolation.
3. Results and discussion 3.1. Performance evaluation of adaptive algorithms in synthetic images 3.1.1. Generation of synthetic images In this subsection, the proposed method is applied to process synthetic images. The results from the proposed method are compared with those obtained from the traditional processing
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
25
Fig. 3. Main procedure of the adaptive processing method.
image. The vectors obtained from the traditional processing method (i.e., uniform sampling and window deformation method) are also shown in Fig. 4, where the Rankin vortex clearly appears in the middle of the images.
Fig. 4. Synthetic image and vector result from the traditional processing method.
method. The synthetic images are generated according to Raffel et al. [29]. The influence of inhomogeneous seeding density in synthetic application is not verified again in this paper. The uniform seeding density is adopted in the verification case, namely, the Rankine vortex, to focus on the performance of the sampling density rate of the velocity. Rankine vortex is a circular flow. The inner part of the vortex is in solid rotation, whereas the outer part is free of vorticity. The distribution of the velocity is defined as follows:
⎧ V ⎪ ⎪ max Vθ ( r ) = ⎨ ⎪V ⎪ ⎩ max
r R R r
3.1.2. Processing of synthetic images The particle image detection algorithm identified approximately 18,221 particles. The average Sd is approximately 0.038/pixels2. Adrian's suggestions indicate that the number of particles within the interrogation windows is set as eight to ensure the robustness of the cross-correlation algorithm [2]. The value of overlap factor is set as 50%. The number of interrogation windows can be computed as 9110 using Eq. (10). The procedure illustrated in Fig. 3 shows that the initial sampling step is conducted. The first distribution of the sampling points instructed by the seeding density can be attained after 10 iterations. The distributions of sampling points in this procedure are also uniform because of the uniform seeding in the synthetic images, see Fig. 5. The gradient and the curl of the velocity can be computed based on the
0≤r≤R r>R
(14)
where Vmax is the max velocity of the vortex and is set as 5 pixels in this verification case, and R is the characteristic distance of the vortex (e.g., the core radius of the vortex) and is set as 80 pixels. The image size is 800 pixels 600 pixels. The vortex center is located in the center of the image. Fig. 4 displays the synthetic
Fig. 5. Distributions of the sampling point based on seeding density.
26
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
Fig. 6. Contour of the sampling density (i.e., the distribution of the ideal distances between the sampling points); the distributions of the sampling points based on the combination of the seeding and velocity information.
preliminary velocity result. The sampling density can be derived accordingly using the above equations. In other words, the distribution of the ideal relative pixel distances of the sampling point can be obtained. Fig. 6 displays the contour of the sampling density, where the value represents the relative distances between the sampling points. The intermediate region should have more sampling points because of the existence of the vortex. Therefore, more measurement points are sampled in Regions A and B. Moreover, the distances between the sampling points in Regions A and B are smaller than the others, and the size of the interrogation windows is also decreased accordingly. Iterative window deformation strategy is implemented after distributing the sampling points to process the images and obtain the final result. Processing the PIV images by using the traditional method requires an appropriate choice of final interrogation area, considering the reliability of the dataset and spatial resolution. A constant sample area of 18 18 pixels2 with a 50% overlap was applied in this case. In most cases of evaluating the performance of the algorithm of PIV, bias error and random error are adopted for the evaluation [30]. The bias error is quantified by the deviation between the theory and the measured result. The random error level is usually quantified as the root-mean-square (RMS) fluctuation of the measured values around the theory value. Here, RMS is regarded as the random error. The formula is defined as follows:
σ=
1 Ns
Ns
∑ (X (i)−XT (i))2 i=1
(15)
where Ns is the total number of the sampling points, X is the local measured velocity, and XT is the local theory result. The two different errors are computed from the U and V components. The average bias error of the V component from the adaptive processing is 0.056 pixels, whereas the one from the traditional method is 0.0604 pixels. The same tendency can be obtained from the U component. The bias error of the adaptive processing method is 7.28% less than that of the traditional method. With regard to the random error, the difference between the two methods becomes relatively small. The error from the adaptive one is lower than that from the traditional method by approximately 4%. At the same time, the contours of the vorticity from the two different methods are computed (Fig. 7). Fig. 7(a) presents the result of the traditional processing method. Fig. 7(b) shows the result of the adaptive method. The theoretical value can be computed based on Eq. (14), the vorticity value in the middle region should be 0.125 pixels/pixels. From Fig. 7, it can be found that the vorticity of the middle region in Fig. 7(b) is more uniform and closer to the theoretical value than that in Fig. 7(a). That also means that the vorticity from the traditional method is fairly disturbed and significantly non-uniform. Besides that, the boundary between the inner and outer flow of Rankine vortex from the adaptive method is much sharper than that from the traditional method. Then, the contour line in Fig. 7(b) is smoother and it should be more accurate to represent the Rankine vortex. Most obviously, it can be observed from Fig. 7 that the fluctuations in the traditional case are much larger than those in the adaptive case. The same phenomena can be verified by the distribution of the vorticity and velocity magnitude in the middle horizontal line from the two cases (Fig. 8). Previous studies showed that the result of cross-correlation could be modeled with a moving average linear filter [31]. Thus, the theory vorticity in Fig. 8(b) is filtered by a moving average filter with a 16-pixels window length for intuitive comparison. The analysis of the bias and random error and the distribution of the vorticity indicate that the adaptive processing method can increase the performance of the PIV algorithm, the reliability of the dataset, and spatial resolution. 3.2. Performance evaluation of adaptive algorithms in experimental images In this subsection, experimental PIV images are processed to demonstrate the good performance of the adaptive processing method qualitatively. Experiments on jet flow are measured by PIV provided by the PIV Challenge 2005 (http://www.pivchallenge. org). The traditional processing method and the adaptive method are applied to obtain the final result for an intuitive comparison.
Fig. 7. (a). Contour of the vorticity from the traditional processing method; (b). Contour of the vorticity from the adaptive processing method.
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
27
Fig. 8. (a). Distributions of the vorticity in the middle horizontal line; (b). Distributions of the velocity magnitude in the middle horizontal line. Red represents the theoretical result, blue represents the traditional result, and green represents the adaptive result.
in the adaptive processing algorithm based on the previously mentioned equations is obtained as 4892. The overlap value is 50% in both methods. The NI is eight according to the suggestions by Adrian [2]. After determining the necessary parameters for the adaptive processing (i.e., the Sd, the sampling density function, and the initial number of sampling points), the sampling points based on the seeding density can be distributed accordingly (Fig. 9). The reliability of the adaptive sampling method based on the seeding density is verified again. More sampling points are placed in the core flow region, and vice versa. The sampling density is then constructed by using the equations combined with the initial velocity information. The final distribution of the measurement points can be obtained (Fig. 10). The white points in Fig. 10 illustrate the location of
Fig. 9. PIV image. The gray level of the image is inversed, and contrast enhancement processing is applied. Blue points represent the sampling points, which are distributed based on the seeding density. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The flow of nitrogen in this experimental application has a velocity of 30 m/s. The viewing area is 29 29 mm2. The resolution of the PIV raw image is 512 512 pixels. The seeding density is fairly inhomogeneous because of inhomogeneous seeding of the flow and the co-flow. With the use of the particle detection algorithm, the number of the detected particle images is approximately 9785, and the seeding density can be obtained as 0.0373 particles/pixels2. Fig. 9 presents the raw PIV images and shows that the Sd in the main flow region is much higher than that in the external flow region. The gray level of the raw image has been inversed. Contrast-limited adaptive histogram equalization method has also been applied to enhance the contrast of the image and the performance of the algorithm. A trade-off choice must be made to select the interrogation window size appropriately. Hence, the final interrogation area under this requirement is set as 32 pixels because of the low Sd in the external flow region in the traditional processing algorithm. The number of sampling points
Fig. 10. Distribution of the sampling points based on the seeding density and the velocity information. White points represent the sampling points. The corresponding vectors are also displayed.
28
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
Fig. 11. (a). Contour of the vorticity from the traditional processing method; (b). Contour of the vorticity from the adaptive processing method.
the measurement points. Figs. 9 and 10 distinguish the significant differences. Fig. 10 samples more points in the boundary of the main flow because of the high vorticity and gradient. Some measurement points were sampled in the external flow under the requirement of sampling density function. Thus, the measurement requirement is more accurately met. The velocity and its derivations are obtained by using iterative window deformation method. Fig. 11 displays the vorticity from the traditional and adaptive method. Due to the augmented measurement points sampled in the main flow, the adaptive method is able to capture more details about the main flow. And, the adaptive method suppresses noise, especially in the left flow region. From the previous paragraph, it can be known that the seeding density in the outer region of the jet is smaller than that in the main flow. Under the limitation of the traditional method, the fixed interrogation areas have been determined. Thus, in the processing of the outer region, the number of particle images in the interrogation domain is decreased, which causes the amplitude of the displacement-correlation peak reduced. Failure to correctly identify these peaks leads to spurious vectors. So some fluctuations exist in the result obtained by the traditional method. On the contrary, adaptive method places few sampling point in the outer region and utilizes large interrogation windows, which meet the requirement of the cross-correlation algorithm and make that the velocity results become reliable. In fact, the conclusion obtained from the synthetic case also can verify this phenomenon. The random and bias errors from the adaptive method are less than that from the traditional method. That is, as to zero displacement, more fluctuations will be shown in the result obtained by the traditional method. In a word, the experimental application can verify the effectiveness and reliability of the proposed adaptive algorithm for PIV processing.
sampling point to the local seeding density and the velocity information unlike the use of the regular mesh in the traditional method. The concept of spring force adopted in previous research is utilized in this study to distribute the sampling points. The proposed method places the sampling points on the basis of the seeding density in the first procedure. After that, a specified sampling density function is constructed on the bias of the vorticity and gradient of the velocity. The final distribution of the points can then be obtained. Two test cases have been chosen to attest the effectiveness of the proposed method: a synthetic test (Rankine vortex) and an experimental test (jet flow). The synthetic example uses the uniform seeding density to verify the operability of the proposed method. The final velocity result indicates that the proposed method can increase the performance of the PIV algorithm by 7% in terms of the bias error and 4% in terms of the random error. The vorticity is smooth in the adaptive processing method, which agrees well with the theory. The actual images are then processed; results verify that the proposed method is reliable and effective. Compared with iterative window deformation method, spatial resolution is locally enhanced, and the dataset is more reliable. Further efforts will be directed toward the consideration of temporal information in PIV processing algorithm.
Acknowledgments We would like to acknowledge the continued support of the NSFC (Natural Science Fund of China) for contract numbers 50876042 and 90916023.
References 4. Conclusions An adaptive PIV processing algorithm is proposed to enhance the performance of the processing algorithm and increase spatial resolution. The proposed method adapts the distributions of the
[1] C.E. Willert, M. Gharib, Digital particle image velocimetry, Exp. Fluids 10 (4) (1991) 181–193. [2] R.J. Adrian, J. Westerweel, Particle Image Velocimetry, Cambridge University Press,, 2011. [3] R.D. Keane, R.J. Adrian, Theory of cross-correlation analysis of PIV images, Appl. Sci. Res. 49 (3) (1992) 191–215.
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
[4] F. Scarano, Iterative image deformation methods in PIV, Meas. Sci. Technol. 13 (1) (2002) R1. [5] T. Astarita, G. Cardone, Analysis of interpolation schemes for image deformation methods in PIV, Exp. Fluids 38 (2) (2005) 233–243. [6] T. Astarita, Analysis of interpolation schemes for image deformation methods in PIV: effect of noise on the accuracy and spatial resolution, Exp. Fluids 40 (6) (2006) 977–987. [7] F.F.J. Schrijer, F. Scarano, Effect of predictor–corrector filtering on the stability and spatial resolution of iterative PIV interrogation, Exp. Fluids 45 (5) (2008) 927–941. [8] R.D. Keane, R.J. Adrian, Y. Zhang, Super-resolution particle imaging velocimetry, Meas. Sci. Technol. 6 (6) (1995) 754. [9] A. Stitou, M.L. Riethmuller, Extension of PIV to super resolution using PTV, Meas. Sci. Technol. 12 (9) (2001) 1398. [10] A. Susset, J.M. Most, D. Honoré, A novel architecture for a super-resolution PIV algorithm developed for the improvement of the resolution of large velocity gradient measurements, Exp. Fluids 40 (1) (2006) 70–79. [11] J. Westerweel, P.F. Geelhoed, R. Lindken, Single-pixel resolution ensemble correlation for micro-PIV applications, Exp. Fluids 37 (3) (2004) 375–384. [12] F. Billy, L. David, G. Pineau, Single pixel resolution correlation applied to unsteady flow measurements[J], Meas. Sci. Technol. 15 (6) (2004) 1039. [13] C.J. Kähler, U. Scholz, J. Ortmanns, Wall-shear-stress and near-wall turbulence measurements up to single pixel resolution by means of long-distance microPIV, Exp. Fluids 41 (2) (2006) 327–341. [14] M. Bitter, S. Scharnowski, R. Hain, C.J. Kähler, High-repetition-rate PIV investigations on a generic rocket model in sub-and supersonic flows, Exp. Fluids 50 (4) (2011) 1019–1030. [15] J. Rohaly, F. Frigerio, D.P. Hart, Reverse hierarchical PIV processing, Meas. Sci. Technol. 13 (7) (2002) 984. [16] D. Di Florio, F. Di Felice, G.P. Romano, Windowing, re-shaping and re-orientation interrogation windows in particle image velocimetry for the investigation of shear flows, Meas. Sci. Technol. 13 (7) (2002) 953. [17] B. Wieneke, K. Pfeiffer, Adaptive PIV with variable interrogation window size
[18]
[19] [20] [21] [22] [23] [24] [25] [26]
[27] [28] [29] [30]
[31]
29
and shape. 15th Int Symp on Applications of Laser Techniques to Fluid Mechanics, Lisbon, Portugal, 2010. F. Becker, B. Wieneke, S. Petra, A. Schröder, C. Schnörr, Variational adaptive correlation method for flow estimation, Image Process. IEEE Trans. 21 (6) (2012) 3053–3065. R. Theunissen, F. Scarano, M.L. Riethmuller, An adaptive sampling and windowing interrogation method in PIV, Meas. Sci. Technol. 18 (1) (2007) 275. R. Theunissen, F. Scarano, M.L. Riethmuller, Spatially adaptive PIV interrogation based on data ensemble, Exp. Fluids 48 (5) (2010) 875–887. K. Yu, J. Xu, L. Tang, J. Mo, New adaptive sampling method in particle image velocimetry, Meas. Sci. Technol. 26 (3) (2015) 037002. M. Novara, A. Ianiro, F. Scarano, Adaptive interrogation for 3D-PIV, Meas. Sci. Technol. 24 (2) (2013) 024012. P.O. Persson, G. Strang, A simple mesh generator in MATLAB, SIAM Rev. 46 (2) (2004) 329–345. L. Kocis, W.J. Whiten, Computational investigations of low-discrepancy sequences, ACM Trans. Math. Softw. (TOMS) 23 (2) (1997) 266–294. K. Zuiderveld, Contrast Limited Adaptive Histogram Equalization (In Graphics gems IV), Academic Press Professional, Inc 1994, pp. 474–485. M. Sambridge, J. Braun, H. McQueen, Geophysical parametrization and interpolation of irregular data using natural neighbours, Geophys. J. Int. 122 (3) (1995) 837–857. D. Garcia, Robust smoothing of gridded data in one and higher dimensions with missing values, Comput. Stat. Data Anal. 54 (4) (2010) 1167–1178. J. Westerweel, F. Scarano, Universal outlier detection for PIV data, Exp. Fluids 39 (6) (2005) 1096–1100. M. Raffel, C.E. Willert, J. Kompenhans, Particle Image Velocimetry: A Practical Guide, Springer, 2007. L. Gui, S.T. Wereley, A correlation-based continuous window-shift technique to reduce the peak-locking effect in digital PIV image evaluation, Exp. Fluids 32 (4) (2002) 506–517. F. Scarano, Theory of non-isotropic spatial resolution in PIV, Exp. Fluids 35 (3) (2003) 268–277.
Contents lists available at ScienceDirect
Flow Measurement and Instrumentation journal homepage: www.elsevier.com/locate/flowmeasinst
Adaptive PIV algorithm based on seeding density and velocity information Kaikai Yu, Jinglei Xu n Jiangsu Province Key Laboratory of Aerospace Power System, Nanjing University of Aeronautics and Astronautics, Nanjing 21006, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 18 January 2016 Received in revised form 14 June 2016 Accepted 7 August 2016 Available online 8 August 2016
An adaptive particle image velocimetry (PIV) processing algorithm to increase local spatial resolution and obtain accurate results is presented. Adaptive sampling criterion is based on seeding density and the analysis of velocity information. The proposed methodology places more measurement points in the region where high seeding density exists or flow parameters vary drastically. This methodology is effective in PIV image processing, especially in non-optimal PIV experimental conditions. The working principle of the adaptive processing method in this paper is to generate spring force on the bias of the deviation between the ideal and actual distances of measurement points. Under artificial force, the sampling points are moved to desirable places to satisfy sampling density function. Sampling density function is important in processing and ingeniously conceived based on the combination of vorticity and velocity gradients. The viability of the method is elaborated through synthetic and experimental tests. The synthetic test verifies the effectiveness of the adaptive processing on the basis of flow information. The bias and random error can be reduced by 7.3% and 4.0%, respectively. The result of the experimental test shows that the sampling points and the interrogation window size can be arranged according to the spatial distribution of seeding density and the actual flow features. Compared with iterative window deformation method, spatial resolution is locally enhanced, and the robustness and reliability of the result are increased, particularly in poor seeding regions. & 2016 Elsevier Ltd. All rights reserved.
Keywords: Particle image velocimetry Adaptive Sampling points Seeding density Velocity information Spring force
1. Introduction Particle image velocimetry (PIV) has become the universal measurement technique to obtain instantaneous velocity and its derivatives, study physical flow mechanism, and design the contour of vehicles. PIV images are analyzed by segmenting the images into small interrogation windows and using the correlation method based on fast Fourier transform to obtain the resulting velocity [1,2]. The most important parameters for the conventional PIV processing are the size of the interrogation windows, the window deformation strategy, and the overlap factor. The result is located in a discrete position, which depends on the PIV processing parameters set by the users (i.e., velocity information can be obtained only at a fixed instance with finite resolution). This strategy, which commonly uses a regular mesh grid in the processing method, is irrespective of the underlying signal [3]. Almost all the flow causes variations of local velocity and inhomogeneity of seeding density. Hence, the traditional PIV algorithm using the
n
Corresponding author. E-mail address: [email protected] (J. Xu).
http://dx.doi.org/10.1016/j.flowmeasinst.2016.08.004 0955-5986/& 2016 Elsevier Ltd. All rights reserved.
uniform interrogation window cannot satisfy the actual requirement completely. PIV has undergone considerable progress on the analysis algorithm over more than a decade, and many researchers have made efforts to improve the performance of the processing algorithm. One of the typical improvements, the iterative strategy, has been implemented in the PIV algorithm. Both the accuracy of the measurement and spatial resolution have been improved [4]. This method has been commonly adopted because of its significantly augmented measurement capabilities and high spatial resolution. Moreover, extensive research has been conducted on the iterative image deformation method to improve the performance of the algorithm. Astarita [5,6] studied the influence of interpolation schemes in the image deformation method and provided useful suggestions on selecting the interpolation scheme. Schrijer [7] researched the stability of the iterative process and concluded that a well-designed spatial filter would improve process stability and increase the reliability of the velocity result. The iterative method improves the quality of the PIV result, but the spatial resolution remains limited because the uniform sampling method ignored the variations of seeding density and velocity dataset. Hence, two approaches can be applied to obtain high resolution in the processing algorithm: Particle Tracking velocimetry (PTV) and its
22
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
developments, and the adaptive method for making PIV correlation on the basis of seeding density and velocity information. The disadvantages of the traditional PIV processing method are caused by the adopted uniform spatial interrogation windows, which resolve velocity. Therefore, the actual situation of particle seeding cannot be considered. PTV has the potential to obtain a much higher spatial resolution because the velocity is based on the seeding particle although it uses a similar method to acquire the experimental images as PIV. However, PTV has inherent disadvantages, such as low particle density; by contrast, particle density is much higher in PIV. A reliable heuristic method called super PIV has been proposed to solve this problem [8–10]. The method has been successfully applied in actual experiments, thereby demonstrating that the method can increase spatial resolution. However, the performance of this method has been limited by the predictor algorithm, particle identification process, and its localization algorithm. Further development along this direction was made by Westerweel [11] by applying a single-pixel ensemble correlation and has been extended further by Billy et al. [12] and Kähler et al. [13]. Recently, the method was applied to compressible flows with high Mach numbers [14]. This innovative method has also increased the spatial resolution significantly for many microscopic applications when the flow is mostly (quasi-) stationary. Another approach is the adaptive method. Rohaly et al. [15] proposed a hierarchical processing scheme, which increased the correlation area gradually until the signal-to-noise ratio satisfied the requirement. This scheme is the basis of the adaptive method that considers seeding density only. A universal method is to apply the 2D Gaussian function on the basis of the velocity information to the interrogation area to reshape the window [16,17]. Becker [18] followed this idea but designed the Gaussian function on the basis of an error model function. The adaptive method is further implemented by resampling the measurement points on the basis of seeding density and velocity information [19–21]. Theunissen et al. [19,20] relocated measurement points by using the 2D PDF method, whereas Yu et al. [21] used artificial force to adjust the locations of the measurement points. Adaptive PIV algorithm is a research topic that has been attracting considerable interest. The concept of adaptive processing method has been recently applied in the tomography PIV by Matteo Novara [22]. The reliability of the adaptive method was verified when more sampling points were dynamically redistributed to better resolve the region where more spatial resolution can be obtained or need to be processed. So an adaptive PIV processing algorithm to increase local spatial resolution and obtain accurate results is presented in this paper. The proposed method utilizes a different approach with high performance to implement the adaptive procedure. The concept of the implementation is based on the authors’ previous work [21]. But the seeding density and velocity information are combined to construct the sampling density function rather than based only on the seeding density in the previous work. The sampling points change gradually under the influence of the artificial spring force until they satisfy the requirement of the local sampling density. Using the proposed sampling method, more details of the flow can be resolved. The first section of this paper describes the detail and the procedure of the adaptive method. The sampling density of the velocity is also introduced in detail. The performance, reliability, and robustness of the proposed method are verified through processing synthetic and experimental images. The advantages of the proposed adaptive sampling method are also presented.
2. Adaptive algorithm 2.1. Construction of adaptive sampling criterion The adaptive algorithm described in this paper samples the measurement points and is designed as a self-contained subroutine. The method can be added to the existing PIV processing algorithm with minimal modification, which is convenient for users to utilize and make some improvements. The intrinsic processing algorithm can be implemented to obtain the final result once the positions of sampling points are obtained. High quality and resolution results can be obtained by concentrating sampling points on the positions where they are most needed at the expense of a reasonable computational cost. Hence, the two main components of the adaptive processing method are an optimal sampling criterion and a strategy to redistribute the points. From this perspective, the adaptive criteria should be designed carefully. The adaptive processing method is described in this paper. It can relocate the sampling points by adapting to the flow features and the particle density. Thus, the sampling criterion must be organized through a special method. The sampling criterion cannot be organized clearly given that researchers do not always know a priori what constitutes the complex flow structures of the problem. The computational errors of the cross-correlation first come to mind because cross-correlation is the principal method to process the PIV images. However, computational errors remain unclear and cannot instruct the sampling procedure directly. The computational errors of cross-correlation tend to increase in regions where seeding density is low or velocity changes rapidly. Thus, seeding density must be the main criterion to instruct the sampling procedure. Yu et al. [21] chose seeding density as the first candidate in their work to implement the sampling procedure, which demonstrated the reliability and robustness of the proposed adaptive method based on spring force model. The details are given in our previous work [21]. Velocity and its derivatives are considered in the construction of the sampling criterion in this paper. Velocity must be adapted after the first sampling procedure based on the seeding density because velocity information remains unclear before sampling the points. The adaptive sampling procedure based on the seeding density is the first step. The sampling procedure based on the velocity information is the second step. The adaptive mesh is used in computational fluid dynamics (CFD) to finely resolve flow feature. A choice of characteristic flow parameters can be detected in CFD, such as velocity, density, pressure, and temperature. Hence, the solution-based adaption method adopted in CFD is usually based on solution features, such as gradients, curvature, or specific features of the flow field (e.g., shock waves and expansion waves). The velocity and its derivatives should be the best candidates for the adaptation given that PIV algorithm obtains only the velocity results. Moreover, a combination of suitable parameters is necessary. Therefore, velocity gradients should be used to detect viscous regions and vorticity to detect shear layers. A new distribution of sampling points can be generated based on the solution obtained from the first step of seeding density-based sampling distribution. This distribution is more suitable for capturing flow features. Hence, the key point of the adaptation procedure is the construction of the sampling density specification function, which specifies the distances between sampling points through the experimental region. The distributions of sampling points can be generated by the influence of the spring force once the function has been established. Prior knowledge of the adaptive mesh in CFD presents that the combination of the flow feature should be considered. Otherwise, the sampling criterion may cause the sampling points to become excessively concentrated or sparse. Hence, an optional sampling
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
23
criterion can be constructed according to the choice of Theunissen [19] to use the standard deviation of the local variations of the velocity.
1 W2
σ k ( x , y) =
i = W /2
j = W /2
∑
∑
( uk ( x + i, y + j) − uk )2 , k = 1, 2
i =−W /2 j =−W /2
(1)
where x and y are the local coordinates and k are the index of the U and V component velocity. uk and uk are the local velocity and the mean velocity of the specified area W, respectively. A combination of the standard deviation is utilized in general.
σ=
1 2
2
∑ σk2
(2)
k=1
This sampling adaptive criterion has been constructed and applied in the experimental processing procedure by Theunissen [19]. Its robustness and reliability have also been verified. A different sampling criterion is adopted and constructed in this paper. Two different parameters which have been found to work well are related to the local gradients and the curl of the velocity.
λ1 (i, j ) = ‖ curl u‖
(3)
2
λ2 (i, j ) = ‖gradient u ‖
u =
∑ uk2
(4)
k=1
where u is the velocity vector and k are the index of the U and V component velocity. The reasons for this choice are the gradients of velocity can measure the local change of the flow and sense shock waves effectively and the curl criterion can use the local rotationality of the flow and performs well in locating shear layers. The standard deviations of the two parameters are computed, and the ratios of the local detected parameter to the relative standard deviation are used to characterize the refinement region.
ζk (i, j ) =
λ k ( i, j ) 1 XY
i=X ∑i = 1
j=Y ∑ j=1
(5)
where i and j are the coordinates and k are the index of the different parameters. X and Y are the matrix sizes of the parameters. The parameters are normalized by intermediate values to distinguish the relative value of the flow feature parameters. The function is defined as follows:
ξk (i, j ) = ζk (i, j )/median (ζk (i, j ))
(6)
After the relative value of the entire region is computed, additional locations where the sampling points should be placed can be determined. Given that the initial distribution of the sampling points has been determined, the function cannot be applied directly and must be transformed to another form. The above functions shows that refinement should be conducted when ξk is larger, which means the distances of the sampling points decrease as the value of ξk increases, and vice versa. Under this requirement, a specified function g should be defined as:
ρk (i, j ) = g (ξk ) k = 1, 2
(7)
The range of the sampling density is limited, which is from 0.5 to 1.5. When ξk → 0, g is 1.5. When ξk → ∞, g is 0.5. Moreover, g is equal 1 when ξk = 1. Hence, a typical function has been constructed.
⎧ − 0.5ξk 0.88 + 1.5 0 < ξk < 1 k = 1, 2 ρk = g (ξk ) = ⎨ ⎩ 0.4ξk−1.10 + 0.6 1 ≤ ξk ⎪
⎪
Fig. 1 displays the distribution of the sampling density function. This figure illustrates that the sampling ratio decreases as the value of ξk increases. Moreover, the max value of the sampling ratio is 1.5, and the minim one is 0.5. This setting can avoid excessive refinement or coarsening of the sampling points. From the foregoing section, two detected flow features can be selected as the adaptive criterion, which causes different distributions of the sampling points. A combined function is used in this paper to make better use of the two features.
(
ρfinal = min ρgradient , ρcurlz
)
(9)
The sampling density of the velocity information has been determined. The main procedure of distributing the sampling points and its applications in the processing synthetic and experimental images are introduced in the next subsections. 2.2. Adaptive sampling procedure
k = 1, 2
(λ k ( i, j ) − λ k )2
Fig. 1. Distribution of the sampling density as a function of ξk .
(8)
Measurement point sampling can be implemented once the sampling criterion is settled (including the seeding density and the criterion of the velocity). This subsection details the procedure of distributing the sampling points. The sampling criterion for the velocity and the seeding density can be obtained according to the method described in Section 2.1. The main procedure of the method is similar to that of Yu et al. [21]. The concept of artificial spring force is utilized. The initial distances of the sampling points are determined after initializing the distributions of the sampling points. The deviations of the initial and the desired distances of the points produce spring force. Spring force pushes the sampling points to the ideal place under the iterative procedure. The concept is first proposed by Persson and Strang [23] in the generation of CFD mesh with high quality. The mesh generated by the spring mode method is verified to be higher than that filtered by the Laplacian method. The detailed procedure is illustrated briefly as follows. 1. The number of the interrogation windows is set according to the number of identified particle images (N), the number of particles within the interrogation window (NI), and the mean overlap value (OF). Hence, the number of the interrogation windows can be defined as
Nw =
N (1 − OF )2NI
(10)
24
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
in Subsection 2.1. 5. The seeding density and the sampling density are attained. The velocity information causes the desired distance to change although the sampling points meet the requirement of the seeding density in the first sampling step. The sampling points move to the new iteration place after the iteration because of the spring force. The final positions and the number of the sampling points are then determined. Then, the final velocity is obtained by using the interrogation method. The procedure of the method can be summarized as follows: The seeding density is computed. The points are placed according to seeding density. The flow feature is computed, and the sampling density function is constructed. The points are placed according to the sampling density based on the position obtained from the first sampling step. The final velocity is computed. Fig. 3 schematically displays the procedure of the adaptive processing method. 2.3. Interrogation method Fig. 2. Typical 2D Delaunay triangulation and generation of the artificial spring force.
The sampling points can then be distributed randomly. The random algorithm in the paper adopts the method proposed by Kocis and Whiten [24]. 2. The sampling points are triangulated by Delaunay triangulation. Fig. 2 displays the typical 2D Delaunay triangulation. The neighboring sampling points of the specified point Pi can be recognized once the triangulation is generated. The desired distance between Pi and Pj, which is associated with the local seeding density, the deviation of the original, and the desired distances, can be obtained. The force is then generated.
⎧ if hj > L ij ⎪ k ( hj − L ij ) Fij = ⎨ ⎪ if hj ≤ L ij ⎩0
(11)
where Fij is the repulsive force to the node Pi from Pj, and hj is the ideal distance between Pi and Pj. The ideal distance is deduced from the seeding density in the first sampling step and from the combination of the distances from the first step and the sampling density in the second sampling step. Lij is the original distance between Pi and Pj. The resultant force of each point is the addition of Fij. Sd is seeding density (in particles per pixel 2)).
hk, j
⎧ NI ⎪ (1 − OF ) k=1 =⎨ Sd ⎪ k=2 ⎩ ρ ⋅L ij
(12)
3. After finishing the computation of the resultant force, the equilibrium positions where all the measurement points meet the requirement of seeding density and the sampling ratio are solved from the system F(p) ¼0. p is the location of sampling points. The details of the solution of the system are not introduced because they are beyond the scope of this work. Under the sophistic hypothesis, the new positions of the sampling points can be obtained by using the forward Euler method .
( )
pin + 1 = pin + Δt⋅F pin
(13)
where n is the iterative number. The sampling points move gradually to the equilibrium positions according to the equation. 4. After approximately 10–20 iterations, the sampling points in the first sampling step arrive at the equilibrium positions. On the basis of these points, the velocity can be obtained by using the cross-correlation method and the sampling density introduced
Fig. 3 schematically illustrates the procedure of the adaptive PIV processing method. The displacements of each sampling points are obtained by cross-correlation based on fast Fourier transform. An iterative procedure is adopted (i.e., iterative window deformation method [4]). Contrast-limited adaptive histogram equalization [25] is utilized before the core processing to enhance the contrast of the original PIV images and increase the reliability of the final result. A high-pass image enhancement method is applied in the computation of the seeding density to decrease the computational cost. For comparison, the traditional method utilized in the following two verification cases adopts the same way to obtain the displacements of the interrogation windows, including iterative window deformation method, the image preprocessing method and the verification method for the final dataset. And the final interrogation areas in iterative window deformation method are determined according to the actual conditions. The predictor needed by the construction of the sampling criterion is acquired by one-pass cross-correlation after the first sampling based on seeding density. Given that the sampling points are not uniformly distributed spatially, an interpolation scheme must be selected to obtain the gridded data. The velocity results are interpolated over the regular grid from the discrete sampling points using natural neighbor interpolation in this paper [26]. The sampling criterion is then computed and filtered by the method proposed by Garcia [27] to increase the stability of the algorithm. The sampling points are relocated using artificial spring force once the final sampling criterion is obtained. Finally, the velocity on the sampling point is computed by using the window deformation method. Moreover, the interrogation areas are decreased to the specified value iteratively. The final dataset is verified by the global standard derivation and the universal outlier detection method proposed by Westerweel [28]. For comparison with the traditional processing method, the scattered data can be interpolated to the uniform grid by natural neighbor interpolation.
3. Results and discussion 3.1. Performance evaluation of adaptive algorithms in synthetic images 3.1.1. Generation of synthetic images In this subsection, the proposed method is applied to process synthetic images. The results from the proposed method are compared with those obtained from the traditional processing
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
25
Fig. 3. Main procedure of the adaptive processing method.
image. The vectors obtained from the traditional processing method (i.e., uniform sampling and window deformation method) are also shown in Fig. 4, where the Rankin vortex clearly appears in the middle of the images.
Fig. 4. Synthetic image and vector result from the traditional processing method.
method. The synthetic images are generated according to Raffel et al. [29]. The influence of inhomogeneous seeding density in synthetic application is not verified again in this paper. The uniform seeding density is adopted in the verification case, namely, the Rankine vortex, to focus on the performance of the sampling density rate of the velocity. Rankine vortex is a circular flow. The inner part of the vortex is in solid rotation, whereas the outer part is free of vorticity. The distribution of the velocity is defined as follows:
⎧ V ⎪ ⎪ max Vθ ( r ) = ⎨ ⎪V ⎪ ⎩ max
r R R r
3.1.2. Processing of synthetic images The particle image detection algorithm identified approximately 18,221 particles. The average Sd is approximately 0.038/pixels2. Adrian's suggestions indicate that the number of particles within the interrogation windows is set as eight to ensure the robustness of the cross-correlation algorithm [2]. The value of overlap factor is set as 50%. The number of interrogation windows can be computed as 9110 using Eq. (10). The procedure illustrated in Fig. 3 shows that the initial sampling step is conducted. The first distribution of the sampling points instructed by the seeding density can be attained after 10 iterations. The distributions of sampling points in this procedure are also uniform because of the uniform seeding in the synthetic images, see Fig. 5. The gradient and the curl of the velocity can be computed based on the
0≤r≤R r>R
(14)
where Vmax is the max velocity of the vortex and is set as 5 pixels in this verification case, and R is the characteristic distance of the vortex (e.g., the core radius of the vortex) and is set as 80 pixels. The image size is 800 pixels 600 pixels. The vortex center is located in the center of the image. Fig. 4 displays the synthetic
Fig. 5. Distributions of the sampling point based on seeding density.
26
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
Fig. 6. Contour of the sampling density (i.e., the distribution of the ideal distances between the sampling points); the distributions of the sampling points based on the combination of the seeding and velocity information.
preliminary velocity result. The sampling density can be derived accordingly using the above equations. In other words, the distribution of the ideal relative pixel distances of the sampling point can be obtained. Fig. 6 displays the contour of the sampling density, where the value represents the relative distances between the sampling points. The intermediate region should have more sampling points because of the existence of the vortex. Therefore, more measurement points are sampled in Regions A and B. Moreover, the distances between the sampling points in Regions A and B are smaller than the others, and the size of the interrogation windows is also decreased accordingly. Iterative window deformation strategy is implemented after distributing the sampling points to process the images and obtain the final result. Processing the PIV images by using the traditional method requires an appropriate choice of final interrogation area, considering the reliability of the dataset and spatial resolution. A constant sample area of 18 18 pixels2 with a 50% overlap was applied in this case. In most cases of evaluating the performance of the algorithm of PIV, bias error and random error are adopted for the evaluation [30]. The bias error is quantified by the deviation between the theory and the measured result. The random error level is usually quantified as the root-mean-square (RMS) fluctuation of the measured values around the theory value. Here, RMS is regarded as the random error. The formula is defined as follows:
σ=
1 Ns
Ns
∑ (X (i)−XT (i))2 i=1
(15)
where Ns is the total number of the sampling points, X is the local measured velocity, and XT is the local theory result. The two different errors are computed from the U and V components. The average bias error of the V component from the adaptive processing is 0.056 pixels, whereas the one from the traditional method is 0.0604 pixels. The same tendency can be obtained from the U component. The bias error of the adaptive processing method is 7.28% less than that of the traditional method. With regard to the random error, the difference between the two methods becomes relatively small. The error from the adaptive one is lower than that from the traditional method by approximately 4%. At the same time, the contours of the vorticity from the two different methods are computed (Fig. 7). Fig. 7(a) presents the result of the traditional processing method. Fig. 7(b) shows the result of the adaptive method. The theoretical value can be computed based on Eq. (14), the vorticity value in the middle region should be 0.125 pixels/pixels. From Fig. 7, it can be found that the vorticity of the middle region in Fig. 7(b) is more uniform and closer to the theoretical value than that in Fig. 7(a). That also means that the vorticity from the traditional method is fairly disturbed and significantly non-uniform. Besides that, the boundary between the inner and outer flow of Rankine vortex from the adaptive method is much sharper than that from the traditional method. Then, the contour line in Fig. 7(b) is smoother and it should be more accurate to represent the Rankine vortex. Most obviously, it can be observed from Fig. 7 that the fluctuations in the traditional case are much larger than those in the adaptive case. The same phenomena can be verified by the distribution of the vorticity and velocity magnitude in the middle horizontal line from the two cases (Fig. 8). Previous studies showed that the result of cross-correlation could be modeled with a moving average linear filter [31]. Thus, the theory vorticity in Fig. 8(b) is filtered by a moving average filter with a 16-pixels window length for intuitive comparison. The analysis of the bias and random error and the distribution of the vorticity indicate that the adaptive processing method can increase the performance of the PIV algorithm, the reliability of the dataset, and spatial resolution. 3.2. Performance evaluation of adaptive algorithms in experimental images In this subsection, experimental PIV images are processed to demonstrate the good performance of the adaptive processing method qualitatively. Experiments on jet flow are measured by PIV provided by the PIV Challenge 2005 (http://www.pivchallenge. org). The traditional processing method and the adaptive method are applied to obtain the final result for an intuitive comparison.
Fig. 7. (a). Contour of the vorticity from the traditional processing method; (b). Contour of the vorticity from the adaptive processing method.
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
27
Fig. 8. (a). Distributions of the vorticity in the middle horizontal line; (b). Distributions of the velocity magnitude in the middle horizontal line. Red represents the theoretical result, blue represents the traditional result, and green represents the adaptive result.
in the adaptive processing algorithm based on the previously mentioned equations is obtained as 4892. The overlap value is 50% in both methods. The NI is eight according to the suggestions by Adrian [2]. After determining the necessary parameters for the adaptive processing (i.e., the Sd, the sampling density function, and the initial number of sampling points), the sampling points based on the seeding density can be distributed accordingly (Fig. 9). The reliability of the adaptive sampling method based on the seeding density is verified again. More sampling points are placed in the core flow region, and vice versa. The sampling density is then constructed by using the equations combined with the initial velocity information. The final distribution of the measurement points can be obtained (Fig. 10). The white points in Fig. 10 illustrate the location of
Fig. 9. PIV image. The gray level of the image is inversed, and contrast enhancement processing is applied. Blue points represent the sampling points, which are distributed based on the seeding density. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The flow of nitrogen in this experimental application has a velocity of 30 m/s. The viewing area is 29 29 mm2. The resolution of the PIV raw image is 512 512 pixels. The seeding density is fairly inhomogeneous because of inhomogeneous seeding of the flow and the co-flow. With the use of the particle detection algorithm, the number of the detected particle images is approximately 9785, and the seeding density can be obtained as 0.0373 particles/pixels2. Fig. 9 presents the raw PIV images and shows that the Sd in the main flow region is much higher than that in the external flow region. The gray level of the raw image has been inversed. Contrast-limited adaptive histogram equalization method has also been applied to enhance the contrast of the image and the performance of the algorithm. A trade-off choice must be made to select the interrogation window size appropriately. Hence, the final interrogation area under this requirement is set as 32 pixels because of the low Sd in the external flow region in the traditional processing algorithm. The number of sampling points
Fig. 10. Distribution of the sampling points based on the seeding density and the velocity information. White points represent the sampling points. The corresponding vectors are also displayed.
28
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
Fig. 11. (a). Contour of the vorticity from the traditional processing method; (b). Contour of the vorticity from the adaptive processing method.
the measurement points. Figs. 9 and 10 distinguish the significant differences. Fig. 10 samples more points in the boundary of the main flow because of the high vorticity and gradient. Some measurement points were sampled in the external flow under the requirement of sampling density function. Thus, the measurement requirement is more accurately met. The velocity and its derivations are obtained by using iterative window deformation method. Fig. 11 displays the vorticity from the traditional and adaptive method. Due to the augmented measurement points sampled in the main flow, the adaptive method is able to capture more details about the main flow. And, the adaptive method suppresses noise, especially in the left flow region. From the previous paragraph, it can be known that the seeding density in the outer region of the jet is smaller than that in the main flow. Under the limitation of the traditional method, the fixed interrogation areas have been determined. Thus, in the processing of the outer region, the number of particle images in the interrogation domain is decreased, which causes the amplitude of the displacement-correlation peak reduced. Failure to correctly identify these peaks leads to spurious vectors. So some fluctuations exist in the result obtained by the traditional method. On the contrary, adaptive method places few sampling point in the outer region and utilizes large interrogation windows, which meet the requirement of the cross-correlation algorithm and make that the velocity results become reliable. In fact, the conclusion obtained from the synthetic case also can verify this phenomenon. The random and bias errors from the adaptive method are less than that from the traditional method. That is, as to zero displacement, more fluctuations will be shown in the result obtained by the traditional method. In a word, the experimental application can verify the effectiveness and reliability of the proposed adaptive algorithm for PIV processing.
sampling point to the local seeding density and the velocity information unlike the use of the regular mesh in the traditional method. The concept of spring force adopted in previous research is utilized in this study to distribute the sampling points. The proposed method places the sampling points on the basis of the seeding density in the first procedure. After that, a specified sampling density function is constructed on the bias of the vorticity and gradient of the velocity. The final distribution of the points can then be obtained. Two test cases have been chosen to attest the effectiveness of the proposed method: a synthetic test (Rankine vortex) and an experimental test (jet flow). The synthetic example uses the uniform seeding density to verify the operability of the proposed method. The final velocity result indicates that the proposed method can increase the performance of the PIV algorithm by 7% in terms of the bias error and 4% in terms of the random error. The vorticity is smooth in the adaptive processing method, which agrees well with the theory. The actual images are then processed; results verify that the proposed method is reliable and effective. Compared with iterative window deformation method, spatial resolution is locally enhanced, and the dataset is more reliable. Further efforts will be directed toward the consideration of temporal information in PIV processing algorithm.
Acknowledgments We would like to acknowledge the continued support of the NSFC (Natural Science Fund of China) for contract numbers 50876042 and 90916023.
References 4. Conclusions An adaptive PIV processing algorithm is proposed to enhance the performance of the processing algorithm and increase spatial resolution. The proposed method adapts the distributions of the
[1] C.E. Willert, M. Gharib, Digital particle image velocimetry, Exp. Fluids 10 (4) (1991) 181–193. [2] R.J. Adrian, J. Westerweel, Particle Image Velocimetry, Cambridge University Press,, 2011. [3] R.D. Keane, R.J. Adrian, Theory of cross-correlation analysis of PIV images, Appl. Sci. Res. 49 (3) (1992) 191–215.
K. Yu, J. Xu / Flow Measurement and Instrumentation 51 (2016) 21–29
[4] F. Scarano, Iterative image deformation methods in PIV, Meas. Sci. Technol. 13 (1) (2002) R1. [5] T. Astarita, G. Cardone, Analysis of interpolation schemes for image deformation methods in PIV, Exp. Fluids 38 (2) (2005) 233–243. [6] T. Astarita, Analysis of interpolation schemes for image deformation methods in PIV: effect of noise on the accuracy and spatial resolution, Exp. Fluids 40 (6) (2006) 977–987. [7] F.F.J. Schrijer, F. Scarano, Effect of predictor–corrector filtering on the stability and spatial resolution of iterative PIV interrogation, Exp. Fluids 45 (5) (2008) 927–941. [8] R.D. Keane, R.J. Adrian, Y. Zhang, Super-resolution particle imaging velocimetry, Meas. Sci. Technol. 6 (6) (1995) 754. [9] A. Stitou, M.L. Riethmuller, Extension of PIV to super resolution using PTV, Meas. Sci. Technol. 12 (9) (2001) 1398. [10] A. Susset, J.M. Most, D. Honoré, A novel architecture for a super-resolution PIV algorithm developed for the improvement of the resolution of large velocity gradient measurements, Exp. Fluids 40 (1) (2006) 70–79. [11] J. Westerweel, P.F. Geelhoed, R. Lindken, Single-pixel resolution ensemble correlation for micro-PIV applications, Exp. Fluids 37 (3) (2004) 375–384. [12] F. Billy, L. David, G. Pineau, Single pixel resolution correlation applied to unsteady flow measurements[J], Meas. Sci. Technol. 15 (6) (2004) 1039. [13] C.J. Kähler, U. Scholz, J. Ortmanns, Wall-shear-stress and near-wall turbulence measurements up to single pixel resolution by means of long-distance microPIV, Exp. Fluids 41 (2) (2006) 327–341. [14] M. Bitter, S. Scharnowski, R. Hain, C.J. Kähler, High-repetition-rate PIV investigations on a generic rocket model in sub-and supersonic flows, Exp. Fluids 50 (4) (2011) 1019–1030. [15] J. Rohaly, F. Frigerio, D.P. Hart, Reverse hierarchical PIV processing, Meas. Sci. Technol. 13 (7) (2002) 984. [16] D. Di Florio, F. Di Felice, G.P. Romano, Windowing, re-shaping and re-orientation interrogation windows in particle image velocimetry for the investigation of shear flows, Meas. Sci. Technol. 13 (7) (2002) 953. [17] B. Wieneke, K. Pfeiffer, Adaptive PIV with variable interrogation window size
[18]
[19] [20] [21] [22] [23] [24] [25] [26]
[27] [28] [29] [30]
[31]
29
and shape. 15th Int Symp on Applications of Laser Techniques to Fluid Mechanics, Lisbon, Portugal, 2010. F. Becker, B. Wieneke, S. Petra, A. Schröder, C. Schnörr, Variational adaptive correlation method for flow estimation, Image Process. IEEE Trans. 21 (6) (2012) 3053–3065. R. Theunissen, F. Scarano, M.L. Riethmuller, An adaptive sampling and windowing interrogation method in PIV, Meas. Sci. Technol. 18 (1) (2007) 275. R. Theunissen, F. Scarano, M.L. Riethmuller, Spatially adaptive PIV interrogation based on data ensemble, Exp. Fluids 48 (5) (2010) 875–887. K. Yu, J. Xu, L. Tang, J. Mo, New adaptive sampling method in particle image velocimetry, Meas. Sci. Technol. 26 (3) (2015) 037002. M. Novara, A. Ianiro, F. Scarano, Adaptive interrogation for 3D-PIV, Meas. Sci. Technol. 24 (2) (2013) 024012. P.O. Persson, G. Strang, A simple mesh generator in MATLAB, SIAM Rev. 46 (2) (2004) 329–345. L. Kocis, W.J. Whiten, Computational investigations of low-discrepancy sequences, ACM Trans. Math. Softw. (TOMS) 23 (2) (1997) 266–294. K. Zuiderveld, Contrast Limited Adaptive Histogram Equalization (In Graphics gems IV), Academic Press Professional, Inc 1994, pp. 474–485. M. Sambridge, J. Braun, H. McQueen, Geophysical parametrization and interpolation of irregular data using natural neighbours, Geophys. J. Int. 122 (3) (1995) 837–857. D. Garcia, Robust smoothing of gridded data in one and higher dimensions with missing values, Comput. Stat. Data Anal. 54 (4) (2010) 1167–1178. J. Westerweel, F. Scarano, Universal outlier detection for PIV data, Exp. Fluids 39 (6) (2005) 1096–1100. M. Raffel, C.E. Willert, J. Kompenhans, Particle Image Velocimetry: A Practical Guide, Springer, 2007. L. Gui, S.T. Wereley, A correlation-based continuous window-shift technique to reduce the peak-locking effect in digital PIV image evaluation, Exp. Fluids 32 (4) (2002) 506–517. F. Scarano, Theory of non-isotropic spatial resolution in PIV, Exp. Fluids 35 (3) (2003) 268–277.