4836
Ind. Eng. Chem. Res. 2001, 40, 4836-4843
Porous Media Characterization Using Mercury Porosimetry Simulation. 2. An Iterative Method for the Determination of the Real Pore Size Distribution and the Mean Coordination Number V. G. Mata, J. C. B. Lopes, and M. M. Dias* Laboratory of Separation and Reaction Engineering, Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, 4200-465 Porto, Portugal
A methodology for the determination of the actual pore size distribution and the mean coordination number of a porous sample using experimental mercury porosimetry data together with a porosimetry simulator is developed. The actual intrusion and extrusion curves are used and optimization by means of an iterative method is performed. The method is validated against computer-generated intrusion and extrusion curves using the porosimetry simulator and applied to an oil-bearing rock sample as a test case. Introduction Mercury porosimetry is a technique frequently used for the determination of the pore size distribution (PSD) of porous materials. Interpretation of its data is generally done using the parallel pore model (p-p-m).1,2 This model ignores the existence of pore interconnectivity, resulting in erroneous results for the derived PSD and the inability to interpret the hysteresis phenomenon for mercury invasion and extrusion curves, sometimes referred as the “ink bottle” effect.3 An early attempt to determine the actual distribution of pore sizes in a porous rock from porosimetry experiments was based on probability concepts.4 To understand the effects of pore interconnectivity, mercury porosimetry simulators based on pore network models have been developed. In the pioneer work of Androutsopoulos and Mann,5 a 2D square network model was developed and used for mercury porosimetry simulations. This simulator was later used for the determination of a more accurate estimate of the PSD for an oilbearing rock and for a catalyst pellet.6,7 From the experimental intrusion curve, a histogram of 5% volume segments was constructed. To identify the real distribution function, a trial-and-error method on a square 2D network with 220 cylindrical pores ranging from 100 to 40 000 Å in size was used. The pore sizes were then grouped in 10 groups, each consisting of a constant number of segments. To each group of segments was allocated a pore size interval within the global size range. The size of each interval was then continuously adjusted until coincident simulated and experimental intrusion curves were achieved. A constant connectivity value of 4 was used in all cases. Large differences between the PSD predicted from the parallel pore model and the PSD obtained from the trial-and-error method were reported, mainly in the larger pores region. This trial-and-error method associated with a mercury porosimetry simulator was used more recently with 3D pore networks to infer the PSD of packing materials used in perfusive chromatography.8 In this work, the PSD was made up of three Gaussian distributions to represent the three types of pores present in the * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: 351-22 508 1661. Fax: 351-22 508 1674.
column: interstitial pores and two kinds of intraparticle poressmacropores (diameters on the order of 1 µm) and micropores (diameters on the order of 0.1 µm). The distribution parameters and the connectivity were systematically changed to adjust the experimental and the simulated intrusion curves. The extrusion curve was not studied. A different method for determining connectivity from experimental and simulated mercury porosimetry curves was developed by Portsmouth and Gladden,9,10 using a 3D spherical network model. This method is based on the generation of secondary mini-hysteresis porosimetry loops, which are generated on both the overall intrusion and extrusion curves. The value of the connectivity is obtained by an analysis of the retained mercury as a function of the network coordination number. The volume of mercury that re-enters the pore structure during the secondary intrusion process is summed to produce a cumulative secondary intrusion curve, providing an improved estimate of the PSD. Payatakes and co-workers have worked extensively on this type of problem and developed and studied a mercury simulator based on a chamber/throat network.11,12 The results were later compared with mercury intrusion and retraction data in model porous media.13 A method for the determination of topological and geometrical parameters was then proposed and used to infer the pore size distribution and other parameters for typical sandstones.14 The porosimetry curves were fitted to analytical functions and estimation of the parameter values was performed by fitting the simulated porosimetry curves to the experimental ones. In part 1 of this work, a mercury porosimetry simulator based on a 3D network of cylindrical pores was described. The influence of the network parameters, such as the network size, the mean coordination number, the accessibility, and the pore size distribution, on the simulated porosimetry curves was studied. It was shown that the most significant parameters are the mean coordination number and the pore size distribution. The effect of the accessibility, described by REI, the ratio between the number of external pores to the number of internal pores, was shown to have a lesser effect on the simulated porosimetry curves. For this
10.1021/ie0101137 CCC: $20.00 © 2001 American Chemical Society Published on Web 10/09/2001
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4837
reason, in the following, the value of REI was kept constant throughout all simulations. The objective of part 2 of this work is the development of a methodology for the determination of the actual PSD and the mean coordination number using experimental mercury porosimetry data together with the porosimetry simulator and performing an optimization by means of an iterative method.15 The method takes the actual intrusion and extrusion curves and uses the raw data to infer the Real-PSD and the mean coordination number. Description of the Iterative Method The iterative method can be divided into four main steps. I. The first step is the determination of the ApparentPSD, that is, the PSD obtained using the parallel pore model on mercury porosimetry intrusion data. II. In the second step, a value for the mean coordination number, C, defined as the average number of pores that are connected to one node, is assumed. III. The third step is itself an iterative methodology. The Apparent-PSD is used as a first iteration on the mercury porosimetry simulator based on a network of interconnected channels with the assumed value of C. A series of iterations is performed by successively modifying the network PSD until the simulated intrusion curve converges to the experimental intrusion curve. IV. In the fourth step, the extrusion curve is simulated using the PSD obtained in the previous step. The simulated extrusion curve is compared with the experimental extrusion curve, and if they agree to within some previously defined criteria, this PSD is taken as the actual pore size distribution (Real-PSD) and the assumed value of C as the actual value of the mean coordination number. If not, step II is repeated with a different value for the mean coordination number, followed by steps III and IV. Here, the notation follows that used in part 1 of this work. For both the diameter probability density function, f(di), and the diameter probability distribution function, F(di), subscript N refers to the number-based function, and subscript V refers to the volume-based function, where di represents the pore diameter. The PSD is described in terms of the number-based diameter probability density function, f(di), or in terms of its discretized form, f(j) for j ) 1, 2, ..., J diameter classes. To validate this method, simulated porosimetry curves were obtained using a known Initial-PSD, f in N (di), and a network with a given value of the mean coordination number, Ci. The simulated intrusion and extrusion curves were used as “experimental” curves. A Gaussian distribution function with a mean value of µ ) 750 µm and a standard deviation of σ ) 250 µm was used for the determination of the pore diameters. The InitialPSD, f in N , is determined from randomly generated pore diameters, truncated at a minimum value of dmin ) 0 and a maximum value of dmax ) 1500 µm. Figure 1a shows the experimental intrusion, FV, and extrusion curves simulated on a network with a mean coordination number of Ci ) 6. The parallel pore model (p-p-m) curve, that is, the intrusion curve that would be obtained if all pores were in direct contact with the external mercury, is also shown. The Initial-PSD, f in N, is shown in Figure 1b.
Figure 1. Simulation with a Gaussian Initial-PSD on a network with Ci ) 6. (a) Experimental intrusion, FV, and extrusion curves. ap ap (b) Initial-PSD, f in N ; Apparent-PSD expressed as f V and as f N .
Next a detailed description of each step is given, together with validation of the method. Step I. Determination of the Apparent-PSD. Experimental porosimetry data are usually given in terms of the cumulative intrusion volume curve from which the volume-based diameter probability distribution function, FV(di), can readily be obtained. Differentiation of FV(di) results in the volume-based diameter probability density, f ap V (di). The Apparent-PSD, expressed in terms of the number-based diameter probability density function, f ap N (di), is obtained using the relationship
fN(di) )
fV(di) di
∫d
2
dmax min
(fV(ξ)/ξ2) dξ
(1)
This function is assumed as the starting point, I ) 0, of the iteration process, that is, fN|I)0 ) f ap N , where I is the iteration number. Figure 1b shows the Apparent-PSD ap expressed as f ap V and as f N obtained from the experimental intrusion curve, FV(di). Step II. Mean Coordination Number. In this step, a value for the mean coordination number, Cf, is assumed. Step III. Iteration on the Intrusion Curve. The objective of this step is to obtain the PSD for which the simulated intrusion curve matches the experimental intrusion curve. This is done by an iterative procedure, with each iteration I involving the following steps: (1) NT random numbers are generated from fN|I-1 and are randomly assigned to the network pore diameters.
4838
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001
(2) The intrusion simulation is executed, and the is determined; the simulated intrusion curve FV|sim I sim and f | are obtained as described in functions fV|sim N I I step I. - fN|I)0, is (3) The difference function, ∆fN|I ) fN|sim I determined. (4) A new density function, fN|I ) fN|I-1 - 1/k ∆fN|I, is obtained, where k is a relaxation constant, usually set to k g 2. The iterative process stops when the simulated intrusion curve, FV|sim I , converges to the experimental intrusion curve, FV(di). The convergence criterion, eV, that quantifies the difference between the experimental and the simulated intrusion curves and is normalized by the maximum deviation that occurs for I ) 1 and defined as J
eV )
2 sim [Fexp ∑ V (j) - FV|I (j)] j)1 J
(2)
2 sim [Fexp ∑ V (j) - FV|1 (j)] j)1
where FV(j) stands for the discretized form of FV(di) as described in part 1 of this work for j ) 1, 2, ..., J classes. In this work, convergence was assumed to be reached when eV < 1%, and the iteration at which this occurs is called iteration N. The new density function obtained in this iteration, fN|I)N, is referred to simply as f N N. At this point, the validation of the iterative procedure on the intrusion curve is performed. As described above, the simulated intrusion and extrusion porosimetry curves obtained using a known Initial-PSD, f in N (di), and a network with a given value of the mean coordination number, Ci, are used as experimental curves. Taking the Apparent-PSD obtained in step I as the starting point, the iterative process is first applied using a network with a known mean coordination number, that is, Cf ) Ci ) 6. Figure 2 shows the evolution of the comparison of the apparent density function, fN|0 ) in f ap N , and the initial density function, f N , with the sim , for simulated and modified functions, fN|I and fN|mod I iterations, I ) 1, 5, and 10. It is clear that the process converges and the modified density function approaches the initial one. This process converged for N ) 25, and the density function for the final iteration, f N N, is shown in Figure 3, together with the initial and the apparent density functions. Figure 3 also shows that N the final intrusion curve, FN V , computed from f N is coincident with the experimental intrusion curve, FV, obtained from f in N . Furthermore, the simulated extrusion curve resulting from f N N is also coincident with the experimental extrusion curve. Just for comparison purposes, the simulated intrusion and extrusion curves obtained directly from the f ap N are also shown as dashed lines. In most practical cases, the connectivity of a porous sample is not know a priori. Figure 4 shows the application of the same iterative procedure using the Apparent-PSD obtained for a network with Ci ) 6, but now using a different value for the mean coordination number, Cf ) 3. The iterative process also converges after N ) 25 iterations to an intrusion curve, FN V , that is coincident with FV, but now the final density function, in fN N, is clearly different from f N . This shows that differ-
Figure 2. Comparison between the apparent density function, in fN|0 ) f ap N (dashed line), and the initial density function, f N (dotted line), and the simulated and modified functions, fN|sim I and fN|mod , for iterations, I ) 1, 5, and 10 (solid lines). I
ent pore sizes distributions can result in similar intrusion curves if the connectivity is different, so the result of the first iterative process is not unique. On the other hand, one observes that the final and experimental extrusion curves are now different, leading to different values of the dimensionless residual mercury Vres Hg . This observation suggested a new methodology for the determination of the mean coordination number, C, as described below. Step IV. Iteration on the Extrusion Curve. After the final iteration on step III, an extrusion simulation is performed. The value of Vres Hg obtained from the final simulation is compared with the experimental value. If N these two values are similar, this means that f re N ) f N and the Real-PSD has been obtained, together with the actual value of the mean coordination number. If not, a new value of Cf must be chosen and the iteration
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4839
Figure 3. Iterative method carried out using the Apparent-PSD of Figure 1 on a network with the same mean coordination number, Cf ) Ci ) 6. (a) Initial-PSD f in N (dotted line); final ap density function, f N N (solid line); Apparent-PSD, f N (dashed line). (b) Intrusion and extrusion curves: experimental (dotted line), final (solid line), obtained from f ap N (dashed line).
procedure on the intrusion curve must be repeated until the two values of residual mercury coincide. To obtain some insight on how to proceed with the choice of Cf, a sensitivity analysis of Vres Hg as a function of Ci and Cf is shown in Figure 5. Here, V*Hg ) Vres Hg |Cf/ res | represents the ratio between V | , the value of Vres C C Hg i Hg f the residual mercury obtained in the final simulation in the iterative method with mean coordination number of Cf, and Vres Hg |Ci, the value of the residual mercury in the experimental extrusion obtained with Ci. Figure 5 shows that the curves present a maximum for 3 < Cf < 4, and thus, the iterative method might result in two possible solutions for the Real-PSD and the mean coordination number, C. In these cases, it might be necessary to obtain independent information on the value of C using other methods such as serial tomography sectioning and 3D reconstruction and analysis. Application to the Test Case Mann and co-workers used a 2D network simulator together with a trial-and-error method to infer the PSD of an oil-bearing rock.6 The porosimetry data for both intrusion and extrusion curves, shown in Figure 6a, were extracted from Mann’s paper and used in this work as a test case for the iterative method. The ApparentPSD expressed in terms of the volume-based probability density function, f ap V , obtained by direct differentiation of the intrusion curve and the corresponding numberbased probability density function, f ap N , obtained using
Figure 4. Iterative method carried out using the Apparent-PSD of Figure 1 on a network with a different mean coordination number, Ci ) 6 and Cf ) 3. (a) Initial-PSD, f in N (dotted line); final ap density function, f N N (solid line); Apparent-PSD, f N (dashed line). (b) Intrusion and extrusion curves: experimental (dotted line), final (solid line).
res Figure 5. Plot of V*Hg ) Vres Hg |Cf/VHg |Ci as a function of the assumed mean coordination number, Cf.
the methodology described in part 1 of this work, are shown in Figure 6b. Because the range of diameters covers about 4 orders of magnitude, these functions were discretized into 40 classes, equally spaced on a logarithmic scale. The f ap V distribution shows three distinct peaks: the largest one at d ≈ 0.02 µm and two smaller ones at d ≈ 0.1 µm and d ≈ 1 µm. In the f ap N distribution, the first peak is also present, but the second peak is dislocated to d ≈ 0.08 µm, and the third peak does not appear. The mean pore diameters for the volumeand number-based distributions are 1.1 and 0.03 µm, respectively.
4840
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001
Figure 6. Experimental porosimetry data for an oil-bearing rock sample.7
Because the coordination number is not known, the iterative method was applied to networks with mean coordination numbers of Cf ) 3, 4, and 6. The iteration results are shown in Figure 7, where large differences between the experimental intrusion curve, FV, and the intrusion curve (iteration 1) obtained by simulation using the Apparent-PSD, f ap N , are clear. The application of the iterative procedure on the intrusion curve (step III) in all three cases results in the convergence of the simulated intrusion curves to a curve close to the experimental intrusion curve after 500 to 600 iterations. From the previous results, Figure 8 shows the final number- and volume-based probability density funcN tions, f N N and f V , respectively, for each value of Cf, as well as the corresponding Apparent-PSD curves, f ap N N ap and f ap V . The major difference between f N and f N in all three cases is the emergence of a peak at d ≈ 1 µm, showing that a number of large-diameter pores had not been accounted for in the Apparent-PSD. The volumebased distributions, f N V , calculated from the correap sponding f N N are very different from f V ; the peak at d ≈ 0.1 µm is still present, but it is much smaller, and a large amplitude peak is now present at d ≈ 0.2 µm. This shows that the presence of even a few large pores produces large differences in the volume-based probability distribution function, as most of the volume is now associated with these diameters. The f N N distributions for the three values of Cf are similar, but some of the differences are important. All three curves show a peak between d ) 0.02 µm and d ) 0.03 µm, but the height of this peak increases with Cf from a value of less than 100 for Cf ) 3 to almost 200 for Cf ) 6. For intermediate values of d, a clear peak at around d ≈ 0.1 µm appears for Cf ) 3 and Cf ) 4, whereas for Cf ) 6, no such peak is present. The peak for larger values of the diameter shifts from a broad peak for Cf ) 3
Figure 7. Experimental (dotted line) and simulated (solid lines) intrusion and extrusion curves for and oil-bearing rock sample obtained for various iterations on the intrusion curve with (a) Cf ) 3, (b) Cf ) 4, and (c) Cf ) 6.
located between d ) 0.2 µm and d ) 3.0 µm to a narrow peak centered around d ) 1.5 µm. The Real-PSD and the actual mean coordination number, C, are chosen by comparison of the extrusion curves from the last simulation, plotted in Figure 7. The simulated extrusion curves are similar for all three values of the coordination number, showing that mercury starts to be extruded at pressures corresponding to d ≈ 0.2 µm, whereas the experimental curve extrusion starts at lower values of pressure for which d ≈ 1.5 µm. The dimensionless residual mercury, Vres Hg , is plotted in Figure 9 as a function of Cf. Judging from these results, it seems clear that the actual value of coordination number is around C ) 4. The mean pore diameters for the volume- and the number-based distributions are 1.9 and 0.38 µm, respectively. These
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4841
Figure 8. Apparent (dashed line) and final (solid line) pore size distribution curves for and oil-bearing rock sample obtained by application of the iteration on the intrusion curve process with (a) Cf ) 3, (b) Cf ) 4, and (c) Cf ) 6.
results are in agreement with the values reported by Mann and co-workers,6 who used a 2D square network (coordination number equal to 4) and concluded that most of the pore volume is in pore diameters larger than 0.2 µm and that the mean network diameter (numberbased) is 0.3 µm. Figures 5 and 9 show that, depending on the type of sample under study, very close values of Vres Hg can be obtained in different ranges of Cf. In these cases, the determination of the actual coordination number and PSD can be solved by using additional techniques, such as serial tomography, to infer an independent value of the mean coordination number. Conclusions A systematic methodology for the determination of the actual PSD and the mean coordination number of a
porous material has been developed. This method uses experimental mercury porosimetry data together with a porosimetry simulator and performs an optimization using an iterative method. This method uses not just the information from the mercury intrusion curves, but also extrusion data to infer the mean coordination number. The method was validated by means of simulated intrusion and extrusion curves obtained from networks with known initial pore size distributions and mean coordination number. The method was applied to a test case: data obtained by Mann and co-workers6 for an oil-bearing rock was treated, and the resulting apparent pore size distribution used to infer the real pore size distribution and the mean coordination number. Good agreement between the experimental and simulated intrusion
4842
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 j ) discretization class J ) number of discretization classes k ) relaxation constant for iterative process NT ) total number of pores REI ) ratio between the number of external pores to the number of internal pores Vres Hg ) dimensionless residual mercury volume V*Hg ) residual mercury volume ratio defined as Vres Hg |Cf/ | Vres C Hg i Greek Letters µ ) distribution mean value ξ ) dummy variable (diameter) σ ) distribution standard deviation
Figure 9. Plot of dimensionless residual mercury, Vres Hg , as a function of the assumed mean coordination number, Cf.
curves was obtained for all values of C, but C ) 4 was the value for which the simulated dimensionless residual mercury better approximated the experimental value. Despite the good agreement in the value of Vres Hg , the experimental and simulated extrusion curves were different; in particular, in the simulated curve, the mercury starts to be extruded at higher pressures than indicated by the experimental data. This result is probably due to a nonhomogeneous spatial distribution of the different sizes of pores, which is not taken into account in the random network model used. If the actual spatial distribution is known from independent experiments, such as serial sectioning, the network can be modified, and the iterative method can be applied to this nonrandom network. Another problem with these simulations results from the fact that the test case sample presents diameters that vary over about 4 orders of magnitude. When the volume-based probability distribution function, f ap V , is used to determine the number-based probability distribution function, f ap N , and then transformed into a limited number of pores by means of a random number generator, the probability of finding large pore diameters is very small, and a nonrepresentative network might be used. This problem can be solved by a suitable variable transformation on the diameter as will be reported in a future publication. List of Symbols C ) mean coordination number Ci ) mean coordination number for Initial-PSD Cf ) assumed mean coordination number d ) pore diameter dmax ) distribution maximum diameter dmin ) distribution minimum diameter eV ) convergence criterion f ) diameter probability density function fN ) number-based diameter probability density function fV ) volume-based diameter probability density function F ) diameter probability distribution function FN ) number-based diameter probability distribution function FV ) volume-based diameter probability distribution function i ) pore number I ) iteration number
Superscripts ap ) apparent exp ) experimental in ) initial re ) real sim ) simulated mod ) modified N ) final iteration
Acknowledgment Financial support for this work was in part provided by the research program PRAXIS XXI (PRAXIS/2/2.1/ QUI/07/94; PRAXIS/2/2.1/CEG/2564/95) for which the authors are thankful. V.G. Mata acknowledges her Ph.D. scholarship by JNICT/PRAXIS. Literature Cited (1) Washburn, E. W. Note on a method of determining the distribution of pore sizes in a porous material. Proc. Natl. Acad. Sci. 1921, 7, 115-116. (2) Ritter, H. L.; Drake, L. C. Pore-size distribution in porous materials. Pressure porosimeter and determination of complete macropore-size distribution. Ind. Eng. Chem. Anal. Ed. 1945, 17, 782-786. (3) Drake, L. C.; Ritter, H. L. Macropore-size distributions in some typical porous substances. Ind. Eng. Chem. 1945, 17, 787791. (4) Meyer, H. I. Pore distribution in porous media. J. Appl. Phys. 1953, 5, 510-513. (5) Androutsopoulos, G. P.; Mann, R. Evaluation of mercury porosimeter experiments using a network pore structure model. Chem. Eng. Sci. 1979, 34, 1203-1212. (6) Mann, R., Androutsopoulos, G. P.; Golshan, H. Application of a stochastic network pore model to oil-bearing rock with observations relevant to oil recovery. Chem. Eng. Sci. 1981, 36, 337-346. (7) Mann, R.; Golshan, H. Application of a stochastic network pore model to a catalyst pellet. Chem. Eng. Commun. 1981, 12, 377-391. (8) Loh K.-C.; Wang, D. I. C. Characterization of pore size distribution of packing materials used in perfusion chromatography using a network model. J. Chromatogr. A 1995, 718, 239255. (9) Portsmouth, R. L.; Gladden, L.F. Determination of pore connectivity by mercury porosimetry. Chem. Eng. Sci. 1991, 46, 3022-2036. (10) Portsmouth, R. L.; Gladden, L. F. Mercury porosimetry as a probe of pore connectivity. Trans. Inst. Chem. Eng. 1992, 70, 63-70. (11) Tsakiroglou, C. D.; Payatakes, A. C. A new simulator of mercury porosimetry for the characterization of porous materials. J. Colloid Interface Sci. 1990, 137, 315-339. (12) Tsakiroglou, C. D.; Payatakes, A. C. Effects of pore-size correlations on mercury porosimetry curves. J. Colloid Interface Sci. 1991, 146, 479-493.
Ind. Eng. Chem. Res., Vol. 40, No. 22, 2001 4843 (13) Tsakiroglou, C. D.; Payatakes, A. C. Mercury intrusion and retraction porosimetry in model porous materials. Adv. Colloid Interface Sci. 1991, 75, 215-253. (14) Tsakiroglou, C. D.; Payatakes, A. C. Analysis of the topological and geometrical characteristics of the pore space of permeable solids using serial tomography, mercury porosimetry, and theoretical simulation, Presented at COPS-IV, Bath, U.K., Sep 15-18, 1996.
(15) Mata, V. G. Characterisation of porous media: Porosimetry, 3D simulation and serial tomography. Application to catalytic supports. Ph.D. Thesis, University of Porto, Porto, Portugal, 1998.
Received for review February 1, 2001 Accepted August 1, 2001 IE0101137
