Preferential concentration of heavy particles : A Voronoï analysis

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Preferential concentration of heavy particles: A Voronoï analysis Romain Monchaux, Mickaël Bourgoin, Alain Cartellier


I. INTRODUCTION
The study of inertial particles laden flows is relevant to many industrial and environmental issues ͑chemical reactors and engine optimization, plankton or pollutant dispersion, clouds formation, etc.͒, but it is also of fundamental interest.4][5][6] Another interesting feature is the enhancement of the settling velocity of particles in turbulent flows.Since an explanation of seeding density dependence of this phenomenon through collective effects has been proposed, 7 different authors have tried to quantify and characterize particles clustering from numerical simulations.Nevertheless an appropriate equation of motion for an accurate modeling of particle dynamics is still lacking, only the point particle limit is generally considered 8,9 and even further simplifications are usually required to achieve numerical simulations. 10,11Experimental investigations are still important to reach a better understanding of the underlying mechanisms.3][14][15] We consider here the many particles case to address preferential concentration and clustering matters including possible collective effects.Do clusters exist as a whole in these flows?How do they form?What is their structure and how does this structure evolve with time?Which effect do they have on the single particle dynamics?Here are some of the questions that still need to be answered.To date, the preferential concentration/clustering problem has been studied with global or Eulerian tools ͑such as box counting methods, pair correlation function estimation, correlation dimension, or topological indicators͒.A dynamical study of the Lagrangian dynamics of particles and of the local concentration field would bring a new insight in these processes.In particular, it would be worthwhile to get access to the concentration along a particle trajectory, a quantity that has been recognized as very important in models. 16,17n that scope, we propose in this article a new approach to analyze particle concentration fields based on Voronoï tesselations that gives a measure of the local concentration field at interparticle length scale.By itself, this data processing technique is particularly enlightening to explore and quantify the preferential concentration phenomenon while, combined to Lagrangian tracking, it also gives access to simultaneous measurements of velocity, acceleration, and local concentration along particles trajectories that are crucial for a better insight of clusters and particle dynamics.As a first step, this article focuses on the preferential concentration problem on a statistical ground.Section II is dedicated to the description of our experimental setup and of the tools used to postprocess the acquired data.In Sec.III, we use Voronoï analysis to quantify the preferential concentration and to identify and characterize clusters.Conclusions and expected further dynamical studies are summarized in Sec.IV.

A. Two-phase flow generation and characterization
Experiments are conducted in a large wind tunnel with a square cross-section of 0.75 m ϫ 0.75 m where an almost ideal isotropic turbulence is generated behind a grid whose mesh size is 7.5 cm ͑Fig.1͒.We can adjust the mean velocity from 3 to 15 m s −1 ͑the turbulence level remains relatively low of the order of 3% at the measurement location and the anisotropy level between the transverse and longitudinal fluctuating velocities is smaller than 10%͒.Inertial particles are water droplets generated by four injectors placed in the convergent part of the wind tunnel, one meter upstream the grid to ensure a homogeneous seeding of the flow.These injectors consist of two tubes carrying water and air that merge into a specifically designed nozzle ͓a Schlick-Dusen ͑Germany͒ 942 two fluid nozzle͔ where pressurized air atomizes the liquid and forms the exiting conic jet.The particle properties relevant to our study can be adjusted within accessible range.Volume loading is directly linked to the water flow rate in the nozzles while the diameter distribution is a subtle function of both water flow rate and air pressure in the nozzles, as detailed below.We always consider regimes of relatively low particles volume loading ͑volume fraction in our experiments covers the range 2 ϫ 10 −6 Ͻ v Ͻ 3 ϫ 10 −5 ͒ so that no significant turbulence modulation by two-way coupling is expected to occur.
We have performed particle image velocimetry ͑PIV͒ and hot-wire measurements in the test section in order to characterize the turbulence underlying our experiment.Hotwire acquisitions have been done for the single phase air flow.They show a slowly decaying nearly homogeneous turbulence.When the mean flow velocity is changed from 3 to 7.5 m s −1 , the associated Taylor microscale Reynolds number R varies from 72 to 114 at the measurement location situated 3 m downstream the grid.The evolution of the turbulence characteristics of the single phase flow at different locations downstream the grid for the explored range of mean velocities is given in Table I.In order to probe the maximum influence on the carrier turbulence of the additional perturbation from the injectors, we have also measured the turbulence fluctuation level when the injectors blow only air ͑hot-wire measurements are not possible in the two-phase configuration͒.It is found to increase typically by 50% com-pared to the case without injectors.Nevertheless, this enhancement of turbulence intensity does not depend much on the imposed air pressure at the injectors as shown on Fig. 2͑a͒ ͑except at very low wind speed͒.It should also be noted that this influence is expected to be reduced when the water flow is turned on since atomized droplets will carry a significant amount of the global additional momentum injected by the nozzle's air flow.PIV was carried out from high speed imaging ͓we use a Phantom V12 camera from Vision Research, Inc. ͑New Jersey, USA͔͒ of the injected water droplet fog.We note that although such inertial particles are not expected to follow the flow, these PIV measurements give a good characterization of the global flow homogeneity and isotropy.Figure 2͑b͒ presents PIV vertical profiles of the axial velocities averaged over x along the whole measurement volume showing the good homogeneity of the flow in the cross-section of the tunnel; moreover, the very slow evolution of these profiles with the axial coordinate x along the measurement volume ͓Fig.2͑c͔͒ evidences the slow decay of turbulence within this volume.On Fig. 2͑d͒, a typical velocity spectrum allows to identify a narrow inertial range as expected at this moderate R value.
Regarding the droplet spray, one of the main goals of our study is to explore the influence of particle diameters d ͑or equivalently of their Stokes number, defined as Stϵ p / = ͑d / ͒ 2 ͑1+2⌫͒ / 36, where p is the particle viscous relaxation time, and are, respectively, the dissipation time scale and length scale of the carrier flow, and ⌫ = water / air is the ratio of particles density to carrier fluid density͒ on the preferential concentration phenomenon.The two-phase nozzles that we use to produce the spray present the advantage to easily allow a systematic variation of the mean diameter of droplets by controlling air and water flow rates in the injectors.The drawback of this versatility is the difficulty to produce a monodisperse spray.The injection process leads indeed to a polydisperse seeding of the flow for which particles diameter distribution has to be estimated.We have used the SprayTech instrument developed by Malvern, Inc. ͑England͒ and based on the diffusion of a laser beam by the spray.Due to the size of the SprayTech apparatus, these granulometry characterizations cannot be done simultaneously with other measurements; therefore, they have been performed once for the entire set of experimental parameters and in the same measurement volume as the main measurements.Figure 2͑e͒ shows typical particle size distributions.The main control parameters governing the diameter distributions are: the air pressure and the water flow rate in the injectors as well as the wind velocity in the tunnel.Diminishing the air pressure while keeping the other parameters constant results in an increase of particle diameter.Oppositely, diminishing the water flow rate produces a decrease of particles diameter, a decrease of the volume loading but an increase of the number of particles seen per image for a monodisperse spray.This nontrivial effect could be interpreted as the result that for a given water flow rate j w , particle number density C 0 is inversely proportional to droplets volume ͑C 0 ϰ d −3 ͒, whereas the particle diameter d is a nontrivial injector dependent growing function of j w ͓let us say d = g͑j w ͔͒, which for our precise injectors must increase more  x ͑cm͒ slowly than j w 3 so that in the end seeding density C 0 ϰ g͑j w ͒ −3 is a decreasing function of water flow rate j w .Since particles are injected in the convergent part of the tunnel, upstream the grid for low wind velocities ͑typically V 0 Յ 4 m s −1 ͒, the velocity in the convergent is not sufficient for the heavier particles to reach the grid and they all settle before entering the fast section of the tunnel.This results in a severe decrease of particles diameters observed in experiments at fixed air pressure and water flow rate for low wind velocities.We have shown that above V 0 Ӎ 4 m s −1 , particles diameter distributions in the measurement region are not much affected by the wind velocity.Once particles diameter distributions are obtained for a given set of injection parameters, we build a corresponding most probable Stokes number defined as St= ͑d max / ͒ 2 ͑1+2⌫͒ / 36, where d max corresponds to the maximum of diameters distribution ͓see Fig. 2͑e͔͒.Figure 2͑f͒ represents the evolution of d max as a function of air pressure and water flow rate in the nozzles for a given carrier flow Reynolds number.We stress that because of the high polydispersity the standard deviation St of Stokes number ͑based on measured diameters distribution͒, which could be interpreted as an error bar for the Stokes number estimation, is large ͑ St / St easily exceeds 50%͒ and will not be displayed on the figures presented in the sequel.
To summarize, each experiment is characterized by the set of parameters ͑R ,St,C 0 ͒ with R , the carrier flow Reynolds number ͑based on Taylor microscale͒, St, the average Stokes number obtained as described above, and C 0 , the global seeding density ͑which we cannot accurately measure, but which we assume to be directly related to the number of particles per image͒.Another usual relevant parameter is the Rouse number ͑defined as the ratio of the terminal velocity of the particles to the turbulent fluctuations intensity͒, which in our experiments, varies in a relatively narrow range from 0.4 to 2. A total of 90 experiments covering a set of about 40 different parameters triplets ͑R ,St,C 0 ͒ have been explored in order to investigate systematic effects on preferential concentration phenomenon.In this study, we focus mainly on the influence of St and C 0 .

B. Acquisitions and postprocessing
Acquisitions are performed using a Phantom V12 high speed camera ͑Vison Research, USA͒ operated at 10 kHz and acquiring 12 bits images at a resolution of 1280 pixelsϫ 488 pixels corresponding to a 125 mm ͑along x͒ ϫ 45 mm͑along y͒ visualization window on the axis of the wind tunnel ͑covering slightly less than an integral scale in the vertical y direction and almost two integral scales in the streamwise x direction͒, located 2.95 m downstream the grid.The camera is mounted with a 105 mm macro Nikon lens opening at f / D = 2.8.An 8 W pulsed 2. ͑Color online͒ ͑a͒ Evolution of the ratio of turbulent fluctuations ͑V rms ͒ to mean wind velocity ͑V 0 ͒ both measured from from hot-wire anemometry as a function of the injector pressure without water flow ͑note that the impact of the injectors on the turbulence level would even be lower with the presence of water͒ for three different mean velocities ͑equivalent turbulent fluctuation rate is constant and of the order of 3.5% when no air is blown in the injectors͒.͑b͒ Mean axial velocity averaged over x in the measurement volume ͑x ͓289-301͔ cm͒ obtained from PIV. ͑c͒ Mean axial velocities at different locations x in the measurement volume ͑also from PIV͒.Different symbols correspond to different location, note that the mean velocity is almost constant within the measurement volume.͑d͒ Velocity spectrum in the measurement volume.͑e͒ Particles diameter PDF evolution with air pressure ͑varying from 2 to 5 bars͒ at fixed water flow rate ͑1.2 l/mn for each injector͒ and fixed wind velocity ͑V 0 = 4.5 m s −1 ͒.Note that PDFs are given from particles volume ͑and not particles number͒.͑f͒ Evolution of the maximum of these PDF with water flow rate and air pressure at the same fixed velocity ͑V 0 = 4.5 m s −1 ͒.

103304-3
Preferential concentration of heavy particles Phys.Fluids 22, 103304 ͑2010͒ copper laser synchronized with the camera is used to generate a 2 mm ͑i.e., 3 -4͒ thick light sheet illuminating the field of view in the stream-wise direction.The camera is orientated with a 50°forward scattering observation angle with respect to the laser sheet to increase the light budget.The resulting deformations are compensated by a Scheimpflug mount.Each experiment consists in a 0.9 s acquisition of 9000 images ͑corresponding to the available on board memory storage of the camera͒ at fixed wind velocity, water rate and air pressure in the injectors.Particles are identified on the recorded images as local maxima with intensity higher than a prescribed threshold.As a consistency test, we have checked that changing the threshold around the selected value does not impact significantly the number of particles found.The subpixel accuracy detection is obtained by locating the particles at the center of mass of the pixels surrounding the local maxima.
As already mentioned, our aim is to study particles concentration fields in order to quantify preferential concentration effects and to identify particles clusters if any.Usual approaches to do this consider the pair correlation function to quantify preferential concentration effects while a box counting method is preferred to access local concentration fields.We propose to use a single tool to tackle simultaneously these two problems: the Voronoï diagrams.Such a Voronoï diagram is the unique decomposition of two-dimensional ͑2D͒ space into independent cells associated to each particle.One Voronoï cell is defined as the ensemble of points that are closer to a particle than to any other.Use of Voronoï diagrams is very classical to study granular systems and has also been used to identify galaxies clusters.The Voronoï diagram computation is very efficient ͑we use MATLAB algorithm͒ with the typical number of particles per image ͑up to several thousands͒ we have to process.Figures 3͑a͒ and 3͑b͒ show a raw acquired image as well as the located particles and the associated Voronoï diagram.

A. Voronoï diagrams: Properties and advantages
Why use Voronoï tessellations?From the definition of the Voronoï diagrams, it appears that the area A of a Voronoï cell is the inverse of the local 2D-concentration of particles; therefore, the investigation of Voronoï area field is strictly equivalent to that of local concentration field.We recall that usually local concentration fields are obtained through box counting methods 7 that shows several disadvantages: they are computationally inefficient and they require to select an arbitrary length scale ͑the box size͒, whereas in Voronoï diagram computation, no length scale is a priori chosen and the resulting local concentration field is obtained at an intrinsic resolution.Similarly, the pair correlation function only gives global ͑nonlocal͒ information and is also associated to the choice of a length scale that spans the whole values of interest increasing dreadfully the computation time.Finally, another interest of Voronoï diagrams is that as each individual cell is associated to a given particle at each time step, thus tracking in a Lagrangian frame the particles directly gives access to the Lagrangian dynamics of the concentration field itself along particles trajectories.Although we will not discuss such Lagrangian aspects in this article, they represent an important opening which will be addressed in forthcoming studies.
Some relevant properties of Voronoï diagrams.Whatever the measurement and data analysis technique used, when one refers to preferential concentration, it is implicitly assumed that one deals with statistical preferential concentration compared to the case where particles would be spatially distributed as a random Poisson process ͑RPP͒.In order to quantify preferential concentration, we therefore compare for each experiment the probability density function ͑PDF͒ of the measured Voronoï areas to that expected for a RPP.The main known properties of Voronoï diagrams associated to RPP can be found in a short review by Ferenc and co-worker 18 and references herein.The first moment of Voronoï area PDF has nothing to do with the spatial organization of the particles since the average Voronoï area A ¯is always identical to the mean particles concentration inverse.Therefore, in the sequel, we will generally focus on the distribution of the normalized Voronoï area V = A / A ¯that is of unit mean.The only known exact result for RPP Voronoï areas statistics concerns the second order moment in the 2D case that is equal to ͗V 2 ͘ RPP = 1.280 corresponding to a standard deviation V RPP = ͱ͗V 2 ͘ RPP −1Ӎ 0.53.Regarding the shape of the PDF of Voronoï areas statistics for a RPP, no analytical solution is known ͑most of the authors fit them with Gamma distribu-tions͒.Ferenc and co-worker proposed a compact analytical expression involving the space dimension as a single parameter.We use this analytical expression as a RPP reference.

B. Experimental Voronoï area PDF
From instantaneous to statistical results.After postprocessing the data as explained in Sec.II B, we obtain for each experiment ͓corresponding to a given set of control parameters ͑R ,St,C 0 ͔͒ a Voronoï diagram at each acquisition time step.Due to the finite size of the field of view, the cells at the border of the diagram have infinite areas.We reject these border cells and keep only the particles whose Voronoï cell is fully included in the visualization window.To improve the statistical convergence of Voronoï area PDF estimation, we compute Voronoï diagram statistics from several uncorrelated images from a given experiment.The unit mean nor- malization ͑V = A / A ¯͒ may then be performed for each image individually or for the whole set of images.We have checked that the impact of the normalization process is negligible.Figure 6͑a͒ shows Voronoï area PDFs associated to ten different statistical sets obtained by gathering 500 uncorrelated instantaneous Voronoï diagrams ͑UIVD͒ from a same experiment.The dispersion between the ten samples is quite small indicating that a good statistical convergence is already reached with 500 UIVD, except maybe in the extreme tails of the PDF.We find that the statistical convergence of the two first moments of these distributions is indeed achieved with one set of 500 UIVD.In the sequel, detailed PDF analysis was usually calculated using one set of 5000 UIVD from each experiment, while statistical error bars on second order moments for a given experiment ͓for instance, for the standard deviation of Voronoï areas V in Fig. 5͑a͔͒ are estimated from the dispersion when these 5000 UIVD are split into ten sets of 500 UIVD analyzed individually.
Experimental Voronoï area distributions.Figure 4͑a͒ displays the PDFs of the dimensional Voronoï cells area A for 40 different experiments.When the dimensional area A is considered, one observes that the maximum of these PDFs spans over 2 decades.This evolution is representative of the average number of particles per image ͑or equivalently of the global seeding concentration C 0 ͒ that for the ensemble of experiments represented goes from 50 to 5000.Note that as the average number of particles per image decreases ͑i.e., as the mean Voronoï area increases͒, the scattering of the right tail on these PDFs increases as a consequence of the lesser statistical convergence.The same PDFs for the normalized Voronoï area V = A / A ¯collapse somehow ͑not shown͒.The tails associated to large Voronoï areas, i.e., to regions of low particles concentration, actually collapse within measurements uncertainty, whereas the tails associated to relatively small Voronoï areas, i.e., to the regions of high particle concentration show a much stronger dependency with the considered experimental conditions, and particularly with the Stokes number as shown in Fig. 4͑b͒.Concerning the Reynolds number dependency, the robustness of PDF right tails ͑associated to depleted regions͒ may be interpreted as the fact that the integral scale of the carrier flow is found not to change much when R is increased ͑on the contrary the dissipation scale does decrease significantly͒.If large depleted regions are mostly associated to large eddies whose size is not affected when R is increased, the distribution of Voronoï areas in these large depleted regions may also be expected to remain robust as Reynolds number changes.It is more difficult to comment on the R dependency of the PDF left tails ͑associated to dense regions͒ for two reasons.The range of Reynolds numbers tested is narrow and St and R are very correlated ͑see I͒ so that sets of experiments who span the explored range of Reynolds numbers at fixed St and C 0 are too scarce to be conclusive.Nevertheless, we observe changes in the PDF left tails with the Reynolds number that deserves a more thorough study we may undertake.The dependency on Stokes number is more subtle and is analyzed in the next paragraph.While Voronoï area PDF of RPP are usually described by Gamma distributions, 18 Fig.4͑c͒ shows that, for the inertial particles we investigate, Voronoï areas statistics are well described by a log-normal distribution.As seen on the figure, the superimposition with a log-normal distribution is almost perfect in the interval Ϯ2 log͑V͒ , where log͑V͒ stands for the standard deviation of log͑V͒ ͑note that extreme PDF tails for statistics of the logarithm of V are beyond experimentally accessible statistical convergence͒.Log normality is further tested in Fig. 4͑d͒ via the well-known log-normal relations between the two first statistical moments of V and log͑V͒ where V stands for the standard deviation of the normalized area V, and where by definition mean normalized area V ¯=1.
In this figure, each point for each plotted relation corresponds to each of the 90 carried experiments evidencing the very good verification of log-normal scaling ͑1͒ and the relatively good verification of relation ͑2͒ that is much more severe considering the difficulty to reach statistical convergence for logarithm of second moments.At the first order, we can therefore reasonably assume normalized Voronoï areas to be log-normally distributed.To date, we do not have any theoretical interpretation of this log-normality, but this result shows that in the limit of experimental convergence, normalized Voronoï area PDFs can be described with one  single scalar quantity that we choose here to be the standard deviation of the normalized Voronoï areas V .Quantifying preferential concentration.It is generally assumed that preferential concentration is primarily governed by particles Stokes number.Most former numerical and experimental studies have evidenced that preferential concentration effects are more significant as the Stokes number approaches unity.In Figs.5͑a͒ and 5͑b͒, we present the evolution of the normalized Voronoï area standard deviation V as a function of Stokes number in our experiments.In Fig. 5͑a͒, each point corresponds to one single experiment.Lines connect different experiments for which the Reynolds number and the number of particles per image are kept constant.In Fig. 5͑b͒, only the Reynolds number is constant along the lines, and each point is the average of points from the previous plot corresponding to experiments with the same Stokes number but possibly different number of particles per image ͑error bars represent the dispersion between such experiments͒.As mentioned earlier, the standard deviation of Voronoï areas for a RPP is analytically known to be V RPP Ӎ 0.53 that defines a reference value to compare with.A standard deviation V significantly exceeding 0.53 reveals the existence of high and low concentration events compared to the RPP case.Oppositely, a standard deviation V below this reference value would evidence the tendency of particles to distribute in a more organized arrangement ͑ V = 0 in the limit of a perfect crystal͒.As seen on Figs.5͑a͒ and 5͑b͒, for the range of explored Stokes numbers ͑spanning from 0.2 to almost 6͒, the standard deviation of the normalized Voronoï areas always exceeds 0.57 and reaches values as high as 0.85, which shows that preferential concentration is always present over this range of Stokes numbers, consistently with former experimental and numerical investigations.Furthermore, both figures peak around St pk Ӎ 2 -3 that is consistent with the generally assumed feature that preferential concentration is maximal for Stokes number around unity corresponding to a better adjustment of particles response time to turbulence forcing time.In order to compare our Voronoï analysis with usual tools, we have also performed the box counting and correlation dimension analysis ͑not shown͒ that are consistently found to give a maximum effect of preferential concentration also for St pk Ӎ 2 -3. Figure 5͑b͒ also suggests a possible Reynolds number effect as the curve for R = 90 seems more peaked than curve at R = 114 and appears to reach its maximum for a slightly lower value of St pk .However, one can hardly be conclusive on such a Reynolds number effect as it is mainly supported by only one point on the R = 90 curve ͓point St= 2.2; ͗ V = 0.6͘ in Fig. 5͑b͔͒ that furthermore presents relatively large statistical error bars.The possible specific influence of Reynolds number on preferential concentration will be addressed more deeply in forthcoming investigations.Moreover, as discussed in the following, the precise shape of trends in Figs.5͑a͒ and 5͑b͒ might be slightly affected by a bias on the estimation of particles typical Stokes number and subsequently on the exact value of St pk as well as the width of the curves.For example, we find St pk տ 1 while other studies ͑mostly nu-merical͒ suggest a maximum effect for StӍ 1.This might be attributed to an overestimation of the typical Stokes number of particles in our experiments as the actual turbulent dissipation scales ͑ or ͒ might be slightly smaller than the values reported in Table I due to the slight turbulence level enhancement from the spray injection process as described in Sec.II A ͑measurements in Table I, used for St estimation, were performed for the carrier flow alone, with particles injectors off, since accurate measurements of turbulence dissipation scales cannot be done at present with the sprays on͒.An increase of the actual turbulence level by 30% ͑as typically observed when the nozzles blow only air but no water is injected͒ would, for instance, reduce by a factor of 2 the Stokes number estimation bringing St pk to a value closer to 1.
Finally, we also note that since our experiments cover a wide range of particles seeding concentrations, it is necessary to avoid any possible statistical bias on V depending on the number of particles per image that are processed.This is a requirement for V to be considered as an actual quantitative indicator of preferential concentration.For instance, experiments with large particles ͓large Stokes numbers in Figs.5͑a͒ and 5͑b͔͒ generally have less particles per image ͑typically less than 1000 ppi͒ than experiments with smaller particles ͑which typically have 3000-4000 ppi͒.To test such a possible bias, we have estimated V from a set of originally highly loaded images from which we randomly removed particles.We have shown that the estimation of V is extremely robust and not biased by this subsampling procedure as it is only reduced by less than 1.5% when up to 80% of the particles are randomly removed from the images.As a partial conclusion, we emphasize that Voronoï analysis allows to robustly quantify preferential concentration with a single scalar quantity ͑the standard deviation of normalized Voronoï areas͒ that is easily accessible and efficiently computed.This analysis confirms the trend of inertial particles to preferentially concentrate with a maximal effect for particles with Stokes number of order unity.However, we note that Voronoï area PDF on their own do not contain any information concerning the spatial structure of particles concentration field that is a key point for the study of clustering and of its dynamics.In Sec.III C, we show how Voronoï area distributions can nevertheless be used to analyze the concentration field structure and dynamics.

C. Toward clusters identification and characterization from Voronoï analysis
The Voronoï area PDF may be used to identify clusters of particles as follows.Consider Fig. 6͑a͒ that presents the superposition of the Voronoï PDF for a typical experiment and for a RPP.These PDFs intersect twice ͓which is more visible on Fig. 6͑b͒ showing the ratio of both PDFs͔.For low and high values of normalized Voronoï area corresponding to high and low values of the local concentration, experimental PDF is above the RPP one while we observe the opposite for intermediate area values.This is consistent with the intuitive image of preferential concentration.Inertial particle concentration field is more intermittent than the RPP with more probable preferred regions where concentration is higher than the Poisson case and subsequently also more probable depleted regions where concentration is lower than in the Poisson case.We consider these intersection points V c and V v as an intrinsic definition of particle clusters and voids.For a given experiment, Voronoï cells whose area is smaller than the first intersection V c are considered to belong to a cluster while those whose area is larger than the second intersection V v are associated to voids.We insist on the fact that these thresholds are intrinsically chosen experiment wise and vary from one experiment to another.In particular, their evolution with the seeding concentration C 0 is affine.Figure 6͑c͒ displays a full Voronoï diagram corresponding to one image taken from one experiment.On this diagram, cells corresponding to clusters ͑resp.to voids͒ have been colored in dark gray ͑resp.light gray͒ while the remaining cells have been patched with white.It appears that dark gray cells ͑resp.light gray cells͒ tend to be connected in groups of various sizes and shapes that we identify as clusters ͑resp.voids͒ whenever they belong to the same connected component.
We then analyze the geometrical structure of the identified clusters and voids.We have computed their area and perimeter from the area and perimeter of the constitutive Voronoï cells.We present the distributions of the clusters area A c in Fig. 7͑a͒ and the scatter plot of their perimeter P c versus A c 1/2 in Fig. 7͑b͒.This last figure shows that for small and compact clusters, P c and A c 1/2 are almost linearly related while large structures exhibit a fractal behavior as also reported in previous experimental and numerical studies ͑see Ref. 7 and references herein͒.As shown in Fig. 7͑a͒, cluster areas A c are found to be algebraically distributed with an exponent Ϫ2 that is found to be independent of both R and St within measurements accuracy ͑see inset͒.This interesting observation implies that clusters do not have a typical size ͑not even an average͒.This result is contrary with the previous study by Aliseda et al. 7 who fitted clusters area distributions ͑measured from a box counting method͒ with a decaying exponential from which they could estimate a typical cluster size.However this estimation is questionable as the clusters area distributions they report are not clearly exponential and generally defined on so few points that an algebraic behavior cannot be ruled out.We have performed the same analysis with the empty regions ͑not shown͒ and we find very similar results than for clusters: large depleted regions have a fractal structure and their areas are algebraically distributed with an exponent close to Ϫ1.8.This last result is in agreement with recent direct numerical simulations 19 where void regions between inertial particles have been analyzed.
We now analyze the concentration within the clusters formerly identified.We define C, the mean concentration inside a cluster, as the inverse of the mean Voronoï area within this cluster and C 0 as the average particles concentration in the whole images during one experiment ͑C 0 is the average number of particles per image divided by the area of the whole visualization domain͒.Figure 8͑a͒ presents the PDF of C / C 0 for nine typical experiments spanning our control parameters space.Though these PDFs are poorly converged due to an evident lack of statistics ͑each identified cluster counts as one statistical sample͒ they are still well described, at first order, by a Gaussian distribution.Their means ͗C / C 0 ͘ define the averaged normalized concentrations within clusters for each experiment; it measures the average relative overloading of particles inside clusters.In our experiments, mean concentration in clusters is generally found to exceed twice the global seeding concentration and it can reach up to eight times the latter depending on experimental conditions ͓see Figs.8͑b͒ and 8͑c͔͒.Figure 8͑b͒ shows the evolution of these mean normalized concentration with Stokes number, a globally growing trend with St can be seen, though scattering of data points is relatively large.On the contrary, as evidenced in Fig. 8͑c͒, the mean normalized concentration in clusters turns to be clearly dependent on the global seeding concentration C 0 following a decreasing power law ͗C / C 0 ͘ ϰ C 0 −0.39 .While one might have expected the concentration inside the clusters to grow linearly with the global concentration ͑͗C / C 0 ͘ remaining constant͒ as a trivial geometric effect, this surprising decreasing power law dependency shows that the average concentration inside clusters grows more slowly than the global seeding concentration.As a consequence, the relative overloading of particles inside the clusters is reduced when the global seeding increases: for the lowest values of C 0 clusters are up to 8 times overloaded while they are only twice as dense as the average for the largest values of C 0 .To rule out any possible bias due, for instance, to a statistical undersampling of low concentration experiments, we have performed the same whole analysis on artificially subsampled data: starting from the highest concentration experiment ͑containing several thousands particles per image͒, we randomly removed up to 90% of the particles and we recomputed Voronoï tesselations and clusters identification.We found that the algebraic distribution of clusters area is always maintained with the same −2 exponent, that the distribution of concentration within clusters is still Gaussian, and that mean normalized concentration in the clusters with the subsampled number of particles per image is weakly affected by the subsampling ͓see dots data on Fig. 8͑c͔͒.This robustness of the analysis to arbitrary subsampling shows that the effect observed in Fig. 8͑c͒ does result from a physical process related to particles loading and not from a statistical bias.͓We have also checked that this behavior was not due to the definition of clusters we introduced based on the distance of Voronoï area PDF to RPP.Indeed, the same analysis was repeated for various fixed concentration thresholds with respect to C 0 ͑which is the way clusters are commonly defined in the literature͒, and the same law ͗C / C 0 ͘ ϰ C 0 −0.39 was recovered.͔Such a nontrivial effect, which has never been reported to our knowledge, is therefore intrinsic to the clustering of dense particles by turbulence.It probably reveals the existence of collective effects inside the clusters, and certainly deserves further investigations as well as cross-analysis with other experiments and numerical simulations.In particular, it should be underlined that this behavior is most probably not due to a steric effect.If so, one should expect both a linear increase of ͗C͘ with C 0 in the dilute limit and a saturation of ͗C͘ ͑irrespective of C 0 ͒ at very large concentrations.None of these trends appear in Fig. 8͑c͒ ͑note that the concentration covers more than one decade in this figure͒.Finally, we further check a possible dependency of ͗C / C 0 ͘ on Stokes number ͓which was not clearly visible at first sight in Fig. 8͑b͔͒ by considering data on Fig. 8͑b͒ divided by the empirical power law ͗C / C 0 ͘ ϰ C 0 −0.39 ͓obtained from Fig. 8͑c͔͒ as a first attempt to separate Stokes number effects from the global seeding concentration effects just described.The result is presented in Fig. 8͑d͒ that shows a possible dependency on Stokes number with a maximum mean concentration in clusters occurring around StӍ 2 reminiscent of the maximum of preferential concentration observed in Figs.5͑a͒ and 5͑b͒.Although further investigations will be needed to be conclusive on this point, as discussed in Sec.IV, such a dependency can be consistently interpreted in the framework of a sweep-stick mechanism. 20,21

IV. CONCLUSIONS AND OUTLOOKS
We have introduced the analysis of preferential concentration in turbulent particles laden flows using Voronoï diagrams as a new tool to get a quantitative insight of the phenomenon.This very computationally efficient tool, not only gives access to the concentration field at an intrinsic local resolution ͑inverse of Voronoï areas directly gives local particles concentration͒ but it also offers a remarkably simple and nonambiguous way to define particles clusters ͑as well as complementary voids͒ and to analyze their structural properties.Several known behaviors ͑previously reported in experimental and numerical works based on classical tools, mainly pair correlations, correlation dimension, and box counting͒, as the maximum of heterogeneity of the concentration field for particles with Stokes numbers around unity, have been successfully recovered and further analyzed at the light of this new tool.By systematically varying the triplet of parameters ͑St, R , C 0 ͒, we have shown that particles Voronoi's area distributions are always reasonably lognormal, so that preferential concentration can be quantitatively measured by a single scalar ͑the standard deviation of these distributions͒.We have characterized clusters ͑and voids͒ geometries and their inner concentration.Geometrical properties ͑mainly the algebraic distribution of mean cluster areas and the fractal structure of clusters͒ may be interpreted as an evidence of the self-similar nature of preferential concentration in particle laden flows.In particular, clusters do not appear to have any characteristic typical scale.The analysis of particle concentration inside the clusters has revealed two new and so far unpredicted results.

͑i͒
The average particle concentration inside the clusters depends on the global particle loading in a nontrivial way.
͑ii͒ After the compensation of this particle loading dependency, the average concentration inside the clusters exhibits a nonmonotonic dependency on the Stokes number with a maximum around unity values.
To summarize Stokes number effects, we find that the overall preferential concentration is enhanced for particles with Stokes number around unity, cluster sizes, and geometries do not depend on Stokes number while mean concentration within clusters is also enhanced ͑at the second order when seeding concentration effects are decoupled͒ around unity Stokes numbers.All together, these results can be consistently interpreted in the framework of recently introduced sweep-stick mechanism by Vassilicos and co-workers. 20,21hey propose a new physical process for the particle preferential concentration phenomena where clusters form around zero acceleration points of the carrier flow to which particles tend to stick and to be advected with.These special points of the carrier flow have been thoroughly investigated in DNS 20,21 and found to have at least two main remarkable properties: ͑i͒ their spatial distribution is not uniform, and they tend to clusterize and ͑ii͒ inertial particles tend to stick to these points with an optimal stickiness for particles with Stokes number around unity.Physically, the sticking mechanism is dominated by the convergence of particles toward zero acceleration points along the maximum compression direction of acceleration gradients tensor of the carrier flow.This convergence is optimal when particles response time p is of the order of the inverse of the corresponding contraction rate 21 ͑if the response time is too small, contraction is too slow, while for too large response times contraction is fast but overshot by the particle inertia͒.In the context of the stick-sweep mechanism, cluster shape and size are therefore expected to be mostly prescribed by the carrier turbulence itself ͑in terms of clustering properties of zero acceleration points themselves͒.For instance, the fractal structure of zero acceleration point distribution has been identified in DNS. 21his is consistent with the self similarity we have observed in our experiment.In the same time, as a result of the Stokes number dependency of the stickiness of zero acceleration points, the concentration of particles within clusters is also expected in the context of this model to depend on Stokes number ͑with an optimum around Stϳ 1͒ and consequently so does the overall preferential concentration.This is also consistent with Stokes number trends observed experimentally.However some observations still remain to be understood.Yet, considering the sweep-stick mechanism alone, it is hard to qualitatively identify the origin of a departure from a linear increase of ͗C͘ with C 0 .Such a quantification remains to be done from available simulations. 20,21One possibility could be the indirect role of gravity ͑gravity was not accounted for in the above mentioned simulations͒.Indeed, and in addition to the Stokes number, the Rouse number, which is the ratio of the terminal velocity to the turbulent intensity, is known to affect the behavior of particles in turbulence and in particular the preferential sweeping ͑see, for example, Ref. 22͒.In addition, particle-turbulence interactions are known to enhance the settling velocity of dense inclusions. 7,8,23One may expect that such a settling velocity would tend to decorrelate somehow the particle locations from those of zero acceleration points, making clustering less efficient than expected.An argument in favor of such a scenario is that the enhanced settling velocity of inclusions in clustered regions is a collective effect that arises from the local concentration alone and thus is not Stokes dependent. 7n this framework, the influence of turbulence intensity on settling velocity enhancement is not completely clarified.If, as identified from simulations, 8 this enhancement is proportional to the turbulence intensity, one may expect an evolution of the relationship between ͗C͘ and ͗C 0 ͘ with Reynolds number.Such a possibility is worth to be addressed in future investigations.
Similarly, a thorough comparative Voronoï analysis of zero acceleration points and particles clusters as well as a parametric Reynolds number study would shed more light on our results.The former is only accessible to numerical works so far, and regarding the latter, a drastic and challenging reduction of the polydispersity of the seeding fog is crucial for any experimental attempt.Numerical studies are naturally monodisperse and could also be of great help to this regard.Incorporating gravity effects in these simulations would also be useful regarding the question of the amount of demixing.To finish, we would like to stress the great potential of Voronoï tesselations that offer a whole range of new openings for further investigations of particle laden flows.In particular, combined to classical Lagrangian particle tracking, it allows to follow the local concentration in a Lagrangian frame ͑see Fig. 9 that illustrates preliminary attempts of such Lagrangian local concentration tracking͒ giving simultaneous access to statistics of particle dynamics ͑velocity and acceleration͒ and local concentration field around particles.Other important aspects concern clusters dynamics, and more specifically clusters lifetime as well as statistics of particles residence time inside clusters ͑Lagrangian evolution of local concentration in Fig. 9 shows, for instance, two particles traveling between clustered and void regions͒.Such conditional and dynamical information are expected to be key ingredients to improve our understanding and modeling capabilities for turbulent particle laden flows.

FIG. 3 .
FIG. 3. Left: a typical raw image.Right: particles located in this image and the associated Voronoï diagram.For clarity, we show only one third of the full acquired image.

FIG. 5 .
FIG. 5. Standard deviation of Voronoï areas as a function of average Stokes number.͑a͒ One point corresponds to one single experiment, lines connect different experiments for which Reynolds number and the number of particles per image is identical.͑b͒ Only the Reynolds number is constant along lines, and each point is estimated as the averaged standard deviation from experiments with same Stokes number ͑but possibly different C 0 ͒.Error bars represent the dispersion between such experiments.

FIG. 6 .
FIG. 6.A way to identify clusters.͑a͒ Superposition of the Voronoï area PDF for a typical experiment ͑R = 85, St= 0.33, and C 0 = 500 particles per image͒; ten continuous lines associated to ten sets of 500 UIVD are represented ͑dispersion is negligible͒ and a RPP ͑dotted line͒.͑b͒ Ratio of the two PDFs presented on the left figure.Vertical dash-dotted lines indicates 2 ͑left͒ and L 2 ͑right͒.͑c͒ Colored visualization of clusters ͑dark gray͒ and voids ͑light gray͒.

FIG. 7 .FIG. 8 .
FIG. 7. ͑a͒ PDF of clusters area.The inset shows the evolution of the fitted power law exponent with Stokes number for the 40 experiments of Fig. 4͑a͒.Vertical dash-dotted lines indicates 2 ͑left͒ and L 2 ͑right͒.͑b͒ Geometrical characterization of clusters for the same 40 experiments.

)
FIG. 4. ͑a͒ PDF of dimensional Voronoï area A for 40 experiments spanning all R , St, and the volume loading explored.͑b͒ PDF of normalized Voronoi area A / A ¯for four experiments at R Ӎ 96 with varying Stokes number ͑each PDF is calculated from 5000 instantaneous fields͒; inset displays a zoom around the maximum.͑c͒ Centered and normalized PDF of the logarithm of Voronoï area for the 40 experiments from upper figure; black dashed line represents a Gaussian distribution.͑d͒ Test of the log-normal scalings ͑1͒ ͑gray dots͒ and ͑2͒ ͑black dots͒ on 80 experiments spanning all R , St, and C 0 ; one point corresponds to one experiment.If V has a log-normal distribution, the dots should gather on the dash-dotted lines.Relation ͑1͒ is very well verified; relation ͑2͒ deviates slightly from the lognormal scaling, but is still relatively well verified considering the severeness of the test that involves second order moments of the variable logarithm for which statistical convergence becomes challenging.