Computational Techniques

Original Microstructure Generation

All of the techniques described in this paper are applied to digital-image based microstructures. Thus, in three dimensions, a microstructure consists of a three-dimensional grid (lattice) in which each site is defined to be either solid or pore. For the systems investigated here, the lattices were always 100*100*100 units for a total of one million sites (pixels). To minimize finite size effects, periodic boundaries were utilized during microstructure generation. If a portion of a solid object, such as a sphere, extended outward through a face of the three-dimensional box, it was moved to the opposite face of the system.

To evaluate the reconstruction techniques employed in this paper, original microstructures based on a penetrable sphere model were selected. This model was selected because its correlation characteristics can be analytically determined [11] and its transport properties had been studied previously [1]. The microstructures were generated by centering solid digitized spheres at random locations (x,y,z) in the lattice using periodic boundary conditions. The spheres were allowed to freely overlap and porosities were varied between about 40 and 15% using two methods for lowering porosity. In method one, the initial system consisting of 500 spheres of diameter 15 pixels (lattice units) was modified by increasing the diameter of each sphere in increments of two up to a value of 19, without changing the centroids of the spheres. In the second method, the same initial system was modified by placing additional spheres into the system. Systems containing 500, 600, and 900 spheres were utilized. Thus, a total of five different original microstructures were generated and analyzed in this study.

Two-Dimensional Slice Selection and Three-Dimensional Reconstruction

For each of the original three-dimensional microstructures, five two-dimensional slices were selected at random locations in the z-plane. Each two-dimensional slice was used as a generating image in reconstructing a separate three-dimensional representation to be compared to the original three-dimensional microstructure.

The reconstruction technique employed for each two-dimensional image was as follows. First, the autocorrelation function for the image was computed. If the two-dimensional image is defined as a discrete valued function, I(x,y), where I(x,y) is equal to one for solid pixels and zero for pore pixels, the two-point correlation function S(x,y) for the image is given by [12]

$$S(x,y)= \sum_{i=1}^{M} \sum_{j=1}^{N} \frac{I(i,j)*I(i+x,y+j)}{M*N}$$ (1)

where M and N=100 for the images used here and periodic boundaries are used to define I(i+x,j+y) when (i+x,j+y) extends beyond the 100*100 two-dimensional image. S(x,y) was determined for values of x and y ranging from 0 to 60 and was then converted to S(r) for distances r in pixels by the equation [12]:

$$S(r)=\frac{1}{2r+1} \sum_{l=0}^{2r} S(r,\frac{\pi l}{4r})$$ (2)

where S(r, θ)=S(rcos θ,rsin θ) was obtained by bilinear interpolation from the values of S(x,y) determined above.

Similar to the approach used by Quiblier [8], the initial reconstructed three-dimensional image consisted of Gaussian distributed noise generated using a uniform random number generator [13] and the Box-Muller method [14] to convert the uniform random deviates to normal deviates. This three-dimensional noise image, N(x,y,z) was directly filtered (or convolved) with the autocorrelation function, F(x,y,z), defined as [15]

$$F(r)=F(x,y,z)=\frac{[S(r=\sqrt{x^2 + y^2 +z^2})-S(0)*S(0)]}{[S(0)-S(0)*S(0)]}.$$ (3)

The resultant image, R(x,y,z) was calculated as

$$R(x,y,z) = \sum_{i=0}^{30} \sum_{j=0}^{30} \sum_{k=0}^{30} N(x+i,y+j,z+k)*F(i,j,k).$$ (4)

This is a simplification of the approach utilized by Quiblier [8], where a matrix of filtering coefficients is computed by solving a huge system of nonlinear equations. With this simplification, the autocorrelation function of the reconstructed porous medium will only approximate that of the original microstructure, and further modification may be required as outlined in section 2.3 below.

R(x,y,z) was converted to a binary (0-pore or 1-solid) image by a thresholding operation. The threshold limits were determined by sampling all of the values of R(x,y,z) and computing a discrete histogram with 500 cells separating the minimum and maximum values of R(x,y,z). This histogram was analyzed to determine the threshold value needed to match the porosity of the computed binary image to that of the original three-dimensional image (not the two-dimensional slice).¹ The porosity of the original three-dimensional microstructure was determined by simply counting the number of pixels in the 100*100*100 system which were porosity. Based on this threshold value, a binary image exhibiting approximately the same porosity as the original three-dimensional microstructure was produced.

Finally, the three-dimensional autocorrelation function of the binary image was calculated for comparison to the autocorrelation function used in the filtering operation.

Hydraulic Radius Matching via Surface Curvature Modification

The final binary images from the image reconstruction process were further modified using an algorithm originally developed to simulate the sintering of ceramic powders [16,17]. Here, the local surface curvature was determined at each pixel in the three-dimensional image by counting how many pore pixels are present in a sphere of some fixed diameter (typically seven pixels) centered at the pixel. Periodic boundary conditions were once again employed during this curvature computation. This measure has been shown to correspond directly to the local curvature [17]. After the curvatures were determined for all pixels, the binary image was modified by interchanging a fixed number (e.g., 200) of solid pixels of highest curvature with the same number of pore pixels of lowest curvature. All curvature values were then updated and the exchange process was repeated. Because the pore volume remains constant and the surface area decreases, this process has the effect of increasing the hydraulic radius (pore volume/surface area) of the microstructure. In every case investigated in this study, the hydraulic radius of the reconstructed binary image was less than that of the original generating structure. The surface curvature modification (exchange) algorithm was used to increase the hydraulic radius of the reconstructed image until it matched that of the two-dimensional generating image, which was somewhat different from the original three-dimensional microstructure due to sampling effects. The effects of this modification were evaluated qualitatively by viewing the original, reconstructed, and modified microstructures and quantitatively by evaluating their transport properties.

Permeability Computation

The permeability computation has been described in detail elsewhere [1,18] but will be briefly outlined here. To compute permeability, the linear Stokes equations are solved for the case of slow incompressible flow. These equations are of the form:

$$\eta \nabla ^2 {\bf v}(r) = \nabla p(r)$$ (5)

and

$$\nabla \cdot {\bf v}(r)=0$$ (6)

where v and p are the local velocity and pressure fields respectively and η is the fluid viscosity. A pressure difference is prescribed across the three-dimensional microstructure and the above equations solved using a finite difference scheme in conjunction with the artificial compressibility relaxation algorithm [18,19]. The digital-image based microstructure is discretized into a marker-and-cell mesh [19] where pressures are defined at the lattice sites and fluid velocity components are defined along the center of bonds connecting sites. Near the solid-pore interfaces, non-centered difference equations are used to improve the accuracy of the solution and force the fluid velocities to zero at each interface. Once the system has sufficiently relaxed, the permeability, k, of the porous medium is calculated by volume averaging the local fluid velocity and applying the Darcy equation [20]

$$<v> = - \frac{k}{\eta} \frac{\Delta P}{L}$$ (7)

where L is the length of the sample (100 lattice units).

Conductivity Computation

The conductivity calculation, basically a solution of the Laplace equation, has been described in detail elsewhere [18,21] but will be briefly presented here. First, the three-dimensional microstructure is converted into a network of resistors by connecting each pair of adjacent pixels by a resistor. A resistor connecting two pore pixels is assigned a conductance of one while a resistor connecting two solid pixels or a pore and a solid pixel is assigned a conductance of zero. A voltage gradient is then applied across the sample and the system is relaxed using a conjugate gradient algorithm [22] to obtain the voltages at the nodes (center of each pixel). Knowing this voltage distribution allows one to calculate the total current and thus the equivalent relative conductivity, σ / σ₀, for the microstructure. Here, σ is the conductivity of the porous medium and σ₀ is the conductivity of the solution in the pore space (taken to be one in these computations).

Critical Diameter Computation

In addition to hydraulic radius, a critical diameter was computed as a characteristic length scale for each microstructure. The computation of this critical diameter is an extension of an algorithm to simulate mercury intrusion porosimetry in two dimensions [23] and is equivalent to the concept of a critical sphere discussed by Thovert et al. [24].

Basically, the goal is to determine the size of the largest digitized sphere which can proceed (percolate) through the pore space from one side of the microstructure to the other. The test sphere is not allowed to overlap any solid regions of the microstructure. The diameter of the test sphere is increased from a value of three in increments of two pixel units. The critical diameter in pixel units is determined as the average of the largest diameter for which percolation occurs and the smallest diameter for which percolation doesn't occur. For cases where the smallest diameter for which percolation doesn't occur is three, a conventional burning algorithm [25] is used to assess the percolation of the pore space for a one-pixel diameter sphere. This length scale, commonly denoted D_c, has been shown, in combination with the relative conductivity, to be an excellent predictor of permeability for both model [1] and real porous media [26].