Image Anal Stereol 2007;26:121-127 Original Research Paper PERCOLATION OF RANDOM CYLINDER AGGREGATES Dominique Jeulin1 and Maxime Moreaud2 1Ecole des Mines de Paris, Centre de Morphologie Mathématique, 35, rue Saint-Honoré, 77305 Fontainebleau Cedex, France; 2NT2I, Pôle Optique et Vision, 20 rue Benoit Lauras, 42000 Saint Etienne, France e-mail: dominique.jeulin@ensmp.fr; maxime.moreaud@mines-paris.org (Accepted November 11, 2007) ABSTRACT The percolation threshold ?c of Boolean models of cylinders with their axis parallel to a given direction is studied by means of simulations. An efficient method of construction of percolating connected components was developed, and is applied to one or two scales Boolean model, in order to simulate the presence of aggregates. The invariance of the percolation threshold with respect to affine transformations in the common direction of the axis of cylinders is approximately satisfied on simulations. The prediction of the model (?c close to 0.16) is consistent with experimental measurements on plasma spray coatings, which motivated this study. Keywords: Boolean model, cylinders, multi-scale, percolation, plasma spray coating, simulation. INTRODUCTION METHODS In many applications in material sciences, thin layers are deposited by means of a plasma spray (Beauvais et al., 2004; Kang, 2005; Sharma et al., 2006). This coating is used to give specific properties to a part, like hardness, thermal or electrical conductivity. When the layer is made of a mixture of various components with a high contrast of properties (like a binary mixture with an insulating and a conductive phase), its physical properties are connected to the possibility for one of the components to percolate through or along the layer. Therefore it is important to know the percolation threshold of the layer in various directions, as a function of the microstructure. In the present work, the toroidal geometry of droplets generating the layer (Beauvais et al., 2004) is simplified by considering cylinders with their axis orthognal to the layer. A fast and efficient numerical method of simulation to generate multiscale 3D complex aggregates is presented, following the approach developed for sphere aggregates (Jeulin and Moreaud, 2005) and for sphero-cylinders (Jeulin and Moreaud, 2006a,b). This method enables us to estimate the percolation threshold of the random microstructure. It is applied to multi-scale random models of cylinders, for which the percolation thresholds of the cylinders and of their complementary sets are estimated. NUMERICAL SIMULATION OF MULTI-SCALE RANDOM AGGREGATES Plasma spray coatings are generated by droplets falling on a substrate. A good representation of the layer is given by a realization of a Boolean model involving random primary grains located on points of a Poisson point process (Matheron, 1967). The initial spherical droplets are spread out on the substrate, and their shape can be approximated by flat cylinders with their axis orthogonal to the layer. In what follows we consider the case of a single primary grain (cylinder with diameter D, height H, and shape factor D/H). The germs are located on the points of a Poisson point process with a constant intensity, or restricted to a random set (multi-scale model) ; in the second case, we obtain a Boolean model based on a Cox point process (Jeulin, 1996). Our method of simulation can be applied to any kind of aggregates made of objects with a bounded diameter (that can be embedded inside a cube with side T, as shown on Fig. 1), with a random implantation on a multi-scale Poisson process (Jeulin, 2004; Moreaud, 2006). 121 Jeulin D et al: Percolation of cylinders number of tests of intersection between objects. This is made by a subdivision of the complete volume into cubes of size T, corresponding to the maximal diameter of the objects (for parallel cylinders, cubes could be replaced by englobing parallelepipeds for a faster simulation). The cubes are explored in a given order (Fig. 2). In each cube, the random number of centres of objects to locate is simulated by a Poisson point process with intensity ?. For each new centre of object to locate is generated: • a uniform location of the centre of the object in the current cube. • a determination of the 13 adjacent cubes of the current cube (according to the exploration order) (Fig. 2). • For all objects already contained in these volumes is performed an intersection test, detailed later. The next step (creation and fusion of aggregates) is explained in the “fusion of aggregates”section. SIMULATION OF A TWO SCALES MODEL We now consider a two-scale model built according to a Cox point process (Jeulin, 1996), where the centres of objects O1 with maximal diameter T1 are randomly located outside exclusion objects O2 with maximal diameter T2. Objects of the small scale O1 are stored only if their centres are outside large size objects O2 (Fig. 3, left). Fig. 2. Left, division of a volume into eight parts ordered along x, y and z. Middle: the 13 adjacent volumes (light grey) to the current volume (dark grey). Right: the 26 adjacent volumes (bright grey) to the volume V2 (dark grey). 122 Fig. 1. Left, cylinder (T = max(diameter, height)); right, any shape (T = maximum caliper diameter of the shape). SIMULATION OF A ONE-SCALE MODEL Let us consider an implantation of centres of objects with a maximal diameter T, according to a Poisson point process with intensity ?. Every volume V contains a random number of objects N following a Poisson distribution with mean ?V. Two properties of the Poisson point process are used: the number of points N1 and N2 contained in disjoint volumes V1 and V2 are independent variables. Conditionally to Ni , the points falling in Vi are uniformly distributed in Vi. The aggregates are obtained by overlapping objects. For the simulation, the following information concerning each object is stored: coordinates of the centre, shape parameters (like heigth and radius for a cylinder), and an aggregate label. To get a fast method of construction, it is useful to limit the Image Anal Stereol 2007;26:121-127 Fig. 3. Examples of two scales simulations. Left: the cylinders corresponding to the small scale objects are kept if their centers are outside of exclusion cylinders corresponding to large size objects. Right: the cylinders, corresponding to the small scale objects, are kept if their centers are inside inclusion cylinders corresponding to large size objects. The simulation is built in two steps: • In a first step, a one-scale simulation is performed for the large exclusion objects O2, without intersection test. To generate this simulation, the full volume is subdivided into cubes V2 with side T2. • In a second step, a one-scale simulation is generated for the small size objects O1, with addition of a new condition. For every new object O1i to be implanted, the cube V2 containing its center and its 26 adjacent cubes are extracted (Fig. 2). Then is performed a test of exclusion of the centre of O1i into each large size object O2 included in these volumes. FUSION OF AGGREGATES The fusion of aggregates is a critical point concerning the computation time processing. To get a fast fusion process, the aggregate labels of objects are updated at the end of the simulation by means of a labeling algorithm (Moreaud, 2006). This method avoids multiple browses of the objects during the simulation, which can take a very long time when several millions of objects are implanted. An intermediate table “AggregateList” is used to link aggregate labels of objects during the simulation with aggregate labels of objects at the end of the simulation. Each object Oi owns an aggregate label Oi.ag. After intersection tests between a new object On to locate and previously located objects Oi, we have to consider several cases: • If On does not intersect objects, a new aggregate is created: On.ag = NewLabel and AggregateList[On.ag] = NewLabel. • If On intersects an object Oi, and if AggregateList[Oi.ag] ^ AggregateList[On.ag], then the two aggregates have to be merged. We update Oi.ag = On.ag = AggregateList[Oiag] = AggregateList[On.ag] = min (AggregateList[Oiag], AggregateList[On.ag]). When the full volume is simulated, we reflect the updates of aggregate labels only in the AggregateList, and then we update aggregate labels of the objects. METHOD OF ESTIMATION OF THE PERCOLATION THRESHOLD A simulation percolates when an aggregate of objects connects two opposite faces of the simulation. To test if a simulation percolates, the highest label N of objects belonging to the cubes of the first plane z is determined, and is compared to the label numbers of objects belonging to the cubes of the last plane z. If all the labels are larger than N, the simulation cannot percolate. If we only look for the percolation threshold, this test, performed for each plane along z, enables us to quickly stop the simulation when it does not percolate. With our method of construction, several realizations of 3D simulations are processed and the presence of at least an aggregate connecting two opposite faces of a field is tested. The percolation threshold pc is obtained when 50% of the realizations percolate for a given volume fraction of objects, and is given by this volume fraction. For two-scale simulations, it is interesting to estimate the percolation threshold of the smaller scale objects, keeping constant the volume fraction of the other objects (exclusion in the present case). A research by dichotomy is used to estimate pc. TIME PROCESSING The time processing for a one scale simulation is depending of the number of tests of intersections to process. With a rough algorithm, the number S1 of tests of intersections is given by: ? n(n -1) 1 2 51 = i =------?— n with n the number of tt 2 2 objects. By using our method which divides the volume of simulation, the number S2 of tests of intersections can be approximated by: T3 T3 2 52 ?n 14 n 3?14 3n with n the number of v v objects, T the diameter of the objects and V the volume of the simulation. For instance, if we process a simulation on a volume V = v3 = 20003 with 123 Jeulin D et al: Percolation of cylinders 5000000 objects with diameter T = 10, we obtain 51 = 12.5×1012 tests of intersections of objects and 52 = 43.75×106. SIMULATION OF RANDOM AGGREGATES OF CYLINDERS AND PERCOLATION THRESHOLD For the study of coatings, we use our method with cylinders with a fixed orientation. By means of large size 3D simulations (typically 5000000 objects by simulation, and 20 realisations by shape factor), the clusters connected to the boundaries of the field are easily extracted thanks to the labelisation of the aggregates, in order to accurately estimate the percolation threshold (Moreaud, 2006). SIMULATIONS WITH CYLINDERS A Boolean model of cylinders (BMC) is used for a one scale structure (cf. Fig. 4a). We also use a two scales model (cf. Fig. 4b): • an exclusion scale modelled by a Boolean model of large cylinders, corresponding to a domain containing no small size cylinders. • a scale of aggregates made of small size cylinders. The two scales model can be used to simulate random aggregates, or to estimate the percolation threshold of the complementary set of the BMC. To implement our method, we have to detail the tests of intersection and of exclusion between two cylinders C1 and C2 (with radii R1 and R2, centers M1(x1,y1,z1) and M2(x2,y2,z2), heights h1 and h2: • Test of intersection: the two cylinders overlap if the two following inequalities are satisfied: (x1 - x2)2 + (y1 - y2)2 < (R1 + R2)2 and 4(z1 - z2)2 < (h1 + h2)2 • Test of exclusion: the center of the cylinder C1 is outside of the cylinder C2 if the following condition is fulfilled: (x1 - x2)2 + (y1 - y2)2 >R22 or 4(z1 - z2)2 > h22 ESTIMATION OF THE PERCOLATION THRESHOLD OF THE BOOLEAN MODEL OF PARALLEL CYLINDERS The percolation threshold of the Boolean model with various grains (spheres, sphero-cylinders...) can be estimated by different methods (Rintoul and Torquato, 1987; Mecke and Stoyan, 2005; Jeulin and Moreaud, 2005; 2006a,b). We used our method of simulation for a Boolean model of parallel cylinders with the same size (a single scale) for different shape factors (ratio diameter/height), and using 20 realizations for each case. The structure being anisotropic, we considered the percolation threshold for propagation along the axis of cylinders (axis 1) and orthogonally to their axis (axis 2) (cf. Fig. 5). The simulations are periodic on the faces of the domain parallel to the direction of propagation. The results are given in table 1. The percolation threshold is slightly larger along axis 1 (close to 0.17), as compared to a propagation along axis 2 (close to 0.15). It is practically invariant with the shape factor, as could be expected: as we already indicated in (Charollais et al., 1997), the connectivity properties of the Boolean model are invariant by affine transformations, which changes the shapes and the sizes, but does not introduce any new connection (this property was used to estimate the 3D connectivity number of a Boolean model of ellipsoids). In the case of parallel cylinders, changing the shape factor is equivalent to performing an affinity of the microstructure in the direction parallel to the axis of cylinders. During this transformation, we keep a Boolean model (the Poisson point process of germs being still a Poisson process) and its connectivity, and therefore its percolation threshold. It is possible that for very large fields of simulations the percolation thresholds in the two directions would converge to the same values, as observed for sphero-cylinders parallel to the same plane in (Jeulin and Moreaud, 2006a), but this point should be checked with very large simulations. The percolation threshold of parallel cylinders is smaller than for spheres (close to 0.29) and higher than for an isotropic distribution of sphero-cylinders (between 0.01 and 0.0003 for aspect ratios h/r in the range 10-3000). It is very close to the percolation threshold of sphero cylinders with their axis parallel to the same plane (up to ± 5 degrees of disorientation): in (Jeulin and Moreaud, 2006a), for 2r/h = 0.2 the percolation thresholds along axis 1 and 2 were 0.1294 and 0.16125 respectively. More interestingly, recent experimental data on the in-plane conductivity of air plasma sprayed molybdenum coatings allowed an empirical estimation of the percolation threshold close to 0.16 (Sharma et al., 2006), in an excellent agreement with our predictions with simulations. Other estimates of the percolation threshold are available: the excluded volume (Balberg et al., 1984), and the zeros of the connectivity number (Bretheau and Jeulin, 1989; Mecke and Stoyan, 2005), being only valid for isotropic Boolean models. 124 Image Anal Stereol 2007;26:121-127 (a) (b) Fig. 4. (a) BMC and the complementary set of the BMC (one scale); (b) two scale Boolean model;an exclusion scale is used to simulate the complementary set of the BMC. Tab. 1. Percolation threshold for the Boolean model of parallel cylinders for different shape factors. Axis 1 is parallel to the axis of cylinders. Axis 2 is orthogonal to axis 1. Shape Factor (2r/h) vol. size ?c axis 1 ?c axis 2 100 4000 0.1763 0.1540 50 2000 0.1739 0.1522 10 400 0.1672 0.1497 2 200 0.1544 0.1537 1 200 0.1524 0.1543 1 100 0.1518 0.1562 1 40 0.1508 0.1531 0.1 40 0.1250 0.1596 ESTIMATION OF THE PERCOLATION THRESHOLD OF THE COMPLEMENTARY SET OF A BOOLEAN MODEL OF PARALLEL CYLINDERS Our method of simulation of multi-scale random aggregates can be used to estimate the percolation threshold of the complementary set of a Boolean model of parallel cylinders (cf. Fig. 4). In that case, we need to use two scales simulations (cylinders for aggregates and for exclusion zones) and to use the method of estimation of the threshold explained in a previous section with a slight change: we estimate by dichotomy the limit volume fraction of exclusion cylinders allowing a percolation of small size aggregate cylinders for 50% of realizations (the small size cylinders are used here to find paths in the complementary set of large cylinders). We get by this way an estimation of the percolation threshold of the complementary set of the Boolean model of parallel cylinders. The results are given in Table 2. The cylinders for aggregates have a diameter and a height equal to 0.1, the cylinders for exclusion zones have a diameter equals to 1 or 2 and a height equal to 1. We can notice that the percolation threshold along axis 1 is higher than for the grains of the Boolean model (compare to Tab. 1), while it is much lower for axis 2: it is much easier to find paths out of cylinders along their faces, than along their axis. To our knowledge, such results are not available in the literature. 125 Jeulin D et al: Percolation of cylinders Fig. 5. Illustration of the two axis for the estimation of the percolation thresholds. Tab. 2. Percolation threshold for the complementary set of a Boolean model of parallel cylinders for different shape factors. Axis 1 is parallel to the axis of cylinders. Axis 2 is orthogonal to axis 1. Shape Factor vol. size ?c axis 1 ?c axis 2 (2r/h) 1 10 0.2167 0.0406 2 20 0.2320 0.0177 ESTIMATION OF THE PERCOLATION THRESHOLD OF A TWO SCALE BOOLEAN MODEL OF PARALLEL CYLINDERS Two scales simulations were performed as described earlier, using large cylinders for inclusion zones, and small size cylinders for the generation of aggregates (cf. Fig. 3, right). The results of simulations are given in Tab.3. As expected, a much lower threshold is obtained, as compared to the one-scale model. The two directions give close results. They should be compared to the lower bound of the percolation threshold, asymptotically expected for a very large ratio between the two scales ?c(2) = (?c (1))2. In the present case, we have ?c(2) = 0.0256. It turns out that the lower bound is not reached since the ratio of scales is limited to 10 in our simulations. Tab. 3. Percolation threshold for two scales simulations. Axis 1 is parallel to the axis of cylinders. Axis 2 is orthogonal to axis 1. cyl. agr cyl. inc vol. ?c axis 1 ?c axis 2 (2r - h) (2r - h) size 1 1 10 10 400 0.0367 0.0382 2 1 20 10 400 0.0377 0.0362 Fig. 6. Example of a two scales simulation (volume: 3003, cylinders for aggregates (r = 2, h=1, Vv=0.2), cylinders for inclusion zones (r = 20, h = 10, Vv = 0.1)). CONCLUSION We presented an original method of simulation of multi scale random aggregates of parallel cylinders based on the Boolean model. This kind of model is well suited to the simulation of coatings deposited by means of a plasma spray. By subdividing the studied volume and using the properties of the Poisson point process, the computation time is significantly reduced to obtain aggregates connected to a boundary of the simulated domain, from which is numerically estimated the percolation threshold. We could check from simulations that the percolation threshold is practically invariant by affinity along an axis with the same orientation as the axis of cylinders, as expected from the invariance of the connectivity by affinity. Furthermore, the predicted percolation threshold is in agreement with experimental data obtained on plasma spray coatings. Using two-scale depositions would significantly lower the percolation threshold of the layer. ACKNOWLEDGEMENT The authors are indebted to V. Guipont (Centre des Matériaux Mines Paris, Paristech) for stimulating discussions. REFERENCES Balberg I, Anderson CH, Alexander S, Wagner N (1984). Excluded volume and its relation to the onset of percolation. Physical review B 30(7):3933-43. Beauvais S, Guipont V, N'Guyen F, Jeandin M, Jeulin D, Robisson A, et al (2004). Study of the porosity in plasma-sprayed alumina through an innovative 3-dimensional simulation of the coating build-up. In: Thermal Spray Solutions – Advanced in Technology and Aplications. Proceedings, International Thermal Spray Conference and Exposition (ITSC 2004), Osaka, Japan, 10-12 May 2004, Pub. ASM International 2004, 782:9. 126 Image Anal Stereol 2007;26:121-127 Bretheau T, Jeulin D (1989). Caractéristiques morphologiques des constituants et comportement a la limite élastique d'un matériau biphasé Fe/Ag. Revue Phys Appl 24: 861-9. Charollais F, Bauer M, Coster M, Jeulin D, Trotabas M (1997). Modelling the structure of a nuclear ceramic obtained by solid phase sintering, Acta Stereol 16(3): 315-21. Jeulin D (2004). N-06-04-MM. Paris Schol of Mines Publication, March 2004. Jeulin D, Moreaud M (2005). Multi-scale simulation of random spheres aggregates – application to nano-composites. In: Chraponski J, Cwajna J, Wojnar L, eds. 9th European Congress on Stereology and Image Analysis 2005 May 10-13; Zakopane, Poland, Vol I, 341-8. Jeulin D, Moreaud M (2006a). Percolation of Multi-Scale Fiber Aggregates, S4G (Stereology, Spatial Statistics and Stochastic Geometry). In: Lechnerova R, Saxl I, Benes V, eds. 6th International Conference, Prague. Jeulin D, Moreaud M (2006b). Percolation d’agrégats multi-échelles de spheres et de fibres: Application aux nanocomposites; Matériaux, Dijon, France, Nov 2006, 13-7. Jeulin D (1996). Modeling heterogenous materials by random structures. Invited Lecture, European Workshop on Application of Statistics and Probabilities in Wood Mechanics, Bordeaux, 22-23 March 1996, N-06/96/ MM, Internal report, Ecole des Mines de Paris. Kang H K (2005). Microstructure and electrical conductivity of high volume Al2O3 reinforced copper matrix composites produced by plasma spray. Surface and Coating Technology 190:448-52. Matheron G (1967). Eléments pour une théorie des milieux poreux. Paris: Masson. Mecke K, Stoyan D (2005). The Boolean Model: from Matheron till today In: Bilodeau M, Meyer F, Schmitt M, eds. Space, Structures, and Randomness, Contributions in Honor of Georges Matheron in the Fields of Geostatistics, Random Sets, and Mathematical Morphology, Lecture Note in Statistics 183, Springer, 151-82. Moreaud M (2006). Propriétés morphologique multiéche-lles et prévisions du comportement diélectrique de nanocomposites. These, Ecole des Mines de Paris. Rintoul M D, Torquato S (1987). Precise determination of the critical threshold and exponents in a three-dimensional continuum percolation model. J Phys A Math Gen 30: L585-92. Sharma A, Gambino R J, Sampath S (2006). Electrical conduction in thermally sprayed thin metallic coatings. Mater Res Symp Proc 890:287-92. 127