Instructions to Authors

11TH INTERNATIONAL SYMPOSIUM ON PARTICLE IMAGE VELOCIMETRY – PIV15 Santa Barbara, California, September 14-16, 2015

Least-square reconstruction of instantaneous pressure field around a body based on a directly acquired material acceleration in time-resolved PIV Young Jin Jeon1, Ludovic Chatellier 1, Anthony Beaudoin1 and Laurent David1 1

Institut P PRIME, UPR3346, CNRS – Université de Poitiers – ISAE-ENSMA, France [email protected]

ABSTRACT This paper summarizes the recently introduced procedure for evaluating an instantaneous pressure field around a body from time-resolved particle image velocimetry (TR-PIV) [1]. Instantaneous velocity and material acceleration fields, from which a pressure gradient field is derived, are provided by the multi-frame PIV method, fluid trajectory evaluation based on ensemble-averaged cross-correlation (FTEE) [2]. The pressure gradient field is then delivered to two processes: discretization of an integration domain and sequential integration. The discretization is introduced to separate the integration domain into several subdomains according to the local reliability of the measurement. And then the pressure field of the most reliable subdomain is reconstructed first and then the remaining subdomains are sequentially reconstructed according to those decreasing reliability. In order to reconstruct the pressure field by spatially integrating the pressure gradient field, a least-square method, which minimizes a sum of squares of the difference between the reconstructed and the measured pressure gradients, is introduced. A set of matrix equations, generated adaptively to the shape of subdomain, yields a unique solution without dependency on integration path. An application to a detached flow, from an experimental TR-PIV measurement of a NACA0015 airfoil with an angle of attack of 30 degrees at Re = 10 5, is proposed. The result is analyzed in terms of both instantaneous and spectral aspects. 1. INTRODUCTION An instantaneous pressure field around a body allows both understanding flow structures in relation to local pressure fluctuations and obtaining force acting on the body. This involves two practical issues: how to integrate the spatial pressure gradients around the body and how to consider the body and its influence on the flow. The present study has proposed a practical solution to these problems by introducing two processes: discretization of the interrogation domain based on the local reliability and least-square reconstruction of the pressure field. In a flow around a body, the body disturbs the flow over the edges and in addition cause a wake region behind the body, particularly for turbulent flows. Although the recent studies in TR-PIV have shown significant improvements [2, 4, 5, 6] in terms of the dynamic range, there is a certain limitation in the disturbed flow regime due to the high order fluid motion. Therefore, this issue should be considered either while reconstructing the pressure field. Recently, Tronchin et al. [7] have proposed a sequential integration with respect to subdomains, into which an entire integration domain is divided, in order to prevent the propagation of error from the less reliable region to the more. We have improved this technique by introducing new criteria based on both recorded body image and Frobenius norm of pressure gradient tensor. As intensively reviewed by van Oudheusden [3], the instantaneous pressure field can be obtained by integrating spatial pressure gradients. The pressure gradient can be derived from a sequence of time-resolved velocity fields and/or a material acceleration field and then spatially integrated by following either of two main strategies: solving a Poisson equation and spatial-marching integration. As an alternative to both strategies, the spatial integration can be considered to be an inversion of differentiation. As an image processing method, Harker and O’Leary [8] have proposed a least-square based twodimensional method, which reconstructs a surface level distribution from vertical- and horizontal-gradient images. We have simplified this approach by introducing vectorizations of planar/volumetric fields, in order to build the matrix equation for both three-dimensional pressure gradients and an arbitrary domain shape. This new least-square based method has two characteristics: resulting a unique solution independent from an integration path (same with the solving a Poisson equation


[9, 10], omni-directional integration [11, 12, 13] and iterative convergence [7]) and directly dealing with pressure gradients and Dirichlet boundary conditions only (same with the spatial-marching integration). In order to validate the present procedure, a high-rate TR-PIV experiment of a flow around an airfoil has been conducted, and FTEE was applied to process the particle image sequence in order to provide the instantaneous fields of velocity and material acceleration, by which the pressure gradient is derived (section 2). The discretization of the integration domain into several subdomains with respect to the local reliability is described (section 3), and the least-square formula for integrating pressure gradients according to the discretized subdomains is follows (section 4). The instantaneous pressure fields and their spectral characteristics are shown (section 5). The further possibility of the present method is discussed and summarized in section 6. 2. EXPERIMENTAL SETUP A TR-PIV experiments of a flow around a NACA0015 airfoil of chord c = 80mm with a fixed angle of attack α = 30° has been carried out in a closed-loop water tunnel with a cross section of 235 × 235 mm2 (Figure 1). The free stream velocity is 1.25 m/s, corresponds to Re = 10 5 based on the chord length. The image sequences were acquired by two 1024 by 1024 pixels cameras operated in tandem at a framerate of 1.5 kHz. A total of 20000 images were continuously recorded by each camera during 14.6 seconds. The test section was illuminated by two laser systems placed top and bottom sides, in order to produce a planar laser sheet without shadowed areas.

Figure 1 Schematic of TR-PIV experiment FTEE [2] provides instantaneous fields of velocity and material acceleration based on a second order polynomial trajectory over 9 images. The final interrogation window is 32 × 32 pixels with 75% overlap with Gaussian weighting [14] while considering a total 243 by 125 grid points. The images are linearly deformed during FTEE based on the sinc interpolation scheme with an 8 × 8 pixel kernel [15]. Here the fluid motion at each grid point, x0, is assumed to follow the polynomial trajectory in the same manner as previous approaches [5, 6, 16, 17, 18], K

x p  n   x 0   ai ni

(1)

i 1

where xp is the position of fluid parcel, K is the given polynomial order, ai is the i-th order polynomial coefficient, and where n = t/Δt is the normalized time step while Δt is the temporal interval between two subsequent recordings. The derivatives at n = 0 then yield the velocity up and material acceleration Du/Dt|p as follows up  0 

dxp dt

 0   a1 and

d2 x Du  0   2p  0   2a2 . Dt p dt

(2)

Since the dynamic ranges of the velocity fields are improved by FTEE, there is another possibility to obtain the material acceleration from the sequential velocity fields by means of the Lagrangian approach [19]. In this study, three sequential velocity fields are considered to iteratively obtain material acceleration as follows,


Du Dt

k 1

 x, t  

u  x p k  t  t  , t  t   u  x p k  t  t  , t  t  2t

1 Du , where x p  t  t   x  u  x, t  t  2 Dt k

k

(3)

 x, t  t

2

here k denotes the iteration number, and the material acceleration (Du/Dt)|k is initially set as 0 at k = 0. Then the material acceleration can be iteratively converged into a certain vector. The velocity vectors at fractional locations within the velocity fields at the adjacent time steps, t ± Δt, are linearly interpolated. However, because this Lagrangian approach has a drawback in the wake region due to the regressed dynamic range of the velocity fields, the present study utilizes two kinds of material accelerations, from FTEE and the Lagrangian approach, with respect to those local reliabilities. 3. DISCRETIZATION OF INTEGRATION DOMAIN In order to consider the body immersed in fluid and its effect on the flow, the integration domain is discretized into five subdomains: masked, near-edge, near-body, wake, and outer regions. Note that the integration domain is identical to the PIV grid scheme. Figure 2 illustrates the stacked particle images while being overlaid with the discretized subdomains.

Figure 2 Schematic of discretization and stacked particle images of 9 continuous recordings: real wing profile and its recorded image are shown as red and blue lines respectively, and interface between the outer region and the others is figured by black dashed line. The grid points are assigned to each subdomain according to the local measurement reliability. Firstly, the masked region Dmask is introduced both to exclude grid points inside the recorded body image and to consider the perspective effect. Since the pressure gradient near the leading and trailing edges are imprecise due to the high shear, grid points inside circles of radius d edge whose centers are the leading or trailing edge are considered as the least reliable domain D4 . In addition, to take into account the possibility of error in a vicinity of the body, the grid points within a distance of dbody from the recorded body image are considered as the near-body region D3. In the TR-PIV experiment, D4 and D3 are set with d edge = 72 pixels and dbody = 44 pixels (= 4.5Δh) respectively. Criteria for the discretization are listed in Table. 1. Table 1 Criteria for the discretization Integration order 1 2 3 4 -

Subdomains Outer region Wake region Near-body region Near-edge region Masked region

Information for consideration Frobenius norm of pressure gradient tensor Recorded body in particle image Physical position of the body Recorded body in particle image


In contrast to the image-based subdomains, Dmask , D4 and D3, the time-varying fluid motion should be considered while assigning the wake region D2. Based on an assumption that the material acceleration suddenly changes at the boundary of the wake region, it is proposed here to distinguish the wake region by comparing a local Frobenius norm of the pressure gradient tensor divided by the density: here the density is introduced to express the formula in terms of a temporal unit only. 2

p



F

   pi    -2      0.005 s  x i j  j 

(4)

Note that the unit of Eq. 4 is 1/Δt2, and equal to a Frobenius norm of the material acceleration tensor. This is only valid when the random error of measured pressure gradient at the outer region is less than a specific criterion that is further detailed in [1]. Unfortunately, this criterion cannot provide a region with smooth boundary but a set of split pieces and scattered grid points, D20, as shown as dashed lines in Figure 3. Therefore, additional morphological processing methods are needed. In order both to obtain the wake region with a smooth boundary and to avoid an isolated or narrow region in D1 and D2 , the reclamation procedure that fills a space between the split regions in D20 is introduced. The reclamation consisted of sequential morphological processing methods: binary dilation, erosion, and dilation, D2 = ((D20 ⊕ S4.5Δh) ⊝ S8.1Δh ) ⊕ S3.6Δh,

(5)

here Sr is the circular structure, Sr = { x | |x| ≤ r }, ⊕ is the binary dilation operator, and ⊝ is the binary erosion operator. Figure 3 demonstrates the sequential morphological processes based on Equation 5.

Figure 3 Schematic of the binary processes As mentioned in the last two paragraphs in section 2, the present study uses two kinds of evaluation schemes for material acceleration, FTEE and the Lagrangian approach. In order to apply better evaluation scheme with respect to each subdomain, the Lagrangian approach is applied to the outer region and FTEE to the other subdomains. Further details are described in [1]. 4. SEQUENTIAL LEAST-SQUARE RECONSTRUCTION OF PRESSURE FIELD The discretized subdomains and the corresponding pressure gradient fields are then delivered to the sequential least-square reconstruction. Note that the equations are expressed in two-dimensional notations, but compatible with three-dimensional notation.


The cost function s is defined as a sum of squares of the difference between the measured and the reconstructed pressure gradients over the certain integration domain D as follows, s   pPIV  p

(6)

xD

Where ∇pPIV is the measured pressure gradient, ∇p is the reconstructed pressure gradient, and x is the position of grid point. The staggered-grid scheme is adopted for both evaluating the pressure at each PIV grid point and considering the pressure gradient between the grid points (Figure 4). Therefore ∇pPIV is given by averaging the pressure gradients of two grid points. To express the right hand side of Eq. 6 as a matrix equation, the pressure and the pressure gradients are vectorised. p  vecD  p  x     p1 p2 p3 ...

T

p x  vecD  p  x  x    px , 1 px , 2 px , 3 ...

T

p y  vecD  p  x  y    p y , 1 p y , 2 p y , 3 ...

(7)

T

Here vecD denotes the vectorization of elements which belong to the domain D. Relations between these vectors can be then expressed as, (8)

p x  G x p and p y  G y p

where Gx and Gy are the discrete gradient matrices. Note that the element indexes of p, px and py can be arbitrarily defined, and the sizes of the vectors depend on the shape of domain. Figure 4 illustrates an example of the vectorization and corresponding matrix equations.  px , 1  1  p  0 1   x, 2     px , 3   1 1      0 1  px , 4     px , 5   1   1 1   x   1 1  px , 6    p    0 1  x, 7     px , 8  1 1       p  1 1    x, 9 

px

Gx

 py, 1  0 1  p   1 0 1   y, 2     py, 3  0 0 0 1      p  1 0 0 1  y, 4   1    p y , 5  y   1 0 0 1     1 0 0 1  py, 6    p   1 0 0 1   y, 7     p y , 8  1 0 0 1  

py

Gy

 p1    pb, 1  p   p   2  b, 2   p3   0       p4    pb, 3   p5   1  0    x    p6   0  p   p   7  b, 4   p8   0      p  0   9

p

p x, b.c.

 p1    pb, 1  p   0   2    p3    pb, 2       p4  0   p5   1    y  0     p6   0  p   0   7    p8   0     p9  p y, b.c. p

Figure 4 Schematic of staggered grid scheme and corresponding vectorizations while integrating the pressure field at Dn: yellow squares denotes the grid points whose pressures are now being considered, red squares indicates the grid points regarded as Dirichlet boundary, and gray squares the grid points which belong to the less reliable subdomain Dn+1 than the current subdomain Dn. The additional vectors, px, b.c. and py, b.c., are introduced to consider the Dirichlet boundary points.


The cost function (Equation 6) can be now written as,



s   p x , PIV  G x p  p y , PIV  G y p xD



(9)

This can be differentiated with respect to the vector p and equalized to 0, s  2  G Tx G x  G Ty G y  p  2  G Tx p x , PIV  G Ty p y , PIV   0 p

(10)

This yields the following matrix equation.

G G T x

x

 GTy G y  p  GTx p x , PIV  GTy p y , PIV

(11)

Note that (GxTGx + GyT Gy) is a symmetric block tri-diagonal matrix. When the vector p has N elements, a rank of (GxT Gx + GyT Gy) is N if at least on Dirichlet boundary point is given, otherwise Equation 11 is singular with the rank of (N - 1). Therefore, one or more reference pressures for the first subdomain D1 should be given. Because we have observed a periodic fluctuation of pressure between the bottom-left and the top-left corner in the experiment, both corner points were considered as the reference points. Firstly, the pressure value at the bottom-left corner is initially set as 0, and the relative pressure field, Δp (x), is reconstructed. Then the reference pressure is modified by considering Bernoulli’s principle and the pressure difference between the two corners,

p (x) 

1 1   U 2  ubottom-left 2   p  xtop-left   p  x  when x  D1 . 2 2

(12)

After the pressure field of D1 is reconstructed, the others are reconstructed by following the integration order. Table 2 shows the number of grid points according to each subdomain for the TR-PIV experiments. Table 2 Number of grid points for each subdomain Integration order 1 2 3 4

Subdomains Outer region Wake region Near-body region Near-edge region

Number of grid points N (mean ± standard deviation) 19847 ± 742 9423 ± 742 380 ± 0 338 ± 0

5. RESULTS and SPECTRAL ANALYSIS Instantaneous pressure fields of 20000 continuous time steps are evaluated, and a total 16384 (= 214 ) instantaneous fields are then spectrally analyzed by means of discrete Fourier transform (DFT). The vortex shedding frequency is observed as 5.68 Hz, corresponding to Strouhal number St = 0.182 based on the characteristic length of 40 mm (= c sin 30°). 5.1 Instantaneous pressure fields Figure 5 shows the instantaneous fields of Frobenius norm of pressure gradient, velocity and pressure while being overlaid with the interfaces between the outer region D1 and the others. The results are sampled on the basis of the vortex shedding period ΔtV.S (1/5.68 sec). Therefore, the result at each time step is somewhat similar to the others. The centers of shed vortices behind the trailing edge are uniformly observed at x / c = 1.35. It is apparent that the boundaries of the wake region are well identified along where the Frobenius norm sharply increases. Although the particle images are recorded with low image quality under the high-rate TR-PIV experiment, the smooth contours of the velocity fields, especially in the outer region, are favorable by FTEE. The detached shear layers are observed just after the leading and trailing edges and exist in the obtained wake region.


The pressure fields seem to be reliable in the outer region unlike the unsteady results in the wake region. Although the material acceleration is available from FTEE in the wake region, the resulting pressure would seem to contain some errors due to the experimental limitations such as the out-of-plane motion in two-dimensional recordings and insufficient spatial and temporal resolution. Nevertheless, the vortex structures inside the wake can be identified.

Figure 5 Instantaneous fields: (left) Frobenius norm of pressure gradient, (center) velocity field, (right) pressure field. Blue dashed lines indicate the interface between the outer region D 1 and the others. ΔtV.S. indicates the vortex shedding period. 5.2 Spectral analysis To investigate the generation and propagation of vortical structures, the amplitude and phase fields with respect to each frequency are considered When the vortical structure, whose pressure is local minimum, is periodically generated at a certain position and moved by the flow, its trace can be found in the phase field of the characteristic frequency in terms of the phase gradient. The propagating speed of the pressure structure corresponding to the frequency f, uf., can be expressed as,

u f  x   2π f

 f  x   f  x 

2

(13)

here ϕ (x) is the phase field in radian, and the unit vector -∇ϕ /|∇ϕ| means the direction of the propagation. Since propagating speeds of the vortical structures cannot exceed the maximum velocity of the flow, approximately 1.5U∞, the


pressure structure which satisfies | uf | < 1.5U∞ can be regarded as a phenomenon caused by the interaction between the flow and body. We have investigated on the basis of the propagation speed in the frequency domain, found several characteristic frequencies, caused by the shed vortices (5.68 and 11.4 Hz) and the leading edge vortices (several peaks at f > 10 Hz) [1]. Figure 6 shows distributions of amplitude and with respect to the sampled characteristic frequencies.

Figure 6 (left) Amplitude field, (center) phase field, and (right) corresponding pressure fluctuation when the phase shift is zero: frequencies caused by (a), (b) the shed vortices and (c), (d), (e) the leading edge vortices Figure 6(a) shows two high amplitude regions, antiphase synchronized with vertically symmetric contour lines. The shed vortices, generated behind the trailing edge, have the maximum amplitude of 0.212 at x/c = 1.38, that is 2.36 times of those behind the leading edge, 0.090. This is because the flow over the trailing edge encounters a suddenly open space behind the inclined airfoil. The shed vortices are altered with each other, thus forming the high amplitude with the doubled frequency, 11.4 Hz, at the centerline of its vortex street (Figure 6(b)). Figures 6(c), 6(d) and 6(e) show the results which emphasise the leading edge vortices. To present clear fields with a small external effect, the frequencies are chosen by checking probable fields one by one. Note that it is difficult to identify


trailing edge vortices because the shed vortices are generated closely behind the trailing edge. The region whose amplitude is relatively high and phase decreases steadily along the flow direction can be regarded as the vortex street, which consists of the leading edge vortices. In this region, the distribution of uf can be assumed to be the same regardless of frequency. Therefore as the frequency, f, increases, the magnitude of the phase gradient, |∇ϕ|, increases. In addition, since vortex structures generally contain a circular local pressure minimum just after the detachment, the width of the vortex street is inversely proportional to the frequency. If the frequency is high than the specific value, about 72 Hz, the specific pattern in phase field is invisible owing to the limitation in spatial resolution. 6. CONCLUSION The practical procedure for evaluating an instantaneous pressure field from TR-PIV on the basis of the least-square method [1] has been summarized. This procedure mainly consists of two processes: discretization of the integration domain and sequential integration based on the least-square method. The integration domain was discretized into subdomains according to local measurement reliabilities: outer region (D 1), wake region (D2), near-body region (D3), near-edge region (D4), and masked region (Dmask). The integration was sequentially conducted from the most reliable region to the most erroneous region; D1  D2  D3  D4. The pressure fields of more reliable subdomains were used as Dirichlet boundary conditions for reconstructing the pressure fields of the latter subdomains. As a result, the erroneous pressure gradients were not retained while considering the more reliable subdomains. Another advantage of the sequential integration based on discretized subdomains was that the scheme for material acceleration can be chosen with respect to each subdomain. In addition, it will be also possible to change parameters, such as polynomial order or number of images in multi-frame PIV or other options, according to each subdomain. The pressure field, which minimizes a sum of squares of whose gradient and the measured pressure gradient, was obtained. The least-square based pressure field reconstruction only required a Dirichlet boundary condition, and yielded a unique solution without dependency of integration path. The relation between the pressure field and its gradient fields were derived into a matrix equation by the vectorizations which convert a field into a column vector. The set of matrix equations is then automatically generated based on each subdomain, and the corresponding example was shown. As the matrix equation is flexible to introduce additional treatments, new approaches can be expected in future works. To explore practical issues of the present procedure while applying it to TR-PIV, a high-rate TR-PIV experiment of flow around an inclined airfoil was carried out. Velocity and material acceleration fields were obtained by means of FTEE. Instantaneous pressure fields were reconstructed by following the present procedure. The possibility of applying the present procedure to three-dimensional measurements was also discussed while deriving the matrix equations. Since the matrix equation, based on the least-square method, took into account the arbitrarily domain shape, it is likely that the present method will be applicable and effectual to flows around a deformable and/or moving body. Finally, the flexibilities of the present method are able to offer further applications, for understanding fluid-structureinteraction phenomena, such as a load evaluation of wing surface. ACKNOWLEDGEMENTS The current work has been conducted as part of the NIOPLEX project, Non-intrusive Optical Pressure and Loads Extraction for Aerodynamic Analysis, funded by the European Commission program FP7, grant n. 605151. REFERENCES [1] Jeon YJ, Chatellier L, Beaudoin A and David L “Least-square reconstruction of instantaneous pressure field around a body from TR-PIV and its application to spectral analysis,” Experiments in Fluids (2015), submitted [2] Jeon YJ, Chatellier L and David L “Fluid trajectory evaluation based on an ensemble-averaged cross-correlation in time-resolved PIV,” Experiments in Fluids 55 (2014) pp.1-16 [3] van Oudheusden BW “PIV-based pressure measurement,” Measurement Science and Technology 24 (2013) 032001


[4] Sciacchitano A, Scarano F and Wieneke B “Multi-frame pyramid correlation for time-resolved PIV,” Experiments in Fluids 53 (2012) pp.1087-1105 [5] Cierpka C and Lütke B and Kähler CJ “Higher order multi-frame particle tracking velocimetry,” Experiments in Fluids 54 (2013) pp.1-12 [6] Lynch KP and Scarano F “A high-order time-accurate interrogation method for time-resolved PIV,” Measurement Science and Technology 24 (2013) 035305 [7] Tronchin T, David L and Farcy A “Loads and pressure evaluation of the flow around a flapping wing from instantaneous 3D velocity measurements,” Experiments in Fluids 56 (2015) pp.1-16 [8] Harker M and O’Leary P “Least squares surface reconstruction from measured gradient fields,” In: IEEE Computer Vision and Pattern Recognition, Anchorage, United States of America (2008) [9] Violato D, Moore P and Scarano F “Lagrangian and Eulerian pressure field evaluation of rod-airfoil flow from time-resolved tomographic PIV,” Experiments in Fluids 50 (2011) pp.1057-1070 [10]de Kat R and van Oudheusden BW “Instantaneous planar pressure determination from PIV in turbulent flow ,” Experiments in Fluids 52 (2012) pp.1089-1106 [11] Liu X and Katz J “Instantaneous pressure and material acceleration measurements using a four-exposure PIV system,” Experiments in Fluids 41 (2006) pp.227-240 [12] Liu X and Katz J “Vortex-corner interactions in a cavity shear layer elucidated by time-resolved measurements of the pressure field,” Journal of Fluid Mechanics 728 (2013) pp.4147-457 [13] Charonko JJ, King CV, Smith BL and Vlachos PP “Assessment of pressure field calculations from particle image velocimetry measurements,” Measurement Science and Technology 21 (2010) 105401 [14] Astarita T “Analysis of weighting windows for image deformation methods in PIV,” Experiments in Fluids 43 (2007) pp.859-872 [15] Scarano F and Riethmuller ML “Advances in iterative multigrid PIV image processing,” Experiments in Fluids 29 (2000) S051-S060 [16] Malik NA, Dracos T an Papantoniou DA “Partic le tracking velocimetry in three-dimens ional flows,” Experiments in Fluids 15 (1993) pp.279-294 [17] Guezennec YG, Brodkey RS, Trigui N and Kent JC “Algorithms for fully automated three-dimens ional particle tracking velocimetry,” Experiments in Fluids 17 (1994) pp.209-219 [18] Li D, Zhang Y, Sun Y and Yan W “A multi-frame particle tracking algorithm robust against input noise,” Measurement Science and Technology 19 (2008105401

[19] Lynch KP and Scarano F “Material acceleration estimation by four-pulse tomo-PIV,” Measurement Science and Technology 25 (2014) 084005