Next: Intrusion Results: Model Up: Main Previous: Introduction

# Description of Algorithm

The MP simulation algorithm has been developed in Fortran, for a supercomputer or workstation, and in C for an imaging computer. On the latter, the mercury flow can be followed visually as it progresses. A better understanding of the manner in which mercury intrudes a porous material can thus be obtained.

A model porous material that is instructive is represented in the computer as a set of solid circles digitized onto a 2-D square pixel array of dimension L x L. The circles are randomly "parked" [8] and are not allowed to overlap, with the additional constraint that each circle must be centered at a pixel. Parking means that once a circle has been placed, it is not allowed to move. Periodic boundary conditions are employed to minimize edge effects. This means that if part of a circle would extend outside the cell boundary, it is wrapped around to the other side. Our simulations were carried out with L = 500 and L = 1000. No significant differences were seen for the different values of L, so that by L=500 finite size errors were very small. Each pixel is classified as being either pore space (black) or solid material (white). Fig. 1 shows a 500 x 500 structure containing 360 21-pixel diameter circles, with a porosity of 0.5.

To begin the mercury intrusion simulation, the L x L box is surrounded with mercury by classifying all pixels within a fixed distance from the box edge as active mercury. These pixel locations are stored in an array as the active locations for mercury flow. This step corresponds to the initial placement of a sample into the sample holder and filling with mercury.

Intrusion is simulated by creating a flow field at each of the active mercury locations. Assuming a contact angle of 180 degrees between mercury and solid material, the flow field can be represented as a circle of a fixed diameter d. This value of the contact angle is not realistic for most porous materials. However, the exact value of the contact angle used in this work is not important, since according to eq. (1), the cos() term is only an overall scale factor when changing intrusion pressure into an equivalent pore size. The focus of this paper is the dependence of the threshold diameter, dc, on microstructural quantities, not on predicting exact pore size distributions, which cannot be done even with the correct contact angle. This diameter d is assumed to be inversely proportional to the applied pressure in a mercury intrusion experiment, via eq. (1). If a circle of diameter d can be centered at an active mercury location without overlapping any solid material, all pixels that it encompasses that are still empty pore space are reclassified as active mercury and their locations stored in a new array. Following this flow propagation step, the original active pixel is changed to deactivated mercury. After all flow fields from the original active mercury array have been processed, the procedure is repeated for the new active array. This iterative process continues until there are no new flow fields to explore, and the new active mercury array is empty. At this point all the area of the pore structure accessible by mercury from the exterior, for flow diameter d, will have been identified.

For a contact angle of 180 degrees, in two dimensions, this method of generating the intruding mercury cluster should give an accurate depiction of what an actual experiment would be like, for the following reasons. First, using circles to build the mercury cluster guarantees that the mercury meniscus will touch solid particles at the correct contact angle, and should give a reasonably accurate shape between particles. Even if the true contact angle was 120-140 degrees, the meniscus would not change that much. Second, all the mercury starts from the outside, and intrudes inwards on all sides, thus duplicating how actual mercury intrusion experiments are carried out.

To check for the continuity of the mercury cluster, from top to bottom or left to right, a burning algorithm is employed [9]. Conceptually, the burning algorithm can be thought of in the following way. To check continuity from left to right, for example, the mercury pixels are regarded as being combustible. Mercury pixels along the left edge are set on fire, and the fire is allowed to propagate to neighboring mercury pixels until no more pixels can be burned. If any mercury pixels on the right edge are found to have been burned, then the mercury cluster must extend completely across the sample. The burning algorithm is simple to program, and requires little execution time.

By executing the mercury intrusion algorithm for a range of values of d, complete simulated mercury porosimetry curves can be generated. The value of dc, the threshold diameter for mercury continuity, can be directly determined using the burning algorithm. This computer simulation algorithm has some interesting differences with experiment. In a mercury intrusion experiment, the pressure is directly measured and the intruded pore diameter interpreted using eq. (1). In our simulation, the pore diameter being intruded is fixed by our choice of the diameter of the mercury circle we try to place. There is no need to even use eq. (1), or worry about what is the appropriate contact angle, since intrusion pressure does not come into the algorithm at all.

It should be noted that the extension of this MP algorithm to three dimensions would be easy computationally, but difficult theoretically. In two dimensions, choosing a circle to build up the flow field of the intruding mercury is appropriate and accurate, as the mercury meniscus must be circular in the completely non-wetting limit, since there is only one radius of curvature to be fixed by the mercury's surface tension and the applied pressure. However, in three dimensions it is the sum of two radii of curvatures that is fixed by the surface tension and applied pressure, so that the meniscus shape does not have to be spherical. So if mercury spheres were used to implement the MP algorithm in three dimensions, the volume intruded at a given pore diameter would only be a lower bound, as a pore that would not admit a sphere of that diameter might very well admit a different shape that also satisfied the meniscus curvature equation. However, the information from such an algorithm might well still be useful, but only in the sense of a lower bound. Allowing for ellipsoidal shapes in the algorithm would be very difficult.

Next: Intrusion Results: Model Up: Main Previous: Introduction