Ψ‐PIV: a novel framework to study unsteady microfluidic flows

In microscopic particle image velocimetry (micro-PIV), correlation averaging over multiple frames is often required, leading to a loss in temporal resolution, therefore limiting the measurement accuracy for unsteady flows. Here, we present a new PIV method suitable to study steady and unsteady laminar flows between parallel plates (i.e., Hele-Shaw flow), which is a common flow configuration in microfluidic applications. Our method reduces the effective seeding density and yields similar if not higher signal-to-noise ratio (SNR) compared to conventional micro-PIV. We call this algorithm Ψ-PIV. Ψ-PIV requires a much smaller number of frames to reach the same SNR compared to the widely used correlation averaging method. This leads to a significant improvement of the temporal resolution. The Ψ-PIV algorithm is used in an experimental investigation of steady and unsteady flows in a Hele-Shaw cell. Our experiment shows that Ψ-PIV reduces the number of required frames by 8 times and 30 times compared to the frames required by conventional PIV for steady and unsteady laminar flow, respectively. In this study, PIV and Ψ-PIV use a single-pass cross-correlation to present the underlying difference between the two approaches.


Introduction
Microscopic particle image velocimetry (micro-PIV) differs from conventional PIV in two main aspects (Adrian and Westerweel 2011;Lindken et al. 2009;Wereley and Meinhart 2010): (1) the image density is usually much lower, N I ≤ 1 and (2) the entire measurement volume is illuminated, rather than a thin light sheet. The low image density is usually compensated by performing correlation averaging (Meinhart et al. 2000), i.e., the spatial correlation for a given interrogation domain over several image pairs. Typically, the correlation averaging over 10-20 frame pairs is required to reach a sufficient effective image density of around 10 (Meinhart et al. 2000;Vennemann et al. 2006). In addition, the use of volume illumination also increases the number of frames N F required. In correlation averaging, N F helps in building up the spatial correlation by the superimposition of the correlation data from the image density of each pair multiplied by the number of frames in the recording sequence. A volume illumination implies that the measurement domain is limited by either the depth of correlation (DOC) (Olsen and Adrian 2000) or the physical edges of the measurement section. When the DOC is small with respect to the depth of the flow domain, there is a little variation in the velocity of the tracer particles (other than the variation due to Brownian motion). However, when the DOC is not small with respect to the depth of the flow domain, or even exceeds the dimension of the flow domain, there is substantial variation in the tracer velocities (Kloosterman et al. 2011). Such is the case for the flow in a Hele-Shaw cell, for which the distance between the two flat plates h is smaller than the DOC. In this case, the displacement correlation peak (see Fig. 1) is substantially broadened because of the Poiseuille flow profile between the plates. In the figure, displacements ΔX are normalized by the maximum centerline displacement ΔX m . In addition, the detection of the displacement peak leads to a bias in the measured velocity (Kloosterman et al. 2011). As a consequence, a much larger number of frame pairs is needed to obtain a sufficient signal-to-noise ratio (SNR) to detect the displacement correlation peak. In the case of a Hele-Shaw flow cell where the measurement volume includes the full height of the cell and the variation of the velocity spans the full parabolic velocity profile, it is necessary to use much more than 100 frame pairs for the correlation averaging (Ehyaei and Kiger 2014). This makes it very difficult to apply micro-PIV to unsteady microfluidic flows.
In conventional micro-channels the channel width is on the order of the channel height, but in this work we focus on microfluidic flows where the width of the channel is much larger than the channel height. In these cases, the Hele-Shaw condition applies. Such test cases are of interest when studying particle manipulation in a microfluidic chip (Shenoy et al. 2016;Akbaridoust et al. 2018), density-driven flows (Green and Foster 1975;Sebestikova et al. 2007;Almarcha et al. 2010;Ehyaei and Kiger 2014;Vreme et al. 2016), flow behavior of power-law fluids (Drost and Westerweel 2013), magnetic field-driven instabilities (Erglis et al. 2013), Marangoni effects (Kim et al. 2015;Hou and Chu Hong 2015) and flow around a cell (Wereley and Meinhart 2001;Rossi et al. 2009). In practice, for micro-PIV measurements with a large field of view, the DOC is also large (Kloosterman et al. 2011). For small DOC, only the particles in the thin measurement plane are in focus and hence more image frames are required to get an effective image density of 10-20 particles. In this paper, we present a novel approach to micro-PIV that makes it possible to obtain reliable results for the estimation of the flow field in a Hele-Shaw celllike geometry, for which the velocity of the tracers varies significantly over the depth of the measurement domain. Correlation average over the channel height h for 50 frames for uniform flow and Poiseuille flow. The synthetic images were created with particle image size, d = 3 pixels, maximum in-plane displacement ΔX m = 8 pixels It was theoretically shown by Ho and Leal (1974) that the particles migrate to an equilibrium position at a migration length, X = 36 ( ∕a)∕Re c and Re c = V ∕ , where is the channel height, a is the particle radius, V is the velocity and is the kinematic viscosity. For a generic microfluidic condition: = 350 μm , a = 10 μm , Re c = 1, the migration length required for the particles to reach equilibrium is approximately 10 m. This makes it challenging to manipulate the particles into their equilibrium position. In our method, we require less than 10 frame pairs, even considering the effect of uniform particle concentration in the Poiseuille flow. The method is explained in Sect. 2, while in Sect. 3 we validate the method using synthetic data. In Sect. 4, we present an experimental validation of our approach for the flow around a 2D cylinder in a Hele-Shaw cell. Sect. 5 details the experimental setup and results from an investigation of unsteady flow around a developing Rankine half-body. The advantages and limitations of our method are discussed in Sect. 6.

Ψ-PIV principle and workflow
Here we consider the laminar flow between two parallel plates, i.e., a Hele-Shaw flow. The two plates of length l and width w are separated by a distance h along the z-direction (see Fig. 2). The Hele-Shaw condition only applies when the in-plane dimension l is much greater than the channel height h, i.e., l ≫ h . Hele-Shaw condition becomes valid at a in-plane distance of ∼ h from the boundary of the flow domain. The velocity profile is parabolic in the z-direction due to the pressure gradient in the (x,y) plane.
where u is the velocity in the x, y-direction, h is the channel height, z is the wall-normal location in the channel, is the viscosity and p is the pressure gradient. From Eq. (1), the velocity field u can be fully determined from the depth-averaged velocity field. We therefore simplify the notation and refer to u = (u, v) as the two-dimensional depth-averaged velocity field in the (x, y)-plane. Our approach is based on the observation that for a Hele-Shaw flow, the flow direction can be established from fewer frames than the averaged flow magnitude. This is because the flow in Hele-Shaw flow is two-dimensional and at a given (x,y) location, the direction of the velocity of a particle is independent of its location along the z-axis, see Eq.
(1). The flow direction can be estimated at each interrogation position to construct the direction field of the velocity field and in turn, used to reconstruct the instantaneous streamlines. For a Hele-Shaw cell, the depth-averaged velocity field u is a potential flow and the velocity potential can be identified as the pressure field, see Eq. (1). Assuming the depth-averaged velocity field to be a potential flow, one can further deduce the depth-averaged velocity field from the direction field as the velocity field.
In the following, we describe the workflow used to reconstruct the depth-averaged velocity field using the Ψ-PIV algorithm. Figure 3 represents the flow chart associated with Ψ-PIV and is compared to conventional micro-PIV. Similar to conventional PIV, single frames are divided into interrogation sections and image pairs are used to compute the cross-correlation function over each interrogation window. Unlike conventional PIV, the correlation map is not used to find the displacement correlation peak but instead to find a direction correlation peak. For the directional correlation map, single frames of PIV images are divided into interrogation sections as shown in Fig. 3 and the image pairs are used to compute the ensemble correlation averaging (Meinhart et al. 2000) of fewer frames. Next, for each angle, we sum up all the values of the averaged correlation map along a line with its origin fixed at the center of the correlation map. The length of the line is identical to the half-length of the correlation map. We do this by retrieving the intensity values of pixels from the correlation averaged map along the line. The summed-up value along the line is calculated for every angle from − 180 • to 180 • . The pixel intensity value along the line for each angle is calculated by interpolation from the nearest pixel values of the correlation average map. The angle corresponding to the highest peak of the directional correlation gives the most plausible direction of the flow in that interrogation window. The measured flow direction field is used to deduce the depth-averaged velocity field. When the velocity becomes zero, the correlation will have a single peak at the origin of the correlation, and it is no longer possible to determine an unambiguous flow direction. In our correlation analysis, this situation is detected and labeled accordingly. In practice, this is mitigated by the fact that zero velocity is only encountered on stagnation points or separation points, where streamlines end (or would Fig. 2 Schematic of a Hele-Shaw cell with the length l, width w, and height h. The flow is uniform in the x-direction and has a parabolic distribution over the height in the z-direction. The depth-averaged velocity in the x-direction is denoted as u originate). These are limited to small areas in the flow and can be identified from the surrounding streamlines. For our method, we noticed that the ambiguity in the determination of flow direction is least if the average particle image in-plane displacement is around 1/4th of the interrogation window size. Previous studies have shown the particle displacement gradient to reduce the amplitude of the correlation peak and broaden its width. At low effective image density N eff , the broadened correlation peak splits into a series of multiple aligned peaks, which is beneficial for our method as the sum of pixel values of more peaks increases the probability of finding the true flow direction. We first compute instantaneous streamlines from the direction field using a fourth-order Runge-Kutta integration scheme. Here we consider the flow at the inlet to be uniform, and determine the directional stream function. The value of the stream function at the center point of each interrogation window can be deduced from the (instantaneous) streamline pattern. Assuming a uniform flow at the inlet u = (U, 0), we determine the values of the stream function for each streamline at the inlet. To assign the inlet stream function values to each point in the image domain, the streamlines at the inflow side of the image should preferably be parallel to each other, i.e., uniform flow. In that case, it is trivial to assign a value to each of the streamlines. In the case of a non-uniform inflow it would generally be possible to determine an appropriate stream function. For example, the flow velocity and the stream function from a source enclosed by the walls can be computed using a panel method. In general, this does not represent a significant complication, because, in most practical microfluidic application of Hele-Shaw flows, the inflow is generally uniform. The stream function ( Ψ ) in the rest of the domain is calculated at the center of each interrogation window by interpolation from the nearest streamline value. Finally, the velocity components in x-and y-direction are deduced from the stream function: We expect this method to reduce the data acquisition time required to measure the displacement field with the desired accuracy. We proceed to validate this approach with simulations using synthetic images and implement our algorithm for velocity measurements in an experimental study of unsteady Hele-Shaw flows.

Generation of synthetic data
Synthetic images were generated to validate the Ψ-PIV algorithm described in the previous section. Image sequences represent a uniform flow at the inlet, along the span-wise y-direction, and a parabolic Poiseuille flow along the z-direction. The main objective is to assess the reliability of Ψ-PIV to estimate the velocity field. We do not take into account the

Fig. 3
The principal difference between the micro-PIV and Ψ-PIV algorithms. In micro-PIV, correlation averaging of a large number of frames is taken to estimate the velocity field which is shown by solid arrows. In Ψ-PIV, the local direction of the flow is estimated from fewer frames. The direction of the flow is shown as the unit length arrows for each interrogation window. The red lines show streamlines that are determined by advecting virtual particles along the measured directions. The stream function values are obtained by nearest neighbor interpolation of streamlines and are assigned to the center location of each interrogation area (marked as black circular symbols). The velocity marked as solid arrows is calculated from the spatial differentiation of the stream function values with respect to the distance between two consecutive interrogation areas Brownian motion of the particles and the background noise from the image sensor. We are interested in demonstrating the effectiveness of the Ψ-PIV algorithm in the Poiseuille flow with a uniform distribution of tracer particles. The tracer particles are uniformly scattered in the channel across the channel height and inserted at random positions at the inlet to keep the image density constant across all frames. We do not model the defocusing effect in size and intensity of the particle images in the synthetic data, because we are interested in the measurements where the DOC is larger than the channel height. All images are generated with an image size of 1024 × 1024 pixels with 8-bit gray-level quantization. The particles in the synthetic images are well resolved and have a mean particle image diameter of 3 pixels. The centerline displacement for the Poiseuille flow profile was kept at 8 pixels. The seeding density N I was kept as six particles on an average in a 32 × 32 pixel interrogation window. Synthetic images are processed using a single-pass FFT crosscorrelation with 32 × 32 pixel interrogation windows and without overlap. We first consider a flow that is uniform at the inlet and in the x-direction. Figure 4 represents the contour plot of the cross-correlation function and for a particular interrogation window the corresponding correlation function and directional correlation function. Figure 4a, d, g, j represent the correlation function for the same interrogation window computed using a decreasing number of frames. Using 1000 frames, the correlation map displays low noise and presents a single strong signal peak, see Fig. 4a. The location of the maximum peak corresponds to the depth-averaged flow velocity because of the convolution of particle image size with the sampled velocity probability density function. As mentioned in the previous section, this correlation peak is broadened because of the displacement gradient effects, see Fig. 4b. Figure 4c shows the directional correlation of the correlation map for correlation averaging of 1000 frames. In this case, the highest peak is at 0 • because of the strong signal peak in the correlation map (see Fig. 4a). The directional correlation peak is not broadened, because the former only determines the direction and is not affected by the displacement gradient of the flow.
As discussed before, a higher image density N I is achieved in micro-PIV by correlation averaging of several image frames N F . Hence, the noise level (ratio of second highest peak to the highest peak of the correlation map) increases with a decrease in N F . In this case, the SNR in the correlation map decreases because of erroneous correlation averaging between the moving particles and particles attached to the wall. The correlation map clearly shows that the noise due to random correlations increases monotonically as the number of frames reduces from 1000 frames (Fig. 4a) to 50 frames (Fig. 4d), 10 frames (Fig. 4g) and 5 frames (Fig. 4j). For a larger number of frames, superposition of the spatial correlation from each frame is large enough to determine the depth-averaged particle image displacement as shown in Fig. 4b. For less than 50 frames (Fig. 4e), the peak splits into a chain of multiple aligned peaks, see Fig. 4d and g, leading to the erroneous measurement of the displacement. This happens when the particle image diameter d is smaller than the amplitude of the displacement gradient in the flow. For 10 and 5 frames, no clear peak is determined because of the strong measured random correlations as shown in Fig. 4h, k, respectively. Since directional correlation is determined by integrating the correlation map values at each angle Θ from − 180 • < Θ ≤ 180 • , it can be robustly computed for a correlation map with higher random noise. For 50 frames (see Fig. 4d) and 10 frames (see Fig. 4g), all the small peaks within the chain of aligned peaks contribute to the peak in the directional correlation, which can be clearly identified, see Fig. 4f, i. For fewer frames than 10, Fig. 4l shows that the directional correlation peak can no longer be determined because of the strong measured random correlations.

Valid detection probability
Next, we compute the valid detection probability for varying effective seeding density, and compare Ψ-PIV to other PIV approaches reported in the literature, see Fig. 5. The valid detection probability is defined as the probability of the highest correlation peak which corresponds to the true mean displacement. Keane and Adrian (1990) demonstrated that the probability of determining the true displacement peak improves as the average number of particle pairs in the interrogation area increases. This average number of particles within the interrogation area is represented by N I F I F O where N I is the image density for a single image pair, F I is the inplane displacement of particles in an image pair and F O is the out-of-plane motion of the particles. For uniform flow, conventional PIV requires a value of N I F I F O ≈ 8-10 particles to reach a valid detection probability greater than 95% (Keane and Adrian 1992;Adrian and Westerweel 2011). However, in our case, N I F I F O is not an accurate measure of the effective image density N eff , because the flow in the channel is two-dimensional and we use volume illumination. Therefore, the loss of correlation is not due to the out-of-plane component and F O = 1. Rather, it is due to the strong velocity gradient along the z-direction of the Poiseuille flow, and the loss of correlation is better represented by the in-plane displacement gradient F Δ (Westerweel 2008). In addition, we also use N F , defined as number of frames used to average the correlation data from each image pair, to improve the spatial correlation. Therefore, we use N eff = N F N I F I F Δ as the effective seeding density to characterize the valid detection probability. Figure 5 represents the valid detection probability for the synthetic data with different methods. Ehyaei and Kiger (2014) reported a similar figure for their method using 1 3 20 Page 6 of 15 2 frames from a sequence of images ( N F = 2 ). Therefore, we also keep N F = 2 in Fig. 5 for consistency in the comparison between methods. For a uniform flow with a uniform flow profile along z-direction (i.e., F Δ = 1), the probability of finding the true displacement correlation peak approaches 0.95 for N I F I F Δ ≈ 8. For a parabolic velocity profile along the z-axis, a much higher effective image density is required to reach a valid detection probability of 0.95 . In Fig. 5, the valid detection probability of the maximum in-plane displacement ΔX m for a parabolic profile shows poor performance. This is due to the particle size biasing effect, where the maximum peak is located at 1 pixel less than the true displacement value. Ehyaei and Kiger (2014) showed that a higher valid detection probability can be reached. This, however, comes at a price of decreasing the spatial resolution. Ehyaei and Kiger (2014) considered peaks at the locations of the true maximum displacement ( ΔX m ), one pixel less than the maximum displacement ( ΔX m − 1 ) and two pixels less than the maximum displacement ( ΔX m − 2 ) to compute the SNR. Using this methodology, they reached a valid detection probability of approximately 0.94 at N I F I F Δ equal to 100, see Fig. 5. This result was obtained in their investigation using synthetic data ( ΔX m = 15 pixels, d = 4 pixels) for a flow profile along the channel height. We applied the algorithm developed by Ehyaei and Kiger (2014) to our synthetic images and computed the associated . We were able to reach similar values for , see Fig. 5. In our case, the valid detection probability curve is slightly higher reaching 0.95 for N eff = 30. This is because the displacement gradient in our synthetic data was not as high as compared with the algorithm developed by Ehyaei and Kiger (2014). Finally we apply Ψ-PIV, where the valid detection probability approaches 1 for a much lower effective seeding density ( N eff of 20 particle image pairs), because only the direction of the flow needs to be determined. The valid detection probability of 0.95 is reached with N eff = 17; see Fig. 5. We further characterize the performance of Ψ-PIV as a function of the interrogation area and the number of frames. Since the interrogation window size was kept constant at 32 × 32 pixels throughout the analysis done in Fig. 5, we studied the valid detection probability of Ψ -PIV as a function of the interrogation area for the same synthetic data, see Fig. 6. As expected, the valid detection probability exceeds 0.95 when there are at least 22 particle image pairs in the interrogation area. The number of frames N F was equal to 2 in the analysis shown in Fig. 5. Hence, we also determined the minimum number of frames required to reach a valid detection probability of 0.95 by Ψ-PIV. The valid detection probability was calculated for synthetic images with similar properties   Keane and Adrian (1992) as used before, but with lower seeding density (less than about 1 particle image on average in a 32 × 32 interrogation window). This would mimic severe particle loss in the flow domain because of particle clogging at the inlet or particles getting attached to the channel wall. Figure 7 shows that if the interrogation area is 32 × 32 pixels then a minimum of 50 image frames are required to have a valid detection probability of 0.95. This translates here to approximately 50 particle images in an interrogation area. In this case, the higher image density requirement could be the result of having less than about 1 particle image on average in an interrogation window.
We characterize the performance of the Ψ-PIV algorithm with separate experimental setups investigating two different flow cases, namely the flow around 2D cylinder and a Rankine half-body. The details of the experiments are described below.

Experimental setup
In the first measurement, the Hele-Shaw cell is built from two PMMA plates of length l = 300 mm and width w = 100 mm, see Fig. 8. The two plates are separated by spacers height of h = 500 μ m and are clamped on the edges to keep h constant throughout the channel length. A cylinder with radius a = 30 mm and thickness h = 500 μ m is placed at the center of the channel at a distance g = 180 mm from the inlet, see Fig. 8. Syringe pumps (Cetoni neMESYS) were used to generate a flow with a centerline velocity of U = (0.06 mm/s, 0) which yields a Reynolds number Re = 12. The working fluid was water, seeded with polystyrene microspheres with a mean particle diameter of 180-200 μ m (Cospheric). Images were recorded with a CCD camera (LaVision Imager Intense) with a 1376 × 1040 pixel image format, 6.45 μ m pixel pitch, and 12-bit gray-level dynamic range. The camera was equipped with a Nikon objective with a 35 mm focal length and a magnification factor of 0.20. An f-stop of 8 was chosen to have a depth of field of 1 mm, which ensures that all particles within the channel are in focus. We chose a DOC of 10 mm such that the DOC is larger than the channel height. This also ensures that all tracer particles are in focus and contribute to the correlation averaging. A camera exposure time of 2 ms was chosen to allow enough light to enter the image sensor for clear identification of the particle images. The active sensor size was cropped to 992 × 992 pixels to coincide with the flow region. An acquisition frequency of 2 Hz was used to capture the average in-plane displacement of about 1 pixel between two consecutive recordings. To get the maximum particle image displacement without correlation loss due to excessive in-plane displacement, a common practice is to keep the particle image displacement smaller than onequarter of the interrogation window size (Keane and Adrian 1990). Hence, for the interrogation analysis, 8 frames were skipped to increase the velocity dynamic range (defined as the ratio of maximum to minimum resolvable velocity) of the recorded data. PIV images were pre-processed to remove the images of the particles, which were attached to the wall. We used a moving average filter in time over 9 frames to find the local minimum intensity value. This was subtracted from each frame, which effectively removes the background.

Directional correlation peak identification
The proposed Ψ-PIV approach is robust regardless of the direction of the flow. For instance, an interrogation window with an arbitrarily chosen flow direction of − 43 • from Fig. 8 Experimental test setup with all the required components. At the top right corner the schematic of the Hele-Shaw cell (l = 300 mm, w = 100 mm, h = 500 μ m) is shown with the location of the cylinder (g = 180 mm) and the inflow direction from three inlet holes to yield an depth-averaged uniform flow U along the transverse direction of the cell. The diameter of the 2D cylinder is 60 mm. The velocity at the surface of the 2D cylinder along the x-axis is denoted as u x the experimental data of flow around a 2D cylinder in a Hele-Shaw cell (described in Sect. 4.1) was considered for the directional analysis. Figure 9a represents the correlation map for 50 averaged images, where the broadened peak is at an angle of −43 • with respect to the x-axis. The corresponding normalized directional correlation function is shown in Fig. 9b. For the same interrogation window, the top view of the correlation averaged map over nine images is shown in Fig. 9c. In this case, multiple strong peaks are present in the direction of the flow, i.e., −43 • . The corresponding normalized directional correlation function is shown in Fig. 9d. For every interrogation window in the image domain, the flow direction is determined in a similar fashion. Figure 10a, b shows the velocity field obtained from the same experiment using a conventional PIV algorithm with correlation averaging over 217 frames and 9 frames, respectively. In both cases, the interrogation windows had a size of 32 × 32 pixel with 0% overlap. Using 9 frames, the conventional PIV yields an estimate of the flow velocity that is significantly lower than the actual velocity. The velocity vectors are not adequately resolved, because the effective image density is insufficient for correlation averaging. On the central lower side of the 2D cylinder (see Fig. 10b), the flow does not increase and the measurement is inaccurate as a result of cross-correlation among non-corresponding particle image pairs. This is a typical problem with micro-PIV measurements at insufficient image density. Figure 10c presents the velocity field constructed from the same measurement data and using our Ψ-PIV algorithm with 9 frames. The measurements are in close agreement with the PIV results obtained from a correlation average over 217 frame pairs (Fig. 10a).

Measurement results
The velocity field measurement can be compared with the analytical solution for the flow around a cylinder expected from potential flow theory. Figure 11 reports the theoretical velocity field in the x-direction at increasing distance from the cylinder surface along the y-axis and the corresponding measured data using PIV and Ψ-PIV. The theoretical solution imposes a no-penetration boundary condition at the surface of the cylinder and allows for a slip velocity.  The velocity on the surface of the cylinder along the x-axis (see Fig. 8) increases to a velocity u x = 2U at the equator of the cylinder. The PIV result for correlation averaging over 217 frames shows that the velocity field is in close agreement with the analytical solution except near the surface of the cylinder. Close to the surface, the velocity decreases because the Hele-Shaw condition is no longer appropriate near the wall and the potential flow solution does not hold because of the no-slip boundary condition. The PIV result for correlation averaging over 9 frames clearly underpredicts the flow field due to the low effective image density. It should be noted that the measured free stream velocity for PIV with 9 frames is less compared to PIV data with 217 frames. This is because, with fewer frames, particles that are attached to the wall significantly contribute to the underestimation of the velocity field. Near the surface of the cylinder, the velocity measured from conventional PIV with correlation averaging over 9 frames is largely erroneous with a difference of up to 37%, owing to the large in-plane displacement and displacement gradient. The Ψ-PIV result for correlation averaging over 9 frames is in agreement with the PIV result for 217 frames. For Ψ-PIV, the difference in measured velocity with respect to the reference velocity is approximately 8% throughout the velocity distribution. Ψ -PIV substantially reduces the requirement in the number of frames necessary for the correlation averaging to reach the desired accuracy. The accuracy of Ψ-PIV is reported in Fig. 12 which displays the root-mean-square error (RMS error) of the measured velocity data using PIV and Ψ-PIV  Fig. 10) from PIV, Ψ-PIV next to the cylinder and compared with the theoretical distribution which is calculated assuming ideal flow conditions. The yellow band represents the region where the Hele-Shaw approximation is not valid ( l ≈ h) in the x-direction at increasing distance from the cylinder surface along the y-axis (marked as the rectangular domain in Fig. 10). For a particle displacement of 8-10 pixels, PIV with 217 frames has a measurement error of 0.05 pixels or 0.5% of the measured displacement. The RMS error for PIV reaches 0.05 pixels for the correlation averaging of 70 frames and remains constant thereafter. With fewer frames, the RMS error of Ψ-PIV is small compared to PIV, because the flow direction can be determined from the ensemble correlation averaging of fewer frames. The Ψ-PIV algorithm reaches a RMS error of 0.05 pixels for the correlation averaging of 9 frames and subsequently remains stable. Thus, for the same PIV dataset, the results from the Ψ-PIV algorithm converges with fewer frames (8 times fewer frames in this experiment) compared to the PIV algorithm. Both the algorithms use a single-pass cross-correlation to present the underlying difference between the two approaches.

Experimental setup
We characterize the performance of the Ψ-PIV algorithm for unsteady flow in a Hele-Shaw cell. In this second experiment, the Hele-Shaw cell has dimensions l = 60 mm and w = 15 mm. The height (h) of the channel is 350 μ m. The flow channel was made of PDMS using standard microfabrication techniques (Xia and Whitesides 1998), and the PDMS layer was bonded to a glass coverslip using plasma oxidation.
The pressure pumps (Fluigent) were used to generate a uniform inflow with a centerline velocity U = (0.2 mm/s, 0). A Rankine half-body is generated in the channel through an additional inlet located at the middle of the channel, see Fig. 13. Standard PTFE Teflon tubing with 1/16 inch outer diameter and 1/32 inch inner diameter were used to connect the Hele-Shaw cell and the pressure pump reservoir. Polystyrene microspheres with a mean diameter of 20-27 μ m (Cospheric) were used as tracer particles. Images were recorded with a CCD camera with similar specifications as mentioned in Sect. 4.1. The camera was equipped with a Nikon objective with a 200 mm focal length, a magnification factor of 0.80 and an f-stop of 8. The computed depth of field was 0.8 mm, which ensured that all particles within the channel were in focus. The DOC is 1.7 mm. The acquisition frequency was 9.8 Hz corresponding to an average in-plane displacement of 2-3 pixels between two consecutive recordings. We generate an unsteady flow in the channel, that corresponds to the following sequence: (1) a uniform flow from t = 0 to t = 74.4 s, (2) a developing Rankine half-body from t = 74.5 s to t = 77.6 s due to an increasing volumetric flow rate through inlet 1, (3) the flow around a developed Rankine half-body generated by a volumetric flow rate Q 1 = 0.4 ml/min at the inlet 1 from t = 77.7 s to t = 240 s. We focus on the time interval between t = 66 s and t = 80.5 s where the flow gradually develops from uniform flow to a fully developed Rankine half-body. The image processing was done using interrogation windows of 32 × 32 pixel with 0% overlap. We used a moving average filter in time over 9 The two additional inlet holes marked as inlet 1 and 2 can be used to generate a source or a sink frames to find the local minimum intensity value. This was subtracted from each frame, which effectively removes the background. Figure 14 compares the streamlines obtained for the same experimental recording and using the velocity field obtained (1) from PIV using correlation averaging over 100 frames (Fig. 14a-c), (2) from PIV using correlation averaging over 3 frames (Fig. 14d-f) and (3) from Ψ-PIV using correlation averaging over 3 frames (Fig. 14g-i). Figure 15 compares the unsteady velocity measured at a given point using conventional PIV and Ψ-PIV. The experimental results are compared with the analytical solution at the location marked as 'x' in the flow domain (Fig. 14). From potential flow theory, the velocity was calculated using a superposition of uniform flow and a single source. Similar to the experimental data, the analytical solution (see dashed line curve in Fig. 15) was kept as a uniform flow for the first 74.4 s (i.e., frames 1-80). The Rankine  half-body was developed by increasing the strength of the source from t = 74.5 to 77.6 s (i.e., frames 81-110), and thereafter the source strength remains constant until the end of the measurement at t = 80.5 s (i.e., frames 111-140). At time t = 66 s (frame 1), the uniform flow is captured well by PIV with correlation averaging over 100 frames (Fig. 14a). At t = 75.5 s (frame 90), although the streamlines are continuous, they do not correspond to the true flow behavior. This is because the temporal resolution is lost due to large number of frames required for the correlation averaging. In our experiments, a full Rankine halfbody develops in 3 s, but the PIV results shown in Fig. 14b take the temporal information of 10.2 s to have enough effective seeding density to find the velocity vectors. Here we see PIV results of a fully developed Rankine half-body, because the PIV, in this case, takes particle image information starting from 75.5 s (frame 90) to 83.7 s (frame 171). This means the correlation averaging also takes 70% of the particle image information from the fully developed Rankine half-body. At t = 78.5 s (frame 120), Fig. 14c shows a fully developed Rankine half-body that also corresponds to the actual flow, because the flow around the developed Rankine half-body is steady for the entire experiment. The velocity field measurement for correlation averaging over 100 frames have low random error but this leads to a notable loss of temporal information and thus PIV is not able to capture the gradual change in the flow rate as the Rankine half-body is developing in the flow field. This can be distinctly seen in Fig. 15 from the significant deviation between the PIV measurements and the theoretical solution. Figure 14d-f show the streamlines using PIV velocity fields computed from correlation averaging over 3 frames. These streamlines are discontinuous and in some cases break midway. This shows that the velocity vectors are erroneous due to lower effective image density. The development of the Rankine half-body becomes visible from the streamlines in Fig. 14e, even though the measured velocity field is inaccurate. For Ψ-PIV with correlation averaging over 3 frames, the streamlines before (at t = 66 s; frame 1) and after (at t = 78.5 s; frame 120) the fully developed Rankine (see Fig. 14g, h) are in agreement with the streamline results of PIV with correlation averaging over 100 frames (see Fig. 14a, c). At t = 75.5 s (frame 90), we can clearly see the evolution of the Rankine half-body from the streamlines in Fig. 14h. In this case, Ψ-PIV has a temporal resolution of 0.3 s with enough effective seeding density to capture the gradual change in the flow field. The velocity estimated using Ψ-PIV is in close agreement with the theoretical solution as the Rankine half-body develops, i.e., from frames 81 to 110, see Fig. 15. For Ψ-PIV, the velocity slightly overshoots the theoretical velocity value for frames 111-130. This may be either due to the control response of the flow rate sensor for the pressure pump or due to a spurious direction vector from the Ψ-PIV algorithm or a combination of both. From frame 131 onward, the results from Ψ-PIV are again in agreement with the theoretical result.

Discussion
The valid detection probability of synthetic data as a function of effective image density N eff shows that a higher SNR can be achieved by Ψ-PIV compared to conventional PIV. Ψ -PIV attains = 0.95 at N eff = 15, see Fig. 5. For ≥ 95%, the minimum effective image density represented by the Ψ -PIV algorithm is less than half compared to the valid detection analysis of tracer particles spread uniformly across the channel height using PIV as shown by Ehyaei and Kiger (2014). Figure 11 displays a comparison between PIV and Ψ -PIV result for the experiment of flow around a 2D cylinder in a Hele-Shaw cell. The Ψ-PIV result for correlation averaging over 9 frames is in close agreement with the PIV result for correlation averaging over 217 frames. Thus, the temporal resolution of Ψ-PIV increases by a factor of 25 compared to conventional PIV. Figure 15 shows that for an unsteady flow, conventional PIV is unable to measure the velocity accurately, because correlation averaging is required over a large number of frames corresponding to a time interval greater than the time scale of the flow. This makes it very challenging to measure the flow field using conventional PIV for unsteady flows. For such flows, Ψ-PIV can obtain higher temporal information compared to conventional PIV for the same number of frames, see Fig. 15. Although Ψ-PIV is an improvement over the current correlation averaging technique, it still requires an effective image density of around 20 particle image pairs. This means that the average effective seeding density needs to be around 20 particle image pairs in 32 × 32 pixels interrogation window to deduce the velocity field from the directional correlation over two consecutive frames. If the image density N I is lower than 20, then the correlation averaging over a few frames is required to reach effective image density ( N eff = N F N I F I F Δ ) of 20.

Conclusion
This paper describes a new algorithm to determine the velocity fields in a Hele-Shaw cell. This method reduces the minimum required effective image density ( N eff = N F N I F I F Δ ) compared to the conventional micro-PIV technique using correlation averaging. This increases the temporal resolution that can be achieved by Ψ-PIV compared to conventional PIV. Ψ-PIV is therefore attractive to measure unsteady flows for microfluidic applications. The major difference lies in the fact that Ψ-PIV requires a lower image density to determine the flow direction for each interrogation window compared to the conventional method. Once the flow direction is determined, the two-dimensional stream function is used to extract and reconstruct the magnitude of the velocity field. Synthetic image evaluation for uniform particle concentration shows that an effective image density of 20 particle image pair is satisfactory to measure the velocity field. This is 5 times lower compared to the required image density in the measurement of the velocity gradient within the correlation depth (Ehyaei and Kiger 2014). Here, we present a steady measurement case of the flow around a 2D cylinder in a Hele-Shaw cell. The Ψ-PIV algorithm using correlation averaging over 9 frames yields similar results as the PIV algorithm using over 70 frames for the single cross-correlation approach. Moreover, measurement results of a developing Rankine half-body for Ψ-PIV can reach significantly higher temporal resolution. This reduction in the number of frames for the correlation averaging in Ψ-PIV would enable a substantial improvement in the temporal resolution in case of time-varying microfluidic flows. In addition, the minimum image density required to determine the direction field could be further improved with the use of advanced PIV processing steps such as multigrid and iterative windows approach (Scarano and Riethmuller 2000).