SSSAJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via ISI Web of Science (19)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Al-Raoush, R.
Right arrow Articles by Willson, C. S.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Al-Raoush, R.
Right arrow Articles by Willson, C. S.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Al-Raoush, R.
Right arrow Articles by Willson, C. S.
Related Collections
Right arrow Tomography
Right arrow Other Approaches
Right arrow Pore-Scale Modeling
Published in Soil Sci. Soc. Am. J. 67:1687-1700 (2003).
© 2003 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

DIVISION S-1—SOIL PHYSICS

Comparison of Network Generation Techniques for Unconsolidated Porous Media

Riyadh Al-Raousha, Karsten Thompsonb and Clinton S. Willson*,a

a Dep. of Civil and Environmental Engineering, Louisiana State Univ., Baton Rouge, LA 70808
b Dep. of Chemical Engineering, Louisiana State Univ., Baton Rouge, LA 70808

* Corresponding author (cwillson{at}lsu.edu).


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 BACKGROUND
 Comparison between Medial Axis...
 CONCLUSIONS
 REFERENCES
 
While network models of porous materials have traditionally been constructed using regular or disordered lattices, recent developments allow the direct modeling of more realistic structures such as sphere packings, microtomographic images, or computer-simulated materials. One of the obstacles in these newer approaches is the generation of network structures that are physically representative of the real systems. In this paper, we present and compare two different algorithms to extract pore network parameters from three-dimensional images of unconsolidated porous media systems. The first approach, which utilizes a pixelized image of the pore space, is an extension to unconsolidated systems of a medial-axis based approach (MA). The second approach uses a modified Delaunay tessellation (MDT) of the grain locations. The two algorithms are validated using theoretical packings with known properties and then the networks generated from random packing are compared. For the regular packings, both methods are able to provide the correct pore network structure, including the number, size, and location of inscribed pore bodies, the number, size, and location of inscribed pore throats, and the connectivity. Despite the good agreement for the regular packings, there were differences in both the spatial mapping and statistical distributions in network properties for the random packings. The discrepancies are attributed to the pixelization at low resolution, non-uniqueness of the inscribed pore-body locations, and differences in merging processes used in the algorithms, and serve to highlight the difficulty in creating a unique network from a complex, continuum pore space.

Abbreviations: MA, medial-axis based approach • MDT, modified Delaunay tessellation


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
BACKGROUND
 Comparison between Medial Axis...
 CONCLUSIONS
 REFERENCES
 
WATER AND SOLUTE TRANSPORT and the distribution of phases in subsurface systems are directly related to the geometry and the topology of the pore space and can have a strong influence on continuum scale properties. While many processes related to site assessment and remediation are usually studied using a continuum approach, pore-scale modeling is considered a powerful tool to better understand the physical processes involved and to determine macroscale constitutive relationships that can be difficult to obtain experimentally. The most critical part of constructing a network model is the identification of the pore structure. This can be done using either using computer-generated porous media systems or through the use of non-destructive imaging techniques (e.g., synchrotron microtomography). The significant challenge lies in the generation of network structures that are physically representative of the real materials of interest.

In this paper we present the development of and compare two improved techniques for generating physically realistic network models of granular porous media. First, we discuss network models and techniques for determining the pore-scale properties and characteristics; second, we describe the regular and random porous media systems used to validate and evaluate the techniques we are proposing; third, we present the improvements made to two network generation techniques based on their known limitations; fourth, we validate these algorithms using regular packings; and finally, we evaluate and compare the network representations of random packings generated by these techniques.


   BACKGROUND
TOP
 ABSTRACT
 INTRODUCTION
 BACKGROUND
Comparison between Medial Axis...
 CONCLUSIONS
 REFERENCES
 
Early models developed to account for pore-scale properties idealized the pore space as collections of capillary tubes and provided simple analytical solutions to predict continuum-scale properties such as permeability. However, these models failed to incorporate the interconnectivity of the pore space. Thus the idea of representing the pore space as a two- or three-dimensional network emerged from the pioneering work by Fatt (1956a)(b,c). Due to the complexity of the pore-space morphology, the pore bodies and throats are usually represented by simplified shapes. Pore bodies have been represented by spheres or cubes, while pore throats have been represented by cylinders or other ducts with non-circular cross-sections. Although network models can be two- or three- dimensional, two-dimensional networks have significant limitations in how they represent the interconnectivity of real three-dimensional media (Chatzis and Dullien, 1977). Common three-dimensional networks are random (e.g., Lowry and Miller, 1995) or cubic lattices with coordination number of six (or smaller if bonds are removed). In actual porous media systems, however, the connectivity is distributed, and can have values larger than six (Kwiecien et al., 1990). One solution to this problem is to use physically representative networks, which do not employ an underlying lattice but instead, map the true pore structure of a medium onto a network (Bryant et al., 1993a). Consequently, they retain more completely the true pore morphology and any inherent spatial correlations.

Network models have been used to study a wide range of single and multiphase flow processes including relative permeability (Blunt and King, 1990, Rajaram et al., 1997, Fischer and Celia, 1999); the effect of pore structure on relative permeability and capillary pressure hysteresis in two phase systems (Jerauld and Salter, 1990); prediction of permeability (Bryant et al., 1993a); investigation of the functional relationship between capillary pressure, saturation, and interfacial areas (Reeves and Celia, 1996); prediction of permeabilities and residence time distributions for mechanical dispersion in packed beds (Thompson and Fogler, 1997); drainage and imbibition (Lowry and Miller, 1995, Hilpert and Miller, 2001); and phase distributions, interfacial areas, and mass transfer (Dillard and Blunt, 2000).

The most critical part of constructing a network model is the identification of the pore-structure. Two general approaches are commonly used. The first attempts to directly map a specific porous medium onto a network structure. The second attempts to create an equivalent network using distributions of basic morphologic parameters. The fundamental difference between the two methods is that direct mappings provide a one-to-one spatial correspondence between the porous medium structure and the network structure, whereas the second type of network is equivalent only in a statistical sense.

Typical parameters used for characterization include pore-body and pore-throat size distributions, throat-length distribution, coordination number, spatial correlation within the system, which includes size correlation between adjacent pore bodies (body-body correlation), size correlation between adjacent pore throats (throat-throat correlation), and size correlation between neighboring pore throats and pore bodies (body-throat correlation), and pore-body/pore-throat aspect ratio. There are different approaches used to identify the pore structure including the use of inverse methods to fit model parameters to an experimentally measured pressure saturation relation and the use of mercury porosimetry to obtain the pore-size distribution of the system.

Direct mapping techniques require some method of obtaining the pore structure before reconstruction of the network. Imaging techniques include the analysis of serial thin sections (e.g., Cousin et al., 1996, Vogel and Roth, 2001) and non-destructive techniques such as X-ray microtomography (e.g., Rintoul et al., 1996, Lindquist et al., 2000). Non-destructive direct imaging approaches (i.e., tomography) are attractive because they provide a detailed and unique description of the pore-space geometry. Equally important is the fact that the connectivity and spatial variation in inscribed pore body and throat sizes are retained. Recent advances in microtomography have allowed researchers to obtain three-dimensional images of porous media systems on the order of 10 µm. The acquisition of these extremely rich data sets has motivated the development of algorithms designed to extract the pore network parameters.

A complementary tool to microtomography is computer reconstruction of various types of porous media (Bakke and Oren, 1997, Nandakumar et al., 1999, Liang et al., 2000). While microtomography is important from an analytical perspective (especially for types of porous media that are not well characterized), the use of computer-simulated media allows one to avoid this experimentally intensive step for packed beds, well-sorted sands, or certain sandstones. Additionally, the use of computer-generated porous media is preferable to directly constructing networks from statistical information because one can work with more fundamental material properties (e.g., grain-size distributions rather then inferred properties of the pore space).

With either approach (microtomography images or computer-generated porous media), a significant challenge lies in the generation of network structures that are physically representative of the real materials of interest. The objective of this paper is to present and compare two different algorithms to extract pore network parameters from three-dimensional images of porous media systems. The first approach is based on a medial-axis analysis of a pixelized image of the pore space. The second method is based on a modified Delaunay tessellation of the grain locations. In this work, the two algorithms are validated using theoretical packings with known properties (e.g., simple cubic and rhombohedral packings). Then, networks generated from random packings are compared. Because various approaches to network generation lead to a variety of structures, it is hoped that this comparison of independent techniques will improve existing algorithms for network reconstruction, will validate the methods, and will give new insight into the fundamental structure of common porous materials.

Porous Media Systems
For this work, we have chosen to use computer-simulated sphere packings for several reasons: (i) to study the MDT; (ii) to compare the MDT (which allows for certain analytical calculations) to the MA-based approach (which relies on the digitized image) and (iii) to allow the creation of digitized images of the system at arbitrarily high resolutions.

The two techniques are validated using cubic and rhombohedral packings. The first packing used has a simple cubic structure having uniform, equally spaced pore bodies and throats. The second packing has a rhombohedral structure, which contains two distinct pore-body sizes that are repeated throughout the regular structure (Kruyer, 1958). The third packing has the same rhombohedral structure, but pixelized at a higher resolution.

  1. The cubic packing consists of 512 spheres, digitized at 25 pixels per sphere diameter. The size of the three-dimensional image is 200 by 200 by 200 pixels. The calculations are based on an interior volume of 175 by 175 by 175 pixels, where a sphere radius is cut out from each side to remove the boundary effects. For this system, there is one pore per sphere, all pores have inscribed sphere radii of 0.7320R and inscribed throat radius of 0.4142R, where R is the sphere radius (Kruyer, 1958). The coordination number is exactly six.
  2. The first rhombohedral packing consists of 512 spheres, digitized at 25 pixels per sphere diameter. The size of the three-dimensional image is 200 by 173 by 163 pixels. The calculations are based on a volume of 175 by 148 by 138 where a sphere radius is cut out from each side to remove the boundary effects. Analytical calculations by Kruyer (1958) indicate two pore sizes for the rhombohedral packing: a smaller one with radius 0.2247R and a larger one with radius 0.4142R; all pore throats have inscribed radii equal to 0.1547R.
  3. A second rhombohedral packing was created at a higher pixel resolution for use with the MA algorithms. (Resolution does not affect the MDT algorithm since it works with the sphere locations directly rather than with digitized packings.) It consists of 512 spheres, digitized at 45 pixels per sphere diameter. The size of the three-dimensional image is 364 by 315 by 297 pixels. The calculations are based on a volume of 319 by 270 by 252 where a sphere radius is cut out from each side to remove the boundary effects. For a granular media system with a d50 of 450 µm, the resolution of this system would be on the order of 10 µm. This is approximately the limit when applying synchrotron microtomography to typical porous media systems. The radius of the spheres in this packing is 22.5 pixels.

Following validation of the algorithms using the regular packings, the two methods are compared using a computer-generated random sphere pack. For random packings, which are generally of more interest than regular packings, the network structure created by these algorithms is not unique. Hence, it is of interest to compare the networks that are generated by each, as well as differences in network geometry that are introduced by parameters in each algorithm. In this paper, we focus on results obtained using a random sphere pack containing 1000 uniformly sized spheres, with a porosity of 38%. The pack is contained in a cubic domain, and is periodic in all directions. The use of periodicity prevents edge effects from influencing the tessellation. However, because the medial axis algorithm does not employ periodicity, the networks were cropped during comparison of the two algorithms. The packing was created using a collective rearrangement technique (Liu and Thompson, 2000). For the medial axis calculations, this system is pixelized using a resolution of 25 pixels per sphere diameter. Thus, the original size of the three-dimensional image is 250 by 250 by 250. After cropping one sphere radius from each side, the image is 200 by 200 by 200.

Network Generation Using Medial Axis Based Approach
The medial axis of an object is the skeleton of the object running along its geometrical middle. For a porous medium, it provides simple and compact basic information about the topology and geometry of the void space. In the case of an image with a range of intensity values, the first step in defining the medial axis is to apply a segmentation algorithm to separate the phases of the image based on the intensity values. A discrete burn algorithm (Kwok, 1988; Sirjani and Cross, 1991; Lindquist and Lee, 1996) is applied to the segmented binary image to obtain a discrete representation of the continuum pore space.

The medial-axis skeleton can be obtained by inverting the image so that the background becomes the new foreground and vice versa and then computing the distance transform of the inverted image. The local maximum of the resulting image will represent points on the object's skeleton. Each skeleton is equally distant from at least two different points on the object border. The border points closest to an interior point are called the projection points of the interior point and can be obtained from a distance transform (Goldak et al., 1991; Brandt, 1992; Lohman, 1998). Figure 1 shows the concept of the medial axis in a simple, two–dimensional cross-section.



View larger version (24K):
[in this window]
[in a new window]
  		Fig. 1. The medial axis concept in two dimensions; circles represent the solid phase and the dashed lines represent the medial axis.

 
Baldwin et al. (1996) developed a morphological thinning algorithm to obtain a representation of the skeleton (medial axis) of the pore space and partition the void space into corresponding pores and throats. No attempt was made to use the medial axis itself as tool to characterize pores and throats; the thinned images were used instead. In their work, the pores are identified by summation of the voxels bounded by the solid matrix and the locations corresponding to local minima are defined as pore necks. The thinning algorithm is based on a segmentation algorithm that leads to misidentification or less accurate partitioning of the void space. Lindquist et al. (2000) developed algorithms to calculate effective pore body and pore throats radii and other geometric properties such as coordination number based on the medial axis of the pore space. Their algorithms primarily have been applied to consolidated systems, and it is not clear whether they are applicable to higher porosity materials (e.g., unconsolidated media).

Medial Axis Generation
The 3DMA software package (http://www.ams.sunysb.edu/~lindquis/3dma/3dma.html; verified 13 June 2003) is used to construct the medial axis on the binary images. The 3DMA medial axis algorithm utilizes the technique developed by Lee et al. (1994) and modified by Lindquist and Lee (1996). A burning algorithm is applied as follows: each voxel in the solid phase is labeled with integer 0, voxels in the exterior region (beyond the boundary of the image) are labeled with integer -1, and voxels in the void space are initially unlabeled. The algorithm proceeds by simultaneously assigning unlabeled voxels in the void space with an integer n+1 for each voxel neighboring a grain-phase voxel. This procedure is iterated at a rate of one layer per iteration until all voxels in the pore space are labeled. There are two methods to define the neighborhood of a voxel. The first is face touching, where each voxel can be neighbored by a maximum of six voxels. The second is boundary touching, where each voxel can be neighbored by a maximum of 26 voxels. For complex three-dimensional systems, the 26-connectivity approach is preferred. Determination of the medial-axis voxels depends on the burning direction; if the burn enters a voxel from two different directions that are not perpendicular, the voxel is considered a medial-axis voxel. If the burn enters two neighboring voxels from different directions, both voxels are considered as medial axis voxels. The result of this burning algorithm is a set of media1 axis voxels each with a unique burn number, which is related to its distance from the closest grain voxels.

The medial axis in its initial form is made up of a network of one-dimensional paths and nodes at path intersections. The spatial location and the burn number of each voxel in the medial axis is recorded and is used as the basis from which to calculate the locations and sizes of inscribed pore body and pore throats as well as the coordination numbers. It should be noted that the medial axis is very sensitive to the noise and the boundary of the solid phase.

Inscribed Pore Radius Calculations
To provide a unique measure of pore location and size, we have chosen to focus on calculating the inscribed pore body radii (i.e., the radius of the largest inscribed sphere in a pore). A node is defined as an intersection of three or more paths on the medial axis. Because of the complex, three-dimensional nature of the medial axis, a node might consist of more than one voxel having the same maximum burn number. In this case, one voxel is chosen randomly as the center of the node. A dilation algorithm is performed to find the radius of the pore as follows: (i) start from the chosen voxel on the node and create a sphere with a radius of one voxel; (ii) insert this sphere in the original segmented data and check if there is any overlap between the inscribed sphere and the surrounding grains. If there is an overlap, the procedure stops and this radius is considered as the radius of the pore. If there is no overlap, the algorithm proceeds by increasing the radius by one voxel (i.e., dilation of the sphere radius), and repeats the previous step until an overlap between the grain-phase and sphere occurs.

In some cases, two or more nodes (and hence pore bodies) occupy what can be considered one void space. Physically, this might not be appropriate because the idea of the pore network is to provide a unique and simple representation of the void space. Therefore, a merging criterion is adopted in which two inscribed pore bodies are merged if they overlap. When merging is implemented, the largest inscribed pore body is chosen to be the inscribed pore body of that void space (Fig. 2) . It is noted, however, that this approach may underrepresent the inscribed pore body sizes for the system.



View larger version (56K):
[in this window]
[in a new window]
  		Fig. 2. Example showing the merging of two inscribed pore bodies that occupy the same void space; if there is any overlap, the largest pore is chosen.

 
Inscribed Throat Radius Calculations
The same dilation algorithm for calculating the pore radii is used to calculate the throat radii. The center of the inscribed-throat sphere is located at the voxel with the minimum burn number along the path connecting two nodes. Because a path is a set of medial-axis voxels that each have exactly two neighboring voxels (on the medial axis) and there might be multiple occurrences of the minimum burn number along the path, the dilation algorithm described above is applied to all the minimum voxels along the path. The smallest radius found is considered as the inscribed radius. Figure 3b presents a schematic of this process.



View larger version (29K):
[in this window]
[in a new window]
  		Fig. 3. (a) The inscribed pore radius-finding algorithm where the burn number is local maxima, and (b) the inscribed throat radius where the burn number is a local minima. Shaded circles represent the solid phase, dashed lines represent the medial axis, and the circles represent the inscribed pore body and throats.

 
Figure 4 graphically illustrates the overall process of creating inscribed pore body and throat data from the medial axis. Figure 4a is a periodic random packing of uniform diameter spheres; Fig. 4b is the medial axis generated from that system and Fig. 4c shows the inscribed pores, with the radius of the sphere equal to the inscribed pore radius.



View larger version (43K):
[in this window]
[in a new window]
  		Fig. 4. (a) Small random packing of uniform spheres, (b) medial axis generated on this packing, and (c) location of inscribed pore bodies on the medial axis.

 
Validation of the Medial Axis-Based Algorithm
Figure 5 shows portions of three-dimensional images of the first two test systems and the location of the pores in these systems, calculated from the medial-axis-based method. (The lighter spheres are the packing; darker spheres are the inscribed spheres denoting pore locations). Table 1 shows the values calculated for the three test systems along with the analytical and MDT-based values (which will be discussed later). For all three systems, the MA algorithm predicted the exact number of pores and for a majority of the pores, the correct coordination number. For the cubic packing, the MA algorithm slightly underpredicted both inscribed pore-body and pore-throat sizes. The values obtained can easily be explained since the system has a finite resolution and the medial axis is not necessarily in the exact center of the pore body and/or throat. Examination of the inscribed pore bodies and throats calculated for the two rhombohedral packings show that the MA algorithm does underpredict some of the throat/body sizes, but that the prediction gets better as the resolution increases. It is to be noted that in the rhombohedral packing with high resolution (Table 1) there are relatively few inscribed pore throats (21 throats) with a radius that is the absolute minimum (i.e., one pixel). These throats should belong to the group of throats that have a radius somewhere between one or two pixels; the error ultimately stems from the pixelization of the original spheres. This observation illustrates one of the significant problems in dealing with a digitized image of a continuum pore space.



View larger version (87K):
[in this window]
[in a new window]
  		Fig. 5. (a) 512–sphere cube packing, (b) location of the inscribed pore bodies in the cubic packing, (c) 512-sphere rhombohedral packing, and (d) locations of the inscribed pore bodies in the rhombohedral packing.

 

View this table:
[in this window]
[in a new window]
  		Table 1. Networks of regular packings as obtained from theoretical calculations (TC), medial axis (MA), and modified Delaunay Tessellation (MDT) based approaches; note that all distances are non-dimensionalized by the sphere diameter (in voxels).

 
Delaunay-Tessellation-Based Approaches
The Delaunay Tessellation
The Voronoi diagram and its dual, the Delaunay tessellation, are techniques for spatial analysis that are used for numerous applications (Okabe, 2000), including the characterization of porous media. Voronoi polyhedra have long been used to analyze local structure in random packings (e.g., Finney, 1970, 1976). Mellor (1989) argued that the simplical cells formed by the Delaunay tessellation are more appropriate than Voronoi cells for characterizing random packings, and this approach has been used widely since. In his original work, Mellor used the Delaunay tessellation to measure geometric properties of the Finney packing (Finney, 1968). He then extended this approach to create a topologic model of the interconnected pore structure, which was used to model capillary-dominated immiscible displacement (Mellor, 1989; Mason and Mellor, 1995).

The reason that the simplical cells in the Delaunay tessellation are a good basis for characterizing packings is their geometry. Specifically, if the centers of the spheres in the packing are used to perform the tessellation, then the largest voids tend to be found in the interior of tetrahedrons, while the tightest constrictions in the packing are necessarily projected onto the faces of the tetrahedrons. Thus, each tetrahedron represents a pore with four throats emanating outward across the four faces (see Fig. 2.17–2.19 from Mellor [1989] or Fig. 5 from Bryant et al. [1993a]). These considerations provide the rationale for characterizing pore and pore-throat geometry from the tessellation. An added benefit is that the tessellation uniquely defines the connectivity of the system because each tetrahedron shares its four faces with four neighboring tetrahedrons in the space-filling array. Hence, the tessellation can be used to map out the interconnected network structure of the packing.

Following Mellor's work, the Delaunay tessellation has been used in a number of applications for pore-space characterization. Vrettos et al. (1989) used it as a basis for modeling various transport processes in a random packing. However, the use of a parking technique to generate their computer-simulated packing resulted in a porosity of nearly 0.6, which is physically unrealistic. Bryant et al. (1993 a, b) again used the tessellation of the Finny packing as a model porous material, but also included compaction and grain growth so as to generate structures more representative of consolidated materials. They performed dynamic modeling, showing quantitative agreement between predicted and measured permeabilities in sandstones, and they demonstrated the effect that pore-scale spatial correlations in the packing have on transport properties. Thompson and Fogler (1997) applied the Delaunay tessellation to computer-simulating sphere packings to create network models of heterogeneous structures, and Fredd and Fogler (1998) used a similar approach to model reactive flows in which the packing structure changes over time.

Limitations of the Delaunay Tessellation
The direct use of the Delaunay tessellation for mapping pore structure has two significant problems. First, the tessellation restricts the network to a fixed coordination number (equal to exactly four in three dimensions). In reality, experimental analyses show that disordered porous media contain a distribution of coordination numbers with a mean value somewhat larger than z = 4. For packed beds of spheres, Yanuka et al. (1986) used serial sectioning to obtain an average coordination number of z = 4.3. Jerauld and Salter (1990) suggest that this value should be amended to z = 5.3 by discarding individual-pore coordination numbers equal to either 1 or 2. The impact of this limitation (i.e., a fixed coordination number of z = 4) depends on the process of interest. For the modeling of permeability or other single-phase processes, its impact may be relatively small. However, it is well known that the coordination number strongly affects percolation properties in a network (Jerauld et al., 1984a, b), which would suggest that altering a network's coordination number would affect the modeling of multiphase flow (Larson et al., 1981; Bryant and Blunt, 1992), the flow of yield-stress fluids (Rossen and Gauglitz, 1990), etc.

The second problem associated with the tessellation, which is of equal or greater concern, is the incorrect identification of pore locations and sizes. Specifically, the Delaunay tessellation tends to subdivide single voids into multiple pores (where these larger voids are defined either intuitively or by other geometric arguments). This problem is most easily conceptualized using a regular packing as illustrated using Fig. 6 . Figure 6a is a schematic of a two-dimensional packing on a nearly square lattice. Clearly, each of the open spaces centered between four particles should be identified as a discrete pore. Similarly, the constrictions between pairs of particles should represent pore throats, thus giving a fixed coordination number equal to four for this particular structure. This intuitive network is shown in Fig. 6b and 6c. In contrast, parts d through f of Fig. 6 illustrate a hypothetical network structure derived from a Delaunay tessellation. In two dimensions, simplical cells in a Delaunay tessellation are triangles, which forces the packing structure to be discretized as shown in Fig. 6d. Using the analogy between the tessellation and the pore structure described above, pores are placed into each of the Delaunay triangles, while the common faces (i.e., sides of the triangles in two dimensions) define the pore throats. Though the sides of all triangles in the tessellation do indeed represent local minima in particle-particle spacing, it can be argued that the diagonal sides are inappropriate choices for pore throats. This point is further emphasized by comparing the two network structures (Fig. 6c and 6f), which clearly share similar overall topology but differ with respect to important details: the Delaunay network contains twice as many pores, it has a fixed coordination number of 3 instead of 4, and it contains pore throats where there should be none. Finally, if the sum of the volumes of all network pores is made to equal the total void volume (which it should for a physically representative network), then individual pores created by the Delaunay tessellation must be split into smaller volumes. This limitation is important because the network statistics will not be representative of the true void sizes in the packing.



View larger version (52K):
[in this window]
[in a new window]
  		Fig. 6. Schematic of proper pore identification by merging simplical cells of the Delaunay tessellation.

 
In three dimensions, an analogous problem occurs in simple cubic packings: pores should clearly be defined as the relatively large voids that exist in the center of each group of eight spheres in the lattice, as shown earlier in Fig. 6b. However, a tessellation will divide each cubic pore into either five or six Delaunay tetrahedrons. The consequences are analogous to what was described in two dimensions: too many pores, a coordination number of 4 instead of 6, and the creation of physically unreasonable pore throats. While similar arguments can be made for other regular packings (e.g., the rhombohedral structures discussed in this paper) definitive arguments regarding random packings are harder to make. However, using visualization, geometric analysis, and coordination-number data, there is certainly evidence that similar problems exist when Delaunay tessellations are used to discretize random packings.

Modified Delaunay Tessellations
Despite the problems mentioned above, the Delaunay tessellation is a very powerful tool for pore-space characterization, and is therefore chosen as a basis for the MDT technique developed as part of this research. In particular, the following important attributes are retained: uniqueness of the spatial map; computational efficiency (Thompson, 2002); ease of defining the network interconnectedness.

We propose a conceptually simple solution to the problems illustrated above, which is to allow for the merging of multiple tetrahedrons into larger polyhedrons in cases where a single void space has been unnecessarily subivided. The major challenge to implementing this approach is the proper identification of sets of tetrahedrons that should be merged.

The algorithm described here relies on the examination of inscribed spheres (contained within the void space) in an effort to identify multiple tetrahedrons that comprise a single pore. The logic for this approach is illustrated in Fig. 7 . Part a of the figure is a schematic of a nearly square pore in two dimensions, divided into two triangles. The inscribed circles are overlapping to a large extent due to the local geometry of the particles. (Note that the inscribed circles are allowed to extend beyond the boundaries of the triangle itself.) Figures 7b and 7c are similar, but the geometry is increasingly elongated, thereby decreasing the amount that the inscribed circles overlap. Once the location of the (largest) inscribed triangles (or tetrahedrons in three dimensions) are known, they can be used to determine whether simplical cells should be merged, based on the amount that they overlap. In the first schematic, intuitive arguments would suggest that the two cells should be merged to form a single quadrilateral pore. In the last case, the pores should clearly remain separate. Arguments for the middle case are more nebulous, emphasizing the need to define a merging criterion based on the degree of overlap. The appropriate choice for this criterion is discussed later.



View larger version (28K):
[in this window]
[in a new window]
  		Fig. 7. Schematic of three potential overlap conditions encountered during pore merging.

 
A related point should be made regarding how the inscribed spheres are defined. In the two-dimensional schematics (Fig. 7), the inscribed circles are uniquely defined as those that contact the set of three vertex circles. In three dimensions, tests showed that this approach works well for packings of uniform-diameter spheres, and is ideal from a performance standpoint because the location and radius of the inscribed spheres can be found analytically by solution of a set of four nonlinear algebraic equations. However, as spheres in the packing become distributed in size, this approach breaks down because of the existence of relatively flat tetrahedrons, which leads to the formation of very large inscribed spheres. This problem is shown schematically (in two dimensions) in Fig. 8 . The unshaded circle is the unique inscribed circle that touches the three gray vertex circles for the triangle on the right. In the three-dimensional packings, the analogous situation generates large spheres that do not lie entirely within the pore space and therefore do not effectively characterize the void structure. Furthermore, due to their large size, they can envelop any number of neighboring spheres in the packing, and the subsequent merging of cells generates unrealistically large pores that have no near neighbors.



View larger version (63K):
[in this window]
[in a new window]
  		Fig. 8. Schematic of unrealistically large inscribed radii. This condition can occur when only the simplical-cell spheres are used for the inscribed sphere calculations.

 
Consequently, the algorithm consists of a local optimization to find the largest inscribed sphere that can fit within the void space immediately adjacent to a tetrahedron. The optimization is performed using a three-dimensional simplex search, modified according to the algorithm of Nelder and Mead (1965): seed points to begin the searches are the centers of mass of the tetrahedrons, and the optimization routine locates the inscribed sphere. Using this approach, pores are defined using a more physically meaningful criterion (equivalent to the definition used by Nolan and Kavanagh [1994]) because it ensures that the inscribed spheres lie entirely within the void phase of the packing. An added benefit in the context of the Delaunay tessellation is that convergence of the local optimization disassociates the resulting inscribed sphere from any particular tetrahedron in the tessellation. This latter point is illustrated in Fig. 9 , which is a schematic of a two-dimensional pore defined by six particles and tessellated into four triangles. For this particular geometry, the optimization could be initiated from any one of the four triangles, but each search would converge on the single inscribed circle shown. An added benefit of this approach is that the final locations of merged pores necessarily correspond to local maxima (in inscribed sphere size), which is not true if pores are tied to tetrahedrons.



View larger version (71K):
[in this window]
[in a new window]
  		Fig. 9. Schematic of single inscribed sphere in a multi-tetrahedron pore. The single sphere is identified regardless of which triangle is used to seed the local optimization.

 
Unfortunately, because many real pore structures contain multiple local maxima in the same void, the critical overlap for merging must still be defined somewhat arbitrarily. A reasonable criterion seems to be whether the center of an inscribed sphere lies within a neighboring inscribed sphere. This choice is discussed briefly below, and probably warrants additional research. Once this criterion is imposed, multiple tetrahedrons become merged together, and the resulting pores are polyhedrons with essentially arbitrary shapes. This modified Delaunay tessellation provides for a better one-to-one correspondence between discrete network pores and the true voids within the packing, and it also allows for a variable coordination number. At the same time, the important attributes of the Delaunay tessellation that were described above are retained.

Figure 10 illustrates the overall process for the periodic 1000-sphere random packing used in this study. The figure shows the packing, the Delaunay tessellation, and the network structure labeled ‘MDT1’ in later results.



View larger version (39K):
[in this window]
[in a new window]
  		Fig. 10. Schematic of network generation algorithm using the modified Delaunay tessellation algorithm.

 
Validation using Regular Packings
As mentioned above, the Delaunay tessellation incorrectly discretizes most regular packings because of the geometric limitations of the simplical cells. Therefore, one important validation of the modified Delaunay tessellation is its ability to produce correct morphologic data for regular packings. As an aside, we note that very small random perturbations (approximately 10-6 x Dsphere in magnitude) were imposed on the spheres in the regular packings because the Delaunay tessellation is unique only for random structures (Bowyer, 1981).

First, a tessellation was performed on the cubic packing, which generated five tetrahedrons per sphere. Subsequent application of the merging algorithm correctly identified that these groups of five should be merged, and that all parameters in the resulting network structure were exactly correct: the number of pores (one pore per sphere), the coordination number (six), and the pore and throat sizes. Computer-aided visualization was also performed, as illustrated in Fig. 11 .



View larger version (51K):
[in this window]
[in a new window]
  		Fig. 11. Visualization of the cubic packing and sample of a pore and its connectivity as obtained using the modified Delaunay tessellation method.

 
Second, a tessellation was performed on the rhombohedral packing, which generated six tetrahedrons per sphere. Subsequent application of the merging algorithm identified two different pore structures, each of which repeats through the regular packing, as shown in Fig. 12 using visualization. The two pore sizes were in exact agreement with the theoretical values given in Table 1 (within the numerical accuracy specified during the run). Further analysis showed that 2/3 of the pores in the packing are four-sphere clusters (which have the smaller value of the inscribed radius), while 1/3 of the pores are six-sphere clusters having the larger radius. All of the throats in the rhombohedral packing are equal in size, with inscribed radii equal to 0.156R, again in agreement with Kruyer (1958). Coordination numbers for the two types of pores are 4 and 8 respectively, differing significantly from the coordination numbers reported by Yanuka et al. (1986). We believe the values presented here are correct, and that simpler empirical or statistical analyses may not be appropriate for all types of pore structures, especially those with extreme porosity values such as the rhombohedral packing (which has the lowest possible porosity of all monodisperse unconsolidated packings).



View larger version (52K):
[in this window]
[in a new window]
  		Fig. 12. Visualization of the rhombohedral packing and sample of the two distinct pores and their connectivity as obtained using the modified Delaunay tessellation method.

 
The Critical Merging Length in the Modified Delaunay Tessellation Algorithm
For regular packings, the network structure created by the modified Delaunay tessellation is independent of the merging criterion chosen (for reasons illustrated in Fig. 9). However, for random packings, the critical merging length has a significant effect on the network structure. While the merging length can be varied over a large range, we compare only three networks for the sake of brevity, which are denoted as follows: (i) DT, a network created directly from the Delaunay tessellation; MDT1, a modified Delaunay network; the merging criterion was set so that two tetrahedrons were merged if the center of one inscribed sphere was inside the other; MDT2, a modified Delaunay network, the merging criterion was set so that two tetrahedrons were merged if the respective inscribed spheres contacted one another. (To help illustrate the differences between MDT1 and MDT2, refer to Fig. 7: The cells in Fig. 7a would be merged in both of the modified networks. The cells in Fig. 7b would remain separate using the MDT1 criterion, but would be merged using the MDT2 criterion.).

Table 2 contains morphologic data for the three networks, and illustrates the significant differences between the three networks. Particularly noteworthy are the following points. The number of pores decreases substantially as the pores become merged. While the minimum void volume increases dramatically, the inscribed sphere radii are essentially unchanged (because inscribed spheres are allowed to bulge out of their tetrahedrons for any of the networks). The average coordination number shows a moderate increase, while the maximum coordination number changes from a value of four (for the Delaunay tessellation) to 14 and 16 respectively for the two merged networks. The last row in the table gives the fraction of throats with a high aspect ratio, defined as Rt,ins/Rp,ins. High-aspect-ratio throats are those that exhibit only a slight constriction (e.g., the diagonals in Fig. 6a), and thus reductions in their number is generally desirable.


View this table:
[in this window]
[in a new window]
  		Table 2. Statistics for Delaunay and modified Delaunay networks.

 

   Comparison between Medial Axis and Modified Delaunay Tessellation-Based Approaches
TOP
 ABSTRACT
 INTRODUCTION
 BACKGROUND
 Comparison between Medial Axis...
CONCLUSIONS
 REFERENCES
 
Comparisons of the network generation techniques are made in this section. For all of these comparisons the MDT2 algorithm was used. When comparing MA versus MDT2 networks, we search for cases where the same pore is identified by both algorithms. These are referred to as ‘matched pores’ in the tables below. The criterion for matching is an overlap of the inscribed spheres found by the two respective methods.

Comparison for Regular Packings
As mentioned previously, the comparison of the two approaches is based on two types of packings: regular and random packings. A comparison of the regular packings is provided in Table 1.

For the cubic packing, there is a good agreement between the two methods. The slight differences in the locations and radii of the pores from both methods are due to the fact that the MDT-based parameters are obtained directly from the actual sphere locations and radii, while the MA-based parameters are obtained from a pixelized representation of the system. There is also a good agreement in the case of the rhombohedral packing for both high and low resolutions. Figure 13 shows the distributions of distance between matched pores for the two resolutions. Clearly, as the resolution increases, the distance between the centers of the matched pores decreases. In the medial-axis-based method, the radius of a pore is very sensitive to the calculated location of that pore because the radius is obtained by a dilation algorithm starting from the center of the pore (i.e., the node location) to the closest surrounding solid space. As mentioned earlier, a node in the medial axis can be made up of several voxels in which the same local maximum was identified. Therefore, the issue of which voxel to choose as the pore center can lead to errors in calculating the radius of the inscribed pore body at that location. A point worth noting is that the minimum throat size predicted by the MA algorithm decreases with increasing resolution. This effect is because of a relatively small number of inscribed pore throats (21) with a radius of one pixel at the higher resolution. In reality, these throats should belong to the group of throats that have a radius of two pixels. The error ultimately stems from the pixelization of the original spheres. This observation illustrates one of the significant problems in dealing with a digitized image of a continuum pore space. Also note however that the average throat radius becomes closer to the true value at the higher resolution, indicating that higher resolutions do indeed improve the description in an average sense.



View larger version (26K):
[in this window]
[in a new window]
  		Fig. 13. Histogram of distance between centers of matched pores normalized by the sphere diameter. (top) Twenty-five pixel diameter sphere resolution and (bottom) 45-pixel diameter sphere resolution.

 
Comparison for Random Packing
Figure 14 shows three-dimensional cutouts of the packing and the inscribed pore bodies obtained from both network generation methods. Table 3 shows the pore network parameters obtained from both methods. We will now discuss several particular examples to highlight our findings. In Fig. 15 and 16 , which are two-dimensional slices through representative pores, the gray spheres are the spheres, the dashed lines are the inscribed pore bodies obtained from the MDT2 algorithm, and the solid lines are the inscribed pore bodies obtained from the MA algorithm. (Note that examining a two-dimensional slice [i] makes the grains appear to have different radii; [ii] may not indicate the true inscribed pore body radii; and [iii] may not show pore body/grain contacts.) In the cases shown in Fig. 15, the inscribed pore bodies are considered to be a good match because of their overlap. Two different situations are illustrated in Fig. 15. On the left is a case where both methods define two pores within the void space shown (pores from each method did not get merged because there is no overlap). On the right, both methods define only a single pore in the void space. On the other hand, Fig. 16 illustrates a situation where there is no match between the inscribed pore bodies defined by the respective methods. While each method does define a pore in the same general void space, the fact that there is no overlap between pores from both methods classifies this case as non-match. This situation explains two features appearing in Table 3. First, approximately 30% of the pores are not matched. Second, a discrepancy exists in the fraction of the void space occupied by inscribed pore bodies when comparing the two methods.



View larger version (55K):
[in this window]
[in a new window]
  		Fig. 14. Uniform random packing system with the dark spheres representing the locations of inscribed pore bodies from the medial axis (MA) and modified Delaunay tessellation (MDT) method. Light color represents grains; dark color represents inscribed pore bodies.

 

View this table:
[in this window]
[in a new window]
  		Table 3. Comparison between Medial Axis (MA) and Modified Delaunay Tessellation (MDT) based approaches on random packing (sphere diameter = 25 pixels).

 


View larger version (20K):
[in this window]
[in a new window]
  		Fig. 15. Two examples showing a match between inscribed pore bodies identified from MA (solid line) and MDT (dashed line) methods.

 


View larger version (79K):
[in this window]
[in a new window]
  		Fig. 16. An example where there is no match between inscribed pore bodies from MA (solid line) and MDT (dashed line) methods.

 
Figure 17 shows the distribution of distances between the centers of matched inscribed pore bodies from the two methods. The largest distances are associated with matches between large pores. The number of inscribed pore bodies from both methods that are separated by a distance larger than 0.2 sphere diameters is 156, or approximately 15%.



View larger version (36K):
[in this window]
[in a new window]
  		Fig. 17. Distribution of the distance between the centers of matched inscribed pore bodies from the two methods for random sphere packing.

 
Figure 18 shows a comparison of the inscribed pore body and throat distributions from the two methods. As expected, the MA method generates inscribed pore bodies with smaller radii than the MDT-based method. Because of this shift in the pore-size distribution, we would expect resulting differences in certain simulations and particularly in imbibition simulations that make use of this parameter to estimate radii of curvature for fluid interfaces.



View larger version (40K):
[in this window]
[in a new window]
  		Fig. 18. Inscribed (top) pore body and (bottom) throat-size distributions (right) MDT and (left) MA for the random sphere packing.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 BACKGROUND
 Comparison between Medial Axis...
 CONCLUSIONS
REFERENCES
 
In this paper, we have presented and compared two different methods for calculating the locations and sizes of inscribed pores and pore throats in three-dimensional porous media systems. Both methods involve significant modifications to existing techniques. The first technique, based on a medial axis, relies on discretization of the three-dimensional image and can be applied to any three-dimensional porous media system. The Delaunay tessellation-based method is traditionally used with sphere packings (in which particle locations are known), but it could easily be used in conjunction with tomographic images from which the centroids of particles had been determined. While the packings analyzed in this work are limited in size, they are representative of systems scanned using high-resolution microtomography; no inherent restrictions on system size exist within the algorithms, within the constraints of available computer memory and/or computation time.

Analysis of the medial axes created from the regular packings showed that, even in these ideal cases, reconstruction of the skeleton is sensitive to the resolution of the image. In addition, the non-uniqueness of the node to pore-body correlation can lead to an overestimate of the number of pores. To provide a more physical representation of the inscribed pore body locations, node-merging algorithms were developed. Results also showed that the pixelization effect (e.g., the finite size of the pixels) plays an important role the reconstruction of the medial axis and extraction of pore network parameters.

As shown in the characterization of the regular packings, both methods are able to provide the correct pore network structure, including the number and location of the inscribed pore bodies, the number and location of the inscribed pore throats, and the connectivity. Being analytically based, the MDT method was able to extract the correct pore body and throat sizes. It was found that, because the MA-based method is resolution-dependent and precise voxel locations of nodes are non-unique, the exact inscribed pore body and throat sizes were underestimated. As the resolution of the image increases, the agreement between the two methods improves. It should be noted, that while inscribed pore body-size distributions generated from the two methods show good agreement and trend, no attempt was made to compare the effectiveness of the two methods to properly capture the spatial correlation of the system. Despite the good agreement of the regular packings, there were differences in the both the spatial mapping and the statistical distributions in network properties for random packings. The discrepancies between the two methods can be attributed to several factors: (i) pixelization at finite resolution; (ii) the arbitrary process of choosing the center of a inscribed pore body (i.e., node) from a group of voxels all having the same burn number; and (iii) the differences in the merging processes. Despite these significant differences, both methods appeared to give reasonable networks as interpreted from three-dimensional visualization and/or examination two-dimensional slices. These observations emphasize the difficulty in creating a unique network from a complex, continuum pore space.

An important consideration that is being addressed in a subsequent study is how the differences in network structure affect relevant transport processes. We expect that the differences in pore sizes and pore-size distributions observed here would have tangible effects on multiphase calculations. In this regard, it should be noted that for the medial-axis networks, one reason that the pore sizes are somewhat smaller than for the MDT networks is because the dilation algorithm proceeds in radial steps of unit-voxel length. This part of the medial-axis algorithm can be refined (at considerable computational expense) to provide more accurate values for the inscribed radii.

Received for publication June 6, 2002.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 BACKGROUND
 Comparison between Medial Axis...
 CONCLUSIONS
 REFERENCES
 





This Article
Right arrow Abstract