Voronoi tesselation analysis of sets of randomly placed finite-size spheres
MM. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020
Vorono¨ı tesselation analysis of sets of randomly placedfinite-size spheres
Markus UhlmannInstitute for HydromechanicsKarlsruhe Institute of Technology76131 Karlsruhe, Germany [email protected]
April 22, 2020
Abstract
The purpose of this note is to clarify the effect of the finite size of spherical particles uponthe characteristics of their spatial distribution through a random Poisson process (RPP).This information is of special interest when using RPP data as a reference for the analysisof the spatial structure of a given (non-RPP) particulate system, in which case ignoringfinite-size effects upon the former may yield misleading conclusions. We perform MonteCarlo simulations in triply-periodic spatial domains, and then analyze the particle-centeredVorono¨ı tesselations for solid volume fractions ranging from 10 − to 0 .
3. We show that thestandard-deviation of these volumes decreases with the solid volume fraction, the deviationfrom the value of point sets being reasonably approximated by an exponential function. Ascan be expected, the domain size for which the random assemblies of finite-size particles aregenerated has a constraining effect if the number of particles per realization is chosen toosmall. This effect is quantified, and recommendations are given. We have also revisited thecase of random point sets (i.e. the limit of vanishing particle diameter), for which we haveconfirmed the accuracy of the earlier data by Tanemura [Forma, 18(4):221247, 2003].
One method to characterize the spatial distribution of particles is with the help of Vorono¨ı tessela-tion. As initially proposed in the context of particulate multi-phase flow by Monchaux et al. (2010)and further elaborated by Monchaux et al. (2012), such a tesselaton can give valuable informationon the tendency of particles to cluster. Apart from partitioning space and therefore allowing foran interpretation as the inverse of a local concentration, the tesselation further generates data onthe connectivity (“neighborhood”) between particles, and it can therefore be directly used for thepurpose of cluster identification and for quantifying Lagrangian aspects of clustering. Anotherbenefit is the intrinsic definition of a clustering threshold (Monchaux et al., 2010) without theneed for a priori selection of a length-scale. A further important advantage is the availability of afast algorithm for performing the tesselation (in the present note we either use the Matlab imple-mentation of the “QHULL” library, Barber et al. 1996, or the “VORO++” library, Rycroft 2009,both of which provide a scaling of the execution time which is linear in the number of particles).For these reasons, various analysis methods based upon Vorono¨ı tesselation of particle positionshave found relatively wide-spread use in the particulate flow community (Obligado et al., 2011;Garc´ıa-Villalba et al., 2012; Fiabane et al., 2012; Tagawa et al., 2012; Dejoan and Monchaux, 2013;Kidanemariam et al., 2013; Uhlmann and Doychev, 2014; Sumbekova et al., 2016; Uhlmann andChouippe, 2017; Monchaux and Dejoan, 2017; Chouippe and Uhlmann, 2018).1 a r X i v : . [ phy s i c s . f l u - dyn ] A p r . Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020The usual procedure in most of these afore-mentioned studies consists in comparing the statis-tics of the actual particle distributions in the respecitve multiphase flow systems (either obtainedfrom experimental measurements or from numerical simulations) with data for particle distribu-tions obtained by a random Poisson process (henceforth denoted as “RPP”), i.e. drawing themrandomly from a distribution which is uniform in space. Statistically significant differences arethen a sign of “structure” in the particle set, and these features can subsequently be interpretedon physical grounds.The spatial distribution of a set of points through an RPP has been investigated in the literature(Tanemura, 2003; Ferenc and Neda, 2007, and references therein). Although no analytical resultsare available, these authors have proposed empirical fits for the probability distribution of theVorono¨ı cell volumes as well as providing reference values for the first few moments which arewidely used. In some applications the shape of the probability density function (pdf) of Vorono¨ıcell volumes does not change much, which means that the data is reasonably described by a singleparameter, i.e. the second moment. It is then common practice to use the difference between themeasured standard deviation and the literature value of an RPP-generated set of points as a proxyfor the tendency to cluster (Monchaux et al., 2010).Due to the finite size of real-world particles, theoretical results derived for the spatial distri-bution of points (such as those given by Ferenc and Neda, 2007) do no longer strictly apply. Thisis due to the fact that there exists a lower bound for the closest packing of finite-size particles asthey cannot overlap. This aspect as well as the impact of a finite domain size (both in numeri-cal simulations and in laboratory-experimental measurements) has so far not been systematicallydocumented. Oger et al. (1996) have observed that the probability distribution of Vorono¨ı cellvolumes changes from a Gamma distribution to a Gaussian when going from one extreme in termsof the solid volume fraction to the other (i.e. from a point set to a close sphere packing).Incidentally, let us note that Tagawa et al. (2012) have reported on the effect of the numberof particle samples upon the standard deviation of Vorono¨ı cell volumes determined for numericaldata-sets of point particles. Monchaux (2012) have analyzed biases in Vorono¨ı-based analysisof experimental data, focusing on: (i) the effect of projecting onto two spatial dimensions frommeasurements with a laser sheet with non-vanishing thickness; (ii) sub-sampling due to “missed”particles; (iii) slight polydispersity.In the present note we first revisit the case of point sets, before addressing the followingquestions pertaining to the statistics of spherically-shaped particles distributed through a randomPoisson process: what is the effect of a non-vanishing particle size? What is the effect of a finitesolid volume fraction? What is the influence of the relative domain size? We consider the case of N p points distributed inside of a cubical domain of side-length L . Forthe sake of generality, we consider N s independently drawn such distributions (i.e. snapshots).Since there is no second length-scale, the problem is then completely determined by the twonon-dimensional numbers N p and N s .The Vorono¨ı tesselation is performed by assuming periodicity of the field over all three spatialdirections. This can be (naively) implemented algorithmically by (redundant) periodic extensionof the fields and later eliminating the duplicate cells on the borders of the fundamental domain.Figure 1 shows the probibility density function obtained by numerical experiments with N tot = N p N s = 5 · samples, by varying the number of particles per snapshot in the range N p =10 . . . and adjusting the number of snapshots N s accordingly. It can be seen that the differ-ences are marginal. The pdf can be fitted to a three-parameter (generalized) Gamma distribution,viz. f ( x ) = c b a/c Γ( a/c ) x a − exp( − b x c ) , (1)2. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020 pd f ( a ) -2 -1 -6 -4 -2 V vor / (cid:104) V (cid:105) pd f ( b ) -6 -4 -2 -6 V vor / (cid:104) V (cid:105) Figure 1: Probability density function of the volume of Vorono¨ı cells for sets of points distributedaccording to an RPP. The number of points per snapshot is indicated by the colors as follows:, N p = 10 ; , N p = 10 ; , N p = 10 ; , N p = 10 ; , N p = 10 ; The total number ofsamples is kept constant ( N tot = N p N s = 5 · ) in the entire series by adjusting the numberof snapshots N s accordingly. Graph ( a ) shows double logarithmic scaling; ( b ) is the same data insemi-logarithmic scaling. The inset shows a zoom of the right tail. σ ( V v o r / (cid:104) V (cid:105) ) N p Figure 2: Standard deviation of the volume of Vorono¨ı cells for sets of points distributed accordingto an RPP, shown as a function of the number of particles N p per snapshot. The total number ofsamples is kept constant ( N tot = N p N s = 5 · ) in the entire series by adjusting the number ofsnapshots N s accordingly.with the following set of coefficients borne out by our data with N p = 10 and N s = 500: a = 4 . , b = 4 . , c = 1 . . (2)These values are very close to those determined by Tanemura (2003) who proposes a = 4 . b = 4 . c = 1 . Figure 2 shows the influence of the number of points per snapshot upon the standard-deviationof the Vorono¨ı cell volumes. It can be seen that the influence becomes insignificant for N p ≥ particles. For smaller numbers, it turns out that the standard-deviation is underpredicted. This We believe that there is a typo in Ferenc and Neda (2007): on their page 524 they list their three-parametergamma fit with coefficient values as a = 3 . b = 3 . c = 1 .
3. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020 σ ( N s ) ( a ) N s ( σ ( N s ) − σ ( )) / σ ( ) ( b ) -10 -8 -6 -4 -2 N s Figure 3: Convergence of the standard deviation of the volume of Vorono¨ı cells for sets of pointsdistributed according to an RPP, shown as a function of the number of snapshots N s , for N p = 10 .( a ) running average; ( b ) relative error with respect to the value obtained for the entire sequence.is due to the small deviations in the right tails (for large cell volumes) which are visible in theinsert of figure 1( b ). The differences are presumably caused by an under-representation of verylarge cell volumes in ensembles with smaller numbers of particles.Let us briefly mention the convergence of the second moment with the number of samples N s . Figure 3 exemplarily shows the data for the case with N p = 10 . Convergence is obviouslynot regular (since there is no order in the sequence), and, therefore, it needs to be monitored inpractice.As a bottom line, we state our best approximation of the standard-deviation of Vorono¨ı cellvolumes for sets of points drawn according to an RPP as follows: σ ( V vor / (cid:104) V (cid:105) ) = 0 . . (3) Now let us turn to the case of a distribution of a set of N p spherical, mono-dispersed particles withdiameter D , again placed by means of an RPP in a cubical box of side-length L . Each draw ofparticle positions is done through a simple Monte Carlo method (Owen, 2013), performed underthe constraint that no two particles should overlap. This is performed for each set by first drawingindividual particle positions sequentially from a uniform and independent probability distribution(as in the case of sets of points treated above), testing each new candidate for geometrical over-lap with all previously determined positions in the set, and, if overlap is found, discarding thatcandidate, then re-drawing, and so on until no overlap is detected (cf. algorithm 1). The presentprocedure is equivalent to random sequential addition (RSA) of hard spheres, as sometimes used inthe context of physical chemistry in order to obtain random sphere arrangements (Widom, 1966).In particular, it has been shown that RSA leads to a so-called jamming limit (i.e. an upper boundfor the solid volume fraction Φ s where no additional spheres can be placed without overlapping)for Φ ( RSA ) s ≈ .
382 (Talbot et al., 1991). Note that the RSA limit is significantly lower than thedensest possible sphere packing (Φ ( densest ) s ≈ . .
3; for values which are even closer to Φ ( RSA ) s it becomes challenging todetermine a sufficiently large number N s of particle assemblies, since the required number of draws N draws in algorithm 1 eventually increases exponentially with Φ s . Ferenc and Neda (2007) state a value of σ = 0 . σ = 0 .
4. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020
Algorithm 1
Determines a single snapshot, i.e. a spatial distribution of N p finite-size particleswhich do not overlap. The number of required draws is stored in the variable N draws . N draws = 0 for i = 1 . . . N p do (cid:46) loop over particle setredraw=true while redraw do draw X ( i ) α from uniform distribution ∀ α = 1 , , N draws ← N draws + 1 if X ( i ) does not lead to overlap with any particle j = 1 . . . i − then redraw=false end ifend whileend for σ ( V v o r / (cid:104) V (cid:105) ) -5 -4 -3 -2 -1 Φ s = 0 (points)Φ s Figure 4: Standard deviation of the volume of Vorono¨ı cells for N p = 10 finite-size particlesdistributed according to an RPP, as a function of the solid volume fraction Φ s . The ratio ofbox-size to particle diameter was varied accordingly, i.e. L/D = ( N p π/ (6Φ s )) / . The number ofsnapshots analyzed in each case was N s = 2000, such that the total number of samples for eachdata point measures N tot = 2 · . Included is the value σ (Φ s = 0) for an RPP distribution ofpoints, as given in (3). The black dashed line indicates the exponential fit (5).Here the relevant variables are N p (number of particles in the set), N s (number of particleassemblies analyzed), D (particle diameter) and L (edge length of the computational cube), whichleads to similarity under three non-dimensional parameters. Those can be chosen as: N p , N s and L/D . Alternatively, one can use the solid volume fractionΦ s = π N p (cid:18) DL (cid:19) , (4)for the characterisation of a parameter point (instead of either N p or L/D ), and/or the totalnumber of samples N tot = N p N s . In a first series we vary the solid volume fraction Φ s , while keeping the number of particles persnapshot constant at N p = 10 , i.e. varying the ratio of box-size to particle diameter accordingly,and systematically using a number of N s = 2000 snapshots ( N tot = 2 · ). Figure 4 shows theprogressive deviation of the standard-deviation of Vorono¨ı cell volumes from the point-set valuewith increasing solid volume fraction. This deviation can be reasonably fitted with the aid of an5. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020 pd f ( a ) -2 -1 -6 -4 -2 V vor / (cid:104) V (cid:105) pd f ( b ) -6 -4 -2 V vor / (cid:104) V (cid:105) Figure 5: Probability density function of the volume of Vorono¨ı cells for sets of spherical, finite-size particles distributed according to an RPP, shown for the data of figure 4. The line-styles areas follows: , Φ s = 10 − ; , Φ s = 10 − ; , Φ s = 10 − ; , Φ s = 5 · − ; , Φ s = 10 − ;, Φ s = 2 · − ; , Φ s = 5 · − . Graph ( a ) shows double logarithmic scaling; ( b ) is the samedata in semi-logarithmic scaling. The solid circles indicate the respective fits to a generalizedgamma distribution (1) for the most dilute and the densest cases.exponential function in Φ s , viz. σ ( V vor / (cid:104) V (cid:105) ) = σ (Φ s = 0) · exp ( − . s ) , (5)which might be useful as an approximation.The corresponding probability density functions of Vorono¨ı cell volumes for these Monte Carloexperiments are shown in figure 5. They can again be fitted reasonably well to three-parameter(generalized) gamma distributions (1), with the shape morphing from the parameter set ( a = 4 . b = 4 . c = 1 .
17) for Φ s = 10 − (i.e. practically identical to the point-set, given in 2) to( a = 16 . b = 31 . c = 0 .
52) for Φ s = 5 · − . Note from (4) that the solid volume fraction isequal to the ratio between the volume occupied by a particle and the mean volume of the Vorono¨ıcells, i.e. Φ s = V p / (cid:104) V vor (cid:105) , where V p = D π/
6. Since for non-overlapping particles a Vorono¨ı cellcannot be smaller than the volume occupied by the particle itself, min ( V vor ) ≤ V p , it followsthat a lower bound on the normalized Vorono¨ı cell volumes is given by the solid volume fractionitself, viz. min ( V vor ) / (cid:104) V vor (cid:105) ≤ Φ s . This lower bound is consistent with what can be observed inthe progressive disappearance of the left tails in figure 5( a ). As a consequence of conservation oftotal volume, the probability of finding large cells correspondingly decreases, i.e. the right tailsare likewise reduced (cf. figure 5 b ). Here we first fix a value for the solid volume fraction Φ s . We then perform simple Monte Carlosimulations for parameter points featuring distinct values of the length-scale ratio L/D . Thismeans that for each parameter pair (Φ s , L/D ) the number of particles per realization is different,i.e. N p = Φ s ( π/ L/D ) . Again we adjust the overall number of samples per parameter point( N tot = N p N s ) by adapting the number of snapshots N s ; for this series we use N tot ≥ · .It can be seen from the pdfs shown in figure 6 that the smallest box-sizes lead to constraintswhich are particularly visible in the far right tails of the distribution (i.e. the probability of findingvery large Vorono¨ı cells is progressively reduced). Figure 7 then shows that the box-size effectupon the standard deviation is rather mild as compared to the effect of the solid volume fractionanalyzed in § pd f ( a ) -2 -1 -6 -4 -2 V vor / (cid:104) V (cid:105) pd f ( b ) -6 -4 -2 V vor / (cid:104) V (cid:105) Figure 6: Probability density function of the volume of Vorono¨ı cells for sets of finite-size particleswith solid volume fraction Φ s = 5 · − , distributed according to an RPP. The relative box-sizeis indicated by the colors as follows: , L/D = 25; ,
L/D = 35; ,
L/D = 150. Graph ( a )shows double logarithmic scaling; ( b ) is the same data in semi-logarithmic scaling.the spatial domains L/D . It can be seen in the figure that the value of
L/D at which the domainbecomes “small” depends upon the solid volume fraction. In figure 8 we have normalized the databy its asymptotic value for large domains (denoted as V vor / (cid:104) V (cid:105) ∗ ). As expected, the curves forthe two different solid volume fractions collapse up to statistical uncertainty; they also collapsewith the normalized curve for the point-set. We can then see that “small” simply means that thenumber of particles is below a certain threshold, say N p ≤ .This threshold can also be argued via a detour, defining a critical length-scale ratio ( L/(cid:96) ) crit below which the domain can be considered as “small” (here “ (cid:96) ” stands for a characteristic lengthscale of the random particle arrangement). Let us assume that the characteristic length can bedefined from the mean cell volume of a Vorono¨ı cell, viz. (cid:96) = (cid:104) V vor (cid:105) / . Then, from the definition (cid:104) V vor (cid:105) = L /N p it directly follows that L/(cid:96) = N / p , indicating that the critical value will onlydepend upon the number of particles, as indeed observed in figure 8. We have investigated random arrangements of point sets and of non-overlapping, finite-size par-ticles in cubical domains assumed to be tri-periodic. The analysis was performed with the aidof Vorono¨ı tesselation, a tool which is frequently employed in the context of particulate flow andother branches of natural sciences and engineering. With this study we have addressed severalaspects that have previously not been documented in a systematic way.First, we have revisited the case of point sets. It was shown that the statistics of a randomarrangement (assigned through a random Poisson process, RPP, which was repeated N s times perparameter point) depend upon the number of particles N p , even if the total number of samples(i.e. the product N p N s ) is sufficiently large. Our present data leads to the following best fit to ageneralized three-parameter Gamma distribution (1): a = 4 . , b = 4 . , c = 1 . , which yields as value for the standard deviation: σ ( V vor / (cid:104) V (cid:105) ) = 0 . . This value is very close to the result stated by Tanemura (2003), while if differs by roughly 3%from the result of Ferenc and Neda (2007). 7. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020 σ ( V v o r / (cid:104) V (cid:105) ) ( a ) Φ s = 0 (points) Φ s = 0 . s = 0 . L/D σ ( V v o r / (cid:104) V (cid:105) ) ( b ) Φ s = 0Φ s = 0 . s = 0 . N p Figure 7: Standard deviation of the volume of Vorono¨ı cells for finite-size particles distributed inspace according to an RPP, shown as: ( a ) a function of the relative box-size L/D ; ( b ) a function ofthe number of particles per box N p . Each line corresponds to a different value of the solid volumefraction, as indicated. The open symbols in blue color in graph ( b ) correspond to the point-particledata of figure 2. σ ( V v o r / (cid:104) V (cid:105) ) / σ ( V v o r / (cid:104) V (cid:105) ) ∗ N p Figure 8: The same data as in figure 7( b ), but the ordinate is normalized with the asymptotic valueof σ ( V vor / (cid:104) V (cid:105) ) for large N p . The line-styles are as follows: , Φ s = 0 (points); , Φ s = 5 · − ;, Φ s = 5 · − .Second, we have analyzed spatial arrangements of finite-size, spherical particles generatedthrough an RPP. The problem depends upon two non-dimensional parameters (apart from thetotal number of samples N p N s ): the solid volume fraction and the relative domain size. Throughnumerical evaluation it is demonstrated that the influence of the solid volume fraction Φ s uponthe standard deviation of the Vorono¨ı cell volumes can be very significant, with the deviation fromthe value for a point-set increasing exponentially with Φ s over the investigated range of 10 − to0 .
3. An approximation which might be useful in practical applications is provided in equation (5).Concerning the box-size dependence, it is found that – as expected – the low-order statisticsof the tesselation are constrained if the box-size is below a certain threshold value. This thresholdis independent of the solid volume fraction if expressed in terms of the number of particles perrealization. If a number N p ≥ is chosen, the relative deviation is kept below 3 · − .8. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020 Acknowledgements
I am indebted to two anonymous referees who have made suggestions which have significantlyimproved the manuscript. The computations were partially performed at SCC Karlsruhe. Thecomputer resources, technical expertise and assistance provided by this center are thankfullyacknowledged.
References
R. Monchaux, M. Bourgoin, and A. Cartellier. Preferential concentration of heavy particles: AVorono¨ı analysis.
Phys. Fluids , 22:103304, 2010. doi:10.1063/1.3489987.R. Monchaux, M. Bourgoin, and A. Cartellier. Analyzing preferential concentration andclustering of inertial particles in turbulence.
Int. J. Multiphase Flow , 40:1–18, 2012.doi:10.1016/j.ijmultiphaseflow.2011.12.001.C.B. Barber, D.P. Dobkin, and H. Huhdanpaa. The quickhull algorithm for convex hulls.
ACMTrans. Math. Softw. , 22(4):469–483, December 1996. doi:10.1145/235815.235821.C.H. Rycroft. Voro++: A three-dimensional Voronoi cell library in C++.
Chaos , 19(4):041111,2009. doi:10.1063/1.3215722.M. Obligado, M. Missaoui, R. Monchaux, A. Cartellier, and M. Bourgoin. Reynolds numberinfluence on preferential concentration of heavy particles in turbulent flows.
J. Phys.: ConferenceSeries , 318(5):052015, dec 2011. doi:10.1088/1742-6596/318/5/052015.M. Garc´ıa-Villalba, A.G. Kidanemariam, and M. Uhlmann. DNS of vertical plane channel flow withfinite-size particles: Voronoi analysis, acceleration statistics and particle-conditioned averaging.
Int. J. Multiphase Flow , 46:54–74, 2012. doi:10.1016/j.ijmultiphaseflow.2012.05.007.L. Fiabane, R. Zimmermann, R. Volk, J.-F. Pinton, and M. Bourgoin. Clustering of finite-sizeparticles in turbulence.
Phys. Rev. E , 86(035301(R)), 2012. doi:10.1103/PhysRevE.86.035301.Y. Tagawa, J Martinez-Mercado, V.N. Prakash, E. Calzavarini, C. Sun, and D. Lohse. Three-dimensional Lagrangian Vorono¨ı analysis for clustering of particles and bubbles in turbulence.
J. Fluid Mech. , 693:201–215, 2012. doi:10.1017/jfm.2011.510.A. Dejoan and R. Monchaux. Preferential concentration and settling of heavy particles in homo-geneous turbulence.
Phys. Fluids , 25(1):013301, 2013. doi:10.1063/1.4774339.A.G. Kidanemariam, C. Chan-Braun, T. Doychev, and M. Uhlmann. DNS of horizontal openchannel flow with finite-size, heavy particles at low solid volume fraction.
New J. Phys. , 15(2):025031, 2013. doi:10.1088/1367-2630/15/2/025031.M. Uhlmann and T. Doychev. Sedimentation of a dilute suspension of rigid spheres at intermediateGalileo numbers: the effect of clustering upon the particle motion.
J. Fluid Mech. , 752:310–348,2014. doi:10.1017/jfm.2014.330.S. Sumbekova, A. Cartellier, A. Aliseda, and M. Bourgoin. Preferential concentration of inertialsub-Kolmogorov particles. The roles of mass loading of particles, Stokes and Reynolds numbers.
Phys. Rev. Fluids , 2(2):024302, 2016. doi:10.1103/PhysRevFluids.2.024302.M. Uhlmann and A. Chouippe. Clustering and preferential concentration of finite-size par-ticles in forced homogeneous-isotropic turbulence.
J. Fluid Mech. , 812:991–1023, 2017.doi:10.1017/jfm.2016.826.R. Monchaux and A. Dejoan. Settling velocity and preferential concentration of heavy particlesunder two-way coupling effects in homogeneous turbulence.
Phys. Rev. Fluids , 2:104302, Oct2017. doi:10.1103/PhysRevFluids.2.104302.A. Chouippe and M. Uhlmann. On the influence of forced homogeneous-isotropic turbulence onthe settling and clustering of finite-size particles.
Acta Mechanica , 2018. doi:10.1007/s00707-018-2271-7.M. Tanemura. Statistical distributions of Poisson Vorono¨ı cells in two and three dimensions.
Forma , 18(4):221–247, 2003. 9. Uhlmann, Vorono¨ı & finite-size spheres Version: April 22, 2020J.-S. Ferenc and Z. Neda. On the size distribution of Poisson Voronoi cells.
Physica A , 385:518–526, 2007. doi:10.1016/j.physa.2007.07.063.L. Oger, A. Gervois, J. P. Troadec, and N. Rivier. Vorono¨ı tessellation of packingsof spheres: Topological correlation and statistics.
Phil. Mag. B , 74(2):177–197, 1996.doi:10.1080/01418639608240335.R. Monchaux. Measuring concentration with Vorono¨ı diagrams: the study of possible biases.
NewJ. Phys. , 14(9):095013, sep 2012. doi:10.1088/1367-2630/14/9/095013.A.B. Owen.
Monte Carlo theory, methods and examples . 2013. URL https://statweb.stanford.edu/~owen/mc/ .B. Widom. Random sequential addition of hard spheres to a volume.
J. Chem. Phys. , 44(10):3888–3894, 1966. doi:10.1063/1.1726548.J. Talbot, P. Schaaf, and G. Tarjus. Random sequential addition of hard spheres.