Modeling Heterogeneous Materials via Two-Point Correlation Functions: II. Algorithmic Details and Applications
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n Modeling Heterogeneous Materials via Two-Point CorrelationFunctions: II. Algorithmic Details and Applications
Y. Jiao
Department of Mechanical and Aerospace Engineering,Princeton University, Princeton New Jersey 08544, USA
F. H. Stillinger
Department of Chemistry, Princeton University,Princeton New Jersey 08544, USA
S. Torquato ∗ Department of Chemistry, Princeton University,Princeton New Jersey 08544, USAPrinceton Institute for the Science and Technology of Materials,Princeton University, Princeton New Jersey 08544, USAProgram in Applied and Computational Mathematics,Princeton University, Princeton New Jersey 08544, USA andPrinceton Center for Theoretical Physics,Princeton University, Princeton New Jersey 08544, USA (Dated: May 31, 2018) bstract In the first part of this series of two papers, we proposed a theoretical formalism that enables oneto model and categorize heterogeneous materials (media) via two-point correlation functions S andintroduced an efficient heterogeneous-medium (re)construction algorithm called the “lattice-point”algorithm. Here we discuss the algorithmic details of the lattice-point procedure and an algorithmmodification using surface optimization to further speed up the (re)construction process. Theimportance of the error tolerance, which indicates to what accuracy the media are (re)constructed,is also emphasized and discussed. We apply the algorithm to generate three-dimensional digitizedrealizations of a Fontainebleau sandstone and a boron carbide/aluminum composite from the two-dimensional tomographic images of their slices through the materials. To ascertain whether theinformation contained in S is sufficient to capture the salient structural features, we computethe two-point cluster functions of the media, which are superior signatures of the microstructurebecause they incorporate topological connectedness information. We also study the reconstructionof a binary laser-speckle pattern in two dimensions, in which the algorithm fails to reproducethe pattern accurately. We conclude that in general reconstructions using S only work well forheterogeneous materials with single-scale structures. However, two-point information via S isnot sufficient to accurately model multi-scale random media. Moreover, we construct realizationsof hypothetical materials with desired structural characteristics obtained by manipulating theirtwo-point correlation functions. PACS numbers: 05.20.-y, 61.43.-j ∗ Electronic address: [email protected] . INTRODUCTION Random heterogeneous multiphase materials or media are ubiquitous. Examples includecomposites, porous media, biological materials as well as cosmological structures, and theirmacroscopic properties are of great interest [1]. In the first part of this series of two papers[2] (henceforth referred to as paper I), we proposed a theoretical formalism to model andcategorize heterogeneous materials via two-point correlation functions S ( r ), which can beinterpreted as the probability of finding two points separated by the the dsiplacement vector r in one of the phases [1]. In particular, we introduced the idea of the two-point correlationfunction space and its basis functions. In general, S of a medium can be expressed bya map ℘ on the associated basis functions, which is composed of convex-combination andproduct operations. We also suggested a set of basis functions by examining certain knownrealizable analytical two-point correlation functions. Moreover, we introduced an efficientisotropy-preserving S -sampling algorithm, namely the lattice-point algorithm, but left thedetails of the algorithm for another paper.Here we will provide algorithmic details of the lattice-point methodology and considerseveral nontrivial applications to illustrate the practical utility of our theoretical formalism.In particular, we will apply the lattice-point algorithm to generate three-dimensional (3D)digitized realizations of a Fontainebleau sandstone and a boron carbide/aluminum compos-ite from two-dimensional (2D) tomographic images of their slices through the materials.To justify whether the reconstructions are successful, one should also measure other sta-tistical descriptors of the media [3]. Here we compute the two-point cluster function C [4] (see definition and discussions in Sec. IV), which contains nontrivial “connectedness”information of the phases of interest. By comparing C of the target and reconstructedmedia, we demonstrate that S is indeed sufficient to capture the salient structural featuresin these cases. We also study the reconstruction of a binary laser-speckle pattern in twodimensions, and find that the algorithm fails to reproduce the target pattern accurately.We conclude that in general reconstructions using S only work well for heterogeneousmaterials with single-scale structures (those having only one characteristic length scale).However, two-point information via S is not sufficient to accurately model multi-scale me-dia (those composed of structural elements associated with multiple characteristic lengthscales). Moreover, we will show how the proper convex combinations of basis functions en-3ble one to obtain two-point correlation functions with desired properties and thus enableone to generate a variety of structures with controllable morphological features. Note thathere we will mainly focus on modeling heterogeneous materials in three dimensions, whichcomplements the two-dimensional examples that we considered in paper I.Yeong and Torquato formulated the (re)construction problem as an energy-minimizationproblem using simulated annealing [5]. This has become a very popular (re)constructiontechnique [6, 7, 8, 9, 10]. In this method, a nonnegative objective function E , called the“energy,” is defined as the sum of squared differences between the target and sampled cor-relation functions. In our case, the energy is given by E = X i [ S ( r i ) − ˆ S ( r i )] , (1)where ˆ S ( r ) and S ( r ) is the target and sampled two-point correlation function, respectively.An important issue in the (re)constructions based on the method of simulated annealing isthe choice of the energy threshold E th , i.e., the error tolerance of discrepancies betweenthe correlation functions of generated medium and the imposed ones. When the energy ofthe (re)constructed medium is below E th , the (re)construction process is terminated. Theenergy threshold is a key indicator of how accurately the medium is (re)constructed. Inour previous work, E th was always chosen to be a very small number, e.g., 10 − ; however,no quantitative analysis on why such a value should be chosen was given. In this paper,we will discuss in detail the significance of the energy threshold for both the orthogonal S -sampling algorithm (only sampling S along convenient orthogonal directions) [3, 5] andthe lattice-point algorithm.In the lattice-point procedure, we consider the digitized medium (pixel system) to bea “lattice-gas” system, in which the black pixels behave like nonoverlapping “particles”moving from one lattice site to another. The two-point correlation function of the mediumis obtained by binning all the distances between the black pixels and dividing the numberin each bin by the total number of distances between all lattice sites within that bin. Togenerate a trial configuration, a randomly selected “particle” (black pixel) is given a randomdisplacement subjected to the nonoverlapping constraint. Only the distances between themoved “particle” and all the others need to be recomputed in order to obtain S of the trialconfiguration. In this way, all directions in the digitized medium are effectively sampled;4oreover, the complexity of the algorithm is linear in N B , i.e., the number of nonoverlappingparticles.For heterogeneous materials containing well-defined inclusions, we will show that thecomputational speed of the (re)construction process can be increased by an algorithm mod-ification using surface optimization, which is essentially a biased pixel-selection procedure.This algorithm modification is based on the fact that in the later stages of (re)constructions,further refinements of the configurations are achieved by the moves of black pixels on thesurfaces of formed pixel-clusters. Thus, only the “surface pixels” are selected and givenrandom moves to generate trial configurations, i.e., the surfaces are optimized. The phys-ical analog of this process is solidification, in the later stage of which the re-arrangementsof solutes only occur on the surfaces of nucleated particles. We will see in the followingthat one can obtain a much lower error (discrepancies between the correlation functions ofthe (re)constructed medium and the imposed ones) when surface optimization is properlyapplied. This improvement enables one to quantitatively study the non-uniqueness problemof the (re)constructions [5, 11, 12, 13, 14, 15].The rest of the paper is organized as follows: In Sec. II, we describe the lattice-pointalgorithm and the algorithm modification using surface optimization in great detail. Thechoice of different pixel-lattices is also discussed. In Sec. III, we discuss the significanceof the energy threshold for both the orthogonal S -sampling method and the lattice-pointalgorithm. In Sec. IV, we suggest a general form of the map ℘ and apply the theoreticalformalism to model a Fontainebleau sandstone, a boron carbide/aluminum composite anda binary laser-speckle pattern. We also show how to construct materials with structuralproperties of interest by manipulating the parameters in the basis functions and choosingproper combination coefficients. In Sec. V, we make concluding remarks. II. ALGORITHMS FOR GENERATING HETEROGENEOUS MATERIALS
In paper I, we derived the exact algebraic equations for the (re)construction problem andshowed that the equations have an infinite number of solutions and cannot be solved rigor-ously. In principle, the Yeong-Torquato scheme enables one to obtain one of the solutionsefficiently. The most time-consuming steps of the scheme are the samplings of the two-pointcorrelation function at every trial configuration. An efficient S -sampling method would5ramatically speed up the (re)construction process. Furthermore, an isotropy-preservingalgorithm is required in (re)constructions whenever a radial two-point correlation function S ( r ) ( r ≡ | r | ) is employed for the case of statistically homogeneous and isotropic media. A. Lattice-Point Algorithm
The lattice-point algorithm is designed to sample the digitized representation of a statisti-cally homogeneous and isotropic medium in all possible directions efficiently. For simplicity,we will illustrate the idea in 2D. Implementation of the algorithm in three dimensions isa straightforward extension. Instead of considering the digitized medium as a collectionof black and white pixels, one can think of the medium as a lattice-gas system: the blackpixels are the “gas molecules” and the white pixels are unoccupied lattice sites, as shown inFig. 1. The “gas molecules” are free to move from one lattice site to another, subject to theimpenetrability condition, i.e., each lattice site can only be occupied by one “gas molecule”.Thus, the black pixels behave like hard particles and the volume fraction of black phase isconserved during the evolution of the system. (a) (b)
FIG. 1: Digitized medium as lattice-gas system: (a) White pixels as unoccupied lattice sites. (b)Black pixels as non-overlapping “gas molecules”.
For a statistically homogeneous and isotropic particle system (e.g., an equilibrium hard-sphere system or an equilibrium lattice-gas system in d -dimensional Euclidean space ℜ d ),the most basic statistical descriptor of the spatial correlations of the particles is the pair6orrelation function g ( r ) [1]. The quantity ρg ( r ) s ( r ) dr is proportional to the conditionalprobability of finding the center of a particle in a d -dimensional spherical shell of volume s ( r ) dr , given that there is another particle at the origin. Here ρ = N/V is the numberdensity of the system and s ( r ) is the surface area of a d -dimensional sphere of radius r ,which is given by s ( r ) = 2 π d/ r d − Γ( d/ , (2)where Γ( x ) is the Euler-Gamma function. Hence, for a finite system, integrating ρg ( r ) overthe volume yields ( N − ρg ( r ) s ( r ) dr is the average number of particles at a radial distance between r and r + dr from a reference particle. Practically, g ( r ) can be obtained from simulations by generatinga histogram for the number of particles n ( r ) contained in a concentric shell of finite thickness(“bin” width) ∆ r at radial distance r from a arbitrarily chosen reference particle.The two-point correlation function S ( r ) of a statistically homogeneous and isotropicmedium can be interpreted as the probability that both of ends of a randomly oriented linesegment with length r fall into the phase of interest, say, the “black” phase. By comparison,we can see that the method of computing g ( r ) of an isotropic particle system implies anatural way of obtaining S ( r ) of a statistically isotropic digitized medium, which efficientlyuses all possible vector-information in the medium. For each black pixel (or “gas molecule”) i , the distances between pixel i and all the other pixels j are computed and binned togenerate a histogram for the number of black pixels separated from each other by distance r . Another histogram for the number of lattice sites separated from a reference site bydistance r is also generated by computing and binning all possible site-separation distances.Suppose the histograms for black pixels and lattice sites are stored in the array BN [ r ]and SN [ r ], respectively. It is easy to show that the two-point correlation function can beobtained from S [ r ] = BN [ r ] /SN [ r ] , (3)(see Fig. 2). In other words, to compute S we first calculate the fraction of occupied latticesites separated by distance r ( r is an integer) from a reference black pixel. Then we averagethis fraction over all black pixels (by choosing every black pixel as reference pixel once) to7btain S ( r ). This procedure is consistent with the geometrical probability interpretationof S . For the digitized media, a natural bin width could be the characteristic size of thepixel of the lattice, e.g., the edge length of a square pixel if square lattice is used. (a) (b) FIG. 2: (color online). Centers of particles contained in the concentric shell of a reference particle:(a) In a equilibrium hard-disk system. (b) In a lattice-gas system (or a digitized medium).
At each step in the simulated annealing process, a trial configuration is generated bymoving a randomly selected black pixel to an unoccupied lattice site. A configuration array can be used to speed up this process. In our 2D implementation, the configuration array isa 2D array with the entries being 1 and 0, corresponding to the occupied and unoccupiedlattice sites, respectively. When a trial move is made, first we need to test whether it violatesthe impenetrability condition by checking whether the underlying lattice site is occupied ornot. This process just requires constant access time to an entry of the configuration array.Note that several trial moves can be attempted before a trial configuration is found for ahigh-density system. After a trial configuration is generated, the old position informationof the selected pixel is stored in the array designated
P O .The next step is to recompute the two-point correlation function for the trial configura-tion. To do this efficiently, a distance matrix is set up when the system is initialized andupdated if a trial configuration is accepted. Note that this matrix is symmetric. For eachtrial configuration, only the position of a randomly selected black pixel is changed, e.g., the k th pixel. Thus, only the k th row and column of the distance matrix need to be updated,which requires N B operations ( N B is the total number of black pixels in the system). The8ld entries in the k th row (or column) of the distance matrix are stored in the array DO .It is also unnecessary to recompute the entire array BN [ r ]. Recall that BN [ r ] contains thebinned number of black pixels separated from each other by distances between r and r + ∆ r (∆ r is the bin width). Once the k th pixel is selected, the distances between pixel k and allthe other pixels are binned and stored in the array BN O [ r ]. After the trial configuration isgenerated, the new distances between pixel k and all the other pixels are binned and storedin the array BN N [ r ]. The array BN [ r ] is updated as follows for each r : BN [ r ] − BN O [ r ] + BN N [ r ] → BN [ r ] . (4) S [ r ] is then recomputed using Eq. (3) and the trial configuration is accepted with theprobability given by Eq. (57) in paper I. If the trial configuration is rejected, all informationof the old configuration can be restored easily using the arrays P O , DO , BN O [ r ] and BN N [ r ]. For example, BN [ r ] of the old configuration can be restored by BN [ r ] − BN N [ r ] + BN O [ r ] → BN [ r ] . (5)The above procedures are repeated until the energy of the (re)constructed medium isbelow the energy threshold [see Eq. (1) in Sec. III] or the total number of evolutions reachesthe prescribed limit value. B. Algorithm Modification Using Surface Optimization
In the (re)construction of media composed of well-defined “particles” or large clusters, wefind out that in the later stages of simulated annealing, the random pixel-selection processis very inefficient. Many trial moves are rejected because a majority of selected pixels areinside the formed “particles” or clusters and it is energetically unfavorable to move themoutside. This requires the use of a biased pixel-selection process, i.e., surface optimization.The idea of surface optimization is analogous to the physical process of solidification:when the nuclei of proper sizes have been formed, they capture more solutes from surroundingsolution to further decrease the total free energy. In the simulated annealing process, once the“nuclei” are formed, the random pixel-selection process is replaced by a biased pixel-selectionprocess, i.e., only those in the surrounding “solution” or on the surface of a “nucleus” areselected to be moved to or along the surface of randomly selected “nuclei”.9t is natural and easy to incorporate the surface optimization with the lattice-point al-gorithm because the black pixels are already considered as “molecules”. Each black pixelis assigned a “free energy”, which is the minus of the number of the nearest neighbors ofthat pixel. If a pixel is inside a “nuclei”, it has the largest number of nearest neighbors andthe lowest “free energy”, which equals − − low-energy subset that contains pixels with lowest “free energy” and the high-energy subset that contains the other pixels. When the trial move is made, only the pixels in thehigh-energy subset are selected to bias the move. If the trial move is accepted, the freeenergy of the moved pixel and its neighbor pixels are recomputed; the subsets are updated. n -10-8-6-4-20 L og ( E ) FIG. 3: Energy E of the constructed medium as a function of stages n . The target autocovariancefunction is the Debye random medium function f D given by Eq. (22), with a = 30 (pixels) and vol-ume fractions φ = φ = 0 .
5. The linear size of the system N = 200 (pixels). Surface optimizationis applied at n = 100, when large clusters have been formed in the constructed medium. Numerical experiments show that applying surface optimization properly will furtherdecrease the final energy of the (re)constructed medium by a factor of 10 − (see Fig. 3), if thesame cooling schedule is used. Thus, surface optimization enables one to efficiently producemore accurate structures associated with the imposed correlation functions. Moreover, by10etting a lower energy threshold, one can address the non-uniqueness issue quantitatively.Progress on this topic will be reported in our future publications. C. Choice of Lattices
A 2D digitized heterogeneous material (medium) can be represented as a 2D array andthe real morphology of the material also depends on the lattice on which the system of pixelsis built [2]. Here, we mainly focus on two commonly used lattices in the literatures, namely,the square lattice and the triangular lattice (see Fig. 4). (a) (b)
FIG. 4: (a) The pixels (squares) and intrinsic directions of a square lattice. (b) The pixels(hexagons) and intrinsic directions of a triangular lattice.
Implementation of the algorithms discussed above on a square lattice is easier than thaton a triangular lattice because the former has orthogonal lattice vectors. However, the pixelsof a square lattice (squares) only have 4-fold symmetry while those of a triangular lattice(hexagons) have 6-fold symmetry, as shown in Fig. 4. We find that for the media withlong-range correlations or large-scale structures, using a square lattice usually introducesundesirable anisotropy in the (re)constructed systems, as shown in Fig. 5; while for themedia without long-range order or large-scale structures, both lattices work well.It is worth pointing out that the triangular lattice is superior to the square lattice in the(re)construction of anisotropic materials. In that case, one cannot count on lattice-pointalgorithm, which only uses radial averaged structural information. Several optimization11 a) (b)
FIG. 5: Constructions of media with large-scale structures. (a) Medium generated on a squarelattice with a square unit cell. (b) Medium generated on a triangular lattice with a rhombus unitcell. directions need to be specified and processed separately. The triangular lattice has sixintrinsic directions due to the higher symmetry of its pixel shape while the square latticeonly has four (see Fig. 4). Methods using more optimization directions have been proposedby Torquato and coworkers [12], i.e., along the 45- and 135-degree directions in the squarelattice; but one needs to pay the cost of complexity of implementation in these cases.
III. SIGNIFICANCE OF THE ENERGY THRESHOLD
An important issue that has not been emphasized in our previous work on the(re)construction algorithms is the choice of the energy threshold E th , i.e., the error tol-erance of discrepancies between the statistical properties of the generated structure and theimposed ones. Recall that in our case, the energy is defined as the sum of squared differ-ences between target ˆ S ( r ) and S ( r ) of the constructed medium [see Eq. (1)] and the energythreshold E th is a prescribed value of E such that when E ≤ E th , the (re)construction isterminated. E th indicates how accurately the medium is (re)constructed. A smaller E th means the two-point correlation function of the generated medium matches the imposed onebetter.In the following, we will estimate the change of the energy of a digitized medium caused12y a perturbation of its original structure (i.e., displacing a randomly selected black pixel).This change of energy can be considered as the energy difference between the target and theconstructed medium. Thus, by imposing a particular value of the energy difference (i.e., E th ),we can in turn estimate the number of perturbed black pixels. In this way, we quantitativelyrelate the energy threshold to the question of how well the medium is (re)constructed. Ingeneral, the energy difference depends on how S is sampled, the linear size of the system N and the number of black pixels N B (volume fraction of black phase φ ), which will bediscussed accordingly. A. E th for the Lattice-Point Algorithm First, we consider the energy threshold of the lattice-point algorithm. The two-pointcorrelation function of a digitized medium is computed from Eq. (3). If the generatedstructure perfectly matches the target structure, the energy given by Eq. (1) should beexactly 0. Suppose the perfect (re)constructed structure is perturbed by moving one of itsblack pixels to an unoccupied lattice site; then S ( r ) is different from ˆ S ( r ) and E becomesa small positive number. Note that we assume the perturbed structure does not have thesame S as the original structure, i.e., there is no structural degeneracy. Let E min denotethe smallest positive value of E . From Eq. (1), we see that the difference between S andˆ S of every bin contributes to E . From Eq. (3) and the fact that the last bin (farthest awayfrom the reference center) contains the largest number of pairs of lattice sites, we have E min = (cid:20) S (cid:18)(cid:20) N (cid:21)(cid:19) − ˆ S (cid:18)(cid:20) N (cid:21)(cid:19)(cid:21) = 4 ω N , (6)where [ N/
2] is the integer part of N/ ω N is the number of elements of the integer setΩ N , which is defined byΩ N = ( ( m, n ) | − (cid:20) N (cid:21) ≤ m, n ≤ (cid:20) N (cid:21) , (cid:18)(cid:20) N (cid:21) − (cid:19) ≤ m + n ≤ (cid:18)(cid:20) N (cid:21)(cid:19) ) . (7)In other words, Eq. (6) estimates E min using the change of E from 0 caused by removing(adding) one pair of black pixels from (to) the last bin (the pixels bounded in pairs due tothe symmetry of the distance matrix). Also note that the pair of pixels are originally in (or13oved to) a “bin” for the pair distance larger than [ N/ S .For a particular value of E th , the maximum number of misplaced pairs of black pixels canbe estimated by N m = E th E min = 14 ω N E th . (8)The ratio of the number of misplaced pairs of black pixels over the total number of pairs ofblack pixels is given by γ = N m N tot = ω N E th φ N . (9)Instead of specifying E th , one can also specify the ratio γ = γ s ; and the threshold can becomputed from Eq. (9): E th = 4 ω N γ s φ N . (10)For example, consider a system composed of 200 ×
200 pixels ( N = 200) and φ = 0 . N S = 40 ,
000 and the total number of black pixelsis N B = 20 , ω N can be obtained numerically, which is ω N ∼ . From Eq. (6),we have E min ∼ − . Suppose we choose the threshold E th = 10 − , from Eq. (8), wehave the number of misplaced pairs of black pixels N m ∼ , which seems to be a largenumber. However, when considering the total number of pairs in the system, we have N tot = N B = 4 × ; from Eq. (9), the ratio is given by γ = N m N tot ≃ × = 2 . × − , (11)which means only one out of a half million pairs is put in the wrong bin. Thus, the mediumis (re)constructed to a very high accuracy. In our simulations, we choose the threshold E th = 10 − . B. E th for the Orthogonal S -Sampling Algorithm For the orthogonal S -sampling algorithm, S ( r ) is sampled line (column) by line (column)by moving a line (column) segment of length r one pixel distance each time and counting the14imes that both ends of the segment are black pixels. This number is then divided by thetotal number of times that one moves the line (column) segment (total number of pixels onthe line (column)) to obtain S of that line (column). Finally, the S sampled from differentlines and columns of the digitized medium are averaged to compute the S of the wholemedium. Similarly, consider the perfect (re)constructed structure is perturbed by movingone of its black pixels to an unoccupied lattice site and there is no structural degeneracy, E min is given by E min = 1 N , (12)where N is the linear size of the system. For a particular E th , the maximum number of misplaced black pixels is given by N m = E th E min = N E th . (13)Thus, the ratio of misplaced black pixels over the total number of black pixels can beobtained by γ = N m N P = 1 φ N E th . (14)If stead, γ = γ s is specified, the required energy threshold E th can be obtained from Eq. (14): E th = γ s φ N . (15)Consider the same system used in the previous section ( N = 200, φ = 0 . γ s = 10 /N = 2 . × − . From Eq. (15),we have E th = 2 . × − × . × ≃ − . (16)In other words, if we choose the threshold E th = 10 − , only 10 black pixels in the systemare misplaced. The medium is (re)constructed to a high accuracy.It is worth noting that in the above discussions, we assume that the perturbed structuredoes not have the same S of the original structure, i.e., there is no structural degeneracy. Ingeneral, however, structural degeneracy does exist (i.e., media with the same S but different15 , S , ...). Results concerning these non-uniqueness issues will be reported in our futurepublications. IV. APPLICATIONS OF THE THEORETICAL FORMALISM
In this section, we illustrate the practical utility of our theoretical formalism in both twoand three dimensions. In particular, we will consider two kinds of applications. Firstly,given a set of basis functions (may not be complete), one can express the scaled autocovari-ance functions (two-point correlation functions) of a statistically homogeneous and isotropicmedium in terms of the specified basis functions with certain accuracy. Realizations of thematerials can be generated using proper (re)construction procedures and subsequent anal-ysis can be performed on the images to obtain effective macroscopic properties of interest;see, e.g., Ref. [10]. Secondly, the map ℘ (see the discussion below) enables one to constructcandidates of realizable two-point correlation functions with properties of interest, which inturn enables one to design and investigate materials with desired structural characteristics.The (re)construction of realizations of 3D medium from the information obtained from a2D micrograph or image is of great value in practice [3]. Therefore, we will apply the lattice-point algorithm to generate three-dimensional digitized realizations of a Fontainebleau sand-stone and a boron carbide/aluminum composite from two-dimensional tomographic imagesof their slices through the materials. In a successful reconstruction, other deemed crucialstructural characteristics besides S obtained from the reconstructed medium should alsoagree closely with that of the target medium. Consequently, in order to judge quantita-tively how well the reconstructions are, we will measure and compare another importantmorphological descriptor, i.e., the two-point cluster function C ( x , x ), defined to be theprobability of finding two points at x and x in the same cluster of the phase of interest [4].For statistically homogeneous and isotropic media, C only depends on the relative scalerdistances between the points, i.e., C ( x , x ) = C ( | x − x | ) = C ( r ). Note that C containsnontrivial topological “connectedness” information. When large clusters are present in themedium, C becomes a long-ranged function and its integral will diverge if the phase ofinterest percolates. The measurement of C for a 3D material sample cannot be made froma 2D cross-section of the material, since it is an intrinsically 3D microstructural function[1]. To sample C from a digitized medium, we associate each pixel with a cluster-index16hat indicates to which cluster the pixel belongs, and only bin the distances of pixel-pairsin the same cluster. Then the number of pair distances in each bin is normalized by thetotal number of distances between the lattice sites within that bin, which is similar to theprocedure of sampling S discussed in Sec. II.We will also study the reconstruction of a binary laser-speckle pattern in two dimensionsand show that the algorithm cannot reproduce the pattern accurately. Moreover, we in-troduce and discuss a classification of heterogeneous materials into multi-scale media and single-scale media . A multi-scale medium (MSM) is the one in which there are multiplecharacteristic length scales associated with different structural elements. Examples of MSMinclude fractal patterns and hierarchical laminate composites [16]. A single-scale medium(SSM) is the one composed of structural elements associated with only one characteristiclength scale, such as a Fontainebleau sandstone and a boron carbide/aluminum composite.We conclude that in general reconstructions using S only work well for heterogeneous ma-terials with single-scale structures. However, two-point information via S is not sufficientto capture the key structural features of multi-scale media. Before presenting the aforemen-tioned applications, we will first consider a general form of ℘ , which includes all possibleconvex combinations of the basis functions. A. A General Form of ℘ In paper I, we introduced the idea of expressing the two-point correlation function of astatistically homogeneous and isotropic medium through a selected set of bases of the two-point correlation function space. In practice, it is convenient to use the scaled autocovariancefunctions that are equivalent to the two-point correlation functions, which are defined asfollows: f ( r ) ≡ S ( i )2 ( r ) − φ i φ φ . (17)Suppose { f i ( r ) } mi =1 is a set of bases of scaled autocovariance functions and ℘ is an map on { f i ( r ) } mi =1 composed of convex combinations and products of f i ( r ) ( i = 1 , ..., m ), thus f ( r ) = ℘ [ { f i ( r ) } mi =1 ] = ℘ [ f ( r ) , f ( r ) , ..., f m ( r )] , (18)17s also a realizable scaled autocovariance function. Note that the choice of basis functions isnot unique.In general, one can consider that ℘ [ { f i ( r ) } mi =1 ] takes the form ℘ [ { f i ( r ) } mi =1 ] = X i α i f i ( r ) + X i,j β ij f i ( r ) f j ( r ) + · · · , (19)where the coefficients satisfy the condition X i α i + X i,j β ij + · · · = 1 . (20)Then one can use standard regression methods to obtain the set of coefficients such that f ( r ) = ℘ [ { f i ( r ) } mi =1 ] is the best approximation of the target function pointwisely.The basis functions we consider here include Debye random medium function f D ( r ), thefamily of polynomial functions f ( n ) P ( r ), damped oscillating function f O ( r ), overlapping-spherefunction f S ( r ) and symmetric-cell material function f C ( r ), which are discussed in paper I. B. Modeling Real Materials
1. Fontainebleau Sandstone
First, we investigate structural properties of a Fontainebleau sandstone from a two-dimensional tomographic image of a slice through the material sample. Sandstone is animportant porous medium in geo-physical and petroleum applications and has been the fo-cus of many studies [3, 17, 18, 19]. A microstructural image of a slice of a Fontainebleausandstone is shown in Fig. 6(a), in which the black areas are solid phases (phase 1) and thewhite areas are void phases (phase 2). The two-point correlation function of the void phase S F S ( r ) is shown in Fig. 6(b).The scaled autocovariance function of the void phase in Fontainebleau sandstone f F S ( r )can be approximated by the convex combination of f D ( r ) and f O ( r ) as follows: f F S ( r ) = α f D ( r ) + α f O ( r ) , (21)where α = 0 . α = 0 .
23 are the combination coefficients and18 a) S FS (r) (b) FIG. 6: (a) A microstructrual image of a slice of a Fontainbleau sandstone [17]. The black phaseis the solid phase with φ = 0 . φ = 0 . S F S ( r ) of the void (white) phase. f D ( r ) = exp( − r/a ) , (22) f O ( r ) = exp( − r/b ) cos( qr + ψ ) , (23)where a = 3 ( pixel ), b = 6 . pixel ) are the effective correlation length; q = 0 . pixel − ) isthe oscillating frequency and ψ = 0 is the phase angle. The two-point correlation function S F S ( r ) is approximated by 19 a) (b) FIG. 7: (color online). (a) The reconstructed 3D realization of the Fontainbleau sandstone from S F S ( r ) of the void phase. The void phase is shown in yellow and the solid phase is transparent foreasy visualization. (b) A 2D slice of the constructed 3D realization. The solid phase is shown inblack and the void phase is shown in white. The linear size of the system N = 160 (pixels). S F S ( r ) = f F S ( r ) φ φ + φ + δS ( r ) , (24)where φ = 0 . φ = 0 .
175 are volume fractions of the solid and void phase, respectively; δS ( r ) is the discrepancy between the sampled two-point correlation function and the basis-function approximation. The average of the absolute values of discrepancies ∆ S , whichindicates how well the sampled two-point correlation function is approximated by the convexcombination, is defined as: ∆ S = 1 N L X r | δS ( r ) | ≃ . × − . (25)where N L = 80 is the sample-length. Note that although the general form of ℘ given byEq. (19) would work well if a complete set of basis functions is given, in the present case,our practical approach is enough.The 3D reconstruction of the Fontainebleau sandstone from S F S ( r ) obtained from thedigitized image of a 2D slice [Fig. 6(a)] is shown in Fig. 7. Visually, the reconstruction pro-vides a good rendition of the true sandstone microstructure, as can be seen by comparing the20
10 20 30 40 50r (pixels)00.050.10.150.2 C (r) C of the target mediumC of the reconstructed medium FIG. 8: The two-point cluster functions C v ( r ) for the void (white) phase of the target and recon-structed 2D slices of the Fontainebleau sandstone.
2D images. As pointed out earlier, to ascertain whether the reconstruction is quantitativelysuccessful, we also measure C of the target and generated media. Since it is an intrinsically3D microstructural function, we only compute and compare the two-point cluster functionsof the 2D slices. The measured C v ( r ) for the void phase of the target and the reconstructedslices of the Fontainebleau sandstone are shown in Fig. 8. The figure reveals that although C v ( r ) of the generated medium is slightly below that of the target medium, the discrep-ancies are acceptable (e.g., the largest discrepancy | δC v | ∼ × − ). Thus, we considerthe reconstruction is successful. Also note that the close agreement of the two-point clusterfunction indicates that C v ( r ) is largely determined by S F S ( r ) of the medium and suggeststhat C v ( r ) might be expressed as a functional of S F S ( r ).
2. Boron Carbide/Aluminum Composite
A 2D digitized image of a boron carbide/aluminum ( B C / Al ) inter-penetrating compositeand the sampled two-point correlation function of the aluminum phase (white phase) S Al ( r )[20] are shown in Fig. 9. As we can see from Fig. 9(b), S Al ( r ) is essentially an exponentiallydecreasing function without any significant short-range correlation. Thus, S Al ( r ) can beapproximated by S Al ( r ) = α ′ φ φ f D ( r ) + α ′ φ φ f O ( r ) + φ + δS ′ ( r ) , (26)21 a) S A l (r) (b) FIG. 9: (a) A digitized image of a boron carbide/aluminum composite [20]. The black phaseis boron carbide with φ = 0 . φ = 0 . S Al ( r ) of the aluminum phase. where α ′ = 0 . α ′ = 0 .
19 are combination coefficients; φ = 0 .
647 and φ = 0 . f D ( r ) and f O ( r ) are given by Eq. (22) and Eq. (23), respectively; the parametersused are a = 3 ( pixel ), b = 10 ( pixel ) and q = 0 .
22 ( pixel − ), ψ = 0. The average of theabsolute values of discrepancies ∆ S ′ is given by∆ S ′ = 1 N L X r | δS ′ ( r ) | ≃ . × − . (27)22 a) (b) FIG. 10: (color online). (a) The reconstructed 3D realization of the boron carbide/aluminum com-posite from S Al ( r ). The aluminum phase is shown in yellow and the ceramic phase is transparentfor easy visualization. (b) A 2D slice of the constructed 3D realization. The ceramic phase isshown in black and the aluminum phase is shown in white. The linear size of the system N = 160(pixels). C A l (r) C of the target mediumC of the reconstructed medium FIG. 11: The two-point cluster functions C Al ( r ) for the aluminum (white) phase of the target andreconstructed 2D slices of the boron carbide/aluminum composite. N L = 80 is the sample-length.The 3D reconstruction of the boron carbide/aluminum composite from S Al ( r ) obtainedfrom the digitized image of a 2D slice [Fig. 9(a)] is shown in Fig. 10. As in the previ-ous section, to verify that the reconstruction is quantitatively successful, we compute andcompare the two-point cluster function C Al ( r ) for the aluminum phase of the target andreconstructed 2D slices of the composite. As shown in Fig. 11, the reconstruction alsoslightly underestimates C Al ( r ) with acceptable discrepancies (e.g., the largest discrepancy | δC Al | ∼ × − ). The medium can be considered as successfully reconstructed. Also notethat C Al ( r ) is long-ranged, indicating large clusters of the aluminum phase present in themedia.
3. Binary Laser-Speckle Pattern
Although in the above two examples our theoretical formalism works well, there aresituations where the microstructural information contained in S ( r ) alone is not sufficientto accurately reconstruct a heterogeneous material. One such example is the multi-scalestructure of a binary laser-speckle pattern (Fig. 12). The figure reveals that there are threestructural elements: “particles”, “stripes” and a background “noise” (individual black pixelsdispersed throughout the white phase). Thus, there are three characteristic length scales inthe medium associated with these structural elements.The reconstruction of the speckle pattern is shown in Fig. 13. Comparing Fig. 13 withFig. 12, we can see that instead of reproducing all the structural elements in the targetmedium, the (re)construction program seems to mix them up to generate a single-scalestructure that has the same (or to a very high accuracy) two-point correlation function asthe target medium. We note that this is a numerical example of structural degeneracy of S ( r ).The laser-speckle pattern we considered is also an example of the multi-scale media (MSM). The separation of length scales in MSM results in the inefficiency of the bulk-based structure characteristics (e.g., n -point correlation functions), since usually they canonly pick up structural information associated with the largest length scale. For example,consider the dislocations in a crystalline solid as the phase of interest; it is clear that thetwo-point correlation function of the “dislocation” phase is identically zero since the “dislo-24 a) S SP (r) (b) FIG. 12: (a) A digitized image of a binary laser-speckle pattern [21]. The volume fraction of theblack phase is φ = 0 . φ = 0 . S SP ( r ) of the black phase. cation” phase has no measure compared with the bulk “crystal” phase (i.e., the volume ofdislocations is zero compared with that of the bulk crystal). In this extreme example, theratio γ of the two characteristic length scales associated with the “dislocation” phase L d andthe “crystal” phase L c is zero, i.e., γ = L d /L c = 0. In digitized media, the ratio of lengthscales is always positive due to the discrete nature of the system; however when the ratiois significantly different from unity, the two-point correlation function alone is not able tocapture the key structural features. On the other hand, extensive experience with successful25 IG. 13: The reconstruction of the binary laser-speckle pattern from S SP ( r ). The linear size ofthe system N = 150 (pixels). reconstructions of single-scale media (SSM) from S shows that the two-point correlationfunctions are indeed sufficient to determine the structures of SSM to a high accuracy. Thus,we conclude that in general reconstructions using S work well for heterogeneous materialswith single-scale structures, while the microstructural information contained in S is notsufficient to accurately model multi-scale media. Note that more morphological informa-tion (e.g., the lineal-path function [22] and the two-point cluster function, etc.) can beused to model MSM more accurately. However, even the cluster-type functions containingconnectedness information of the media such as C can not completely characterize MSMstatistically. C. Constructing Materials with Desired Structural Characteristics
From Eq. (18), one can construct candidates of realizable two-point correlation functionsusing the basis functions. Given a set of basis functions with diverse and interesting proper-ties, one would be able to construct a two-point correlation function that exhibits all the use-ful properties of the basis functions to some extend and generate an “optimal” structure thatrealizes all the desired structural features. Thus, the theoretical formalism enables one todesign materials with structural properties of interest. Given an accurate structure-propertyrelation, one could even design materials with physical properties of interest by manipulat-ing their two-point correlation functions. For example, Adams, Gao and Kalidindi recently26eveloped a methodology to obtain finite approximations of the second-order properties clo-sure in single phase polycrystalline materials, from which the second-order microstructuredesign can proceed [23].As pointed out in paper I, we know very little about the basis function set { f i ( r ) } mi =1 atthis stage. Our choice of { f i ( r ) } mi =1 is based on the criteria that f i ( r )’s should have simpleanalytical forms and they present typical features of certain known two-point correlationfunctions. Important features exhibited by most realizable two-point correlation functionsare monotonically deceasing or damped oscillating, which corresponds to materials withoutor with significant short-range order, respectively. Another feature of two-point correlationfunctions that may affect the structures of the corresponding materials is the smoothnessof the function, which is not emphasized in our previous work. In the following, we willsee through several examples that the properties of basis functions can be observed in thegenerated structures; also the media with desired structural properties can be obtainedby manipulating the combination coefficients and the parameters (e.g., effective correlationlength) in the basis functions.In the paper I, we already investigated hypothetical correlation functions combining theexponentially decreasing and damped-oscillating features. In this section, we provide anexample of non-smooth correlation functions (with discontinuous derivatives), i.e., the poly-nomial function of order two f (2) P ( r ), which is defined as f (2) P ( r ) = (1 − r/c ) ≤ r ≤ c, r > c, (28)where the parameter c is the effective correlation length. Other correlation functions hav-ing this non-smooth feature include the functions of known constructions (e.g., f S ( r ), f C ( r ),etc.) and other polynomial functions [2]. A typical medium associated with this type of func-tions is composed of dispersions of fully penetrable “particles”. The monotonically decayingpart of the functions determines the size and shape of the “particles” and the penetrabilityis consistent with the flat “tail” of the functions, which implies no spatial correlation be-tween the “particles”. When these functions are use as dominant basis functions, one canexpect the generated media also contain well-defined overlapping “particles”; thus, surfaceoptimization could be applied in the (re)constructions. However, when functions like Debye27andom medium function are dominant in the convex combination, this algorithm modi-fication should be used with care because the media could contain “clusters” of all sizesand shapes, which would significantly reduce the efficiency of the modified algorithm usingsurface optimization.Besides f (2) P ( r ), we also use Debye random medium function f D ( r ) and damped-oscillatingfunction f O ( r ) as the basis functions. Consider the simplest form of ℘ for these three basisfunctions, i.e., f ( r ) = α f D ( r ) + α f O ( r ) + α f (2) P ( r ) , (29)where α i ( i = 1 , ,
3) satisfies 0 ≤ α i ≤ P i α i = 1; and f D ( r ), f O ( r ) and f (2) P ( r )are given by Eq. (22), Eq. (23) and Eq. (28), respectively. The characteristic length scalesof the generated structures are determined by the parameters in the basis functions; andthe ratios of the characteristic lengths are chosen to be close to unity, i.e., the hypotheticalmedia belong to SSM. f(r) FIG. 14: The constructed autocovariance function f ( r ) given by Eq. (29). The combinationcoefficients are ( α , α , α ) = (0 . , . , . a = 8 ( pixel ); b = 5 ( pixel ), q = 1 ( pixel ) − , ψ = 0; c = 5 ( pixel ). Suppose we choose the combination coefficients ( α , α , α ) = (0 . , . , . a = 8 ( pixel ); b = 5 ( pixel ), q = 1 ( pixel ) − , ψ = 0; c = 5 ( pixel ). The combined autocovariance function f ( r ) isshown in Fig. 14 and the constructed structure with volume fraction of black phase φ = 0 . a) (b) FIG. 15: (color online). (a) The constructed 3D structure associated with f ( r ) shown in Fig. 14with volume fractions φ = φ = 0 .
5. The “black phase” is shown in yellow and the “white phase”is transparent for easy visualization. (b) A 2D slice of the constructed 3D realization. The linearsize of the system N = 150 (pixels). is shown in Fig. 15. From the figures, we can see the effects of each basis function on thegenerated structure. The short-range correlations in the structure are determined by f (2) P ( r ),which allows spatially uncorrelated “particles” with diameter c to form. These “particles”would form clusters of different sizes as the volume fraction of the black phase increases.The middle-range correlations in the structure are dominated by the oscillation part in f ( r ), which is the contribution of f O ( r ). As we have shown in paper I, the parameter q is manifested as a characteristic repulsion among different structure elements (appearing indifferent forms at different volume fraction) with diameter of order q ; while b controls theoverall exponential damping, and thus, the effective range of the repulsion. The large-scalecorrelation is dominated by the long “tail” of f D ( r ), which allows clusters of all shapesand sizes to form. In Fig. 15, we can clearly identify the structure elements associatedwith each basis function. For example, the stripe-like structures are associated with theoscillating feature of f ( r ), the width of which is approximately q . Several “particles” withdiameter a can be found dispersed in the white phase, which correspond to the contributionof f (2) P ( r ). Note that most “particles” form clusters together with the “stripes” at this29elatively high volume fraction φ = 0 .
5. Finally, the spatial distribution of large-scaleclusters is determined by the “tail” of f D ( r ). f(r) FIG. 16: The constructed autocovariance function f ( r ) given by Eq. (29). The combinationcoefficients are ( α , α , α ) = (0 . , . , . a = 8 ( pixel ); b = 10 ( pixel ), q = 0 . pixel ) − , ψ = 0; c = 15 ( pixel ). Suppose now we would like to generate a similar structure as the one in Fig. 15 but withlarger “particles” and clusters. To “grow” the “particles”, we need to increase the param-eter c in f (2) P ( r ), which controls the effective diameter of the “particles”. To form largerclusters, we need to reduce the “repulsion” between the structural elements; or equivalently,to suppress the oscillation introduced by f O ( r ). This can be done by either increasing theparameter b in f O ( r ) or decreasing the coefficient α . We first adopt the former method.Thus, to construct a new structure with the required properties, we use the following pa-rameters: a = 8 ( pixel ); b = 10 ( pixel ), q = 0 . pixel ) − , ψ = 0; c = 15 ( pixel ); and usethe same α i ( i = 1 , ,
3) as the last example. The combined autocovariance function f ( r ) isshown in Fig. 16 and the constructed structure with volume fraction of black phase φ = 0 . f ( r ).Now we would like to further increase the size of clusters in the constructed medium whilekeeping all the parameters in the basis functions unchanged. This time we need to adjustthe combination coefficients to obtain the required structure. To generate larger clusters,30 a) (b) FIG. 17: (color online). (a) The constructed 3D structure associated with f ( r ) shown in Fig. 16with volume fractions φ = φ = 0 .
5. The “black phase” is shown in yellow and the “white phase”is transparent for easy visualization. (b) A 2D slice of the constructed 3D realization. The linearsize of the system N = 150 (pixels). the repulsion effect need to be further suppressed; thus, we need to reduce α . Since thesum of the coefficients must equal 1, we also need to increase the other two α i ( i = 1 , α makes f (2) P ( r ) dominant in the combination; however, aswe discussed before, f (2) P ( r ) corresponds to spatial uncorrelated “particles” and the sizeof clusters formed by these “particles” is determined by the volume fraction the “particle”phase, which cannot be changed here. So we choose to increase α of f D ( r ) based on the factthat the large correlation length and long “tail” of f D ( r ) significantly stimulate the formingof large clusters. Thus, in the construction we use the following combination coefficients:( α , α , α ) = (0 . , . , . f ( r ) is shown in Fig. 18 and the constructedstructure with volume fraction of black phase φ = 0 . ℘ ) can provide variety of candidates of realizable two-pointcorrelation functions. We can also obtain desired structures by manipulating the combina-tion coefficients and the parameters in the basis functions. Thus, we argue that in general31
20 40 60 80r (pixels)00.20.40.60.81 f(r)
FIG. 18: The constructed autocovariance function f ( r ) given by Eq. (29). The combination coef-ficients are ( α , α , α ) = (0 . , . , . a = 8 ( pixel ); b = 10 ( pixel ), q = 0 . pixel ) − , ψ = 0; c = 15 ( pixel ). given a complete set of basis functions, Eq. (19) enables one to construct candidates of re-alizable functions with desired properties and to design materials with structural propertiesof interest. V. CONCLUSIONS
In this paper, we described in great detail the lattice-point algorithm, which has beenshown to be both efficient and isotropy-preserving for (re)constructing statistically homo-geneous and isotropic media. The digitized medium (pixel system) can be considered as alattice-gas system, in which the black pixels behave like nonoverlapping “particles” movingfrom one lattice site to another. The two-point correlation function of the medium is sampledby binning all distances between black pixels and dividing the number in each bin by the totalnumber of distances between all lattice sites in that bin. The energy threshold indicating towhat accuracy the medium is (re)constructed has been discussed in detail for the first time.This quantity is directly related to the issue of non-uniqueness of (re)constructions. Wehave also described an algorithm modification using surface optimization to further speedup the (re)construction process and discussed its implementation based on our lattice-gasversion of the digitized media. Numerical experiments have shown that by applying the32 a) (b)
FIG. 19: (color online). (a) The constructed 3D structure associated with f ( r ) shown in Fig. 18with volume fractions φ = φ = 0 .
5. The “black phase” is shown in yellow and the “white phase”is transparent for easy visualization. (b) A 2D slice of the constructed 3D realization. The linearsize of the system N = 150 (pixels). surface optimization algorithm properly (i.e., in the (re)constructions of media composedwell-defined “particles”), the final energy can be reduced by a factor of 10 − if the samecooling schedule is used. The choice of different pixel-lattices has also been discussed.We have applied the theoretical formalism proposed in paper I to model several exam-ples of real materials. In particular, we used the lattice-point algorithm to generate 3Drealizations of a Fontainebleau sandstone and a boron carbide/aluminum composite from2D images of their slices through the materials. The two-point cluster functions in the re-constructions in these two different cases matched the corresponding ones in the originalmaterials, thus demonstrating the efficiency and accuracy of the (re)construction algorithmin these cases. We also studied the reconstruction of a binary laser-speckle pattern in 2D,in which the algorithm fails to reproduce the target pattern accurately. We conclude thatin general reconstructions using S only work well for heterogeneous materials with single-scale structures. However, two-point information via S is not sufficient to accurately modelmulti-scale media. In paper I, we pointed out that given a complete set of basis functions,one can obtain the basis-function approximations of the two-point correlation functions for33he materials of interest to any accuracy. Here we only used a simple form of Eq. (19) dueto our limited knowledge of the basis functions set. We also constructed realizations of ma-terials with desired structural characteristics by manipulating the combination coefficientsand the parameters in the basis functions. Thus, Eq. (19) enables one to construct candi-dates of realizable functions with desired properties and to design materials with structuralproperties of interest.We are now developing efficient (re)construction procedures that take into account ad-ditional microstructural information that characterize the media. For example, using ourlattice-gas version of the digitized media, the two-point cluster functions C can be directlyincorporated into the Yeong-Torquato scheme [5], which enables one to obtain more accu-rate (re)constructions of the target media. Also, a quantitative structure-property relationis necessary for the further applications of the theoretical formalism. Finite element anal-ysis on the generated media will be performed to establish the relation. Such work will bereported in our future publications. Acknowledgments
Acknowledgment is made to the Donors of the American Chemical Society PetroleumResearch Fund for support of this research. [1] S. Torquato,
Random Heterogeneous Materials: Microstructure and Macroscopic Properties (Springer-Verlag, New York, 2002).[2] Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. E , 031110 (2007).[3] C. L. Y. Yeong and S. Torquato, Phys. Rev. E , 224 (1998).[4] S. Torquato, J. D. Deasley, and Y. C. Chiew, J. Chem. Phys. , 6540 (1988).[5] C. L. Y. Yeong and S. Torquato, Phys. Rev. E , 495 (1998).[6] K. Wu, M. I. J. Dijke, G. D. Couples, Z. Jiang, J. Ma, K. S. Sorbie, J. Crawford, I. Young,and X. Zhang, Trans. Porous Media , 443 (2006).[7] M. A. Ansari and F. Stepanek, AlChE Journal , 3762 (2006).[8] R. Hilfer and C. Manwart, Phys. Rev. E , 021304 (2001).
9] D. Basanta, M. A. Miodownik, E. A. Holm, and P. J. Bentley, Metall. Mater. Trans. A ,1643 (2005).[10] H. Kumar, C. L. Briant, and W. A. Curtin, Mech. Mater. , 818 (2006).[11] D. Cule and S. Torquato, J. Appl. Phys. , 3428 (1999).[12] N. Sheehan and S. Torquato, J. Appl. Phys. , 53 (2001).[13] M. G. Rozman and M. Utz, Phys. Rev. Lett. , 135501 (2002).[14] R. Hosemann and S. N. Bagchi, Acta Cryst. , 237 (1954).[15] C. Chubb and J. I. Yellott, Vision Research , 485 (2002).[16] J. Quintanilla and S. Torquato, Phys. Rev. E , 4368 (1996).[17] D. A. Coker, S. Torquato, and J. Dunsmuir, J. Geophys. Res. , 17497 (1996).[18] L. M. Schwartz, F. Auzerais, J. Dunsmuir, N. Martys, D. P. Bentz, and S. Torquato, PhysicaA , 28 (1994).[19] C. Manwart, S. Torquato, and H. R., Phys. Rev. E , 893 (2000).[20] S. Torquato, C. L. Y. Yeong, M. D. Rintoul, D. L. Milius, and I. A. Aksay, J. Am. Ceram.Soc. ∼ ecarpenter2.[22] B. Lu and S. Torquato, Phys. Rev. A , 922 (1992).[23] B. Adams, X. Gao, and S. Kalidindi, Acta Mater. , 3563 (2005)., 3563 (2005).