Holographic PIV/PTV for nano flow rates–A study in the 70 to 200 nL/min range

Experimental setup

The microfluidic channel was milled into a polymethyl methacrylate (PMMA) block measuring 40 × 8 × 77 mm (W × H × L). The channel itself was designed to be 20 mm long with cross-section dimensions of 1,000 × 30 µm (W × H), which corresponds to an area of 30,000 µm². The height of the channel is relatively shallow in order to keep the particles close to the sensor to increase their visibility. After fabrication, the channel was polished with plastic polish paste (Xerapol) to remove milling ripples. The ripples, if not removed, diffract too much light, worsening the particle visibility (Figure 3). The precision of our method relies on accurately knowing the cross-sectional area of the channel, therefore, before running the experiments, we optically measured its geometry with focus-variation microscopy (Confovis TOOLinspect, magnification 20x, 275 nm lateral resolution, 300 nm depth resolution). By averaging the area of 2000 cross-sectional profiles, we saw that the channel’s area was in fact 27,893 µm², 7% smaller than expected (Figure 4). The channel was sealed by applying a microfluidic adhesive tape (3 M Diagnostic Tape 9795R) on the milled side of the PMMA block. The tape was 100 µm thick and consisted of a clear polypropylene film coated with pressure sensitive silicone-based adhesive.

Figure 3:

The effects of surface polishing on channel visibility. A) A channel before polishing. B) The same channel after polishing (Images acquired on a lens-based microscope. The channel shown here is not the one used in our measurements).

Figure 4:

Measurement of the channel’s real geometry using focus-variation microscopy. A) 3D point cloud obtained by the microscope. B) Average channel cross-section (an average of 2000 parallel profiles spread over a channel distance of 700 µm). The spikes near at the border of the channel are not physical, they are artefacts of the imaging method due to reflection on the edges. Also note that the axes are not to scale.

The image sequences were recorded at 10 frames per second using a CMOS sensor (iDS UI-3592LE-C Rev.2). The sensor has 4,912 × 3,684 pixels, where each pixel is a square with a side length of 1.25 µm. Therefore, it was possible to image the full width of the channel and a length of 6.14 mm. The sensor provides RGB images, but they were converted to gray-scale because LDHM is intrinsically monochromatic, therefore the holograms should have a single light intensity channel. The sensor is positioned roughly at the middle of the channel’s length and directly under it. The air gap between the sensor and the channel is less than one mm wide. An 80 mW blue laser diode (PL450B, λ = 450 nm) was used as the light source.

The flow was generated and controlled by a syringe pump (Nemensys Low Pressure, NEM-B101-02 B, 14:1 gear). The syringe used (SETonic P/N 3,010,307) had a total volume of 1 mL for a 6 cm stroke allowing a minimum flow rate of 70 nL/min. The pumped liquid was purified water, seeded with polystyrene (PS) beads measuring 5 µm in diameter, which is large enough to avoid Brownian motion. Before each experiment, the mixture was agitated to make sure that the beads were uniformly distributed in the water. The beads had a density of 1.05 g/cm³, which gave them the tendency to slowly sediment. As shown in Figure 2, the channel was mounted such that the flow is in direction of gravity, therefore this sedimentation has to be accounted for when calculating the flow velocity.

A total of 50 image sequences were recorded, each 12 s long. The sequences were roughly equally distributed among seven different flow rates: 0, 70, 100, 125, 150, 175, and 200 nL/min. To minimize any transient effects, we waited 30 min before recording a new sequence, whenever the pump’s flow rate was altered. The sequences were stored on a hard drive and the data analysis was performed offline.

Simulation

Besides collecting experimental data, we also generated synthetic data by simulation. The product of the simulations were particle trajectories inside a rectangular flow channel. Many sets of trajectories were generated by imposing different total flow rates. The trajectories were then converted to synthetic holographic videos. The objective here is to quantify how well our measurement method is able to estimate the imposed flow rate when processing the synthetic holographic videos. More precisely, the synthetic data was used to estimate the combined standard uncertainty of our measurement method. In comparison to the experimental data, the synthetic data is a reference standard with negligible value uncertainty (u _ref ). Which implies that fluid velocity is constant and its value equal to the one imposed to the simulation. The same is not true for the experimental data, for which the real flow rate value can diverge significantly from the flow rate setpoint of the syringe pump for reasons such as the pump’s dosing precision or flaws in the setup (for instance, elasticity in the tubings and leaks in the connections). By using synthetic data, we factor out the measurement error caused by the instabilities of the flow rate during a measurement.

Another advantage of synthetic data is that, once the simulator is built, it is easy to generate high volumes of data. In total 600 videos were generated, 100 for each of the investigated flow rates (70, 100, 125, 150, 175, and 200 nL/min).

To generate the synthetic data, we combined a Computer Fluid Dynamics (CFD) simulator, with a hologram generator. Using the CFD simulator we calculated the flow profile inside the channel for a given flow rate (boundary condition), with which we calculated the trajectory of particles inside the channel. Finally, using the hologram generator, we rendered these trajectories as holographic videos.

The CFD simulations were done using OpenFOAM’s icoFoam [11], which is a single-phase, incompressible, laminar, transient solver. The simulated geometry corresponded to our channel design (1,000 × 30 µm cross-section). The mesh was orthogonal, made only of rectangular cuboids cells measuring 2 × 2 × 20 µm, where the largest dimension is parallel to the flow direction. The walls of the channel received no-slip boundary conditions. The simulation was conducted for 20 ms, which was enough to bring the flow to a steady state. The result of the simulations are velocity profiles as the one shown in Figure 5.

Figure 5:

Velocity profile on the cross-section of the channel for a flow rate 100 nL/min. x and y axes are not to scale.

The next step is to calculate the particle trajectories. Approximately 230 particles are visible in each frame, which is similar to the amount in the experimental data. For each new particle that entered the channel, x and y coordinates were randomly sampled from uniform distributions. Using the x-y coordinates, we can get the velocity v_z from the simulated flow profile. The other velocity components (v_x and v_y) are constant and equal to zero, due to the fact that the particles are perfectly following the uniform laminar flow, as observed experimentally at low Reynolds numbers [12]. The simulation recalculates the position of each particle every 0.1 s, which matches the frame rate used in our experiments. From these sequences of positions the particle trajectory is determined.

To transform the trajectories into a sequence of holograms, we used the light scattering routines provided by the Python package HoloPy [13]. The particles were modeled as spheres with a diameter of 5 µm and refractive index of 1.6 (Polystyrene). The holograms have a pixel size of 1.25 µm and dimensions 800 × 2,187 pixels. These specifications also match the conditions of the experimental setup in terms of tracer particle and image sensor. To increase realism, the holograms are then corrupted with random speckle noise and blurred with a Gaussian filter. The result is shown in Figure 6.

Figure 6:

An example of synthetically generated hologram. For visualization purposes, just a piece (800 × 800 pixel) of the frame is shown here. The real frame is larger (800 × 2,184 pixel).

Analysis

We used Python 3.9 to develop an analysis tool that can process the videos acquired by the sensor and estimate the flow rate in the channel. This tool can be broken down into five sequential modules (Figure 7).

Figure 7:

An example of a video frame after being processed by each module of the analysis tool. For visualization purposes, just a piece (800 × 800 pixel) of a frame is shown here. The real frame is larger (800 × 2,184 pixel). The modules are: A) Conversion to gray-scale, B) crop and background subtraction, C) holographic reconstruction, D) particles segmentation, E1) PIV. Due to the particle sparsity, it is not possible to calculate the velocity for all image patches. Those appear black in the image, E2) PTV. Some trajectories contain zigzags. This may happen when two particles move side-by-side very close to each other. In this situation, the algorithm may fail to differentiate them. Fortunately, this has no strong impact on the final flow rate calculation, because it is relatively rare, and it has minimum effect on the velocity aligned with the flow direction, which is the information that is used to calculate the flow rate.

The first module reads the video frames, converts them from RGB to gray-scale, and, if necessary, rotates them so that the flow is in the z-direction.

As shown in Figure 7A, the images are quite cluttered by the channel’s imperfections and impurities attached to the sealing tape. In the second module of the analysis tool, those artifacts are removed by subtracting from the frames an image where no particles are present. This background image is generated by calculating the median between all frames in the video. After the background subtraction, only the moving particles (and a certain level of noise) are visible on the images, as shown in Figure 7B.

After the background subtraction, the interference patterns are sufficiently visible on the frames, and we can proceed with the holographic reconstruction. For the reconstruction, we relied on the Python package Holopy [13]. The result of a reconstruction is an image where each pixel is a vector with magnitude and phase. However, we discarded the phase, because just the magnitude (light intensity) is relevant to the velocimetry methods. For each hologram, we performed six reconstructions with different focal distances. The final image was the pixel-wise minimum of the six reconstructions. This ensures that all the particles, regardless of how far they are from the image sensor, will be depicted as sharp black dots as shown in Figure 7C.

The next module is responsible for segmenting the images that are the result of the holographic reconstruction. Here, segmentation means to generate a binary image where particle pixels have value one, while all the others have value zero, as shown in Figure 7D. For that, we first applied a median filter with a 3 × 3 kernel to minimize speckle noise, then used Yen’s thresholding method [14] to generate the binary image.

Once the frames have been segmented, we can run a particle velocimetry technique (PIV or PTV) to determine the flow rate. PIV (Figure 7E1) operates on the principle that particles that are close together move with similar velocities. Therefore, it divides each frame into many small patches and calculates the cross-correlation of each patch with its corresponding patch in the subsequent frame. The argmax of the cross-correlation provides an estimate of the particle displacement of that specific frame patch. By averaging each of the patch displacements, across all the frames, it is possible to calculate average particle velocity in an image sequence. In a video with n frames, where each frame F is divided into k rows and i columns, the average particle velocity v can be estimated as:

Equation (1):

Where Δt is the time interval in seconds (1/frame rate), and pixel_size is the sensor’s pixel size in µm. S(t, x, z) is the displacement estimated at frame t for the patch located at row z and column x:

Equation (2):

The cross-correlation (represented by the symbol ∗ ) is generally performed on patches of different sizes; the patch at frame t, often called search area, being larger than the patch t-1, often called interrogation window. This means that the output of the cross-correlation will generally be a 2D function. However, in our PIV implementation, we constrained the width of both patches to be the same; therefore, the output of the cross-correlation is simply a 1D function. In this way, our PIV method only measures displacements along the z-direction (the flow direction). Such a constraint reduces computational costs and reduces the occurrence of spurious displacement estimates due to particles being wrongly correlated to their lateral neighbors.

PTV (Figure 7E2), on the other hand, relies on a very different principle. PTV tries to find the white spots in subsequent frames that represent the same particle. Once this process is done for all frames, we can trace the trajectory of the individual particles. From the length and time stamps of these trajectories, one can calculate the particle velocity. Finally, the average flow velocity v in z-direction can be estimated as the mean particle velocity. Our PTV analysis relies on the package TrackPy [15], which is a Python implementation of the algorithm by Crocker and Grier [16].

For both velocimetry methods, the average flow velocity is calculated as a simple mean of all measured particles (PTV) or patches (PIV) velocities. The simple mean gives a valid estimate because the particles are fairly uniformly distributed in space. Such a distribution has also been observed in other studies working with low Reynolds flows [2, 12].

Finally, the flow rate in nL/min is calculated as:

Equation (3):

Where v is the average flow velocity in µm/s, v _sink is the average particle sedimentation velocity in µm/s (in z-direction, since the channel is mounted vertically), and ch_area is the channel area in µm². We estimated v _sink by evaluating Equation (2) for image sequences with no flow (syringe pump turned off). On average, particles sedimented at 15.6 μm/s.

Measurement uncertainty

Using the top-down approach described by G. White [17], we calculated out method’s expanded relative uncertainty (U _% ):

Equation 5:

where:

1. k: coverage factor. We adopted k = 2 (95.4% confidence interval).

2. q _m : average measured flow rate.

3. u _c : combined standard uncertainty.

4. u _imp : imprecision due to random effects. It is given by the standard deviation of the measurements (σ_q).

5. u _bias : uncertainty in the measurement bias.

6. u _ref : uncertainty in the reference flow rate. Since this data was generated by simulation, this uncertainty is very low and was neglected in the calculations (u_ref = 0).

7. u _rep : uncertainty in the average measurement (q _m ) value. Calculated as the standard error of the mean (σ _q /√n, where n is the number of measurements).

The uncertainty numbers are shown in detail in Table 1 and summarized in Figure 9. Both velocimetry methods performed similarly well. Their relative uncertainty remained under 6% across the investigated flow rate range. However, PTV performed slightly better. Its uncertainty stayed between 3.9 and 4.9%, while for PIV the values ranged from 4.2 to 5.6%.

Figure 9:

Expanded relative uncertainty for a coverage factor of 2.

Experimental data

Comparing the predictions of the velocimetry methods on the experimental data (Figure 10) one can see that the predictions of both methods only differ by a few fractions of a pixel; 0.082 pixel/Δt on average. Such a “cross-validation” between the methods is indicative that the velocimetry modules are functioning as required.

Figure 10:

Results on the experimental data. A) PTV vs PIV estimates. Each blue dot corresponds to one of the 50 measurements in the dataset. The red-shaded region marks where the deviation is less than one pixel/Δt. The two methods predict very similar values. The deviation is always in the sub-pixel range. B) Flow rates obtained with PIV vs. syringe pump setpoint. C) Flow rates obtained with PTV vs. syringe pump setpoint. Figures B and C use standard boxplots.

Figure 10B and C compare the estimated flow rate with the setpoint used by the syringe pump. We observe a good average agreement between these quantities; the MRE was 9.1% for PIV, and 8.8% for PTV. Still, two outliers are present and the errors for the lowest reference flow rate (70 nL/min) were considerably higher. Possible causes of these errors are presented in the next section.