Next: Effective Medium Theory
The problem of determining the effective linear elastic properties of random, multi-phase materials is an old, difficult, and important problem, with applications in almost every area of materials science and engineering [8,9]. Much can be done analytically in the case where the composite is made up of inclusions having simple geometries that are randomly or regularly embedded in a matrix. However, many materials, like polycrystalline ceramics and metals, polymer blends, sandstone and carbonate rocks, and cement- based materials, must be considered as random composites at the micrometer scale or lower. Their microstructures cannot be described by the picture of inclusions of simple geometric shape embedded in a matrix. In fact, usually the only direct microstructural information one has to work with in these kinds of random materials are actual digital images. These can be 2-D images acquired using electron or light microscopy, 3-D images obtained using x-ray microtomography , serial sectioning, or magnetic resonance imaging, or 2-D and 3-D simulated images from digital-image-based microstructural models [4,5]. In all these cases, a numerical algorithm to compute elastic properties must be able to work on an arbitrary digital image in two or three dimensions.
The elastic algorithm operates directly on a digital image by treating each pixel, either a square in 2-D or a cube in 3-D, as a linear finite element. Standard finite element techniques  are used, in combination with a conjugate gradient solver  that is able to handle the many finite elements that are necessary for adequate resolution of a microstructure. Usually periodic boundary conditions are used, but this is not necessary. Once individual phase properties are supplied, the composite material properties can be computed by applying a strain and computing the appropriate stress and/or energy averages . Using the digital image as the finite element mesh simplifies the finite element code enormously, since mesh generation is often the most difficult and time consuming step in using the finite element method . This procedure also uses the maximum resolution possible from an image (one pixel equals one finite element) in the computation of elastic properties. An important feature of this algorithm is that it enables images to be made of the stress or strain fields, which can then be compared with microstructural features.
It should be emphasized that the above algorithm uses only the simplest version of the general finite element method. Higher order interpolations could be used for each pixel, improving the accuracy of the method by using quadratic or higher order interpolations for the displacement field. This would require, however, that additional nodes be assigned to each pixel, say at the middle or at the faces, because of the additional unknown coefficients in such an interpolation scheme. We believe it to be preferable, at least in a general algorithm for arbitrary complex random materials, to stay with a simple cubic lattice and use, if possible, higher resolution (more pixels per microstructural feature) with the same linear interpolation scheme to make the approximate displacement field more accurate. This algorithm has also been used to solve complicated 3-D electrical conductivity problems .
This algorithm in particular, and the finite element method in general, are closely related to spring network problems. By an appropriate choice of force constants, a spring network on a square lattice with spring stretching and bending forces may be mapped exactly onto our finite element scheme for a 2-D array of square pixels .
We use a recent digital-image-based algorithm for generating the random microstructures studied in this paper [6,7]. The algorithm does not repeatedly embed inclusions in a matrix, but rather uses a convolution and thresholding scheme.
Consider an intensity function defined on a simple cubic lattice, Io(i,j,k). Io(i,j,k) is initially a white noise random field, with the values of Io evenly and randomly distributed between 0 and 1. Now define a kernel function Ke(i'-i, j'-j, k'-k), and define a new intensity J(i,j,k) via:
A threshold Jc is then chosen, so that the final two-phase microstructure I(i,j,k) is generated by:
Many choices of the kernel Ke are possible. A particularly simple one, that was used previously to model the pore structure of carbonate rocks , is a Gaussian:
The value of w sets the length scale of the correlations in the microstructure. The resulting microstructure, as represented by I(i,j,k) and for the kernel in eq. (3), has been shown, using different values of Jc, to give microstructures where one phase can be considered as the matrix and the other phase as inclusions, and a broad region, between about 20% volume fraction limits for each phase, where each phase forms a fully connected network . In this bi-continuous region, the usual matrix-inclusion ideas used to understand composites are not necessarily useful. Figure 1 shows a cross-section of one of these 3-D composites, at a volume fraction of 50%. The image is 128 x 128 pixels, sliced from a 128 x 128 x 128 system, using the kernel of eq. (3) with w = 5. In 2-D, of course, only one phase percolates at a time, and the common percolation threshold for systems generated with the 2-D version of this kernel is at 50% area fraction, which agrees with a previous analysis . We next discuss the effective medium theory used to analyze our numerical results.
Figure 1: Showing a 128 x 128 slice through a 128 x 128 x 128 Gaussian
kernel-based microstructure (w = 5, c1 = c2 = 0.50).