Next: Microstructure and properties Up: Main Previous: Introduction

Computational and Experimental Techniques

The microstructure reconstruction algorithm has previously been fully described [10], so that only a brief treatment is given here. The 3-D digital image structure generated is limited to 100³ or 200³ pixels in volume, with periodic boundaries employed to reduce surface errors. Larger images can be generated, but the programs that compute quantities like elastic moduli, fluid permeability, and electrical and thermal conductivities, are much more computer memory-intensive, and so 100³ - 200³ are the largest images used at present. We use a filtering function, applied to a 3-D random noise image, that is based on correlation analysis of a real 2-D image of the porous glass microstructure. Algorithms with some similarities to ours have been previously used for porous Vycor [11], which employed the filtering of a random white-noise image with the convolution of a Laplacian and a Gaussian function. Alternate generation techniques for porous Vycor based on spinoidal decomposition [12], an "intersecting growth mechanism" [13], and phase separation in a three-dimensional Ising model [14] have also been previously presented, although not with the extensive testing of properties given in this paper.

Starting with a 2-D binary image of porous Vycor [15], the 2-D autocorrelation function for the pore phase is determined. If the M x N pixel 2-D image is defined as a discrete valued function, I(x,y), where I(x,y) is equal to one for pore pixels and zero for solid pixels, the two-point correlation function for the porosity is given by [16]

The function S(x,y) is then converted to its polar form and averaged over all orientations to obtain S(r). Finally, the autocorrelation function, F(r), for the porosity is defined by [17]

By definition, F(0) = 1 and F(r) 0 as r , since S(r) [S(0)]² as r . For an isotropic, random material like porous Vycor, F(r) is nearly zero for r on the order of a few times the largest pore size.

Since for a large enough image of an isotropic material, this autocorrelation function will be the same in two and three dimensions, it can be used to filter a 3-D image of white noise [18]. For a finite size image, there will of course be some slice-to-slice fluctuation. For the size image considered here, these fluctuations were negligible. After the filtering process, a gray scale image is obtained. A threshold gray scale is then chosen, such that all pixels with gray level above and below this threshold are set to either solid or pore. This threshold gray scale is such that a 3-D binary image is created in which the porosity matches that of the original 2-D image.

Finally, a rearrangement algorithm based on analysis of local curvature [19,20,21] is used to alter the hydraulic radius (the ratio of pore volume to pore surface area) of the computed 3-D structure. Usually the hydraulic radius is adjusted so as to match that measured for the original 2-D image [10], providing a 3-D computer representation of the material microstructure with the correct porosity and pore surface area. However, it is known that porous Vycor glass is a variable material, with a variable pore size and surface area. In the literature, values for the nitrogen BET surface area range from 100 m²/g [15] all the way up to 173 m²/g [5]. We have had to generate 3-D images with different surface areas in order to be able to accurately compare to experiments, as many (but not all) properties are sensitively dependent on surface area, or, equivalently, since the porosity is fixed, on hydraulic radius and pore size. This may be easily done using the surface rearrangement algorithm. A portion of a computed 3-D microstructure of the porous Vycor glass is provided in Fig. 1, with a surface area of 70 m²/g. It is important to note that the scale of the original TEM picture used [15], and therefore the scale of the resulting generated 3-D image, is approximately 1 nm per pixel, so that features smaller than this have been smoothed out and lost [8]. To convert the measured specific surface area in the model from m^-1 to m²/g, we are assuming that the density of the backbone is approximately that of fused silica, 2.18 10⁶ g/m³ [8].

Figure 1: Reconstructed 3-D image of porous Vycor glass. Pores are shown in white or grey while solids are transparent.

Once a three-dimensional representation of the solid and pore phases in the material has been generated, the next step is to determine the location of the capillary-condensed water as a function of relative humidity, given that there is an adsorbed water layer on all solid surfaces. To analyze the capillary water content, a geometrical algorithm [1] is employed that is based on the Kelvin-Laplace equation:

For desorption, the above geometrical constraint must be augmented by a connectivity analysis [1]. During desorption, a capillary pore can empty only when its size equals or exceeds the Kelvin-Laplace radius and when it is connected to the exterior of the sample by a pathway consisting exclusively of the same size or larger "empty" pores. Including this connectivity analysis in the simulation of sorption isotherms produces a natural hysteresis in the desorption-adsorption curves similar to that observed experimentally [1]. This geometrical analysis is applied to the porous Vycor microstructures to produce a series of three-phase composites consisting of solid, air-filled pores, and water-filled pores.

As was mentioned above, in order to carry out the sorption computation properly, especially for the smaller pores, it is necessary to first determine the thickness, t, of the adsorbed water layer as a function of relative humidity. Here, we use an empirical equation developed for cement-based materials [25]

Two other t(RH) functions in the literature are actually quite similar to this one [26,27], so that the choice of function would not make a large difference in the results. In eq. (4), the thickness of the water layer is given in nanometers. Computationally, for a given relative humidity, a sorbed water layer n pixels thick (n=0 or 1) was first placed on the solid surfaces (using a dilation algorithm) before executing the geometrical algorithms for sorption described above. Non-integer pixel adsorption thicknesses, x, determined using eq. (4), were approximated in a two step approach. First, the integer value of x, int(x), dilations were executed, followed by the geometrical sorption algorithm. After sorption, the total amount of water present was incremented by adding the value obtained by multiplying the exposed (to air) surface area remaining after the capillary condensation analysis by the factor [x-int(x)] [1].

There are some difficulties in using the Kelvin-Laplace equation as described above. Using molecular dynamics and density-functional theory studies of adsorption in a single pore, it was found that eq. (3) starts to break down significantly at pore sizes of about 5-10 nm [28,29, 30]. This is the typical pore size in porous Vycor. However, by first considering a water layer on the pore surfaces, before estimating capillary condensation with eq. (3), we are in effect using a modified Kelvin-Laplace equation, which has been shown to work somewhat better at smaller pore sizes than the unmodified Kelvin-Laplace equation [29]. Also, we desire to compute adsorption and desorption on a general pore structure represented by a digital image, not on a single model shape pore. The computational burden would be enormous if Monte Carlo [28] or density functional theory [31] were used to more correctly handle the actual fluid molecular state in pores less than 10 nm. in width. In the real pore space that we want to consider, there are no nicely-shaped (slit, cylinder, etc.) pores on which to operate these more powerful but more complicated techniques [32]. The simple but approximate geometric technique outlined above can equally well be applied to any type of pore space, and is particularly suited to digital images, which is usually what one has to deal with in looking at real systems. There is also evidence that the main cause of the hysteresis loop between adsorption and desorption is percolation effects due to the pore network, not the size and shape of individual pores [30,33, 34]. In that case, the topology of the pore network will tend to wash out individual pore details that are treated somewhat incorrectly by our algorithm. More discussion of this question will be delayed until Section 4.

The sorption isotherms were measured using a hygrogravimetric method [35]. The mass of the sample was continuously recorded in a chamber in which the temperature and the relative humidity were controlled. The required partial pressure of water vapor was produced by mixing dry air and water-saturated air. For adsorption, starting from a relative humidity of around 3%, the relative humidity was increased in a stepwise manner. As soon as the mass of the sample had reached a constant value, the relative humidity was increased (or decreased for desorption). The relative humidity range of the apparatus was 3 to 98%. It should be noted that differences which exist between experimental curves published in the literature are probably at least partially due to the variability of the porosity and pore structure of porous Vycor glass from one sample to another, and not only differences in experimental techniques. Computationally, we have noticed large differences in sorption/desorption for different surface area models, as the hydraulic radius, for a fixed porosity, acts as a length scale for the pore sizes in the pore space, which of course will greatly affect the amount of water present at a given RH value. The surface area of the experimental sample used was not measured, but was probably somewhere between 100 and 200 m²/g, typical for porous Vycor.

To compute the shrinkage due to capillary condensed water, it is necessary to specify the stresses present in this water in the microstructure. The hydrostatic pressure or capillary tension in the capillary-condensed water is given by the simple expression

where for a given relative humidity, K is calculated using eq. (3).

The shrinkage due to the surface adsorbed water (Bangham shrinkage [36]) is based on the specific surface free energy (J/m²) change which occurs as the thickness of this adsorbed layer changes. This change is given by the Gibbs adsorption equation in the form [37]

(6)

where P is less than P_o. The zero point for surface free energy can be taken to be at 100% relative humidity, in order to convert γ_s to γ_s. The zero point is arbitrary, as only the change in surface energy is physically meaningful. The expression for t given in eq. (4) was substituted into eq. (6) and the integral then computed to give:

(7)

Fig. 2 provides a plot of γ_s vs. relative humidity determined using eq. (7). In Fig. 2, at zero RH, the surface free energy is about 0.2 J/m². In Iler [38], the surface energy of a fully hydroxylated silica surface is given as 0.13 J/m², while for a fully dehydroxylated surface the surface energy is about 0.26 J/m^2.. For the porous Vycor material being considered, the value of 0.2 J/m² falls in between these values, which is a reasonable result, as the surfaces in Vycor are probably not at either of these two extreme hydroxylation points.

Figure 2: Computed change in surface energy as a function of relative humidity from eq. (7) in text.

Once a model microstructure has been generated, it can then be analyzed using a finite element technique especially adapted to digital images. In this technique, each cubic pixel in the 3-D image is taken to be a tri-linear finite element, which results in a finite element problem with as many nodes as there were pixels in the original image. This is a large computational problem, too large for a direct solver like a Gaussian elimination method, and so requires a relaxation solution method like the conjugate gradient algorithm [39].

To find the elastic moduli of a model microstructure, periodic boundary conditions are applied that incorporate a fixed applied strain. The numerical solution (minimum elastic energy) for the displacements is found [40], and the elastic moduli are simply given by the average stress found in the problem, in the usual way of composite theory [41].

Computing the shrinkage due to internal forces is somewhat more complicated. In this case, the overall size and shape of the system must be made variables, so that the energy of the system depends not only on the displacements in each element, but also explicitly on the size and shape of the unit cell of the system. In response to internal forces that tend to make the material shrink (or expand), the size and shape of the unit cell can then adapt itself so as to minimize the elastic energy and find the correct solution for the overall shrinkage. The internal forces in a porous medium come from the negative pressure in the capillary water, which acts on the solid surfaces with which it is in contact, and the changes in the specific surface energy in the adsorbed water layers.

The effect of the negative pressure in the capillary water is computed by applying the appropriate stress directly to the nodes on the capillary water:solid surface. To handle the effect of the specific surface energy on the system size, a term is added to the system energy of the form

where S_A is the surface area in square meters of the exposed surfaces, and γ_s is the specific surface free energy (J/m²). As before, an exposed surface is a surface that is not covered with capillary water. The surface area of the digital structure can be written out in terms of a power series in the displacements of the nodes. Since only small strains are being considered, the series is truncated at linear order. The result is an additional term in the finite element elastic energy that is linear in the nodal displacements. The system shape and size is then modified to minimize this overall energy. For example, for an unstrained system, with no capillary water, the system will shrink, thereby reducing the surface elastic energy by reducing the surface area at the expense of increasing the volume elastic stored energy by the same amount [8].

An important test case for the treatment of surface energy is that of the dilute limit of spherical holes in a solid frame, with a specific surface energy term applied to the surface. This problem is easily solved using spherical polar coordinates [42,43,44]. One applies a uniform radial displacement U, which would arise due to the application of the surface energy, to the surface of a spherical hole embedded in the center of a large solid sphere of elastic moduli K_s and G_s. The resulting volume elastic and surface free energy is minimized with respect to U, and the new size of the system computed. For the dilute case, the result for the linear strain ε that arises from the application of the specific surface free energy γ_s is

where c is the porosity and a is the radius of the hole. In this case, we are starting from γ_s = 0, so actually γ_s = γ_s. For a 15 pixel diameter spherical hole in a 40³ unit cell, the agreement with the exact result is 4.7%.

Another test case can be found in Scherer's paper [8], the shrinkage per unit surface energy in an intersecting cylinder model. This consists of a simple periodic cubic structure, with equal diameter circular cylinders oriented along the x,y, and z directions, which pass through each other. The radii of the cylinders is a, and size of the unit cell is L. The quantity x=a/L is determined by the porosity of the model, and the specific surface area is determined by the size of L in real units. Detailed formulas can be found in Ref. [8]. A test case was set up using a 60³ pixel unit cell, with 41 pixel diameter cylinders. The solid had K_s=20.44 GPa, and G_s=16.30 GPa. The theoretical prediction for the shrinkage per unit surface energy was -1.24 · 10⁻³, while the finite element computation gave -1.17 · 10⁻³, a difference of only 5.6%.

Mackenzie [6] gives a formula for the dimensional change of a fully saturated porous medium, which has a pressure p_cap in the pore fluid. This equation is exact when the solid frame is linearly elastic, and the pores are fully saturated. If K is the bulk modulus of the porous solid, and K _s is the bulk modulus of the material making up the solid frame of the porous solid, then the linear strain ε is given by

The value of K is for empty pores, which assumes that the bulk modulus of the water has no effect. This will be the case when the bulk modulus of the water is much less than that of the solid frame, and/or if the water is free to flow in and out of the porous material as the size of the material changes. For negative pressure (hydrostatic tension), the strain is negative, implying shrinkage. For a partially saturated porous solid, the average of the capillary pressure over the pore space can be used instead of p_cap, which brings in a factor of S, the saturation. The saturation is defined as the fraction of the pore space volume that is filled with capillary water. The resulting approximate shrinkage equation for partially saturated porous solids is:

determined from the dilute elastic moduli limit of a spherical hole [45],

where c is the porosity. The algorithm gives agreement with the exact result within 1-2%. In eq. (10), only the effective bulk modulus needs to be known in order to know the strain that arises due to a pressure in all the pores. An analogous result is that of Hashin and Rosen [46], where for a two-phase random material with arbitrary thermal expansion coefficients and moduli, only the effective bulk modulus needs to be known in order to know the effective thermal expansion coefficient exactly.

We have found that eq. (11) holds very well for saturations of 80% or more, but decreases in accuracy for lower saturations. To actually compute the capillary-driven shrinkage for partial saturations, it was necessary to surround the system with a one half pixel thick wall, as shown schematically in 2-D in Fig. 3. Periodic boundary conditions were impossible to accurately maintain without this wall. Without the wall, numerical instability seemed to arise when the system was only partially saturated and periodic boundary conditions were used. The wall was not necessary for computing shrinkage at saturation. The elastic moduli of this wall were chosen so that when the pore space was fully saturated, the exact result according to eq. (10) was found. A value of 1/5 the effective moduli of the entire system was used, which was independent of the size of the system in pixels. The choice of the wall moduli compensated for the small increase in surface area caused by the wall touching the pore space, therefore preserving the accuracy of the computational method.

Figure 3: Schematic picture of wall construction (2-D) for partially saturated systems. Dashed lines show periodic cell boundaries, black is wall, gray is material.

So for both shrinkage mechanisms, negative capillary pressure and changes in specific surface free energy, the finite element codes give accurate evaluations of exactly known dilute limits, and so should therefore give reasonably trustworthy results for the more complicated porous Vycor glass microstructures. Finite difference codes were also applied to compute electrical and thermal conductivity, and fluid permeability. The details of these algorithms have been described elsewhere [10,44].