Particle image based simultaneous velocity and particle concentration measurement

The aim of this study is the expansion of the application of particle image velocimetry (PIV) to include the determination of particle concentration within the visualized area, in addition to velocity analysis. The assessment of particle concentration is valuable in various lab-scale experiments involving particle dispersion. Additionally, it plays a crucial role in evaluating the quality of PIV images. The research investigates two particle image-based concentration techniques: the exponential averaging-based sliding method and the Voronoi cell-based method on the particle images. The exponential averaging method provides a straightforward approach, utilizing a constant length scale for sliding average application to particle images. However, this method may result in broadened interfaces or a ‘marker-shot’ effect at low concentrations, making it less suitable for scenarios involving highly non-uniform particle distributions, such as concentrated jet emissions into ambient environments. Consequently, detecting interfaces in such cases requires additional effort for reliable results. In contrast, the Voronoi cell-based technique offers the advantage of spatially adaptive resolution, making it well-suited for variable concentration distributions and situations where interface detection is crucial. To comprehensively evaluate the performance of these techniques, a synthetic test case was generated to simulate a diffusion problem featuring an initial step in concentration distribution. Both the exponential averaging and Voronoi cell-based methods were applied and compared using this synthetic test case. Additionally, the effect of particle–particle overlap is analyzed theoretically and experimentally with uniform concentration and comparison with particle counter measurements. A modified Voronoi method is introduced, providing flexibility in capturing a wide range of concentration regions and features. An example experimental scenario involving a turbulent puff was considered demonstrating the application of the developed methods. The results demonstrate that the Voronoi method effectively captures small structures with high concentrations while providing reliable results in regions with low concentrations.


Introduction
In many fluid mechanical systems, the concentrations of the species are time and/or space dependent.Dispersed phase analysis is important in the context of solid/liquid particles in liquid/gaseous flows for analyzing the flow fields as well as heat and mass exchange rates in a variety of systems (Nishino et al 2000).Concentration field measurements are useful to understand the spatial distribution of the particle concentration in various scenarios including detection of boundaries of the species or mixing effectiveness, to name a few (Gnirß and Tropea 2008).Typical methods involve a point measurement sensor system like particle counter.Point measurement techniques, however, encounter difficulties when there is need for instantaneous spatial information on physical properties such as velocity and concentration.Additionally, most point measurement systems like particle counters sample part of the air to count the particles.Therefore, it is intrusive to the system being measured and can only provide an average in longer duration of the test.
For many scenarios involving aerosol dispersion analysis like air pollution (Saïd et al 2005, Zhang et al 2020, Newland and Woods 2023), bio-aerosol transmission studies (Shah et al 2021, Kähler et al 2023), indoor air quality (Tang et al 2020, Kähler et al 2023), a thorough understanding of flow features, direction, and fluctuations in species concentration is imperative to assess air quality.Mean point-wise properties would not be sufficient; instead, comprehending the spatial and temporal variations in concentration becomes essential.Moreover, it is crucial to determine if the sampled air is representative of the sensed region (Vincent 2007).Achieving this requires a quantitative understanding of aerosol dispersion characteristics in flows, enabling a better grasp of the chances of successful aerosol sampling.Remote sensing of aerosols also necessitates this knowledge to optimize detection capabilities.By gaining insight into the dispersion behavior of aerosols, one can effectively maximize their ability to sample aerosols and gather valuable data for further analysis and decision-making.
Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are instantaneous whole-field fluid velocity measurement techniques based on imaging of particle tracers.These methods make it possible to detect spatial flow structures and provide information of the spatial differential quantities of turbulence.In this work, the goal is to determine particle concentration within the field of view in addition to the velocity information.This leverages data from preexisting flow tracers, as opposed to employing tracer gases.Furthermore, in addition to concentration measurements as detailed, the particle concentration within the visualization window is also important for assessing the quality of the PIV images and the resultant velocity field.
An image-based method proposed previously is based on intensity of the scattered laser light from particles (Gnirß and Tropea 2008).However, the method involves many corrections to obtain physical concentrations of which uneven illumination and stray reflections pose difficulty to correct reliably.Another approach utilizing image intensity information is the autocorrelation of the image acquired (Nguyen et al 2012, Warner andSmith 2014).The peak of the autocorrelation is related to the intensity within the interrogation window and can be utilized to derive particle concentration within the window.While these methods are simple to apply, they still depend on reliable background subtraction and appropriate interrogation window selection.
In the present work, methods based on individual particle images are discussed.One of the methods is the exponential sliding average method for concentration determination.It needs to be emphasized that the length scale is fixed once chosen and cannot be varied within the image even if the concentration of the particle images is drastically different.Another method introduced here is based on Voronoi cells, which in principle is an adaptive grid procedure to calculate the local particle concentration.Such method has been demonstrated for application in the cosmological application (Schaap and Van De Weygaert 2000).However, for fluid mechanical measurements only few such considerations are available.Monchaux et al (2010) analyzed clustering of particles in turbulent flows utilizing Voronoi cells.Weber et al (2023) applied Voronoi cells for evaluation of bubbly flows.One method for interface determination in multiphase flows is presented by Li et al (2021) where a complementary to Voronoi cellsthe Delaunay meshing technique-is applied.Some methods discussed above (Schaap andVan De Weygaert 2000, Li et al 2021) focus in getting the sharp features present in the respective applications: cosmological and multiphase flows.However, when one is interested in a particle concentration field encompassing vastly different scales, it is not straightforward to apply the above methods.This consideration, which has not been explored previously, holds significance for practical inquiries related to particle cloud dispersion or emitted particle plumes.In this work, methods based on PIV particle images for particle concentration determination for fluid mechanical applications are presented and analyzed in detail.Particular focus has been applied to Voronoi cell method, which could be applied for variety of real scenarios involving drastically different concentrations like pollution studies, indoor air quality to name a few.

Particle concentration measurement
This study primarily focuses on the particle image based concentration measurement based on two techniques: (a) sliding average method; (b) Voronoi cell based method.

Sliding average method
The sliding average method is a technique that averages the dataset by considering a moving window of data points.Here, an exponential moving average type filter, which applies a local filter to compute the average particle per pixel as below is utilized (Lukin 2007),  Υ(i) is the intensity at a certain pixel i and Υ avg is the computed intensity value at pixel i and i − 1, accordingly.l is the length scale in pixel chosen for the average computation.The same method is applied in all four directions one after the other to compute the sliding average value of the intensity values (Lukin 2007).For simplicity, the technique is referred to as sliding average method from this point forward.
Figure 1 illustrates the results obtained from the sliding average method with a single particle located at position (0, 0).Prior to applying the sliding average, the value at this location is 1, indicating one particle at the center pixel (1 particle per pixel, 1 ppp).It is assumed here that the particle images are available as binary representations (ones and zeros) to focus on the technique.Further details on how these binary images might be obtained from PIV images are detailed later in the experimental section.The method effectively redistributes the particle into all four directions after passes in all directions as given by equation (1).For example at length scale l = 4 px, the peak value is C I ∼ 0.02 ppp at the center and progressively decreases to 0 as the distance from the defined length scale increases.Notably, the peak value diminishes with an increase in the averaging length scale, while the influence of the particle extends further as the length scale grows.This still maintains the sum of the values in the domain to be 1.
Figure 1 reveals the importance of choosing an appropriate length scale for the calculation of the sliding average method.On the one hand, there must be enough particles within the chosen length scale to reveal the average particles in the window.On the other hand, the difference in the value must signify the change in concentration rather than the effect of a single particle.

Voronoi cell based method
The Voronoi cell surrounding a point (or particle) defines a polygonal region encompassing the points closest to it and forming a closed network with neighboring cells without any overlaps or gaps (Okabe et al 2000, Schaap andVan De Weygaert 2000).An illustrative example of Voronoi cells with uniformly spaced particles with different spacing (concentration) in either side of origin is shown in figure 2(a).Figure 2(b) depicts the Voronoi cells for the case with macroscopically same concentration; however, the particles are randomly distributed as expected in practical case with Brownian motion.
The inverse of the cell area technically represents the number of particles per unit area, which can be interpreted as particle density or concentration.This can be visualized as true for the perfect case as shown in figure 2(a) which is similar to calculation of the unit cell as done in material science (Callister and Rethwisch 2018).However, in the context of fluid mechanics, the distribution of particles is less ordered as shown in figure 2(b) with Brownian diffusion.Hence, the above measure has to be rewritten to necessarily reflect concentration in the traditional continuum sense and not focus on individual particles.
To address this, one can approximate the continuum by either considering a sufficiently large number of particles or selecting an appropriate examination window size.Here, we consider K number of nearest neighbor cells to calculate the average area of the K cells and the inverse of this average area is the local concentration centered at the particle considered.
In addition, one is also interested in approximating the local concentration when there are few particles in the domain.This is particularly important when one needs to distinguish the low concentration and high concentration regions either in time or space and not averaging the entire set for sufficiently high number of particles.
An alternate approach to Voronoi method is to utilize the Delaunay triangles for the calculation of the concentration.However, it was found that the size distribution of the Delaunay triangles for practically uniform cases varied significantly, resulting in a broad distribution compared to the Voronoi cell sizes.This results in the need of a larger number of triangles for convergence to macroscopic concentration.The distribution obtained is presented in supplementary information.In this work, only the Voronoi method is explored further.

Test case
For comparison of the different methods, a synthetic case of a diffusion problem is considered.A domain extending from −10 000 px < x < +10 000 px and 0 < y < 1000 px is considered.An example multiplication factor considered here is L = 1 mm = 10 px representative of realistic scenarios.The domain is closed by walls in x-direction.At x = 0, there is a removable wall, which is removed at t = 0 s.A periodic boundary condition is considered along the y-axis.Before removing the wall, uniformly placed particles with C 1 = 10 −4 ppp and C 2 = 0.01 ppp are considered at left (x < 0) and right (x > 0) of origin, respectively (figure 2(a)).Such concentrations are representative concentrations of cases when one can neglect particle-particle overlap, the effect of which is discussed later.Considering a concentration scale C 0 = C 2 = 0.01 ppp, the initial concentrations can be rewritten as C 1 /C 0 = 0.01 and C 2 /C 0 = 1.The diffusion process is modeled with applying diffusion jumps to each particle in the domain.
The diffusion jump ∆r of a particle in a time step ∆t is given by (Daune 1999, Sankaran et al 2020): where R is a random number between 0 and 1 defining the probability of the diffusion jump of length ∆r.This jump ∆r occurs in random direction θ with 0 ⩽ θ ⩽ 2π .A diffusion coefficient of D = 10 −5 m 2 s −1 is assumed here and a time step of 1 s is considered (where the physical scale 1 mm = 10 px is considered as before).
The diffusion process is run for extended time before time t = 0 s resulting in a practically uniform concentration on both sides of the origin, as depicted in figure 2(b).At t = 0 s, the wall separating the two concentration side is removed.
For time scales t ≪ L 0-w 2 /D (where L 0-w = 1000 mm = 10 000 px is the distance from center to the side walls), the concentration distribution for t > 0 (i.e. after the wall separating region is removed), can be compared to the analytical one-dimensional solution with an initial step in concentration.The analytical solution for the diffusion equation starting at t = 0 s is given by (Balluffi et al 2005) One can test the capturing of the initial step in concentration with the methods detailed above.This gives the information on the characteristics of the interface determination by the method employed.The practical initial step distribution before the interaction between the two different concentration obtained at t = 0 s is shown in figure 2(b) (which also depicts the Voronoi cells constructed around the points).Importantly, a reflective boundary condition is applied for all the sides and the corners for Voronoi cell construction, resulting in eight repetitions of the window of interest.This ensures that the cells within the region of interest have appropriate sizes, as the end points lead to open Voronoi cells.
The results for the instance t = 0 s for the distribution shown in figure 2(b) following the sliding average method and Voronoi method are presented in figures 3 and 4, respectively.In the case of the sliding average method with length scale l = 20 px, one can observe that the higher concentration is captured in detail with the local fluctuations and also identifies the step change in concentration.However, at low concentration, one can see a marker-shot like noise with values falling near zero and peaking intermittently.This phenomenon arises from the insufficient presence of particles within the imposed length scale.Therefore, although the selected length scale is effective in capturing the step change and high concentrations, it is not that suitable at low concentrations.At larger length scale of sliding average of l = 80 px, the noise at lower concentration is present but considerably reduced.Further, the jump is broadened towards both the higher and lower concentrations.
In case of the Voronoi method, the effect of different numbers of nearest neighbors (K) averaged is shown in figure 4. With lower number of cells averaged, i.e.K = 30, concentration profile reveals all the features at small scales leading to the fluctuations at high concentrations.It is important to highlight that the marker-shot noise effect is eliminated with accounting for only K = 30 cells (cf figure 4(b)).Further, the jump in concentration is also captured reliably.When the number of averaged cells is increased, there is marginal improvement at low concentrations and some smoothening at high concentration (cf figure 4).However, the jump in concentration is notably diffused to larger length scales.It is interesting that the smoothing effect on the jump is more pronounced at lower  concentrations.In contrast, at higher concentrations, the jump remains relatively unsmoothed because the data from 100 or 500 points still falls within a relatively small physical scale.

Choosing the number of nearest neighbors K
A possible solution which can be easily implemented for the Voronoi method is to incorporate a criterion to choose the number of cells to be averaged depending on the local fluctuations.
Fluctuations due to Brownian motion: To choose a natural number of cells to average, one needs to understand the natural variation present within the uniform concentration regime.In this regard, let us consider the case of the practically uniform concentration in a closed space with Brownian diffusion.A large domain (x × y = 30 000 px × 30 000 px) with the overall non-dimensional concentration (C/C 0 ) of 1 is considered (C 0 = 0.01 ppp).Specifically, 9 × 10 6 randomly distributed particles are considered within this domain.
Firstly, the Voronoi method is applied to the entire field to cover the entire space of interest.1000 random points away from the sides are chosen and the effect of averaging K number of nearest cells is observed.The resulting average concentration and spread of the data with the 1000 points chosen is shown in figures 5(a) and (b).
The interesting part is that this distribution can be combined together to a single distribution by applying the central limit theorem (Billingsley 1995).The obtained distribution is shown in figure 5(c).The obtained histogram follows the expected Gaussian profile as observed in figure 5(c).This is useful to understand the expected deviation in the concentration with K cells considered within the uniform concentration regions.
Consequently, this insight can be extended to nonuniformly distributed points, enabling us to anticipate the uncertainty associated with averaging K cells.This knowledge becomes a crucial factor in determining the appropriate value for K. Specifically, if the deviation within K cells surpasses a predefined threshold, reflecting a certain level of stringency, there may be limited benefit in further increasing the number of averaged cells.In fact, doing so might inadvertently obscure underlying features or concentration gradients by averaging the cells covering the entire variation.Therefore, opting for a smaller value of K could be more pertinent, as it has the potential to unveil these intricate details more effectively.
Cut-off criterion: By referring to the information given in figure 5, one can observe the expected deviation with averaging K number of cells.Here, we consider moving-mean of chosen K C number of cells.For instance, if we take K C = 30, the deviation (δ) is ∼0.2 with respect to the expected value (cf figure 5(b)).
For the real case, however, the expected value is not known a priori and one needs to rely on the local mean.The available estimate for C 0 is the locally computed C Kc .The variation of C Kc itself has to be accounted for, which has highest probability to be within 0.8-1.2times the expected value (cf figure 5(b)).The computed value C Kc has highest probability to be within 0.8-1.2times the expected value (cf figure 5(b)).Hence, the total relative deviation, relative to the mean at first K C , expected could be considered as a criticality δ crit ∼ 0.5 (2 times the deviation with respect to the expected, divided by the lower of the local value).
Following this reasoning, prior to averaging the nearest neighbor cells, one can monitor the variation of the moving mean.If the moving mean lies above/below the set δ crit of the first K C cell average, then there is significant change in the concentration when farther cells are averaged.As a result, one can truncate the cells to be averaged.In other words, if the condition is fulfilled, then truncate the averaging at K = i number of nearest neighbors corresponding to the point considered.
Here, C mm (i) represents the moving mean computed over K C number of cells where i runs from the nearest to farthest cells from the cell in consideration.When K C is odd, the window is centered around the element in the current cell.When K C is even, the window is centered around the current and previous cells.
The above consideration is effectively truncating the number of averaged cells depending on the deviation when K is increased.The criticality can be generalized by following the scaling shown in figure 5(c) for any K C : The value 2.74 is the product 30 1/2 × 0.5 as chosen above and the denominator generalizes for any K C from the central limit theorem as expected from the figure 5(c).
The results for the synthetic diffusion case obtained by applying the cut off criterion (equation ( 4)) with K C = 30 and K C = 20 is shown in figures 6(a) and (b) respectively.In addition, if the length scale of truncated K cells is smaller than the minimum length scale (λ) chosen (selected spatial resolution) then continue with averaging following the same strategy before while ignoring the condition for truncation.This essentially acts like a filter for the fluctuations at very small scales as seen in figure 6.The minimum length does not change the results at the low concentration (cf insets in figure 6) as expected due to the farther spacing of the particles.These results indicate that the cut-off criterion effectively captures the jump in the concentration within small physical scale.
To compare the different methods at different time instances of the diffusion process, two time instances t = 100 s and t = 500 s are shown in figure 7. It can be seen that the concentration profiles obtained by both methods follow similar trend and fluctuate around the analytical solution.The sliding average method with l = 80 px is smoothened in comparison with the Voronoi method with cut-off criterion.

Theoretical consideration of particle image overlap
In the synthetic test case examined previously, it is supposed that all the particles in the visualized volume are detectable in the projected image as visualized by the cameras.However, in real-world scenarios, especially with an increase in concentration, particle image overlaps will occur.This means that as the physical concentration increases, the number of particle images will not reflect the actual number of the particles in the visualized image.The theoretical consideration of this effect is undertaken here.
Consider the total image area of A with N p number of particles in the visualized thin volume.Assume a critical circular disc of area A crit around the particle image in which another particle image center falls for overlap to occur in the image.The total probable number of overlapping particle images, N o , obtained from Poisson distribution of particles can be expressed as (Maas 1992, Cierpka et al 2013): To obtain the effective information/detectable particle images visible with overlap, one needs to account for the probability of exactly m number of particles at the area considered A crit (Maas 1992): The relative overlap probability (P o ) with j number of overlaps (i.e.j = 1 means there is exactly one other particle within the same A crit ) is calculated by where P (0) and P (1) represent the probabilities of having no particle and just one particle within the area, respectively, both of which do not result in overlap.P o (j ) indicates the probability of j overlaps occurring.The effective count of particles, or the number of detectable particle images within the image, is obtained by where ξ is the total number of overlaps to be taken into account for all practically possible overlaps that occur.Here ξ = 30 is considered (where the probability of overlap is practically negligible with P o < 10 −5 ).
One can define the density/concentration in the image plane when all the particles in the visualized volume are considered as ϕ a = N p /A, and the detectable number density of particle images is given by ϕ d = N d /A.At very low concentrations, when the particle images are far from each other, the overlap probability is practically zero and one can expect ϕ d ≈ ϕ a .Now, the actual density of the particles in reality (C a ) in units of particles/cm 3 is the density of the particles in the measurement volume.If one plots the graph of C a vs. ϕ d relationship from experimental data, the initial slope serves as the multiplication factor for ϕ a , yielding theoretical concentration in physical space at very low physical concentrations.Further, at higher concentrations, depending on the value of A crit the detected density of particle images ϕ d deviates from the initial slope depending on the overlaps (equation ( 9)).The resulting theoretical curves and comparisons with experimental results are presented in the next section.

Calibration
In real experiments, one can visualize the particles within the laser-illuminated sheet of certain thickness as done in PIV/PTV experiments.The physical concentration information from the visualized windows requires either the determination of the exact volume of the visualized window or calibration with a particle counter at constant concentration.Here, we take the approach with calibration with particle counter at uniform concentration in the room.More details are given in this section.
The measurement equipment is a stereoscopic PIV system consisting of two cameras, a laser, and a computer.The cameras are set at forward scattering mode w.r.t. the laser sheet (Prasad 2000).This results in substantial signal to noise ratio for accurate particle image detection.Two frames are taken which allows the calculation of the three velocity components in the measurement plane with stereoscopic PIV.Here, for concentration determination only information from one of the camera and one frame is utilized.However, in principle the data from the other frame and other camera could also be considered in order to increase the accuracy of the concentration measurement.
For concentration calibration, the closed room was nebulized with DEHS particles (mean diameter ≈ 0.4 µm) and simultaneous measurements using a particle counter AQ Guard (Palas GmbH, Germany) were performed as reference.The obtained raw PIV images at different particle concentration are shown in figure 8(a).It is essential to underscore the importance of sharp particle images for accurate particle counting and is a critical consideration during experimental setup.The average particle image diameter varied between 1.5 and 2 px for the visualized images.
To count the particle images within the PIV snapshots, it is assumed that the particle images have 2D Gaussian intensity profiles.The local maxima (intensity peaks) are counted as a particle.Before applying this, a sliding local minimum is subtracted to remove background noise and an intensity threshold for the particle cutoff is selected.DaVis (LaVision 2021) is used for performing these steps.This method does lead to undercounting at high particle concentrations due to particleparticle overlap as explained earlier.However, the subsequent calibration, elaborated in this section, addresses this concern, provided that the image particle counts do not reach saturation.
The calibration of the particle image data to the physical concentration is needed to account for the spatial variation of the laser intensity, the non-constant light sheet thickness and shape, and the varying magnification factor.The distribution of the detected particle image concentration (ϕ d ) within the visualized image is shown in figure 8(b).Further, the resultant calibration curve at certain location in space obtained with Voronoi method is shown in figure 8(c).This revealed nonlinearity in dependence of the observed particles per pixel to the physical concentration.A smoothing spline fit is obtained for all points in the x-y-plane as a function of concentration as shown in figure 8(c).
It can be observed that the particle image density ϕ d increases linearly with physical concentration of particles at low concentrations and deviates around 500 cm −3 .However, with further increase in physical concentration, the particle image based density increase is more gradual.This indicates the effect of the particle image overlaps in the observed images.Consequently, the data obtained from visualization is directly proportional to the results from the particle counter, but only at low concentrations when particle overlap is negligible.As the physical concentration increases, the likelihood of particle overlap also rises.Next, we delve into a comparison of these results with theoretical description of the overlap effect described before.
The results obtained from the particle counter and that measured from the visualized window processed with different intensity threshold are depicted in figure 9.For comparison, the results obtained from the theory presented in section 3.3 are overlaid.The theoretical curves generated are at different A crit corresponding to 2.0-3.2px (overlapping) disc diameter.The disc diameter is chosen in order to encompass the experimental results observed.It needs to be noted that the particle diameter in the images varied between 1.5 and 2 px for the visualized images.It is reasonable to assume that the overlapping disc diameter is of the order of the diameter of the particle images.The theoretical curves closely follow the experimental curve indicating the effect of the particle overlap.The theoretical curve corresponding to 3.2 px disc diameter aligns more closely with the experimental curve at a lower intensity threshold of 400 (cf figure 9(a)).Conversely, the curve corresponding to 2.4/2.8px disc diameter aligns more closely with the experimental curve at a higher intensity threshold of 1200 (cf figure 9(c)).The typical intensity of the particle image is expected to follow a Gaussian profile, so reducing the intensity threshold would lead to a larger effective particle image diameter, aligning with the observed trend.Further, the initial slope which represents the curve without overlap effect decreases as the intensity threshold increases.This indicates that the total particles considered decreases as the intensity threshold increases.Consequently, this reduction leads to reduced difference between the overlap curves and the curve without overlap at higher concentrations.At lower intensity thresholds (cf figure 9(a)), particle overlap results in a near saturation effect at higher concentrations.Hence, the overlap effect could be reduced by choosing higher intensity threshold.
It is important to note that the overlap effect is similar, regardless of whether one considers the sliding average method or the Voronoi method.

Measurement uncertainty analysis
The uncertainty of the measurements was assessed by applying the calibration function to one instance of each uniform concentration case.The calculated uncertainty ranges between 11% and 15% of particle counter measurements, considering maximum averaging of 100 cells for the Voronoi method.This variation aligns with the predicted ∼10% expected variation within the uniform case with Brownian diffusion (figure 5).The additional uncertainty is due to the fitting function which also includes the particle-particle overlap effect.For the sliding average method, the calibration function was determined similarly with length scale l = 16 pixel.The uncertainty for this method is approximately ∼11% at a concentration of 1283.1 cm −3 .While the uncertainty decreases slightly at higher concentrations, it significantly increases as the concentration decreases, attributed to marker-shot-like noise, as previously explained.Further, it needs to be emphasized that this approximation is for uniform concentration, and the uncertainty would increase with local gradients.

Experimental results with a pulsatile jet
A pulsatile jet of air nebulized with DEHS particles (mean diameter ≈ 0.4 µm) into a closed room is visualized as an example case with concentration variation.The air is heated through a water bath maintained at 37 • C and led to the source of circular cross section of diameter 3 cm.The air is driven by displacement method which is controlled to provide a defined volume of air at a defined frequency.The measurement equipment and the setup is the same as that utilized for the uniform cases such that the concentration calibration could be utilized in determining the concentration distribution with the case presented here.
The PIV images are captured at specific phases of the cycle to facilitate phase averaging.A total of 30 instances are acquired within the cycle time period.For illustration, instances at time t ′ and t ′ + 0.133 (non-dimensionalized by the cycle time) are considered, where the latter instance is the moment when the turbulent puff is still within the visualized window.
The obtained velocity field V (absolute velocity magnitude) for the case with a total volume of 1000 ml and 16 cycles/minute is shown in figure 10.
The resulting concentration profile for the instantaneous PIV images for two time instances is shown in figure 11.The concentration calibration shown in figure 8 is applied for the Voronoi method.For the sliding average method, similar calibration curve is obtained and utilized here to find the concentration contour.For a closer look, the comparison between the methods along the y = 0 line is plotted in figure 12.It can be observed from figures 11 and 12 that qualitatively both the methods reveal similar structures on the larger scales.However, the Voronoi method reveals finer details and sharper variation in the concentration contour, particularly at high concentration.At low concentrations, the Voronoi method does not reveal the oscillations associated with the marker-shot noise.

Conclusion and summary
In this study, we focus on the determination of particle concentration using images acquired from PIV.We investigate and analyze two distinct methods for concentration determination: the sliding average method and the Voronoi cell-based method.Synthetic data of a diffusion problem, featuring an initial step in concentration, is utilized to explore the characteristics and limitations of each approach.In addition to the aforementioned methods, we introduce a modified Voronoi method with a criterion to select a certain number of cells to average.This method effectively captures large range of concentration regions and features including initial jump conditions in a diffusion problem.This new method offers flexibility by allowing adjustment of the minimum physical length scale and variation with respect to natural variations due to Brownian motion.
Furthermore, calibration of the visualized images with respect to particle counter revealed particle-particle overlap at higher concentration which is described by theoretical consideration.The described methods were applied to a real experimental case of turbulent puff.It is shown that the Voronoi method is able to capture small structures with high concentrations, while providing reliable results also in regions with low concentrations.
The advantages of the sliding average method are: (a) the relatively simple computations, and (b) the integral value of the concentration is the sum of the particle images.The benefits with Voronoi method are: (a) the adaptive length scale, (b) the control over feature length scale, while still maintaining minimum number of points for reliable local concentration.The computational effort for the Voronoi method could mean that it is more logical to apply when one is interested in understanding the features in instantaneous images.When one is interested in overall average field, sliding average method does give the necessary information when averaged over many snapshots.

Figure 1 .
Figure 1.The contour of the variation of the local C I value with one particle at (0, 0) for different length scales for the sliding average: (a) l = 4 pixels; (b) l = 8 pixels; (c) l = 16 pixels.The axis dimensions are in pixels.

Figure 2 .
Figure 2. (a) The Voronoi cells around uniformly spaced particles with concentrations C 1 = 10 −4 ppp and C 2 = 0.01 ppp to the left and right of origin, respectively.(b) Voronoi cells around practically uniform distribution with Brownian diffusion with the same concentrations C 1 = 10 −4 ppp and C 2 = 0.01 ppp to the left and right of origin, respectively.

Figure 3 .
Figure 3. (a) Particle density distribution at t = 0 along the line y = 500 px obtained by the sliding average method, comparing the effect of the length scale l, and the analytical solution.(b) The zoomed-in view at the lower concentration region.

Figure 4 .
Figure 4. (a) Particle density distribution at t = 0 along the line y = 500 px obtained by the Voronoi method, comparing the effect of K and the analytical solution.(b) The zoomed-in view at the lower concentration region.

Figure 5 .
Figure 5.The concentration average and the standard deviation distribution with number of Voronoi cells averaged (K) from 10 to 1000.(b) Zoomed in view for K = 10 to K = 50.(c) The collapsed histogram of the normalized counts with scaling K 1/2 (C − Cexp)/C 0 with the overlaid Gaussian profile.The Cexp (=1) is the expected concentration.

Figure 6 .
Figure 6.Particle density distribution at t = 0 s along the line y = 500 px obtained by the Voronoi method, comparing the effect of the cut-off criterion and minimum length scale λ = 60 px for (a) K C = 30 and (b) K C = 20.The blue line is the analytical solution and the inset shows the zoomed-in view at the transition to the lower concentration region.

Figure 7 .
Figure 7. Particle density distribution obtained at time instances: (a) t = 100 s and (b) t = 500 s.The comparison is obtained along the line y = 500 px obtained by the Voronoi method, with cut-off criterion (K C = 30) and minimum length scale λ = 6 mm and the sliding average method with l = 80 px.The blue line is the analytical solution (equation (3)).

Figure 8 .
Figure 8.(a) The raw particle images at different physical concentrations, where I is the intensity in counts.The axis dimensions are in pixels.(b) Contour of the measured variation of ϕ d (ppp) at a constant physical particle concentration of Ca = 1283.1 cm −3 .The Voronoi method with intensity threshold of 800 counts is applied.(c) The experimental points of ϕ d (ppp) at different physical particle concentrations Ca and the calibration function obtained at point location marked by × in (b) as a smoothing spline fit.The error bars indicate the standard deviation of the particle counter measurements.

Figure 9 .
Figure 9.Comparison of the experimentally measured particle number density Ca from particle counter versus detected in the image (ϕ d ) and the theoretical values at different intensity threshold: (a) experimental intensity threshold of 400 counts; (b) experimental intensity threshold of 800 counts; (c) experimental intensity threshold of 1200 counts.The theoretical critical area (A crit ) corresponds to diameter of 2.0 (cyan), 2.4 (blue), 2.8 (green) and 3.2 (red) px.Further, the dashed line depicts the initial slope of the experimental data points which represents the theoretical curve if there were no overlaps.

Figure 11 .
Figure 11.The concentration contour obtained by (a) and (b) the sliding average method with l = 16 pixel and intensity threshold of 800 counts; (c) and (d) Voronoi cell method with cut off criterion (K C = 30) and minimum length scale λ = 6 mm.The results are for same instances as in velocity field in figure 10, i.e. instances t ′ (a) and (c) and t ′ + 0.133 (b) and (d).

Figure 12 .
Figure 12.The concentration profile comparison with sliding average method and Voronoi method (as presented in figure 11) along the y = 0 line at time instances: (a) t ′ and (b) t ′ + 0.133.