A microstructure defined by a digital image is already naturally discretized and lends itself immediately to numerical computation of many useful quantities. We use a finite element method (FEM) to estimate the elastic properties of the model system. FEM uses a variational formulation of the linear elastic equations and finds the solution by minimizing the elastic energy using a fast conjugate gradient method. Each voxel is taken to be a tri-linear finite element. A strain is applied, with the average stress or the average elastic energy giving the effective elastic modulus. The digital image is assumed to have periodic boundary conditions. Most of the simulations were performed on a modern workstation (500MHz Compaq XP1000). Each FEM simulation at 1203 required ~hour of CPU time and MB of memory. Further details of the theory and the programs can be found elsewhere [Garboczi and Day 1995, Garboczi 1998, Arns 2001]. We assign to the rock skeleton values of the elastic properties of quartz given by [Mavko 1998]: bulk modulus K = 37.0 GPa, shear modulus µ = 44.0 Gpa and mineral density = 2.65 g/cm3. We model the water saturated case at T =200ºC and 40 MPa pressure [KWater = 2.2 GPa, µWater = 0 GPa, [Han 1986, Castagna 1993]], and the oil saturated case (30 API oil) at T = 200ºC and 25 MPa pressure ( KOil = 0.5 GPa, µOil = 0 GPa) using temperature adjusted moduli [Castagna 1993]. The choice of the water saturated condition is made to allow for comparison with the experimental data from [Han 1986].
In order to obtain accurate numerical results it is necessary to estimate and minimize three sources of error: finite size effects, statistical fluctuations and discretization errors. These error analyses have not been carried out in the past due to limitations in computational speed. However, if one wishes to compare computations directly to experiment such error analyses must be carried out [Roberts and Garboczi 2000]. We discuss the three sources of error separately.
In the previous section, we argued that at the length scale of voxels the averaging of the porosity is acceptable. To test if this cell size provides good averaging of the elastic properties we compare data measured at 1203 , 1603, and 2403 cells (see Figure 4(a)). The finite size errors appear to be small for the 1203 cell size. To further quantify these errors we discuss statistical uncertainties for the various cell sizes. To estimate statistical uncertainties we bin the data as a function of porosity, , with bin sizes = 0.025. The error bars reflect twice the standard error (S.E. = / ), where is the standard deviation. There is a 95% chance that the "true" value lies between the indicated standard error bars shown in Figure 4(b). The results are accurate to within 2-3% for most data points. The binned data is presented for different cell sizes (1203, 1603, and 2403). Differences between the predictions at the different cell sizes are well within statistical error, indicating that finite size effects are acceptably small.
Finally, we consider the discretization error; the error due to the use of discrete voxels to represent continuum objects. To do this we coarse-grain the original images by generating realisations of the tomographic data sets at integer multiples of the resolution of the original image (i.e., n x 5.7 µm = 11.4 µm, 17.1 µm, ...). To generate the images at poorer resolution we bin voxel clusters of sizes n3 (n = 2, 3, 4) using a simple majority rule. It has been previously shown [Roberts and Garboczi 2000] that the variation of an elastic property, P, follows
where is the voxel size, a is a fitting parameter, and P0 is the "continuum" value corresponding to infinite resolution. In Figure 4(c) we show the scaling behavior of the discretization error with for different porosities. Figure 4(d) shows the change in elasticity with resolution along with the continuum values given by Eqn. 1. It is important to note that even at a resolution of = 5.7 µm, the predictions P() differ from the "continuum" P0 value by up to 15%. The errors increase approximately linearly with porosity. The statistical error in the continuum limit P0, is estimated from Eqn. 1 using the standard deviations at the different lattice sizes.