- Email: [email protected]
Microstructure characterization of granular materials
ARTICLE IN PRESS
Physica A 377 (2007) 545–558 www.elsevier.com/locate/physa
Microstructure characterization of granular materials Riyadh Al-Raoush Department of Civil and Environmental Engineering, Southern University and A and M College, Baton Rouge, LA 70813, USA Received 12 June 2006; received in revised form 2 November 2006 Available online 3 January 2007
Abstract This paper presents techniques and algorithms to compute microstructure properties of irregular-shaped granulate assemblies utilizing 3D images. The techniques are capable of extracting microstructure properties of particles such as centeroid, particle size distribution, shape indices (i.e., sphericiy and angularity), number of contacts (i.e., distribution of coordination numbers), contact network, packing efficiency, distribution of local void ratio and radial distribution function. Such properties are critical parameters for micromechanical-based numerical models to capture micro- and macromechanical behavior of geomaterials. X-ray microtomography was used to reconstruct high-resolution 3D image of a natural sand system to represent granular materials. Microstructure properties of the sand system were computed and compared with properties of a computer-simulated image of periodic random spheres. Findings indicate that the use of simplified systems of idealized spheres to model micro- and macromechanical behavior of granular systems can lead to inaccurate results due to the differences in the microstructure between both systems. Methods presented in this paper enabled capturing a more realistic microstructure that can be incorporated in micromechanical models to better simulate, understand, or explain macroscale behavior of granular materials based on their actual microstructure. r 2007 Published by Elsevier B.V. Keywords: Microstructure; Granular materials; Micromechanical properties; Microtomography; Random packing
1. Introduction Granular material, one type of frictional geomaterials, is often composed of individual grains that vary in shape, size and surface texture. Such parameters can significantly affect the final packing of granulates and their contact distribution characteristics. Microstructure characterization of irregular particles such as granular materials, pharmaceutical materials and powders, is of a significant importance to many engineering applications, particularly to models that implement micromechanical theories to describe the macroscale behavior [1,2]. It has been recognized that geometry and topology of a granular system control its electrical, thermal, fluid-like flow, stress transfer mechanisms and other mechanical properties. For instance, mechanical conditions of a granular system such as shearing, compaction, transmission of stress and statistical properties of contact force distribution are greatly influenced by its microstructure characteristics (e.g., the distribution of the number of contacts of particles) [3–5]. Considering two adjacent particles within the assembly, one needs Tel.: +1 225 771 5870; fax: +1 225 771 4320.
E-mail address: [email protected]. 0378-4371/$ - see front matter r 2007 Published by Elsevier B.V. doi:10.1016/j.physa.2006.11.090
ARTICLE IN PRESS 546
R. Al-Raoush / Physica A 377 (2007) 545–558
to account for the exact interaction characteristics so that the interparticle forces and relative motions can be estimated. The interparticle interaction is essentially controlled by size and shape of the contact area. Furthermore, the macro- or global behavior of granulate assembly such as fluid-like flow, stress–strain response and multi-phase flows can be described by microbased models. Such models require information about microstructure properties such as local void ratio distribution, contact number distributions and particle size distribution. These properties require robust techniques for characterization. Development of such techniques has been a real challenging issue for many researchers, moreover, there are complications associated with the implementation of the actual microstructural properties in numerical models. To overcome the complexity associated with real granular systems, granular materials are often represented and modeled as ideal packing assemblies of rigid spherical particles, elliptical particles and/or clusters of spheres [1,6]. These models have been found to serve as a good starting point to account, to some extent, for the particle shape effects on the global behavior of granulates. However, fundamental problems and limitations of their applicability to the real granular materials are obvious. Therefore, the actual spatial distributions of the microstructure properties are not quite incorporated and thus, fundamental physics can be missing. While a few studies have been conducted to characterize microstructure properties of irregular systems, these studies have been limited to determination of basic microstructure properties such as local void ratio, surface area, and density [7,8]. The objectives of this paper are: (1) to present algorithms to compute microstructure properties of granular materials from 3D images, and (2) to compare microstructure properties of a computersimulated pack of random spheres to those of natural systems (i.e., sand systems). Synchrotron X-ray microtomography is used to obtain the 3D image of the natural sand system. 2. Materials In this paper, I present two different systems: computer-simulated 3D image of random spheres and X-ray microtomography image of a sand system. The first system is periodic random packing consisting of 1000 spheres uniformly sized and digitized at 25 voxels/sphere diameter. The size of the 3D image is 200 200 200 voxels and has a porosity value of 0.354. The packing was generated by the algorithm developed by Liu and Thompson [9] based on a collective rearrangement method. In these type of models, sphere packs are generated through two steps; random generation and placement of spheres in a given domain allowing them to overlap. Spheres are usually generated from a given probability density function. An iterative relaxation process, such as Monte Carlo or molecular dynamics simulations, is then applied to separate the overlapped spheres to reach a stable position with a specific overlap tolerance. Sphere sizes are reduced if the overlap between spheres is larger than the specified tolerance. The relaxation process is repeated until a stable packing with acceptable sphere overlap is obtained. More details about the algorithm used to generate the sphere packing in this paper can be found in Ref. [9]. The second system is a natural sand system [10,11]. The system was imaged by synchrotron X-ray tomography at the GeoSoilEnviroCARS (GSECARS) BM-13D beamline at the Advanced Photon Source, Argonne National Laboratory. Image resolution is 11.5 mm. Image reconstruction was performed using algorithms developed by GSECARS to convert CT attenuation data to 3D volumetric data. The size of the image is 375 400 340 voxels. The porosity of the system is 0.363 and the mean particle size, D50, is 0.325 mm. Note that porosity values of both systems are similar; this in turn provides reasonable base of comparison between microstructure properties of natural and computer-simulated systems of similar porosity values. 3. Recognition of individual particles Recognition of individual particles in a given system is the first step required to determine its microstructure properties. Particle identification of 3D irregular particles such as granular materials (e.g., sand) is a challenging task. In this paper, I modified the algorithm developed by Al-Raoush that labels individual particles and improved its accuracy and speed [12]. The algorithm is based on watershed transform and proceeds by constructing a distance map of the voxels of the particles in the image. Fig. 1 illustrates the procedure to distinguish particles in contact at an area. The procedure is best explained through the steps shown in Fig. 1 as follows: Fig. 1(a) shows a binary image of 2 touched particles, the image has 2 voxel values;
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
547
Fig. 1. Separation of touched particles; (a) binary image of two touched particles, (b) distance transform of the two particles, (c) flooding simulation process starts from local maxima to define watershed line that separates the particles, (d) watershed line (dividing line) and (e) a voxel is considered part of the dividing line if the filling enters a voxel from the indicated directions [12].
values of ‘‘1’’ depict the foreground (i.e., particles) and values of ‘‘0’’ represent the background. Fig. 1(b) shows the distance map of the particles calculated by distance transform. Each foreground voxel in the image is labeled based on its distance from the nearest background voxel. The distance transform is applied to the foreground voxels of the image starting from the boundaries of the particles towards their centeroids to obtain local maxima. Local maxima are markers that distinguish distinct (individual) particles in the image. In Fig. 1(b), local maxima are the voxels labeled by ‘‘11’’ and ‘‘6’’ in the large and small particles, respectively. Distance map of a 3D image is computed by calculating the distance X ab Euclidean between two voxels aðia ; j a ; ka Þ and bðib ; j b ; kb Þ using the Euclidean metric as 2 2 2 1=2 X ab , Euclidean ¼ ðia ib Þ þ ðj a j b Þ þ ðka k b Þ
(1)
where ia, ja and ka are the row, column, and depth indices of a given voxel (i.e., a, b, y), respectively, in the 3D image. It is worth noting that the row, column and depth indices are equivalent to the x, y and z distances in the Cartesian coordinate system, respectively. The distance map algorithm commences by labeling boundary voxels of the particles in the image by assigning label ‘‘1’’ to each voxel indicating that they are one voxel away from the closest background voxel. A boundary voxel is the one that connects a background voxel through a face, edge or corner. The algorithm searches the 26 neighboring voxels of each foreground voxel and counts their voxel values. If the summation is 26, then the selected voxel is an interior voxel (i.e., surrounded by 26 foreground voxels); otherwise the selected voxel touches at least one background voxel and is considered a boundary voxel. The algorithm proceeds until all foreground voxels are inspected and all boundary voxels are identified and labeled. The algorithm proceeds and searches foreground voxels to identify the neighboring voxels of the voxels that was labeled in the previous step; voxel distance (i.e., voxel labeling) is calculated between these voxels and the closest background voxel using Eq. (1). The algorithm proceeds until all foreground voxels are labeled indicating their voxel distance from the closest background voxel (i.e., distance map of foreground voxels). Fig. 1(c) shows the flooding simulation that starts from local maxima shown in Fig. 1(b). The algorithm simultaneously labels voxels that have the highest values in the distance map and assigns a label ‘‘1’’ to them (i.e., voxels that are labeled by ‘‘11’’ in Fig. 1(b)). The second iteration starts by simultaneously labeling voxels that have the second highest values in the distance map and assigns a label ‘‘2’’ to them (i.e., voxels that are labeled by ‘‘10’’ in Fig. 1(b)). The same labeling procedure used in the previous iteration is repeated until all voxels are labeled. In the flooding simulation, each marker represents a distinct particle (catchment basin) where the simulation moves from the markers with a constant rate of labeling. Dividing lines (watershed lines)
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
548
a
c
b
Local maxima (markers)
Catchment basins
Watershed lines
Direction of flood simulation Local maxima (markers)
d
e Watershed line
Catchment basins Fig. 2. Identification of individual particles; (a) image of two touched particles, (b) distance transform of the two particles, (c) flooding simulation process to define watershed line that separates the particles, (d) watershed line and catchment basins (particles) and (e) labeled particles.
between distinct particles are constructed at voxels where labeling comes from different catchmnet basins (i.e., particles). Fig. 1(e) illustrates the directions that decide voxels of dividing lines. Fig. 2 shows an overview of the algorithm of particle identification and Fig. 3 shows particle labeling of random spheres and sand systems. Image resolution and quality impact the results of particle identification algorithm and other algorithms presented in this paper. Image quality can be improved by applying specific types of filters based on the artifact type of the image. The impact of the resolution on a given parameter can be determined by carrying out a resolution sensitivity analysis where tolerance values can be determined.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
549
Fig. 3. Identification of individual particles, (left) raw image, (right) labeled image; (a) random spheres and (b) sand system.
4. Particle centers In order to obtain microstructure properties of the system, the coordinates of the center of gravity of a particle is computed as [13] P if ði; j; kÞ i , (2) cgi ¼ volume of particle P jf ði; j; kÞ cgj ¼
cgk ¼
j
volume of particle P kf ði; j; kÞ k
volume of particle
,
(3)
,
(4)
where cgi, cgj and cgk are the global coordinates of the center of gravity of the particle and f(i, j, k) is the corresponding voxel value of the particle, f(i, j, k) ¼ 1 in a binary image. The volume of the particle is the total number of voxels that have the same label of the particle as obtained from the particle identification algorithm. 5. Particle size distribution Two different approaches are proposed to measure the sphere-equivalent diameter of an irregular particle. The first approach is based on the principal lengths of the particle, which can be determined using the concept
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
550
of moment of inertia. The algorithm proceeds in five steps: (1) the moment of inertia matrix of each particle is evaluated using Eqs. (5) and (6) below; (2) the eigenvectors of the inertia matrix are computed and used as the direction cosines or the rotation matrix; (3) the matrix that contains the positions of each voxel of each particle is rotated using the rotation matrix to coincide with the axis of the 3D image; (4) major, minor and intermediate lengths of each particle are computed by searching for the maximum and minimum positions of the rotated particle along the axis of the image; (5) the diameter of the particle is calculated as the average value of the major, intermediate, and minor lengths. The moment of inertia matrix is calculated as 2 3 I xx I xy I xz 6 7 I ¼ 4 I yx I yy I yz 5, (5) I zx I zy I zz where I xx ¼
X
j cgj
2
þ
X
j
I yy ¼
X X
(6a)
2 k cgk ,
(6b)
2 j cgj ,
(6c)
k
i cgi
2
þ
X
i
I zz ¼
2 k cgk ,
k
i cgi
2
þ
X
i
I xy ¼ I yx ¼
j
X
i cgi
i
I xz ¼ I zx ¼
X X j
j cgj ,
(6d)
k cgk ,
(6e)
j
i cgi
X
i
I yz ¼ I zy ¼
X
k
j cgj
X
k cgk .
(6f)
k
The second approach for determination the diameter of a particle is based on the distance between the center of gravity of the particle and its boundary voxels. Boundary voxels are determined according to the algorithm described above. The diameter of a particle is calculated as 1=2 Y D¼2 ðcgi V Bi Þ2 þ ðcgj V Bj Þ2 þ ðcgk V Bk Þ2 , (7) B
Boundary voxel, VB
r2
r3 r4
r1 Centeroid
Fig. 4. 2D illustration of the computation of the diameter of a particle.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
551
6 D50=26.53 voxels = 0.305 mm
Frequency (%)
5
D90=37.71 voxels = 0.434 mm
4 D50=28.59 voxels = 0.329 mm D90=40.89 voxels = 0.470 mm
3
2
1
0 0
20
40
60
80
100
D (voxels)
Fig. 5. Particle size distribution of the sand system calculated by two different algorithms; the mean distance between center to boundary voxels (dashed line) and the mean of the principal lengths of particles (solid line).
where V Bi ; V Bj and V Bk are the row, column and depth indices of the boundary voxel VB. Fig. 4 is a simplified 2D illustration of the approach. Fig. 5 shows the particle size distribution as obtained by the two different algorithms. The shape of the distribution is similar in both cases; however, particle diameters obtained based on the principle lengths are larger than those obtained based on the average distance between the center of the particle and its boundary voxels. In both approaches, the value of D50 is close to the value reported in Ref. [10]. 6. Particle shape indices Particle shape can be characterized by two microstructure parameters; sphericity and angularity. Simulations using micromechanical models indicate that shape indices have significant impact on macromechanical behavior of granular materials [3]. Due to the complexity of obtaining physically based characterization of shape parameters, most micromechanical models assume either constant values or arbitrary distributions of shape parameters [14]. In this paper, I present physically representative distributions of sphericity and angularity of granular materials from high resolution 3D images. The sphericity index I PSph of a given particle, p, is computed as I PSph
N Pvox 1 X DPi DPi ¼ P , N vox i¼1 DPmin DPmax
(8)
where N Pvox is the total number of boundary voxels of particle p; DPi the diameter that originates from boundary voxel i and passes through the centeroid of particle p; DPmin and DPmax are the minimum and maximum lengths of particle p, respectively. The roundness index (the inverse of the angularity), I PR , of a given particle, p is calculated as I PR ¼
Ap P 2 , D þDP 4p min 4 max
(9)
where Ap is the surface area of particle p. Surface area is computed by the summation the areas of exterior faces of boundary voxels. Fig. 6 shows distribution of sphericity index of the sand system; indices range between 0.04 and 1.49. I followed the classification suggested in Ref. [15] to classify particles based on their sphericity indices.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
552 3.5
min: 0.04 avg: 0.30 max: 1.49
3
Frequency (%)
2.5 2 1.5 1 0.5 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
Sphericity index
Fig. 6. Sphericity index distribution of the sand system.
7 min: 1.12
6
avg:1.25 max:1.38
Frequency (%)
5 4 3 2 1 0 1.1
1.15
1.2
1.25 1.3 Roundness index
1.35
1.4
1.45
Fig. 7. Roundness index distribution of the sand system.
Particles can be classified as: prismoidal (I PSph 41), sub-prismoidal (I PSph ¼ 1:0 0:5), spherical (I PSph ¼ 0:5 0:2), sub-discoidal (I PSph ¼ 0:2 0:1) or discoidal (I PSph ¼ 0:1 0:0). Fig. 7 shows distribution of roundness index (inverse of angularity) of the sand system; indices range between 1.12 and 1.38. Particles can be classified as: very angular (I PR 41), angular (I PR ¼ 1:5 1:4), subangular (I PR ¼ 1:4 1:3), sub-rounded (I PR ¼ 1:3 1:2), rounded (I PR ¼ 1:2 1:1) and well rounded (I PR ¼ 1:1 1:0) [15]. 7. Number of contacts The average number of particles in contact with a given particle is a critical parameter in the analysis of granular packing; it gives information about several properties in the system and is used to simulate the behavior of granular material [16,17]. Calculating the contact number geometrically is difficult due to the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558 24
sand random spheres
22 20
sand min: 2.00 avg: 7.30 max: 16.00
18 Frequency (%)
553
16 14
spheres 2.00 6.83 11.00
12 10 8 6 4 2 0 0
2
4
6
8
10
12
14
16
18
Number of contacts
Fig. 8. Distribution of the average number of particles in contact with a given particle in sand and random spheres systems.
uncertainty of calculating the accurate position of the center of gravity of the particle. Moreover, calculating the exact number of contacts geometrically becomes impossible in the case of irregular particles such as granular materials. In this paper, I developed an algorithm that gives an accurate estimate of the number of particles in contact with a given particle. The algorithm proceeds according to the following steps, (1) each particle is identified and labeled with a distinct value using particle identification algorithm described previously; (2) boundary voxels of each particle are identified; (3) the 26 neighboring voxels of each boundary voxel and their voxel values are identified. The number of contacts of a given particle is determined as the number of distinct voxel values found in step 3 excluding the ‘‘0’’ value, where ‘‘0’’ values depict the void phase; (4) steps 3 and 4 are repeated for all particles in the system. Fig. 8 shows the distribution of contacts for random spheres and sand systems. The random sphere packing is used to verify the algorithm and to compare the distribution of number of contacts of an ideal packing to the distribution of natural systems. Statistical values of the random packing are similar to the values reported for different monosized glass beads [16]. While the average number of contacts are slightly different in both systems (6.83 in the random spheres packing and 7.30 in the sand system), the distributions are different in shape and statistics. Both systems have the same minimum number of contacts (i.e., 2), whereas they have different maximum values. Approximately, 2% of the sand particles are contacted by 16 particles, which is the maximum number of contacts observed in the sand system. In the random packing system, only 1% of the spheres are contacted by 11 spheres (i.e., the maximum number of contacts). 8. Contact network The topological structure of the number of particles that are l topological distance away from a given particle is very important in modeling force propagation between particles when subjected to external forces. To determine the contact network of a system of particles, I followed the following algorithm: (1) the algorithm starts with a given particle and determines the particles that are in contact with it, this determines the particles that are at topological distance l ¼ 1. This step is achieved using the algorithm described in the previous section (number of contacts algorithm); (2) the contact network at the next topological distance (i.e., l ¼ 2) consists of all particles that are in contact with particles at previous topological distance (i.e., l ¼ 1), excluding the particles that have been counted previously; (3) the algorithm proceeds until all particles are counted; (4) the steps above are repeated for all particles in the system. Fig. 9 shows the average number of contacts of particles as a function of topological distance l. In both systems, the number of particles increases with topological distance until it reaches a certain maximum value above which the number of particles starts to decrease. The minimum value represents the average number of contacts of a given particle whereas the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
554 104
sand random spheres fitted curve fitted curve
Number of contacts at l
103
101
101
100 100
101 Topological distance l
Fig. 9. Distribution of the average number of particles in contact with a given particle at different topological distances (contact network).
maximum value indicates that the topological distance reaches the boundary of the sample. Contact values between minimum and maximum topological distances represent the valid range of contact network. The maximum topological distance is 6 and 8 for the random spheres packing and sand system, respectively. I fit the number of contacts nl at a given topological distance l to the equation suggested by Refs. [18,19]: n l ¼ b 1 l 2 þ b 2 l þ b3 ,
(10)
where b1, b2 and b3 are coefficients. b1 is a measure of the topological density of the system and it is a tool to characterize its structural stability and rigidity. I observed b1 values of 10.99 and 3.22 for random spheres and sand, respectively. b1 of the random spheres packing agrees with the crystalline arrangement classifications of disordered systems [20]. Note that this fit is only valid between the minimum and maximum topological distances as shown in Fig. 9. The effect of angularity and nonuniformity of particles can be observed by the number of contacts at a given topological distance. In the random spheres packing, more contact particles can be found at a given topological distance compared with the sand system. This is mainly due to the uniformity and angularity of the particles in the sand system.
9. Packing efficiency Packing efficiency is a measure of arrangement of particles in a given system [21,22]. The higher the packing efficiency, the more particles can be packed within a certain volume. Packing efficiency of a given system controls its bulk density; therefore, it significantly affects the micro- and macromechanical behavior of the system. Packing efficiency parameter can be determined by the number of particles within a radial distance from a given particle. Fig. 10 shows the average number of particles within a given radial distance r from a particle of a diameter D. It is observed that in the random spheres system, the number of particles within r ¼ D is very close to the average number of contacts, whereas it is much larger than the average number of contacts in the sand system. This is mainly due to the characteristics of particle size distribution of each system. It is also observed that the sand system has better packing efficiency than the random spheres packing; as the radial distance r increases the packing efficiency of the random spheres system approaches that of the sand system. This finding is expected due to the fact that size nonuniformity and angularity of the particles simplify the densification process of a certain system. The angular and different-sized particles have more chance to slide between other particles and fill the void, if it exists.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
555
103
Number of particles
sand random spheres
102
101
100 1
1.5
2
2.5
3 r/D
3.5
4
4.5
5
Fig. 10. Packing efficiency: the average number of particle centers within a radial distance r from a given particle of a diameter D.
Fig. 11. 3D visualization of boundaries of local void regions of the random spheres system; the green color depicts particles and the red color depicts boundaries of local void regions.
10. Local void ratio Local void ratio of a particle is defined as the ratio of the volume of the void space that surrounds the particle divided by the volume of the particle. The challenging problem is how to partition the total void space in local regions associated with each particle. The algorithm developed in Ref. [11] to calculate local void ratios of irregular particles from 3D images is adopted in this paper. Fig. 11 shows a 3D visualization of the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
556 35
sand random spheres
30
sand Frequency (%)
25
spheres
min: 0.30 avg: 0.69 max: 2.10
20
0.23 0.52 1.76
15 10 5 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Local void ratio
Fig. 12. Distribution of local void ratios of the random spheres and sand systems.
random spheres systems; the green color depicts particles and the red color depicts boundaries of local void regions. Fig. 12 shows frequency distributions of local void ratios for the random spheres packing and the sand system; statistics of distributions are shown in the figure. Distribution statistics and 3D visualization show that variation in the distribution of local void ratios is larger in the sand system; this observation is in good agreement with the observations of other researchers [23]. For the random spheres packing, the mean local void ratio is 0.520 and the global void ratio is 0.548; whereas for the sand system, the mean local and global void ratios are 0.690 and 0.569, respectively. The larger difference between local and global void ratios observed in the sand system is mainly attributed to the size and shape distributions of sand particles. Local void ratios range between 0.30 and 2.10 in the sand system compared to a range between 0.23 and 1.76 in the random spheres system. As expected in the sand system, there are more local configurations with large local void ratio (i.e., larger that one). This is mainly due to the orientation and size distribution of sand particles. 11. Radial distribution function The radial distribution function, g(r), is the probability of finding a particle center in a spherical shell of a radius r, given that there is another particle at the origin of the spherical shell [24]. The radial distribution function is related to the average number of particles by Z r2 pðr2 Þ pðr1 Þ ¼ gðrÞ4pr2 dr, (11) r1
where p(r) is the average number of particles in a spherical shell of radius r. Fig. 13 shows the normalized radial distribution of the average number of particle centers in a spherical shell for both systems as a function of the dimensionless distance r/D. The normalized probability distribution function is different for both systems; the sand system shows less oscillations and peaks due to large size distribution. The function for the packing of random spheres starts from a peak at r ¼ D representing maximum probability of the contact particles, then the function decreases to minimum probability at r ¼ 1.25D indicating a minimum likelihood of finding particles in contact. The probability function then increases to another local peak at r ¼ 1.65D. The function continues to oscillate until it approaches a constant value of one at r ¼ 4.5D. Radial probability function of the sand system shows different behavior; while it starts at initial peak at r ¼ 1, the probability decreases to a local minimum at r ¼ 1.25D. The curve then slightly oscillates to reach a local peak at r ¼ 1.55D. The function shows lesser oscillation compared with the random spheres system. The curve continues its damping behavior until it approaches unity.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
557
10 sand random spheres
8
g(r)
6
4
2
0 1
1.5
2
2.5
3
3.5
4
4.5
r/D
Fig. 13. Normalized radial distribution functions for particle centers in the random spheres and sand systems.
12. Summary and conclusion The paper presented methodologies and algorithms to obtain microstructure characteristics of assemblies of irregular particles that represent granular materials. I obtained fundamental microstructure properties such as centeroid, particle size distribution, distribution of number of contacts (i.e., coordination number), shape indices (i.e., sphericity and angularity), contact network, packing efficiency, distribution of local void ratio, and radial distribution function. X-ray microtomography was used to obtain high resolution 3D digitization of a natural sand system representing granular material. Microstructure characteristics of the sand system were extracted and compared with the characteristics of a computer-simulated packing of random spheres. Findings indicate that the use of simplified systems of idealized spheres to model macroscale behavior of granular systems may lead to inaccurate numerical predictions. This is mainly due to the differences in microstructure characteristics between both systems, which will be reflected in the interparticle interactions. Methods presented in this paper can be incorporated in micromechanical models to simulate, understand, or explain macroscale behavior of granular materials based on their real microstructure characteristics. Acknowledgments The image of the sand system shown in Fig. 3(b) was acquired at GeoSoilEnviroCARS (sector 13), Advanced photon source (APS), Argonne National Laboratory. The author greatly acknowledges the support of the staff of GeoSoilEnviroCARS. GeoSoilEnviroCARS is supported by the National Science Foundation-Earth Sciences (EAR-0217473), Department of Energy–Geosciences (DE-FG02-94ER14466) and the state of Illinois. Use of the APS was supported by the US Department of Energy, Basic Energy Sciences, Office of Science, under Contract no. W31-109-Eng-38. References [1] [2] [3] [4] [5] [6]
M. Satake, in: H.H. Shen et al. (Eds.), Advances in Micromechanics of Granular Materials, 1992. M.I. Alsaleh, K.A. Alshibli, G. Voyiadjis, ASCE Int. J. Geomech. 6 (2006) 248. M. Nicodemi, A. Coniglio, H.J. Hermann, Phys. Rev. E 55 (1997) 3962. S.F. Edwards, D.V. Grinev, Chaos 9 (1999) 551. P. Richard, P. Philippe, F. Barbe, S. Bourles, X. Thibault, D. Bideau, Phys. Rev. E 68 (2003) 020301. G.T. Seidler, G. Martinez, L.H. Seeley, K.H. Kim, E.A. Behne, S. Zaranek, B.D. Chapman, S.M. Heald, Phys. Rev. E 62 (2000) 8175.
ARTICLE IN PRESS 558 [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
R. Al-Raoush / Physica A 377 (2007) 545–558 M.A. Taylor, E.J. Garboczi, S.T. Erdogan, D.W. Fowler, Powder Technol. 162 (2006) 1. L.B. Wang, J.D. Frost, J.S. Lai, J. Comput. Civ. Eng. 18 (2004) 28. G. Liu, K.E. Thompson, Powder Technol. 113 (2000) 185. R. Al-Raoush, C.S. Willson, J. Hydrol. 300 (2005) 44. A.H. Reed, R. Al-Raoush, C.S. Willson, J.E. Cable, in: Proceedings of the AGU meeting, Washington, DC, May 28–31, 2002. R. Al-Raoush, K.A. Alshibli, Physica A 361 (2006) 441. G. Lohman, Volumetric Image Analysis, Wiley, New York, 1998. M.I. Alsaleh, Numerical modeling for strain localization in granular materials using Cosserat theory enhanced with microfabric properties. Ph.D Thesis, Department of Civil and Environmental Engineering, Louisiana State University, 2004. K.A. Alshibli, M.I. Alsaleh, ASCE J. Comput. Civ. Eng. 18 (2004) 36. T. Aste, M. Saadatfar, T.J. Senden, Phys. Rev. E 71 (2005) 061302. S.F. Edwards, D.V. Grinev, Physica A 263 (1999) 545. G.O. Brunner, J. Solid State Chem. 29 (1979) 41. J. Conway, N. Sloane, Proc. R. Soc. London Ser. A 453 (1997) 2369. M. O’Keeffe, Z. Kristallogr. 196 (1991) 21. T. Aste, M. Saadatfar, A. Sakellariou, T. Senden, Physica A 339 (2004) 16. T. Aste, T. Di Matteo, S.T. Hyde, Physica A 346 (2005) 20. E.R. Nowak, J.B. Knight, E. BenNaim, H.M. Jeager, S.R. Nagel, Phys. Rev. E 57 (1998) 1971. S. Torquato, Random Heterogonous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2001.
Physica A 377 (2007) 545–558 www.elsevier.com/locate/physa
Microstructure characterization of granular materials Riyadh Al-Raoush Department of Civil and Environmental Engineering, Southern University and A and M College, Baton Rouge, LA 70813, USA Received 12 June 2006; received in revised form 2 November 2006 Available online 3 January 2007
Abstract This paper presents techniques and algorithms to compute microstructure properties of irregular-shaped granulate assemblies utilizing 3D images. The techniques are capable of extracting microstructure properties of particles such as centeroid, particle size distribution, shape indices (i.e., sphericiy and angularity), number of contacts (i.e., distribution of coordination numbers), contact network, packing efficiency, distribution of local void ratio and radial distribution function. Such properties are critical parameters for micromechanical-based numerical models to capture micro- and macromechanical behavior of geomaterials. X-ray microtomography was used to reconstruct high-resolution 3D image of a natural sand system to represent granular materials. Microstructure properties of the sand system were computed and compared with properties of a computer-simulated image of periodic random spheres. Findings indicate that the use of simplified systems of idealized spheres to model micro- and macromechanical behavior of granular systems can lead to inaccurate results due to the differences in the microstructure between both systems. Methods presented in this paper enabled capturing a more realistic microstructure that can be incorporated in micromechanical models to better simulate, understand, or explain macroscale behavior of granular materials based on their actual microstructure. r 2007 Published by Elsevier B.V. Keywords: Microstructure; Granular materials; Micromechanical properties; Microtomography; Random packing
1. Introduction Granular material, one type of frictional geomaterials, is often composed of individual grains that vary in shape, size and surface texture. Such parameters can significantly affect the final packing of granulates and their contact distribution characteristics. Microstructure characterization of irregular particles such as granular materials, pharmaceutical materials and powders, is of a significant importance to many engineering applications, particularly to models that implement micromechanical theories to describe the macroscale behavior [1,2]. It has been recognized that geometry and topology of a granular system control its electrical, thermal, fluid-like flow, stress transfer mechanisms and other mechanical properties. For instance, mechanical conditions of a granular system such as shearing, compaction, transmission of stress and statistical properties of contact force distribution are greatly influenced by its microstructure characteristics (e.g., the distribution of the number of contacts of particles) [3–5]. Considering two adjacent particles within the assembly, one needs Tel.: +1 225 771 5870; fax: +1 225 771 4320.
E-mail address: [email protected]. 0378-4371/$ - see front matter r 2007 Published by Elsevier B.V. doi:10.1016/j.physa.2006.11.090
ARTICLE IN PRESS 546
R. Al-Raoush / Physica A 377 (2007) 545–558
to account for the exact interaction characteristics so that the interparticle forces and relative motions can be estimated. The interparticle interaction is essentially controlled by size and shape of the contact area. Furthermore, the macro- or global behavior of granulate assembly such as fluid-like flow, stress–strain response and multi-phase flows can be described by microbased models. Such models require information about microstructure properties such as local void ratio distribution, contact number distributions and particle size distribution. These properties require robust techniques for characterization. Development of such techniques has been a real challenging issue for many researchers, moreover, there are complications associated with the implementation of the actual microstructural properties in numerical models. To overcome the complexity associated with real granular systems, granular materials are often represented and modeled as ideal packing assemblies of rigid spherical particles, elliptical particles and/or clusters of spheres [1,6]. These models have been found to serve as a good starting point to account, to some extent, for the particle shape effects on the global behavior of granulates. However, fundamental problems and limitations of their applicability to the real granular materials are obvious. Therefore, the actual spatial distributions of the microstructure properties are not quite incorporated and thus, fundamental physics can be missing. While a few studies have been conducted to characterize microstructure properties of irregular systems, these studies have been limited to determination of basic microstructure properties such as local void ratio, surface area, and density [7,8]. The objectives of this paper are: (1) to present algorithms to compute microstructure properties of granular materials from 3D images, and (2) to compare microstructure properties of a computersimulated pack of random spheres to those of natural systems (i.e., sand systems). Synchrotron X-ray microtomography is used to obtain the 3D image of the natural sand system. 2. Materials In this paper, I present two different systems: computer-simulated 3D image of random spheres and X-ray microtomography image of a sand system. The first system is periodic random packing consisting of 1000 spheres uniformly sized and digitized at 25 voxels/sphere diameter. The size of the 3D image is 200 200 200 voxels and has a porosity value of 0.354. The packing was generated by the algorithm developed by Liu and Thompson [9] based on a collective rearrangement method. In these type of models, sphere packs are generated through two steps; random generation and placement of spheres in a given domain allowing them to overlap. Spheres are usually generated from a given probability density function. An iterative relaxation process, such as Monte Carlo or molecular dynamics simulations, is then applied to separate the overlapped spheres to reach a stable position with a specific overlap tolerance. Sphere sizes are reduced if the overlap between spheres is larger than the specified tolerance. The relaxation process is repeated until a stable packing with acceptable sphere overlap is obtained. More details about the algorithm used to generate the sphere packing in this paper can be found in Ref. [9]. The second system is a natural sand system [10,11]. The system was imaged by synchrotron X-ray tomography at the GeoSoilEnviroCARS (GSECARS) BM-13D beamline at the Advanced Photon Source, Argonne National Laboratory. Image resolution is 11.5 mm. Image reconstruction was performed using algorithms developed by GSECARS to convert CT attenuation data to 3D volumetric data. The size of the image is 375 400 340 voxels. The porosity of the system is 0.363 and the mean particle size, D50, is 0.325 mm. Note that porosity values of both systems are similar; this in turn provides reasonable base of comparison between microstructure properties of natural and computer-simulated systems of similar porosity values. 3. Recognition of individual particles Recognition of individual particles in a given system is the first step required to determine its microstructure properties. Particle identification of 3D irregular particles such as granular materials (e.g., sand) is a challenging task. In this paper, I modified the algorithm developed by Al-Raoush that labels individual particles and improved its accuracy and speed [12]. The algorithm is based on watershed transform and proceeds by constructing a distance map of the voxels of the particles in the image. Fig. 1 illustrates the procedure to distinguish particles in contact at an area. The procedure is best explained through the steps shown in Fig. 1 as follows: Fig. 1(a) shows a binary image of 2 touched particles, the image has 2 voxel values;
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
547
Fig. 1. Separation of touched particles; (a) binary image of two touched particles, (b) distance transform of the two particles, (c) flooding simulation process starts from local maxima to define watershed line that separates the particles, (d) watershed line (dividing line) and (e) a voxel is considered part of the dividing line if the filling enters a voxel from the indicated directions [12].
values of ‘‘1’’ depict the foreground (i.e., particles) and values of ‘‘0’’ represent the background. Fig. 1(b) shows the distance map of the particles calculated by distance transform. Each foreground voxel in the image is labeled based on its distance from the nearest background voxel. The distance transform is applied to the foreground voxels of the image starting from the boundaries of the particles towards their centeroids to obtain local maxima. Local maxima are markers that distinguish distinct (individual) particles in the image. In Fig. 1(b), local maxima are the voxels labeled by ‘‘11’’ and ‘‘6’’ in the large and small particles, respectively. Distance map of a 3D image is computed by calculating the distance X ab Euclidean between two voxels aðia ; j a ; ka Þ and bðib ; j b ; kb Þ using the Euclidean metric as 2 2 2 1=2 X ab , Euclidean ¼ ðia ib Þ þ ðj a j b Þ þ ðka k b Þ
(1)
where ia, ja and ka are the row, column, and depth indices of a given voxel (i.e., a, b, y), respectively, in the 3D image. It is worth noting that the row, column and depth indices are equivalent to the x, y and z distances in the Cartesian coordinate system, respectively. The distance map algorithm commences by labeling boundary voxels of the particles in the image by assigning label ‘‘1’’ to each voxel indicating that they are one voxel away from the closest background voxel. A boundary voxel is the one that connects a background voxel through a face, edge or corner. The algorithm searches the 26 neighboring voxels of each foreground voxel and counts their voxel values. If the summation is 26, then the selected voxel is an interior voxel (i.e., surrounded by 26 foreground voxels); otherwise the selected voxel touches at least one background voxel and is considered a boundary voxel. The algorithm proceeds until all foreground voxels are inspected and all boundary voxels are identified and labeled. The algorithm proceeds and searches foreground voxels to identify the neighboring voxels of the voxels that was labeled in the previous step; voxel distance (i.e., voxel labeling) is calculated between these voxels and the closest background voxel using Eq. (1). The algorithm proceeds until all foreground voxels are labeled indicating their voxel distance from the closest background voxel (i.e., distance map of foreground voxels). Fig. 1(c) shows the flooding simulation that starts from local maxima shown in Fig. 1(b). The algorithm simultaneously labels voxels that have the highest values in the distance map and assigns a label ‘‘1’’ to them (i.e., voxels that are labeled by ‘‘11’’ in Fig. 1(b)). The second iteration starts by simultaneously labeling voxels that have the second highest values in the distance map and assigns a label ‘‘2’’ to them (i.e., voxels that are labeled by ‘‘10’’ in Fig. 1(b)). The same labeling procedure used in the previous iteration is repeated until all voxels are labeled. In the flooding simulation, each marker represents a distinct particle (catchment basin) where the simulation moves from the markers with a constant rate of labeling. Dividing lines (watershed lines)
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
548
a
c
b
Local maxima (markers)
Catchment basins
Watershed lines
Direction of flood simulation Local maxima (markers)
d
e Watershed line
Catchment basins Fig. 2. Identification of individual particles; (a) image of two touched particles, (b) distance transform of the two particles, (c) flooding simulation process to define watershed line that separates the particles, (d) watershed line and catchment basins (particles) and (e) labeled particles.
between distinct particles are constructed at voxels where labeling comes from different catchmnet basins (i.e., particles). Fig. 1(e) illustrates the directions that decide voxels of dividing lines. Fig. 2 shows an overview of the algorithm of particle identification and Fig. 3 shows particle labeling of random spheres and sand systems. Image resolution and quality impact the results of particle identification algorithm and other algorithms presented in this paper. Image quality can be improved by applying specific types of filters based on the artifact type of the image. The impact of the resolution on a given parameter can be determined by carrying out a resolution sensitivity analysis where tolerance values can be determined.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
549
Fig. 3. Identification of individual particles, (left) raw image, (right) labeled image; (a) random spheres and (b) sand system.
4. Particle centers In order to obtain microstructure properties of the system, the coordinates of the center of gravity of a particle is computed as [13] P if ði; j; kÞ i , (2) cgi ¼ volume of particle P jf ði; j; kÞ cgj ¼
cgk ¼
j
volume of particle P kf ði; j; kÞ k
volume of particle
,
(3)
,
(4)
where cgi, cgj and cgk are the global coordinates of the center of gravity of the particle and f(i, j, k) is the corresponding voxel value of the particle, f(i, j, k) ¼ 1 in a binary image. The volume of the particle is the total number of voxels that have the same label of the particle as obtained from the particle identification algorithm. 5. Particle size distribution Two different approaches are proposed to measure the sphere-equivalent diameter of an irregular particle. The first approach is based on the principal lengths of the particle, which can be determined using the concept
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
550
of moment of inertia. The algorithm proceeds in five steps: (1) the moment of inertia matrix of each particle is evaluated using Eqs. (5) and (6) below; (2) the eigenvectors of the inertia matrix are computed and used as the direction cosines or the rotation matrix; (3) the matrix that contains the positions of each voxel of each particle is rotated using the rotation matrix to coincide with the axis of the 3D image; (4) major, minor and intermediate lengths of each particle are computed by searching for the maximum and minimum positions of the rotated particle along the axis of the image; (5) the diameter of the particle is calculated as the average value of the major, intermediate, and minor lengths. The moment of inertia matrix is calculated as 2 3 I xx I xy I xz 6 7 I ¼ 4 I yx I yy I yz 5, (5) I zx I zy I zz where I xx ¼
X
j cgj
2
þ
X
j
I yy ¼
X X
(6a)
2 k cgk ,
(6b)
2 j cgj ,
(6c)
k
i cgi
2
þ
X
i
I zz ¼
2 k cgk ,
k
i cgi
2
þ
X
i
I xy ¼ I yx ¼
j
X
i cgi
i
I xz ¼ I zx ¼
X X j
j cgj ,
(6d)
k cgk ,
(6e)
j
i cgi
X
i
I yz ¼ I zy ¼
X
k
j cgj
X
k cgk .
(6f)
k
The second approach for determination the diameter of a particle is based on the distance between the center of gravity of the particle and its boundary voxels. Boundary voxels are determined according to the algorithm described above. The diameter of a particle is calculated as 1=2 Y D¼2 ðcgi V Bi Þ2 þ ðcgj V Bj Þ2 þ ðcgk V Bk Þ2 , (7) B
Boundary voxel, VB
r2
r3 r4
r1 Centeroid
Fig. 4. 2D illustration of the computation of the diameter of a particle.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
551
6 D50=26.53 voxels = 0.305 mm
Frequency (%)
5
D90=37.71 voxels = 0.434 mm
4 D50=28.59 voxels = 0.329 mm D90=40.89 voxels = 0.470 mm
3
2
1
0 0
20
40
60
80
100
D (voxels)
Fig. 5. Particle size distribution of the sand system calculated by two different algorithms; the mean distance between center to boundary voxels (dashed line) and the mean of the principal lengths of particles (solid line).
where V Bi ; V Bj and V Bk are the row, column and depth indices of the boundary voxel VB. Fig. 4 is a simplified 2D illustration of the approach. Fig. 5 shows the particle size distribution as obtained by the two different algorithms. The shape of the distribution is similar in both cases; however, particle diameters obtained based on the principle lengths are larger than those obtained based on the average distance between the center of the particle and its boundary voxels. In both approaches, the value of D50 is close to the value reported in Ref. [10]. 6. Particle shape indices Particle shape can be characterized by two microstructure parameters; sphericity and angularity. Simulations using micromechanical models indicate that shape indices have significant impact on macromechanical behavior of granular materials [3]. Due to the complexity of obtaining physically based characterization of shape parameters, most micromechanical models assume either constant values or arbitrary distributions of shape parameters [14]. In this paper, I present physically representative distributions of sphericity and angularity of granular materials from high resolution 3D images. The sphericity index I PSph of a given particle, p, is computed as I PSph
N Pvox 1 X DPi DPi ¼ P , N vox i¼1 DPmin DPmax
(8)
where N Pvox is the total number of boundary voxels of particle p; DPi the diameter that originates from boundary voxel i and passes through the centeroid of particle p; DPmin and DPmax are the minimum and maximum lengths of particle p, respectively. The roundness index (the inverse of the angularity), I PR , of a given particle, p is calculated as I PR ¼
Ap P 2 , D þDP 4p min 4 max
(9)
where Ap is the surface area of particle p. Surface area is computed by the summation the areas of exterior faces of boundary voxels. Fig. 6 shows distribution of sphericity index of the sand system; indices range between 0.04 and 1.49. I followed the classification suggested in Ref. [15] to classify particles based on their sphericity indices.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
552 3.5
min: 0.04 avg: 0.30 max: 1.49
3
Frequency (%)
2.5 2 1.5 1 0.5 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
Sphericity index
Fig. 6. Sphericity index distribution of the sand system.
7 min: 1.12
6
avg:1.25 max:1.38
Frequency (%)
5 4 3 2 1 0 1.1
1.15
1.2
1.25 1.3 Roundness index
1.35
1.4
1.45
Fig. 7. Roundness index distribution of the sand system.
Particles can be classified as: prismoidal (I PSph 41), sub-prismoidal (I PSph ¼ 1:0 0:5), spherical (I PSph ¼ 0:5 0:2), sub-discoidal (I PSph ¼ 0:2 0:1) or discoidal (I PSph ¼ 0:1 0:0). Fig. 7 shows distribution of roundness index (inverse of angularity) of the sand system; indices range between 1.12 and 1.38. Particles can be classified as: very angular (I PR 41), angular (I PR ¼ 1:5 1:4), subangular (I PR ¼ 1:4 1:3), sub-rounded (I PR ¼ 1:3 1:2), rounded (I PR ¼ 1:2 1:1) and well rounded (I PR ¼ 1:1 1:0) [15]. 7. Number of contacts The average number of particles in contact with a given particle is a critical parameter in the analysis of granular packing; it gives information about several properties in the system and is used to simulate the behavior of granular material [16,17]. Calculating the contact number geometrically is difficult due to the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558 24
sand random spheres
22 20
sand min: 2.00 avg: 7.30 max: 16.00
18 Frequency (%)
553
16 14
spheres 2.00 6.83 11.00
12 10 8 6 4 2 0 0
2
4
6
8
10
12
14
16
18
Number of contacts
Fig. 8. Distribution of the average number of particles in contact with a given particle in sand and random spheres systems.
uncertainty of calculating the accurate position of the center of gravity of the particle. Moreover, calculating the exact number of contacts geometrically becomes impossible in the case of irregular particles such as granular materials. In this paper, I developed an algorithm that gives an accurate estimate of the number of particles in contact with a given particle. The algorithm proceeds according to the following steps, (1) each particle is identified and labeled with a distinct value using particle identification algorithm described previously; (2) boundary voxels of each particle are identified; (3) the 26 neighboring voxels of each boundary voxel and their voxel values are identified. The number of contacts of a given particle is determined as the number of distinct voxel values found in step 3 excluding the ‘‘0’’ value, where ‘‘0’’ values depict the void phase; (4) steps 3 and 4 are repeated for all particles in the system. Fig. 8 shows the distribution of contacts for random spheres and sand systems. The random sphere packing is used to verify the algorithm and to compare the distribution of number of contacts of an ideal packing to the distribution of natural systems. Statistical values of the random packing are similar to the values reported for different monosized glass beads [16]. While the average number of contacts are slightly different in both systems (6.83 in the random spheres packing and 7.30 in the sand system), the distributions are different in shape and statistics. Both systems have the same minimum number of contacts (i.e., 2), whereas they have different maximum values. Approximately, 2% of the sand particles are contacted by 16 particles, which is the maximum number of contacts observed in the sand system. In the random packing system, only 1% of the spheres are contacted by 11 spheres (i.e., the maximum number of contacts). 8. Contact network The topological structure of the number of particles that are l topological distance away from a given particle is very important in modeling force propagation between particles when subjected to external forces. To determine the contact network of a system of particles, I followed the following algorithm: (1) the algorithm starts with a given particle and determines the particles that are in contact with it, this determines the particles that are at topological distance l ¼ 1. This step is achieved using the algorithm described in the previous section (number of contacts algorithm); (2) the contact network at the next topological distance (i.e., l ¼ 2) consists of all particles that are in contact with particles at previous topological distance (i.e., l ¼ 1), excluding the particles that have been counted previously; (3) the algorithm proceeds until all particles are counted; (4) the steps above are repeated for all particles in the system. Fig. 9 shows the average number of contacts of particles as a function of topological distance l. In both systems, the number of particles increases with topological distance until it reaches a certain maximum value above which the number of particles starts to decrease. The minimum value represents the average number of contacts of a given particle whereas the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
554 104
sand random spheres fitted curve fitted curve
Number of contacts at l
103
101
101
100 100
101 Topological distance l
Fig. 9. Distribution of the average number of particles in contact with a given particle at different topological distances (contact network).
maximum value indicates that the topological distance reaches the boundary of the sample. Contact values between minimum and maximum topological distances represent the valid range of contact network. The maximum topological distance is 6 and 8 for the random spheres packing and sand system, respectively. I fit the number of contacts nl at a given topological distance l to the equation suggested by Refs. [18,19]: n l ¼ b 1 l 2 þ b 2 l þ b3 ,
(10)
where b1, b2 and b3 are coefficients. b1 is a measure of the topological density of the system and it is a tool to characterize its structural stability and rigidity. I observed b1 values of 10.99 and 3.22 for random spheres and sand, respectively. b1 of the random spheres packing agrees with the crystalline arrangement classifications of disordered systems [20]. Note that this fit is only valid between the minimum and maximum topological distances as shown in Fig. 9. The effect of angularity and nonuniformity of particles can be observed by the number of contacts at a given topological distance. In the random spheres packing, more contact particles can be found at a given topological distance compared with the sand system. This is mainly due to the uniformity and angularity of the particles in the sand system.
9. Packing efficiency Packing efficiency is a measure of arrangement of particles in a given system [21,22]. The higher the packing efficiency, the more particles can be packed within a certain volume. Packing efficiency of a given system controls its bulk density; therefore, it significantly affects the micro- and macromechanical behavior of the system. Packing efficiency parameter can be determined by the number of particles within a radial distance from a given particle. Fig. 10 shows the average number of particles within a given radial distance r from a particle of a diameter D. It is observed that in the random spheres system, the number of particles within r ¼ D is very close to the average number of contacts, whereas it is much larger than the average number of contacts in the sand system. This is mainly due to the characteristics of particle size distribution of each system. It is also observed that the sand system has better packing efficiency than the random spheres packing; as the radial distance r increases the packing efficiency of the random spheres system approaches that of the sand system. This finding is expected due to the fact that size nonuniformity and angularity of the particles simplify the densification process of a certain system. The angular and different-sized particles have more chance to slide between other particles and fill the void, if it exists.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
555
103
Number of particles
sand random spheres
102
101
100 1
1.5
2
2.5
3 r/D
3.5
4
4.5
5
Fig. 10. Packing efficiency: the average number of particle centers within a radial distance r from a given particle of a diameter D.
Fig. 11. 3D visualization of boundaries of local void regions of the random spheres system; the green color depicts particles and the red color depicts boundaries of local void regions.
10. Local void ratio Local void ratio of a particle is defined as the ratio of the volume of the void space that surrounds the particle divided by the volume of the particle. The challenging problem is how to partition the total void space in local regions associated with each particle. The algorithm developed in Ref. [11] to calculate local void ratios of irregular particles from 3D images is adopted in this paper. Fig. 11 shows a 3D visualization of the
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
556 35
sand random spheres
30
sand Frequency (%)
25
spheres
min: 0.30 avg: 0.69 max: 2.10
20
0.23 0.52 1.76
15 10 5 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Local void ratio
Fig. 12. Distribution of local void ratios of the random spheres and sand systems.
random spheres systems; the green color depicts particles and the red color depicts boundaries of local void regions. Fig. 12 shows frequency distributions of local void ratios for the random spheres packing and the sand system; statistics of distributions are shown in the figure. Distribution statistics and 3D visualization show that variation in the distribution of local void ratios is larger in the sand system; this observation is in good agreement with the observations of other researchers [23]. For the random spheres packing, the mean local void ratio is 0.520 and the global void ratio is 0.548; whereas for the sand system, the mean local and global void ratios are 0.690 and 0.569, respectively. The larger difference between local and global void ratios observed in the sand system is mainly attributed to the size and shape distributions of sand particles. Local void ratios range between 0.30 and 2.10 in the sand system compared to a range between 0.23 and 1.76 in the random spheres system. As expected in the sand system, there are more local configurations with large local void ratio (i.e., larger that one). This is mainly due to the orientation and size distribution of sand particles. 11. Radial distribution function The radial distribution function, g(r), is the probability of finding a particle center in a spherical shell of a radius r, given that there is another particle at the origin of the spherical shell [24]. The radial distribution function is related to the average number of particles by Z r2 pðr2 Þ pðr1 Þ ¼ gðrÞ4pr2 dr, (11) r1
where p(r) is the average number of particles in a spherical shell of radius r. Fig. 13 shows the normalized radial distribution of the average number of particle centers in a spherical shell for both systems as a function of the dimensionless distance r/D. The normalized probability distribution function is different for both systems; the sand system shows less oscillations and peaks due to large size distribution. The function for the packing of random spheres starts from a peak at r ¼ D representing maximum probability of the contact particles, then the function decreases to minimum probability at r ¼ 1.25D indicating a minimum likelihood of finding particles in contact. The probability function then increases to another local peak at r ¼ 1.65D. The function continues to oscillate until it approaches a constant value of one at r ¼ 4.5D. Radial probability function of the sand system shows different behavior; while it starts at initial peak at r ¼ 1, the probability decreases to a local minimum at r ¼ 1.25D. The curve then slightly oscillates to reach a local peak at r ¼ 1.55D. The function shows lesser oscillation compared with the random spheres system. The curve continues its damping behavior until it approaches unity.
ARTICLE IN PRESS R. Al-Raoush / Physica A 377 (2007) 545–558
557
10 sand random spheres
8
g(r)
6
4
2
0 1
1.5
2
2.5
3
3.5
4
4.5
r/D
Fig. 13. Normalized radial distribution functions for particle centers in the random spheres and sand systems.
12. Summary and conclusion The paper presented methodologies and algorithms to obtain microstructure characteristics of assemblies of irregular particles that represent granular materials. I obtained fundamental microstructure properties such as centeroid, particle size distribution, distribution of number of contacts (i.e., coordination number), shape indices (i.e., sphericity and angularity), contact network, packing efficiency, distribution of local void ratio, and radial distribution function. X-ray microtomography was used to obtain high resolution 3D digitization of a natural sand system representing granular material. Microstructure characteristics of the sand system were extracted and compared with the characteristics of a computer-simulated packing of random spheres. Findings indicate that the use of simplified systems of idealized spheres to model macroscale behavior of granular systems may lead to inaccurate numerical predictions. This is mainly due to the differences in microstructure characteristics between both systems, which will be reflected in the interparticle interactions. Methods presented in this paper can be incorporated in micromechanical models to simulate, understand, or explain macroscale behavior of granular materials based on their real microstructure characteristics. Acknowledgments The image of the sand system shown in Fig. 3(b) was acquired at GeoSoilEnviroCARS (sector 13), Advanced photon source (APS), Argonne National Laboratory. The author greatly acknowledges the support of the staff of GeoSoilEnviroCARS. GeoSoilEnviroCARS is supported by the National Science Foundation-Earth Sciences (EAR-0217473), Department of Energy–Geosciences (DE-FG02-94ER14466) and the state of Illinois. Use of the APS was supported by the US Department of Energy, Basic Energy Sciences, Office of Science, under Contract no. W31-109-Eng-38. References [1] [2] [3] [4] [5] [6]
M. Satake, in: H.H. Shen et al. (Eds.), Advances in Micromechanics of Granular Materials, 1992. M.I. Alsaleh, K.A. Alshibli, G. Voyiadjis, ASCE Int. J. Geomech. 6 (2006) 248. M. Nicodemi, A. Coniglio, H.J. Hermann, Phys. Rev. E 55 (1997) 3962. S.F. Edwards, D.V. Grinev, Chaos 9 (1999) 551. P. Richard, P. Philippe, F. Barbe, S. Bourles, X. Thibault, D. Bideau, Phys. Rev. E 68 (2003) 020301. G.T. Seidler, G. Martinez, L.H. Seeley, K.H. Kim, E.A. Behne, S. Zaranek, B.D. Chapman, S.M. Heald, Phys. Rev. E 62 (2000) 8175.
ARTICLE IN PRESS 558 [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
R. Al-Raoush / Physica A 377 (2007) 545–558 M.A. Taylor, E.J. Garboczi, S.T. Erdogan, D.W. Fowler, Powder Technol. 162 (2006) 1. L.B. Wang, J.D. Frost, J.S. Lai, J. Comput. Civ. Eng. 18 (2004) 28. G. Liu, K.E. Thompson, Powder Technol. 113 (2000) 185. R. Al-Raoush, C.S. Willson, J. Hydrol. 300 (2005) 44. A.H. Reed, R. Al-Raoush, C.S. Willson, J.E. Cable, in: Proceedings of the AGU meeting, Washington, DC, May 28–31, 2002. R. Al-Raoush, K.A. Alshibli, Physica A 361 (2006) 441. G. Lohman, Volumetric Image Analysis, Wiley, New York, 1998. M.I. Alsaleh, Numerical modeling for strain localization in granular materials using Cosserat theory enhanced with microfabric properties. Ph.D Thesis, Department of Civil and Environmental Engineering, Louisiana State University, 2004. K.A. Alshibli, M.I. Alsaleh, ASCE J. Comput. Civ. Eng. 18 (2004) 36. T. Aste, M. Saadatfar, T.J. Senden, Phys. Rev. E 71 (2005) 061302. S.F. Edwards, D.V. Grinev, Physica A 263 (1999) 545. G.O. Brunner, J. Solid State Chem. 29 (1979) 41. J. Conway, N. Sloane, Proc. R. Soc. London Ser. A 453 (1997) 2369. M. O’Keeffe, Z. Kristallogr. 196 (1991) 21. T. Aste, M. Saadatfar, A. Sakellariou, T. Senden, Physica A 339 (2004) 16. T. Aste, T. Di Matteo, S.T. Hyde, Physica A 346 (2005) 20. E.R. Nowak, J.B. Knight, E. BenNaim, H.M. Jeager, S.R. Nagel, Phys. Rev. E 57 (1998) 1971. S. Torquato, Random Heterogonous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2001.