Sie sind hier: ICP » R. Hilfer » Publikationen

# Article

[page 1, §1]
[1.1.1.1] Accurate prediction and understanding of material parameters for disordered systems such as rocks [1], soils [2], papers [3], clays [4], ceramics [5], composites [6], microemulsions [7] or complex fluids require geometrical microstructures as a starting point as emphasized by Landauer [8] and numerous authors [9, 10, 11, 12]. [1.1.1.2] Digital three dimensional images of unprecedented size and accuracy have been prepared for the case of Fontainebleau sandstone, and are being made available to the scientific community in this brief report.

[1.1.2.1] Multiscale modelling of disordered media has recently become a research focus in mathematics and physics of complex materials and porous media [13, 14, 15, 16, 17, 18]. [1.1.2.2] Accurate prediction of physical observables for multiscale heterogeneous media is a perennial problem [9, 8]. [1.1.2.3] It requires knowledge of the three dimensional disordered microstructure [11]. [1.1.2.4] Our objective in this brief report is to provide to the scientific public a sequence of fully three dimensional digital images with a realistic strongly correlated microstructure typical for sandstone. [1.1.2.5] Resolutions from to are available for download [19]. [1.1.2.6] Experimental computed microtomographic images of comparable size, resolution or data quality are, to the best of our knowledge, not available at present. [1.1.2.7] More importantly, experimental images of similar size and quality are not expected to become available to the scientific community in the near future.

[1.1.3.1] Despite the impressive progress in fully three dimensional high resolution X-ray and synchrotron computed tomography of porous media in recent years [20, 21] acquisition times for 1500 radiograms needed for a -sample of average quality at the ID19 beamline of the European Synchrostron Radiation Facility are on the order of 30 minutes [22]. [1.1.3.2] Extrapolating to the number of 32768 such blocks, that we provide in this report, would thus require on the order of 16384 hours or roughly 2 years of uninterrupted beamtime. [1.1.3.3] It is unlikely that this amount of beamtime will ever be spent.

[1.1.4.1] The continuum multiscale modeling technology for carbonates developed in [23, 24, 25] was applied to Fontainebleau sandstone in [26], to create a synthetic, non-experimental image at very high resolution. [1.2.0.1] A laboratory sized cubic sample of sidelength cm was generated containing roughly polyhedral quartz grains. [1.2.0.2] For simplicity there are a total of 99 grain types each defined by eighteen intersecting planes. [1.2.0.3] The grains are rescaled, randomly oriented and have a prescribed overlap with each other (see [26] for more modeling details). [1.2.0.4] The model was geometrically calibrated against a well studied experimental microtomogram at resolution [27]. [1.2.0.5] The geometric calibration was based on matching porosity, specific surface, integrated mean curvature, Gaussian curvature, correlation function, local porosity distributions (with and measurement cell size), and local percolation probabilities at the same measurement cell sizes. [1.2.0.6] Comparison of these quantities at was carried out in [26].

[1.2.1.1] The continuum sample generated and characterized in [26] is the starting point for the work reported here. [1.2.1.2] To eliminate boundary effects a centered cubic sample, denoted by , of sidelength cm was cropped from the original deposit. [1.2.1.3] The pore space inside is denoted by , the matrix region is denoted .

[1.2.2.1] The sample region was discretized into cubic voxels ), whose sidelength is a multiple of the base resolution . [1.2.2.2] The discretization employs a cubic array of collocation points placed centrally and distributed uniformly inside each voxel . [1.2.2.3] It yields an integer gray scale value with for each of the voxels (discrete [page 2, §0]    volume elements). [2.1.0.1] The gray value equals the number of collocation points inside voxel , that fall into the matrix region [26]. The gray values approximate the integral

 (1)

which is the porosity inside the voxel at resolution . [2.1.0.2] Thus, approximates , while approximates .

 [h] 128 256 512 1024 2048 4096 8192 16384 32768 1 1 1 1 8 64 512 4096 32768 1.6 3.2 6.5 11 88 704 5632 45056 360448
Table 1: Overview of data sets available from http:www.icp.uni-stuttgart.demicroct for download. is the resolution in units of . is the sidelength of the cubic sample in units of . The total number of voxels is , and is the number of cubic blocks used to represent the sample. One has whenever . gives an estimate of the total CPU time in hours, that was used to generate the discretizations. Calculations were typically performed in parallel on 512 Xeon processors requiring a total wall time of approximately in hours due to an essentially linear speedup. Total wall time was around 800 hours or 34 days.

[2.1.1.1] At the lowest resolution the sample consists of discrete volume elements (voxels), and it can easily be stored in a single file. [2.1.1.2] At higher resolution this becomes technically inconvenient, and we have limited the file size to blocks with voxels corresponding to roughly 1 Gigabyte. [2.1.1.3] The highest resolution at which the sample is still stored in a single file is . [2.1.1.4] For resolutions the data are stored in several blocks. [2.1.1.5] This results in a maximum of 32768 blocks for the highest resolution of . [2.1.1.6] The number of voxels is at this resolution. [2.1.1.7] Table 1 gives a summary of the available resolutions , sample sidelength (in units of ), number of voxels , number of -blocks , and an estimate of the CPU time and wall time expended for the computation.

[2.1.2.1] All computations were performed on the HLRS’s bwGRID cluster at the Universität Stuttgart consisting of 498 compute nodes each holding two Intel Xeon CPU’s capable of 11.32 GFLOPs. [2.1.2.2] The peak performance (Linpack) is 38 TFLOP [28]. For the calculations performed in this work 256 nodes (=512 CPU’s) were used in parallel. [2.1.2.3] The discretization algorithm was parallelized, but not optimized. [2.1.2.4] Every discretized volume element (voxel) requires one byte of storage. [2.1.2.5] The storage requirements without compression amount to roughly 40 Terabytes.

[2.1.3.1] The results can be used to calculate resolution dependent geometric and physical properties. [2.1.3.2] Because the underlying continuum microstructure is available with floating point precision, the resolution can be changed over many decades. [2.1.3.3] Resolution dependent geometrical or physical properties can be compared with, or extrapolated to, the continuum result. [2.1.3.4] This is illustrated with the porosity (volume fraction of pore space) and specific internal surface (surface area per unit volume).

[2.1.4.1] The exact values of and are not known for the continuum model. [2.1.4.2] Depending on the geometrical or physical quantity of interest computations on the continuum model as well as computations on discretizations with more than discrete volume elements are (as a rule) challenging at present.

[2.1.5.1] The discretized samples approximate the exact geometry of the continuum model. [2.1.5.2] Let for and for . [2.2.0.1] Then

 (2)

is the number of voxels with grey value at resolution , and

 (3)

is an estimate for the porosity based on a segmentation of its discretization at resolution with segmentation threshold . [2.2.0.2] The porosity of the continuum model at infinite resolution is expected to obey

 (4)

for all and all . [2.2.0.3] The validity of these bounds depends on the continuum microstructure. [2.2.0.4] They are expected to hold approximately for the sandstone microstructure. [2.2.0.5] Grey values and do not guarantee, that the pore boundary does not intersect the voxel region, although it will be true for the majority of voxels in the limit .

[2.2.1.1] Figure 1 shows the porosity as a function of the segmentation threshold for resolutions .

[2.2.2.1] Figure 2 shows the porosity bounds and for resolutions . [2.2.2.2] It also shows eight extrapolations defined as the linear functions

 (5)

where and the slopes and intercepts are determined such that

 (6)

for and . [2.2.2.3] Table 2 lists the slopes and intercepts of all eight extrapolations. [2.2.2.4] In between the upper and lower bounds Figure 2 also shows (with symbols connected by a solid line) the porosity obtained by maximizing the variance [page 3, §0]    of the grey scale histogram, a standard technique in image processing [29].

-0.004944104133756 0.131677153054625 0.006128869055829 0.133806445926894
-0.005437994277600 0.133652713629999 0.006038194581379 0.13416914382469
-0.005680456774371 0.134137638623542 0.005981708264244 0.134282116458962
-0.005843877976029 0.134301059825200 0.005942813391024 0.13432101133218
Table 2: Slopes and intercepts for extrapolations in Figure 2 according to equation (5).

[3.1.1.1] Figure 3 shows the mean values

 (7)

plotted against the differences

 (8)

as data points. [3.1.1.2] It also shows a linear fit to the data as a dashed line. [3.1.1.3] The linear fit extrapolates to the value

 (9)

where the uncertainty is from the residues of the fit.

8 4 2 1 0.5
[mm] 20.039 20.530 20.774 20.904 20.946
[mm] 9.6371 9.8774 9.9978 10.061 10.065
[mm] 6.0407 6.1915 6.2673 6.3068 6.3096
[mm] 9.1680.02 9.8250.06 10.1590.16 10.355
Table 3: and are the upper and lower bounds computed from eqs. (16) resp. (17) and eq. (18) for specific internal surfaces as functions of resolution . is the mean field estimate computed from eqs. (18)–(23). is calculated by applying the direct algorithm from [30] to thresholded cubic subsamples, and averaging the results as described in the last paragraph of the main text.

[3.1.2.1] We now turn to an estimate for the specific internal surface . [3.1.2.2] Assume, that every voxel at resolution with is cut by a single plane. [3.1.2.3] Let denote a ball with radius centered at the center of the voxel . [3.1.2.4] Then is the inscribed ball with radius , and denotes the circumscribed ball with radius . [3.1.2.5] If the voxel is cut by a single plane, then the surface area of the cut is bounded above by the area of the circle forming the base of the spherical cap cut from the circumscribed ball . [3.1.2.6] Let denote the radius of the circle forming the base, let denote the height of the cap, and let be the radius of the sphere from which the cap is cut. [3.2.0.1] Then

 (10)

and

 (11)

is the volume of the spherical cap. [3.2.0.2] Introducing the voxel porosity

 (12)

and solving the last equation for and gives

 (13)

as a function of voxel porosity, where . [3.2.0.3] The inequality

 (14)

selects the solution with . [3.2.0.4] The area of the circle at the base of the spherical cap is so that

 (15)

is the relation between the circular base area and the volume fraction of a spherical cap cut from a sphere of radius . [page 4, §0]    [4.1.0.1] Combined with eq. (1) this yields the upper bound

 (16)

for the specific internal surface inside a voxel with resolution having grey value under the assumption that the voxel is cut by a single flat plane.

[4.1.1.1] For the lower bound we use the circular base area obtained by the intersection with the inscribed ball of radius . [4.1.1.2] The minimum is attained when the intersection plane is oriented perpendicular to one of the space diagonals of the cubic voxel. [4.1.1.3] The inscribed ball is intersected by the plane only when its distance from the nearest voxel corner exceeds . [4.1.1.4] This occurs at and at . [4.1.1.5] Therefore, one finds the lower bound for and , while

 (17)

holds for . [4.1.1.6] Upper and lower bounds for the specific surface are obtained from these results by summation as

 (18)

where b=min resp. b=max. [4.1.1.7] Table 3 lists the results for and for .

[4.1.2.1] Table 3 shows that the bounds for the specific surface are rather wide. [4.1.2.2] It is then of interest to investigate the following simple geometric estimate: In the limit most voxels are intersected by a single plane. [4.1.2.3] The intersecting planes can have varying orientations. [4.1.2.4] As a simple approximation it may be assumed, that the effect of different orientations is averaged out, and that a single orientation suffices to estimate the specific surface inside a voxel. [4.1.2.5] Here we assume this single direction to be a space diagonal of the voxel. [4.1.2.6] All four space diagonals are equivalent by symmetry so that it suffices to consider one direction. [4.1.2.7] A simple meanfield like appproximation is to assume, that every voxel with porosity given in eq. (1) is intersected by a single plane oriented perpendicular to one of the four space diagonals of the voxel.

[4.1.3.1] Let denote the distance from the corner of the voxel along the diagonal to the intersection point, where the plane intersects the space diagonal. [4.1.3.2] By symmetry it suffices to consider . For the intersection forms an equilateral triangle. [4.2.0.1] In this case is related to the voxel porosity by

 (19)

and the area of the triangle is

 (20)

[4.2.0.2] For the range the intersection forms a hexagon. [4.2.0.3] In this case the analogous results are

 (21)

and

 (22)

[4.2.0.4] The specific surface area of a voxel with grey value is then

 (23)

computed from eqs. (19) and (20) for or from eqs. (21) and (22) for with . [4.2.0.5] For , i.e. smaller than roughly , equations (19), (20) are used for , while (21), (22) are used for , but now in both cases with . [4.2.0.6] The second line in Table 3 is then obtained by inserting eq.(23) into eq. (18) with .

[4.2.1.1] The last line in Table 3 gives the specific internal surface calculated directly using the algorithm from [30]. [4.2.1.2] The data are obtained by averaging segmented -blocks. [4.2.1.3] The segmentation threshold was chosen by maximizing the variance of the grey scale histogram [29]. [4.2.1.4] For all blocks were included into the average, for half of the blocks, and for ten percent randomly chosen blocks were averaged in the measurement. [4.2.1.5] The calculation of took approximately 1100 hours CPU time (without segmentation), while the mean field estimate can be obtained (e.g. within Matlab) directly from the histogram at virtually no computational cost.