Pratical determination of the coordination number in granular media

In a random stacking of particles, the coordination number is defined as the mean number of particles touching a given one. A classical measurement by image analysis is directly based on this definition. It can be applied in the theoretical case of a global analysis but some problems occur for practical applications. We propose here a new measurement method based on the contribution to the Euler-Poincaré characteristic. A single measurement is performed both on the initial structure (particles linked by contacts) and on the segmented structure (isolated particles). The proposed method is robust w.r.t. the shape of the particles and their contacts.


INTRODUCTION
The concept of coordination number arises in various contexts, sometimes along with different acceptations: -in a monoatomic structure, the coordination number is nothing but the number of neighbours of each atom (cf.Fig. 1a) -in a random stacking of particles, different particles may have different numbers of neighbours; in this case, it is more appropriate to call «coordination number» the average number of neighbours per particle.
The coordination number plays an important role in the description of granular media.It can be related to the microstructural evolution of a powder during sintering (cf.Fig. 1b) and to several physical or mechanical properties of sintered materials (Artz, 1982;Jernot, 1991;Tancret et al., 1997).It has been estimated in the past by using several methods (German, 1989;Guyon and Troadec, 1994).We will focus here on its direct measurement by image analysis without any assumption about the particle shape.

GLOBAL MEASUREMENTS
The coordination number cannot be obtained using stereological relationships because it is a topological parameter.As such, it must be measured directly in the space of interest (e.g., mostly 3D space in materials science).The coordination number is defined as the mean number of particles touching a given one.

Sequential analysis
This method directly derives from the definition of neighbours.It starts with a segmentation by watershed of the population of particles (Lantuéjoul and Beucher, 1981).Then, a sequential analysis is performed on this segmented set.At first, one particle is selected and reconstructed by geodesic dilation (Fig. 2a).Then, after dilation, this particle is intersected with the segmented set and its number of neighbours is counted (Fig. 2b).This procedure is repeated for each particle up to a complete examination of the initial set (note that it is time-consuming).The coordination number is the mean number of neighbours measured for all the particles (Fig. 2c).

Parallel analysis
This method is based on the measurement of the Euler-Poincaré characteristic (Hadwiger, 1957;Serra, 1982;Stoyan et al., 1995;Ohser and Nagel, 1996) further denoted by EPC.In 2D space, this topological parameter is simply the number of connected components minus the number of holes within them.In 3D space, it is equal to the number of distinct surfaces minus their genus.
Consider a finite population of particles {P i , i∈I} satisfying the following two conditions: a) each particle is convex; b) pairs of particles can touch each other but triplets cannot.
Assume that one can construct (for instance by erosion or segmentation of the population of particles) another population of particles {Q i , i∈I} such that: c) each Q i is non empty, convex and contained within P i ; therefore the Q i 's can be referred to as subparticles; d) distinct subparticles are pairwise disjoint.
From both populations, one defines the sets . The EPC is denoted by N.
Let N C (P i ) be the coordination number of particle P i : On the other hand, because of d), we have This formula is illustrated on Fig. 3 with the two measurements on the initial connected set X and on the disconnected set Y.

LOCAL MEASUREMENTS
The two previous methods provide exact results only when the structure is globally explored.They must be adapted when the structure has to be analysed locally, which arises when the set of interest cannot be examined at one go.

Sequential analysis
In order to be applicable, this method requires the neighbourhood of each particle to be totally contained within the measurement field.Equivalently, any particle that has its neighbourhood hitting the edges of the field cannot be properly measured and therefore must be discarded (cf.Figs. 4 and 5).Moreover, the biggest particles have the largest propensity of hitting the edges of the frame.The «bias» incurred can be compensated by the Miles-Lantuéjoul correction (Miles, 1974;Lantuéjoul, 1980).It consists of assigning each particle a weight that is proportional to the chance it has to be contained within the measurement field.Parallel analysis As will be seen in the following sections, the parallel analysis can be easily adapted to the local case.

Local contribution to the Euler-Poincaré characteristic
In a recent paper (Jernot et al., 2004), it has been established that the EPC of a bounded set X of R d can be obtained by superimposing arbitrarily a tessellation on X (i.e., a partition of R d into convex polyhedral cells), measuring the contribution of each cell to the EPC of X and summing all contributions.
The local contribution of a cell Z to the EPC of X is a weighted sum of the EPC of X within each facet F § of the tessellation contained in Z.The weights are chosen to compensate for the fact that facets belong to several cells.It is explicitly defined by the formula: ( ) ( ) in which dimF denotes the dimension of F (0, 1, 2, ..., d) and ordF stands for its multiplicity order, i.e., the number of cells sharing F as a facet.In the case of a 2D square tessellation of the space the orders of the facets are 4 (vertices), 2 (edges) and 1 (face).The measurement is illustrated in Fig. 6 on a 2D example: in the left field of this figure, the EPC values 1, 2 and 2 come respectively from the non-empty intersections of X with one vertex, two edges and one face.

Local measurement of the coordination number for a finite family of particles
In terms of local contributions, since , the coordination number can be expressed as: ( ) or, equivalently, as a ratio of expectations obtained by selecting one cell at random among all the cells hitting X: The set presented in Fig. 4 can be analysed using this method: the contributions are equal in both measurement fields, C(X;Z) = 1/2 and C(Y;Z) = 3/2, which gives 3 4 N C = by applying Eq. 2. In Fig. 7, the population of particles presented in Figs. 2 and 3 is analysed through four adjacent fields Z 1 , ... , Z 4 .Because this population is bounded, both local and global analyses lead exactly to the same result:

Local measurement of the coordination number for an infinite family of particles
At this stage, only finite populations of deterministic particles have been considered.Now we may wonder whether the definition of the coordination number can be extended to random, stationary and locally finite populations of particles.
Let X & (resp.Y & ) be the random set made up of the union of all particles (resp.all subparticles).To define the coordination number, a natural idea is to replace in Eq. 1 the connectivity numbers by their corresponding specific values.Namely we put , whatever the convex domain D with nonempty interior (Stoyan et al., 1995).
In practice, Eq. 4 is not directly applicable because both specific connectivity numbers are unknown.Various estimators are available for -if the domain consists of two parallel section planes separated by a known distance, the disector technique can be used (Sterio, 1984).
-if the domain is a cell of a regular tessellation (i.e., all cells are identical up to a shift to a reference cell Z) the contribution technique can also be effective.It yields the intuitive but nontrivial fact that the proof of which is left in a forthcoming paper.
Starting from (5), unbiased estimators for ( ) can be obtained, either by averaging the contributions of different cells or by assigning them different weights provided by a Kriging exercise (Chilès and Delfiner, 1999), which makes it possible to integrate the spatial structure of the random set under study.
It should be pointed out that the knowledge of one realization X of X & in two parallel section planes is insufficient to properly determine the corresponding realization Y of Y & .A close look at Fig. 9 will easily convince the reader that the determination of Y requires X to be known in a full-dimensional domain.For the sake of consistency, it is deemed preferable to estimate both connectivity numbers also in a fulldimensional domain, thus resorting to the contribution technique rather than to the disector technique.
From Eqs. 4 and 5, it comes This is basically the same formula as Eq. 3, except that the expectations have different interpre-tations.They are taken over all cells hitting X in Eq.
3 and over all realizations of X & in Eq. 6.
As a consequence, suppose that the contributions of the realizations X and Y have been measured on the cells J) j , Z ( j ∈ , then we propose to use ( ) as an estimate of the coordination number.In this formula, both sets of coefficients J) j , ( j ∈ λ and J) j , ( j ∈ µ are either constant weights (equal to 1/Card J) or Kriging weights.They add up to one.

A PRACTICAL EXAMPLE
The population of particles displayed in Fig. 8 is a simulation of a disc packing (Bideau, 1983).It has been chosen because its coordination number is explicitly known to be equal to 4 as a consequence of the underlying physical process (gravity) that governs its construction.On this population of particles both analyses, sequential and parallel, were performed.As can be seen in Fig. 8, both results (4.095 and 4.109) are quite similar even if 28 traces of particles have been discarded during the sequential analysis.
It should be pointed out that sequential analysis gives valuable results only when many particles have their whole neighbourhood accessible.For example, examining the population of particles through the fields of Fig. 9 cannot yield any relevant result.In contrast to this, the parallel analysis can cope with this limitation.On each field the contributions for X and Y were determined.The sums of all contributions, gives an estimated value of the coordination number equal to 4.109.On the other hand, it is also interesting to estimate the coordination number on each field (see Fig. 9).The average of these estimates is equal to 4.113 which is not very far from the theoretical value 4.This probably comes from the fact that all the particles have the same size and the distribution of the number of neighbours per particle is tight.This 2D example may look somewhat academic but it stresses the limitations of the sequential analysis.It turns out that these limitations are even more severe in 3D space.At the outset, there is a problem in the choice of the resolution that creates a conflict between the size of the particles and that of the measurement field.A high resolution leads to a good definition of contacts but to a small number of particles, whereas a low resolution allows a large number of particles to be analysed but the contacts are poorly defined.

CONCLUSION
In order to estimate the coordination number of particles in a powder stacking, a new method, based on the Euler-Poincaré characteristic, has been proposed.This method is independent from the workspace dimension.It is easy to implement and fast.It has been extended to the case of a local analysis using the local contribution to the Euler-Poincaré characteristic.This method is robust w.r.t. the nature of contacts (punctual or not) and the particle shape.
Fig. 1.Contacts in granular media: a) punctual contacts in a periodical arrangement of particles, b) necks between sintered bronze particles.

Fig. 4 .
Fig. 4. Local measurement by sequential analysis: the global value of the coordination number (4/3) is out of reach with this method.

Fig. 5 .
Fig. 5. Local measurements on the previous set of Fig. 2. The coordination number calculated by summing up the results from each measurement field is 91 . 2 11 32 N C ≈ = .
numbers of X & and Y & .These are defined as domain where the realizations of X & and Y & are known:

CFig. 8 .
Fig. 8. Local estimation of the coordination number for an infinite structure: comparison between sequential and parallel analyses.
Touching particles are called neighbours.In what follows, two global measurement methods are presented.