7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Inertial particles clustering in turbulent flows : a Voronod analysis
R. Monchaux* and M. Bourgoint and A. Cartelliert
ENSTA ParisTechUnit6 de M6canique, Palaiseau, 91761, FRANCE
t Laboratoire des 6coulements gdophysiques et industries, Grenoble, 38000, France
romain.monchaux@tensta.fr
Keywords : Grid turbulence, Lagrangian Particles Tracking, particleladen gas, preferential concentration, Voronoi diagrams
23 mars 2010
Study of inertial particles laden flows is relevant to many industrial and environmental issues, but it is also of fun
damental interest. Since the obtention of the equations governing the particles dynamics have only been approached
in some limited situations (Maxey and Riley (1983)), different experimental investigations focus on the single particle
problem by measuring the Lagrangian acceleration in order to get insight on the relevant forces acting on the particle.
In the many particle case, a striking feature is the tendency to preferential concentration that has been observed for
long (Eaton and Fessler (1994)) and is still studied experimentally and numerically. An other interesting feature is
the enhancement of the settling velocity of the particles in turbulent flows. Since the explanation of this phenomenon
through the existence of clusters has been proposed by Aliseda et al. (2002), different authors tried to quantify and
characterise this clustering numerically. Nevertheless, the "kinematic" simulations used neglect most of the forces
acting on the particles that may play an important role in the particles dynamics and experimental investigations are
required.
Do clusters exist as a whole in these flows ? How do they form? What is there structure and how does this structure
evolve with time ? Which effect do they have on the single particle dynamics ? Here are some questions that have to
be answered. To date, the preferential concentration/cluster problem has been studied with global tools such as box
counting methods, pair correlation function estimation or topological indicators. Nevertheless, such tools do not allow
the required simultaneous dynamical study of the Lagrangian dynamics of the particles and of the local concentration
field.
We propose to use Voronoi tesselations that give a measure of the local concentration field at the interparticles length
scale : through Lagrangian measurements, we compute simultaneously the velocity, the acceleration and the local
concentration field associated to each particle. We focus on the preferential concentration problem on a statistical
ground and we show that previously obtained results can be recovered with the Voronoi diagrams as a single tool :
preferential concentration is quantified and reaches a maximum for Stokes numbers around unity, clusters show not
to have any typical length scale and their geometry exhibits fractal structure. This analysis is a first step before stu
dying the dynamics of the clusters/particles with simultaneous measurements of Lagrangian velocity, acceleration and
concentration fields allowed by the Voronoi analysis. We work in a wind tunnel grid turbulence experiment where the
flow is seeded with a dense water fog made of droplets whose diameters may range from 15 to 150 pm. The Reynolds
number R, is of order 100. Measurements are performed by Lagragian tracking of particles with a high speed camera
sampling the flow at 10 kHz.
1 Experimental setup and postprocessing
1.1 Two phases flow generation and characterization
Illustrations of this paragraph are given in figure 1.
Experiments are conducted in a large wind tunnel with a square crosssection of 0.75 m x 0.75 m where an almost
ideal isotropic turbulence is generated behind a grid whose mesh size is 7.5 cm. We can adjust the mean velocity from
3 to 15 m.s the turbulence level remaining relatively low, of the order of :'. The inertial particles are a water drop
plets fog generated by four industrial injectors placed in the convergent part of the wind tunnel, one meter upstream
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
x (cm) Vo (m.s 1) 1 (m .s3 L (mm) A (mm) r (mm) ,, (ms) v, (cm.s1 R
290 2.92 0.0084 61.4 12 0.79 42 1.90 72
290 4.51 0.0309 61.4 9.7 0.57 22 2.59 90
290 6.33 0.0854 61.4 8.2 0.45 13.3 3.38 114
TAB. 1: Evolution of the turbulence characteristics for the single phase flow for three different velocities and at the
two downstream locations limiting our measurement volume. x : distance downstream the grid, Vo : mean longitudinal
velocity, e : dissipation rate, L, A, j : integral, Taylor and Kolmogorov length scales, T,,, v, : Kolmogorov time and
velocity scales, last column display the Taylor scale based Reynolds number.
the grid to insure an homogeneous seeding of the flow. These injectors consist of two tubes carrying water and air
that merges into a specifically designed nozzle (a SchlickDiisen 942 two fluids nozzle) where the air under pressure
atomizes the liquid and form the exiting conic jet. The particles properties relevant to our problem can be more or less
adjusted : mass loading is directly linked to the water flow rate while the diameter distribution is a subtle function of
both water flow rate and air pressure.
We have performed PIV and hotwire measurements in the test section in order to characterise the turbulence under
lying our experiment. Hotwire acquisitions have been done for the single phase air flow while PIV was carried out
from high speed imaging (we use a phantom V12 camera) of inertial particles laden flows. Hot wire measurements
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 varies from 72 to 114 at the measurement location situated 3 m
downstream the grid. Evolution of the turbulence characteristics of the single phase flow at different locations downs
tream the grid and for the explored range of mean velocities is given in table 1.
We present PIV vertical profiles of the axial velocities averaged over x along the whole measurement volume showing
the relative homogeneity of the flow along the cross section of the tunnel, the very slow evolution of these profiles
with the axial coordinate x along the measurement volume evidences the slow decay of turbulence within this volume.
When the injectors blow only air, turbulence level increases by about .1' compared to the case without injectors, ne
vertheless, this enhancement of turbulence intensity does not depend much on the air pressure imposed except maybe
at very low speed.
The polydispersed seeding of our flow has been estimated using the SprayTech instrument developed by Malvern
and based on the diffusion of a laser beam by the spray. These granulometry estimations have been performed in
dependently of any other acquisition, but within the same measurement volume. The control parameters governing
the diameter distributions are : the air pressure, the water flow rate and the wind velocity in the tunnel. From these
measurements, we build an average Stokes number defined as St (d/r)2 (1 + 27)/36, where d is the diameter cor
responding to the maximum of the particles diameter distributions similar to those of fig. 1, rj is the Kolmogorov length
scale and 7 is the ratio of the particles density to the fluid density. An error bar on this average Stokes number built
from the particles diameters root mean square is of order ".~' and will not be displayed on the figures presented in
the sequel. One experiment is characterized by a Reynolds number which is properly estimated from velocity measu
rements, an average Stokes number obtained as described above, the mass loading (p, < 10 5 in all our experiments)
and the Rouse number defined as the ratio between the terminal velocity of the particles and the turbulent fluctuations
intensity (in our experiment, Rouse number ranges from 0.4 to 2). At fixed RA and St, we assume that mass loading
is directly related to the number of particles per image.
1.2 Acquisitions and postprocessing
Acquisitions are performed using a phantomvl2 high speed camera operated at 10 kHz and acquiring 12 bits images
at a resolution of 488 pixels x 1280 pixels corresponding to a 45 mmx 125 mm visualisation window centered with
respect to the wind tunnel section and situated 2.95 m downstream the grid. The camera is mounted by a 105 mm
macro Nikon lens opening at 2.8. A 8 W pulsed copper Laser synchronised with the camera is used to generate a
2 mm thick light sheet illuminating the field of view in the streamwise direction. The camera aims at the measurement
volume with an angle of 50 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 at wind ve
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
4 2 0
Vx (m.s 1)
U
mean flow
direction
0.08
2.3 m/s
V 3.9 m/s
0.07 A 7 m/s
V
V V
A A A
3 4 5
injector air pressure
C)
0.06
0.05
0.04
4 2 0 0.03L
Vx (m.s1 ) 2
10
0
x 8
a)
E 6
> 4
UL
'i ~.
I
I
I
I
I
I
ii
0 1
10 10 1
droplet diameter
P=5
P=4
....... P=3
P=2
.I
2'3
i
02 103
(g m)
S1.8
E
1.6
I 1.4
0
I1.2
3 4
air pressure (bar)
FIG. 1: Top : experimental setup. Middle, from left to right, for three different mean velocities : mean axial velocity
averaged over x in the measurement volume (x e [289 301] mm); mean axial velocities at different locations x in the
measurement volume; evolution of the ratio of turbulent fluctuations (V,,, ) to velocity average (< V >) as a function
of the injector pressure without water flow in (note that the impact of the injectors on the turbulence level would even
be lower with the presence of water). Bottom left : particles diameter Probability Density Function evolution with the
air pressure (P varies from 2 to 5 bars) at fixed water flow rate (1.2 1/mn for each injector) and fixed wind velocity
(Vo = 4.5 m.s 1); right : evolution of the maximum of these PDF with the water flow rate and the air pressure at the
same fixed velocity.
44a.
Laser.
sheet
30
20
10

Phantom V12
highspeed X
camera
VO=6 3ms
VVO=45ms
VO=2 9 ms
0
>,
30L
I t . ....
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
700
7..70,* .* .
.._. 600
60
E
x (pixels) x (mm)
FIG. 2: Left: a typical raw image. Right : particles located in this image and the associated Voronoi' diagram. For
clarity, we show only one third of the full acquired image.
locity, water rate and air pressure fixed. Particles are identified to local maxima of the pictures intensity higher than a
prescribed threshold. The impact of this threshold on the number of particles found is weak around the used threshold.
Subpixel accuracy detection is obtained by locating the particles at the center of mass of the pixels surrounding the
local maxima.
We aim at studying the particles concentration fields in order to quantify the preferential concentration effects and
to identify the clusters if any. Usually the pair correlation function is computed to quantify the preferential concentra
tion effects while a box counting method is preferred to gain access in a local concentration field. In this paper, we
propose the use of a single tool to tackle these two problems : the Voronoi' diagrams. Such a Voronoi' diagram is the
unique decomposition of the 2D space into independent cells associated to each particle. It is defined as the ensemble
of points that are closer to a particle than to any other. Use of Voronoi' diagrams is very classical in the study of granu
lar systems and has also been used to identify galaxies clusters. To compute the Voronoi' diagrams, we have used the
Matlab algorithm that shows to be very fast with the typical number of particles per image we had to process. Figure
2 shows a raw acquired image as well as the located particles and the associated Voronoi diagram.
2 Preferential concentration evidence and quantification
2.1 Voronil diagrams
From the definition of the Voronoi diagrams, it appears that the area of a Voronod cell is the inverse of the local
concentration in particles, thus the studies of the local concentration field obtained by any mean or from the Voronoi
area field are equivalent. As preferential concentration implicitly means preferential concentration with respect to a
uniform random Poisson process (RPP), we aim at quantifying the former as the departure from a uniform distribution
of particles. Ferenc and Neda (2007) give a review of the few known properties of the Voronoi diagrams associated
to RPP : the probality density functions (PDF) of the Vorono areas are well fitted by Gamma distributions involving
the space dimension as a single parameter (empirical result), the second moment of these distributions in the 2D case
is equal to 1.280 (exact result), the first moment has nothing to do with the organisation of the particles since the
average Voronoi area is equal to the area of the considered plane divided by the number of particles present (in the
sequel, PDF of Voronoie areas are normalized to be of unit mean).
2.2 Experimental Voronig area PDF
Postproc Poio oessing of the data described in 1.2 gives for each experiment and at each time step a Vorono diagram. To
increase size of statistics sets, we gather Voronoif diagrams associated to several images. Left figures in fig. 3 display
the PDF of the Voronoi cells area associated to 40 different experiments. These PDF are obtained by analysing 5000
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
images per experiment. Considering the PDF of the dimensional Voronoi areas, one sees the maxima of these PDF
spanning over two decades associated to the number of particles per image that evolves from 50 to 5000, note that as
the average number of particles per image is getting lower, the scattering of the right tail on these PDF increases. The
same PDF normalized to be of unit mean collapse somehow (not shown). The tails associated to large Voronoi areas,
i.e. to regions of low particles density, actually collapse within measurements uncertainty whereas the tails associated
to relatively small Voronoi area, i.e. to regions of high particles concentration, show a much stronger dependency with
the considered experimental conditions, the Stokes number for example as shown on fig. 3 top right. The sizes of the
large almost empty regions of the flow are thus probably associated to large eddies whose sizes depend only on the
integral scale which does not vary when the Reynolds number increases in the wind tunnel while the dissipation scale
is getting smaller.
Fig. 3 presents the centeredreduced PDF of the logarithm of the Voronoi areas measured in our experiment. These
PDF correspond to lognormal distributions in the interval 2a, where a stands for the standard deviation of these
distributions. This assertion is also tested in fig. 3 bottom right where we plot simultaneously the mean of the logarithm
of the Voronoi area as well as the expected lognormal scaling of the two first moments of the distributions as a function
of the logarithm of the mean Voronoi area (see caption for details). On this figure, one experiment is represented by one
point, evidencing the almost perfect verification of the lognormal scaling law. We can not give so far any theoritical
interpretation of this lognormal scaling, but this result shows that at first order, the normalized Voronoi area PDF can
be described with a single scalar quantity : the standard deviation of the normalized area distributions.
Figure 4 shows the evolution of the normalized Voronoi area standard deviation as a function of the Stokes number.
We present two figures : on the left one, one point corresponds to one experiment, lines connect different experiments
for which Reynolds number and number of particles per image are identical. On the right one, along lines, only
Reynolds number is constant, one point averages experiments with same Stokes number and same number of particles
per image. Error bars are the resulting dispersion between such experiments. A standard deviation of the normalized
areas of 0.52 corresponds to the RPP analytical value of 1.280 described earlier and defines a reference basis for
our measurements. If the measured standard deviation significantly exceeds this 0.52 value, events of high and low
concentration compared to the RPP case are present and oppositely, standard deviations below this value evidence a
tendency to a more organised particles distribution, this standard deviation equaling zero in the limit of a perfect cristal.
From St 0.2 to St 6, the standard deviations of the Voronoi areas in our experiments vary from 0.57 to 0.84
evidencing preferential concentration over this entire range with a maximum found for St 2 3. We have performed
box counting and correlation dimension methods on the same data and we report with these tools a maximum effect of
preferential concentration around the same values when people usually report maximum of preferential concentration
for St 1. This discrepancy can not thus be attributed to our postprocessing investigation but rather to a small
decrease of r due to the injection process. Note that for large Stokes numbers, the number of particles per image is
quite low (typically < 1000) when it usually reaches 3000 4000. This could question the relevance of the standard
deviation calculation performed as the statistical set is shrinking, but we have shown that this calculation is robust :
when particles are randomly removed from a given set of images, the standard deviation is reduced by 1.5% when the
number of particles used for the calculation is divided by 5. It has to be noted that the Voronoi PDF do not contain
any information about the structure of the concentration fields which is a key point of the dynamical study we want
to conduct in the future. Next paragraph shows how one can however use these distributions to get an insight into the
concentration field structure and dynamics.
2.3 Clusters identification and characterisation
One open question is to know wether clusters exist as a whole in turbulent flows laden with inertial particles. To
address this question, one has to study dynamically the behaviour of areas of high concentration in particles, but former
to this study, one has to identify possible candidates to clusters. Voronoi PDF may be used in this context. Top left
figure on fig. 5 presents the superimposition of the Voronoi PDF associated to a typical experiment and to a RPP, these
PDFs intersect twice (which is more visible on the top right figure that shows their ratio) : for low (resp. high) values
of the area corresponding to high (resp. low) values of the concentration, the experimental PDF is above the RPP one
while we observe the opposite for intermediate values which is in good agreement with the usual picture of preferential
concentration. We define the intersection points as intrinsic definitions of the clustering limits : for a given experiment,
Voronoi cells whose area is smaller than the first intersection are considered to belong to a cluster and those whose area
is larger than the second intersection are associated to "holes". Figure in the middle of fig. 5 displays a full Voronoi
diagram corresponding to one image taken from one experiment. On this diagram, cells corresponding to clusters
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10
LL 10:
10
3
)4
35
6
101 10 101
Voronoi Area
102
103
10
LL 2
S10
103
10
i4
4 2 0 2
log(Normalized Voronoi Area)
05
0
Y)4
*..
0 1
og(
1 2
log(<
2 3
Voronoi area>)
FIG. 3: PDF of Voronoi areas (from top left to bottom right) : dimensional area for 40 experiments spanning all Rx,
St and mass loading explored; normalized Voronoi area for 4 experiments at Rx 96 with varying Stokes number
(each PDF is calculated from 5000 instantaneous fields), inset displays a zoom around the maximum; logarithm of
normalized area (40 experiments of figure top left) and a gaussian. Bottom right: test of the lognormal scaling on 80
experiments spanning all Rx, St and mass loading explored, one point correspond to one experiment: black dots :
log(V) vs log(V), grey dots : log(V) vs log(V) or(V)2/2, where V stands for the Voronoi areas and on(V) for the
standard deviation of V. If V has a lognormal distribution, the grey dots may gather on the dashdotted line.
100
101
102
LL
o10
10
10
10
r4
4 5
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.9 0.9
V  R=72
CU CU, R =90
20.8 2 0.8
CU CU
F _R = 114
0.6 l:* i 0.6
E E
0 0
0.5 0.5
0 2 4 6 0 2 4 6
averaged Stokes numbers averaged Stokes numbers
FIG. 4: Standard deviation of Voronoi areas as a function of average Stokes number. Left: one point corresponds to
one experiment, lines connect different experiments for which Reynolds number and number of particles per image are
identical. Right : along lines, only Reynolds number is constant, one points averages experiments with same Stokes
number and same number of particles per image. Error bars are the resulting dispersion between such experiments.
(resp. to holes) have been coloured in black (resp. light grey) while the remaining cells have been patched with white.
It appears that black cells tend to form groups of various sizes and shapes that we identify to clusters whenever they
belong to the same connected component. The same algorithm is used for the grey cells to isolate depleted regions.
We now turn to the statistical properties of the identified clusters and empty regions. On each image, we have identified
the different connected component of low/high concentration cells to identify the empty regions/clusters (note that
this procedure applied to synthetic RPP fields does not lead to any cluster identification since the largest connected
component never contains more than 10 particles). Then, we have computed their area (A) and perimeter (P) from the
area and perimeter of the Voronoi cells that constitute them. In fig. 5, we present the distributions of the clusters area
as well as P as a function A1/2 for the identified clusters. This last figure shows that for small and compact clusters, P
and A1/2 are almost linearly related while large structures exhibit a fractal behaviour which has already been reported
experimentally and numerically (see Aliseda et al. (2002) and references herein). The clusters sizes (observed from
the A distribution) are algebrically distributed with an exponent 2 which seems independent of both RA and St
within our measurement accuracy (see inset). This implies that clusters do not have typical nor average size which
would be contradictory with the similar study conducted earlier by Aliseda et al. (2002) who reported an exponential
distribution. The same analysis has been conducted on the empty regions (no figure shown) and leads to very similar
results : large empty regions have fractal structure and their sizes are algebrically distributed with an exponent close to
1.8. This last result is in agreement with the one obtained by Yoshimoto and Goto (2007) in a DNS study of exactly
void regions formed by inertial particles.
3 Conclusions and outlooks
This article proposes the Voronoi diagrams as a new tool to study preferential concentration in particles laden flows.
We argue that this method allows a quantitative insight into formerly known results. We have quantified the preferen
tial concentration with a single scalar quantity evidencing the maximum heterogeneity of the concentration field for
averaged Stokes numbers of order unity. We have proposed a way to identify clusters and depleted regions and we have
characterized their geometries : both of them have sizes algebraically distributed and exhibit fractal structure which
testify of the selfsimilatity of preferential concentration in particles laden turbulent flows. All these results obtained
with a single tool are consistent with former experimental and numerical works based on more classical diagnosis
tools. This pushes to consider the Voronoi tesselations as an efficient tool in this context and incites to go further and
use the great potential of this analysis to follow the concentration field in a lagrangian frame. As Voronoi diagrams
allow to study simultaneously the concentration field and the particles dynamics, we will address in the future the rele
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
UL 101
D
S10o
101 100 1C
normalized Voronoi area
101
102
101 10 101 102
normalized Voronoi area
50 x(mm) 100
103
()
C
a2
U)
S101
10
100
101 102 103
10 101
(cluster area)1/2
FIG. 5: A way to identify clusters. From top to bottom. Left: superimposition of the Voronoi areas PDF for a typical
experiment (10 continuous lines associated to ten sets of 500 UIVD : dispersion is negligible) and a RPP (dotted line).
Right : ratio of the two PDF presented on the left figure. Middle : colour visualization of clusters (black) and holes
(light gray). Middle left : PDF of clusters area, inset shows the evolution of the fitted slope with the Stokes number
for the 40 experiments of fig. 3. Middle right : geometrical characterisation of clusters for the same 40 experiments.
Dashdotted lines on topright and bottomleft figures indicates r2 and L2.
100
10
10
L_ 02
S10
13
104
50
E
E
0
101
102
UL 
S10
10
150
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
vant matters of clusters lifetime, clusters dynamics and statistics of particles velocity and acceleration conditioned
to the local concentration. It would also be of great interest to reanalyse former numerical and experimental studies
in the Voronoi framework.
This research is supported by ANRDSPET. We would like to thanks other members of the ANR for fruitfull dis
cussions and Raphael Candelier for providing an efficient way to identify connected components.
References
A. Aliseda, A. Cartellier, F. Hainaux, and J. C. Lasheras. Effect of preferential concentration on the settling velocity
of heavy particles in homogeneous isotropic turbulence. Journal ofFluidMechanics, 468 :77105, October 2002.
J. K. Eaton and J. R. Fessler. Preferential concentration of particles by turbulence. Int. J Multiphase Flow, 20
169209, 1994.
M. R. Maxey and J. J. Riley. Equation of motion for a small rigid sphere in a nonuniform flow. Physics ofFluids, 26:
883889, April 1983.
J. S. Ferenc and Z. Neda. On the size distribution of Poisson Voronoi cells. Physica A Statistical Mechanics and its
Applications, 385 :518526, Nov. 2007.
H. Yoshimoto and S. Goto. Selfsimilar clustering of inertial particles in homogeneous turbulence. Journal of Fluid
Mechanics, 577 :275+, April 2007.
