Porous Media Characterization Using Mercury Porosimetry Simulation

A porosimetry simulator associated with a tridimensional pore network model was developed. This simulator, presented in two joint papers, is used to i...
0 downloads 0 Views 161KB Size
Download PDF
Ind. Eng. Chem. Res. 2001, 40, 3511-3522

3511

Porous Media Characterization Using Mercury Porosimetry Simulation. 1. Description of the Simulator and Its Sensitivity to Model Parameters Vera G. Mata, Jose´ Carlos B. Lopes, and Madalena M. Dias* Laboratory of Separation and Reaction Engineering, Departamento de Engenharia Quı´mica, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

A porosimetry simulator associated with a tridimensional pore network model was developed. This simulator, presented in two joint papers, is used to infer the actual pore size distribution of a catalyst particle by means of an iterative method. In this first paper, the network model and porosimetry simulator are introduced and a sensitivity study of the simulator to the model parameters is performed. These parameters were divided in two main groups. The first relates to the spatial arrangement of the three-dimensional pore network and includes the number of pores, the mean coordination number, the accessibility, and the network geometry. The second group is related to the characteristics of the pore size distribution, namely, its spreading, symmetry, and number of modes. Introduction Diffusion, convection, and dispersion phenomena in porous materials are strongly dependent on the pore structure of the medium. Immiscible displacement of nonwetting fluids by a wetting one,1 chromatographic separations,2 and heterogeneous reactors with catalytic or noncatalytic reactions3 are some examples of processes that depend on the transport of fluid mixtures through porous materials. For heterogeneous porous catalysts, the connectedness and geometry of the pores play a role in transport and reaction that is often more important than other relevant physical data.4 In heterogeneous catalysts, the activity, selectivity, and related deactivation-regeneration phenomena can be drastically influenced by the pore structure properties of the catalyst pellet. The distribution of the active material throughout the pellet and the accessibility of the intraparticle surface area are, among other factors, strongly dependent upon the pore size distribution (PSD) and pore connectivity.5,6 Descriptions of porous media using pore network models that incorporate geometrical and topological properties of the pore space have been used since the pioneering work of Fatt.7-9 Two-dimensional (2D) and three-dimensional (3D) network models have been used to simulate and predict transport and flow processes.10-17 These models are simplified mathematical representations of real porous media, and their objective is to provide a reasonable idealization of the complex structure of the porous medium in such a way that transport processes can be treated mathematically. To this end, these models must incorporate the most relevant characteristics of the solid material, while their complexity should be kept at a manageable level.12 However, the full potential of network models for simulating capillary and flow phenomena in porous media has not been fully exploited. The customary way of using network models * To whom correspondence should be addressed. E-mail: [email protected]. Phone: 351-22 508 1661. Fax: 351-22 508 1674.

consists of selecting a 2D or a 3D network, assuming most of its geometrical and topological properties, and proceeding with the simulation of the process under study, disregarding the actual properties of the real porous medium and their relationships to the network model parameters. Most of the works mentioned above have shown that the PSD and the pore connectivity are two of the most important properties for the characterization of any porous material and for the simulation of transport and flow phenomena. For decades, mercury porosimetry has been one of the techniques used for the determination of the PSD. This technique includes an experimental partsmercury intrusion and extrusionsfollowed by its interpretation by means of a model.18 The conventional parallel pore model generally used19,20 leads to results that often are not very meaningful, as will be shown, because this model ignores the existence of a network of interconnected void spaces within the porous solid. Porosimetry simulators associated with 2D and 3D network pore models that account for the connectivity of the porous space have been previously developed by several authors. Androutsopoulos and Mann21 developed a 2D square network of cylindrical pores. Each pore was considered to be of fixed length and of variable random diameter. The simulations were carried out by a combination of trial-and-error and Monte Carlo techniques, and they were able to obtain theoretical and experimental intrusion and extrusion curves within reasonable qualitative agreement.10 Cox22 studied the effect of various types of crystalline lattices (face centered cubic, body centered cubic, and simple cubic) and by assuming a normal distribution for the pore diameters obtained a rapid estimation of the appropriate parameters for the PSD. Portsmouth and Gladden23,24 developed a spherical, 3D, random network of cylindrical pores in order to investigate the effects of PSD and pore connectivity on mercury porosimetry experiments. For a Gaussian PSD, a characteristic entrapment value of mercury after extrusion was obtained and found to be

10.1021/ie010112e CCC: $20.00 © 2001 American Chemical Society Published on Web 07/10/2001

3512

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001

dependent only on the distribution parameters (mean and standard deviation). A sequence of mini-hysteresis loops provided information on the medium connectivity. Tsakiroglou and Payatakes25,26 developed a simulator based on their earlier throat-and-chamber model.27 This simulator included the mechanisms by which mercury menisci move in pores and stop at entrances to throats or chambers, and it took into account the mechanism of snap-off during extrusion. They studied the effects of throat-and-chamber size distributions, the coordination number, and pore size correlations on the pressure curves. Their main conclusions were that the intrusion curve depends mainly on the throat size distribution and on the mean coordination number and that the extrusion curve depends strongly on the ratio between the average chamber diameter and the average throat diameter. In previous works, mercury porosimetry network simulators were mostly used to study the effect of the various model parameters on the intrusion and extrusion curves. However, the reverse process of deconvoluting the information obtained from simulations in order to obtain the actual porous medium properties was not done in a systematic way. In addition, the extrusion curves and the information contained in them, strongly related to the pore connectivity, are ignored. The aim of this work is the development of a methodology for obtaining the real PSD and pore connectivity of a porous sample, based on experimental mercury porosimetry data and on a porosimetry simulator. In this first paper (part 1), a mercury porosimetry simulator based on a 3D pore network model developed by Mata28 is described. The influence of the network parameters, such as the network size, the mean coordination number, the accessibility, and the PSD, on the simulated porosimetry curves is studied. In the second joint paper (part 2), an iterative methodology will be described to obtain the PSD and the mean coordination number of a porous sample using experimental mercury porosimetry data and the porosimetry simulator. For the sake of simplicity in implementation and calculation, a lattice of cylindrical pores together with a simple set of rules for the intrusion and extrusion phenomena was used. However, more complex network geometries of the throat-and-chamber type and more elaborate rules could be implemented and used with this iterative methodology. Pore Network Model The pore network model used in this work is a cubic lattice of interconnected pores. The lattice total number of nodes is M ) mxmymz where mx, my, and mz are the number of nodes in the x, y, and z spatial directions. The nodes are connected by cylindrical pores of constant length, l, and with diameters di, i ) 1, ..., NT, following a given PSD, fN. There are two types of pores: external pores that have direct access to the network exterior and internal pores within the network. The total number of pores, NT, is given by the sum of the number of internal pores, NInt, and the number of external pores, NExt. The model connectivity is set by the mean coordination number, C, defined as the average number of pores that are connected to each node. The cubic lattice has a maximum value of Cmax ) 6, in which case each internal pore is connected to 10 neighbors and each external pore to 5 neighbors. Networks with C < 6 are obtained by

randomly removing a number of internal pores, ΘInt, and a number of external pores, ΘExt, always preserving at least two pores connected to each node. The actual value of C is related to Cmax by C ) (NInt + NExt)Cmax/ Nmax, where Nmax ) 3mxmymz + mxmy + mymz + mxmz. The accessibility to the internal pores is determined by the parameter REI ) NExt/NInt. In the study of different network sizes and geometries, a constant value of REI ensures that the access of an external fluid to an internal pore has the same probability to occur. For a given set of parameters, C and REI, a sufficiently high value of the total number of pores, NT, needs to be selected so that the results are not influenced by the network size, as will be discussed later in this paper. PSD: Some Definitions The PSD defined between a minimum and a maximum diameter, dmin and dmax, is described in terms of the number diameter probability density function

fN(ξ) )

( )

d N(ξ) dmax > ξ > dmin dξ NT

(1)

where N(ξ)/NT represents the number fraction of pores with diameters between ξ and ξ + dξ, and subjected to max fN(ξ) dξ ) 1. To describe the intruthe restriction ∫ddmin sion curves, a number diameter probability distribution function is introduced

FN(di) )

∫dd

max

fN(ξ) dξ

i

(2)

where FN(di) represents the fraction of pores with diameter greater or equal to di and is subject to FN(dmax ) ) 0 and FN(dmin ) ) 1. [The probability distribution function is usually defined as FN(di) ) i fN(ξ) dξ, but eq 2 is a more convenient definition ∫ddmin commonly used in porosimetry.] Because experimental intrusion curves are obtained in terms of volume, the PSD can also be described by the volume diameter probability density function

fV(ξ) )

( )

d V(ξ) d F (ξ) ) dξ VT dξ V

(3)

where FV(di) is the volume diameter probability distribution function. With the simplifying assumption of constant pore length in the network model, it can be easily shown that the number- and volume-based probability density functions are related by

fV(di) )

di2fN(di)

∫dd

ξ2fN(ξ) dξ

max

min

(4)

The number-based function required by the simulator is determined from the volume-based function obtained directly from porosimetry experiments, solving eq 4 in order to obtain fN(di).

fN(di) )

fV(di)

∫d

di2

dmax min

[fV(ξ)/ξ2] dξ

(5)

For practical purposes and ease of graphical representation and for the calculations in part 2 of this work,

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 3513

discretized forms of the diameter probability density functions, fN(di) and fV(di), and of the diameter probability distribution functions, FN(di) and FV(di), must be used. In the following simulations and its results, these functions have been discretized in j ) 1, 2, ..., J ) 40 equally spaced classes. In the fN(di) and fV(di) plots, each class is represented as a single point connected by a straight line to the neighboring classes. The use of column histograms was avoided so that more than one curve may be represented and compared in a single plot. Mercury Porosimetry Simulation Rules The mercury porosimetry simulator described below is based on the above pore network model and on a set of heuristic rules representative of the real porosimetry process, which has two distinct parts: intrusion of mercury into the pore network and extrusion of mercury from the pore network. Mercury Intrusion. During the intrusion process, the mercury pressure, in contact with the surface of the porous sample, is increased in a stepwise way, leading to the filling of the porous space with mercury. For an empty pore, the conditions necessary for mercury penetration are as follows: (i) Contact between the pore and the bulk mercury has been established; i.e., a continuum mercury path between the external mercury and the pore exists. (ii) The pore size is such that the Washburn19 criterion for penetration is met

di g -

4σmv cos θ P

(6)

where di is the pore diameter and σmv, θ, and P are the mercury surface tension, contact angle, and pressure, respectively. The intrusion simulation progresses with a stepwise pressure increase. Initially, only the external pores are in contact with the bulk mercury. The pressure is raised, and each pore that is exposed to mercury is tested to verify if it fulfills the Washburn criterion. Pores that meet this requirement are filled with mercury, and the neighboring pores are stored in a list of potential mercury intrusion pores. This list is scanned and updated until the mercury front is incapable of advancing any further at the given pressure. The volume of mercury that has intruded into the porous space at this pressure is then recorded. This procedure is successively repeated for higher pressures, until the entire network has been filled with mercury. Mercury Extrusion. The extrusion process reverses the intrusion process; i.e., the pressure is lowered, and the mercury is allowed to escape from the porous sample. Three simultaneous criteria must be met before mercury can leave a pore: (i) The mercury within a pore has a free surface, i.e., a meniscus. (ii) A percolation (continuous) path from the pore to the bulk mercury exists. (iii) The Washburn criterion for withdrawal is verified.

di < -

4σmv cos θ P

(7)

At the end of the invasion process, a pair of interfaces within the smallest pore of the network is supposed to

remain uncoalesced, and it is from these two interfaces that mercury extrusion begins upon sufficient decrease of pressure. The mercury leaves this pore at the pressure dictated by the Washburn criterion. Pores connected to this one belong to a list of potential pores for mercury drainage. For each decreasing pressure value, this list is scanned and the pores that fulfill all of the three criteria become empty. The volume of mercury that has left the porous space at this pressure is then recorded. At each pressure step, the whole network is checked for the existence of possible clusters of mercury that become isolated from the bulk liquid. These clusters are trapped and remain in the network up to the end of the porosimetry process. The pressure is further decreased, and this procedure is repeated until it is impossible to drain any more mercury from the network. The volume of mercury retained in the network at this point is referred to as the residual mercury volume and is expressed by the dimensionless variable Vres Hg representing the residual mercury volume fraction of the total intruded volume, VT. Mercury Porosimetry Simulator The mercury porosimetry simulator described below is based on the above pore network model and simulation rules. The global setup of a porosimetry simulation is done in two steps: (i) generation of the network and (ii) intrusion and extrusion simulations. For the network generation, a priori knowledge of the following parameters is necessary: the network size, mx × my × mz; the mean coordination number, C; and the accessibility parameter, REI. The mean coordination number is modified by starting with a network with C ) 6 and randomly removing a number of internal pores, ΘInt, and a number of external pores, ΘExt. A given Initial-PSD, in terms of the number diameter probability density function, f in N , is chosen, and NT pore diameter values are generated and randomly assigned to each network pore. Using eq 4, the corresponding volume diameter probability density function, f in V , is obtained. The porosimetry simulation progresses as described, and the simulated intrusion and extrusion curves are obtained and computed in terms of the volume diameter probability distribution function, FV(di). Assuming the parallel pore model,19,20 an Apparent-PSD can then be directly derived from the simulated intrusion curve. The derivative of FV(di) (see eq 3) leads to the apparent volume diameter probability density function, f ap V , and then using eq 4 and the proposed iterative scheme, the apparent number diameter probability density function, f ap N , is obtained. The word apparent used here means that the obtained PSD from differentiation of the simulated intrusion curves will be different from the one ideally obtained if the network consisted solely of a bundle of parallel pores. In reality, the porosimetry curves and the corresponding derivative curves are affected by the pore network interconnectivity and by the inaccessibility of the internal pores to the external mercury. This simulator was implemented on Fortran 77 code, and the simulations were run on an IBM Power RISC/ 3CT computer with a performance of a 120 MFlop Linpack benchmark. A typical simulation on a network with C ) 6 and 20 × 20 × 20 nodes takes approximately 15 min of computing time, 5 min for intrusion, and 10 min for extrusion.

3514

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 Table 1. Pore Network Parameters for the Various Simulations figure C mx × my × mz

Sensitivity Analysis The final objective of this work is the determination of the actual PSD and the mean coordination number of a given porous sample by means of the iterative method described in part 2, using experimental porosimetry data and the porosimetry simulator described above. Before advancement in this direction, a preliminary study, of the simulator and its sensitivity to the various network parameters, becomes necessary, so that the most important parameters are clearly singled out. First, it is necessary to address the question of the network size, expressed by the total number of pores, NT. On the one hand, the pore network model is assumed to be a representative part of a given porous sample, so it should be as large as possible. On the other hand, computational limitations require that the network be finite and sufficiently small not to exceed reasonable computing time and memory requirements. The question that must be addressed is then, what is the minimum network size to be used to obtain simulation results that are statistically representative, or in other words, the results are independent of the network size. Second, a sensitivity analysis of the simulated porosimetry curves to the network connectivity, expressed in terms of the mean coordination number, C, and the external mercury accessibility, expressed in terms of the accessibility parameter, REI, is presented. Finally, to study the effect of different PSDs on the simulated intrusion and extrusion curves, results with different Initial-PSD curves are analyzed. To this end, a characteristic reference Gaussian distribution (see the Appendix) with mean pore size, µ ) 750 µm, and standard deviation, σ ) 250 µm, was used for the generation of the pore diameter distribution. A discrete Initial-PSD, f in N (di), was determined from randomly generated NT pore diameters, properly truncated from a minimum value of dmin ) 0 and a maximum value of dmax ) 1500 µm. The random pore diameters were generated by the mapping of a sequence of uniformly distributed numbers in the range of 0-1, on the corresponding number diameter probability distribution function Fin N (di). Figure 1 shows the resulting derived number-based discrete function with µ ) 750 µm and σ ) 247 µm and the corresponding volume-based discrete function with µ ) 896 µm and σ ) 224 µm. Also, the network parameters used in all of the subsequent sensitivity studies are summarized in Table

NInt

NExt

ΘInt

ΘExt

REI

0 0 0 0 0 0 0 0 0 0 0 0

0.50 0.22 0.14 0.11 0.083 0.069 2.00 0.57 0.33 0.24 0.18 0.15

6 6 6 6 6 6 3 3 3 3 3 3

5×5×5 10 × 10 × 10 15 × 15 × 15 20 × 20 × 20 25 × 25 × 25 30 × 30 × 30 5×5×5 10 × 10 × 10 15× 15 × 15 20 × 20 × 20 25 × 25 × 25 30 × 30 × 30

450 3 300 10 800 25 200 48 750 83 700 225 1 650 5 400 12 600 24 375 41 850

300 2 700 9 450 22 800 45 000 78 300 75 1 050 4 050 10 200 20 625 36 450

150 0 600 0 1350 0 2400 0 3750 0 5400 0 150 225 600 1 650 1350 5 400 2400 12 600 3750 24 375 5400 41 850

5, 9, 10 6 5 4 3

20 × 20 × 20 20 × 20 × 20 25 × 25 × 25 30 × 30 × 30

25 200 21 000 32 500 41 850

22 800 19 000 29 405 37 864

2400 0 0 2000 3 800 400 3095 15 595 655 3986 40 436 1414

6

5 5 5 5

20 × 20 × 20 20 × 20 × 20 20 × 20 × 20 20 × 20 × 20

21 000 21 000 21 000 21 000

18 600 2400 20 792 208 20 979 21 20 998 2

7a

6 20 × 20 × 20 25 200 22 800 2400 6 10 × 10 × 80 25 700 22 300 3400 6 5 × 10 × 160 26 450 21 550 4900

0 0 0

7b

6 20 × 20 × 20 25 200 22 800 2400 6 10 × 10 × 80 24 647 22 300 2347 6 5 × 10 × 160 23 818 21 550 2268

0 0 0.11 0 1053 0.11 0 2632 0.11

8

6 20 × 20 × 20 25 200 22 800 2400

0

2-4

Figure 1. Gaussian-type initial PSD showing the number diameter probability density function, f in N , and the corresponding volume diameter probability density function, f in V.

NT

0.11 0.11 0.11 0.11

4 200 0 0.13 2 008 2192 0.01 1 821 2379 0.001 1 802 2398 0.0001 0 0.11 0 0.15 0 0.23

0 0.11

1, where for simplicity the data are grouped by figure number. For each case, the total number of pores, NT, the number of internal and external pores, NInt and NExt, the number of internal and external pores removed, ΘInt and ΘExt, and the resulting value of REI are listed. Effect of the Network Size. For this study, six different network sizes (5 × 5 × 5, 10 × 10 × 10, 15 × 15 × 15, 20 × 20 × 20, 25 × 25 × 25, and 30 × 30 × 30) were studied for the maximum and minimum values of the mean coordination number (C ) 3 and 6). The number of nodes ranged from 75 to 27 000 and the total number of pores from 225 to 41 850 for C ) 3 and from 450 to 83 700 for C ) 6. In each case, NT pore diameters were generated following the chosen characteristic theoretical probability density function using a random number generator function. The function arguments seedsused on the random number generator, often designated as the random number initiator, can have an arbitrary value, and for each seed value, it determines a particular random number sequence. It was then necessary to find the minimum network size for which the pore diameter function and the simulated curves are not dependent on the seed value. The effect of different seed values in the generation of f in N was analyzed for the different network sizes. Figure 2 shows, for the six network sizes, the comparison between the values of the mean and the standard deviation for the reference distribution function and the corresponding computed values obtained with four discrete functions generated using different seeds. It can be seen that, as the total number of pores increases, the computed mean approaches the reference value, µ ) 750 µm, and the standard deviation approaches the reference truncated function value, σ ) 247 µm. It was observed that for larger network sizes, with NT g12 600, the difference between the discrete and the reference values is within 0.1%. Therefore, it was accepted that the discrete Initial-PSD should always be based on a

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 3515

Figure 2. Plots of (a) mean, µ, and (b) standard deviation, σ, as functions of the total number of pores, NT, for different size networks, mean coordination number (C ) 3, solid symbols; C ) 6, open symbols), and four different initial seeds.

network size with NT g12 600 or roughly a minimum 3D network size of at least 20 × 20 × 20. Once the network size and the Initial-PSD were chosen, the indexation process of assigning a diameter value to each pore in the network was implemented, by using a uniformly random number generation for the pore-to-diameter indexation. Then, intrusion and extrusion simulations were carried out for various network sizes, using the same initial PSD and four different poreto-diameter seed values for the indexation. Figure 3 shows the corresponding simulated porosimetry curves for C ) 3 and 6, where each curve represents the average curve over four different seed simulations. For comparison, 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 (dashed line). Analysis of Figure 3 shows that all intrusion curves lie below the p-p-m curve, whereas the extrusion curves lie above it; as the network size increases, the simulated intrusion curves get further and further apart from the p-p-m curve, approaching a limiting curve. This indicates that there is a network size above which the simulations can be representative of experiments in a large sample. To characterize the effect of seeding in the intrusion curves, Figure 4 shows the square root of the sum of the square of the residuals of the intrusion curve for each seed value, FV(di) relative to the mean value of all seed values, F h V(di), that is, S ) x∑i[FV(di)-F h V(di)]2, for different network sizes. From Figure 4 it is observed that S decreases logarithmically with NT, and for small networks, of say NT < 10 000, the values of S show a considerable spread of the order of 0.1 for the different

Figure 3. Intrusion and extrusion curves averaged over four different pore-to-diameter indexation seeds, for various network sizes and (a) C ) 6 and (b) C ) 3. The dashed line shows the corresponding p-p-m curve.

Figure 4. Plot of the sum of the squares of residuals, S, as a function of the total number of pores, NT, for different size networks, mean coordination number (C ) 3, solid symbols; C ) 6, open symbols), and four different pore-to diameter indexation seeds.

seed simulations. Analysis of Figure 3a shows that for C ) 6 (open symbols) a network of 20 × 20 × 20 (NT ) 25200) gives in essence the same result as that of a network of 30 × 30 × 30 (NT ) 83 700). Figure 4 shows that the overall deviation for different indexation seeds in the case of the 20 × 20 × 20 network is already around 0.1. In the case of C ) 3 (solid symbols), only for larger network sizes (30 × 30 × 30 and NT ) 41 850), the simulated intrusion curves seem to approach the limiting curve (Figure 3b) and S takes values on the order of magnitude of 0.1. The same type of study was conducted for other intermediate values of the mean

3516

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001

coordination number. In the remainder of this work, networks with a size of 20 × 20 × 20 (8000 nodes, NT ) 25 200) are used for C ) 6 and 5, while for C ) 4 and 3, larger networks with sizes of 25 × 25 × 25 (15 625 nodes, NT ) 32500) and 30 × 30 × 30 (27 000 nodes, NT ) 41 850), respectively, are necessary. Tsakiroglou and Payatakes,25 who also used a cubic network, did similar studies of the minimum required network size and concluded that a 20 × 20 × 20 network was adequate for all of their simulations. In addition, Portsmouth and Gladden24 used 3107 nodes on a spherical network model, independently of the value of the coordination number, giving a number of nodes between 18 642 for C ) 6 and 9321 for C ) 3, but no sensitivity study was reported. The network size also affected the extrusion process, and this effect may be characterized and quantified by the calculation of the residual mercury volume, Vres Hg . Figure 3a shows that, as the network size increases, Vres Hg increases, because the probability of finding a percolation path from an interior pore to the bulk mercury decreases. This decrease in Vres Hg is also related to the decrease in the accessibility of the internal pores to the external mercury, as will be seen later in this work. This observation is in agreement with the results of Tsakiroglou and Payatakes25 and Portsmouth and Gladden.24 Effect of the Mean Coordination Number. The effect of the pore interconnectivity on the simulated intrusion and extrusion curves was studied by means of the mean coordination number, taking the reference Initial-PSD and obtaining the simulated intrusion and extrusion curves for various values of the mean coordination number, namely, C ) 6, 5, 4, and 3. The mean coordination number is modified by starting with a network of C ) 6 and randomly removing a number of internal, ΘInt, and external, ΘExt, pores, as reported in Table 1, taking care to keep constant the ratio between the internal and external pores at REI ) 0.11. The influence of C on the simulated intrusion and extrusion curves is shown in Figure 5a, where one observes that, as the mean coordination number decreases, both the intrusion and the extrusion curves get further apart from the p-p-m curve. This is in agreement with the previous results reported by Tsakiroglou and Payatakes25 and may be explained by the fact that, as the coordination number decreases, the number of percolation paths decreases. Hence, a given volume of mercury is reached at higher pressures (smaller diameters) during intrusion and at lower pressures (larger diameters) during extrusion. Parts b and c of Figure 5 show the corresponding volume diameter probability density function, f ap V , and the number diameter probability density function, f ap N, for the different values of C (solid lines), compared with the initial-PSD (dashed line). It is observed that the Apparent-PSD does not describe correctly the actual network PSD. While for small values of di the two curves are relatively similar, for large pore diameters, the Apparent-PSD underestimates the Initial-PSD, and this is done by overestimation in the region of intermediate values of di. This observation is valid for all values of the mean coordination number, but as C decreases, these deviations become increasingly important. Conner and Lane5 and Tsakiroglou and Payatakes25 also observed a similar behavior, expressing it in terms of “the PSD moves to smaller sizes as the mean coordination

Figure 5. Simulations with a Gaussian Initial-PSD. Effect of the mean coordination number, C, on (a) the intrusion and extrusion curves, (b) the volume diameter probability density function, f ap V , and (c) the number diameter probability density function, f ap N. Dotted lines show the corresponding p-p-m or Initial-PSD curves. Table 2. Apparent PSD Mean and Standard Deviation and Residual Mercury for Various Values of C f ap V REI ) 0.11, Figure 5

f ap N

C

µ (µm)

σ (µm)

µ (µm)

σ (µm)

Vres Hg

6 5 4 3

840 823 809 775

164 154 153 165

739 733 721 679

208 197 192 195

0.39 0.46 0.51 0.55

number decreases”. This is clearly due to the pore interconnectivity and the “hiding” of large pores by the smaller pores. This effect is also stressed by the calculation of the mean, µ, and the standard deviation, σ, for ap f ap V and f N , given in Table 2 and showing that both values are always underestimated. Figure 5a shows that the extrusion curves are also

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 3517

Figure 6. Simulations with a Gaussian Initial-PSD. Effect of the ratio of external to internal pores, REI, on (a) the intrusion and extrusion curves and (b) the number diameter probability density function, f ap N . Dotted lines show the corresponding p-p-m or Initial-PSD curves.

affected by C, and this effect may be quantified by the calculation of residual mercury, Vres Hg , that is the mercury retained in the network after extrusion is completed. Data in Table 2 show that Vres Hg increases as the mean coordination number decreases. This may be justified by the observation that, as C increases, the probability of finding adjacent pores filled with mercury from which the mercury may retract increases. Similar qualitative trends were obtained by Conner and Lane5 and by Tsakiroglou and Payatakes25 where in both cases cubic networks were used. Quantitative results are harder to compare because the networks used in these two works were of the throat-and-pore type, potentially giving rise to high mercury retentions (between 75 and 85% for C between 3 and 6) as reported by Tsakiroglou and Payatakes.25 Conner and Lane5 report much lower values of Vres Hg (between 13 and 23% for C between 4 and 6), but this result may be due to the relatively small networks (5 × 5 × 5 and 9 × 9 × 9) used. Portsmouth and Gladden24 using a spherical network reported values of Vres Hg in the range of 25-40% for C between 3 and 6. They also noticed the “striking feature” of the existence of a maximum in the value Vres Hg for values of C lower than 4, which was not observed in the present work. Effect of the Accessibility. The accessibility of the internal pores to the external mercury, that is, to the external pores, is measured by REI, the ratio between the number of external pores, NExt, and the number of internal pores, NInt (see Table 1). Figure 6 shows the results of simulations performed in networks with C )

Figure 7. Simulations with a Gaussian initial PSD. Effect of network geometry on the intrusion and extrusion curves with (a) varying REI and (b) constant REI. Dotted lines show the corresponding p-p-m or Initial-PSD curves. Table 3. Apparent PSD Mean and Standard Deviation and Residual Mercury for Various Values of REI f ap V REI C ) 5, Figure 6 0.13 0.01 0.001 0.0001

f ap N

µ (µm) σ (µm) µ (µm) σ (µm) 840 811 804 780

164 141 135 125

733 728 725 716

208 189 185 172

Vres Hg 0.45 0.47 0.48 0.49

5 and at various REI values. As REI decreases, the intrusion curves get further apart from the p-p-m curve because of interconnectivity effects; that is, it is more difficult for the mercury to reach the internal pores. This is clearly reflected in the shape of the apparent PSD, where the fraction of intermediate size diameters becomes larger and the mean, µ, decreases as REI decreases. The value of the residual mercury (see Table 3) increases as REI decreases, again because of the interconnectivity effects. The problem of accessibility was first addressed by Portsmouth and Gladden,24 who ensured that the access to internal pores was held constant as the size or the connectivity of the spherical network was altered, but they did not study the effect of changing the accessibility on the simulations results. Nevertheless, the effect of accessibility seems to be of less importance than some of the other parameters; for example, while varying the mean coordination number from 6 to 3 results in an increase of Vres Hg of ca. 40%, a decrease of 4 orders of magnitude in REI resulted in only a 10% increase in Vres Hg (see Tables 2 and 3). Effect of the Network Geometry. In the previous simulations, an equal number of nodes in each direction

3518

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001

was used, that is, mx ) my ) mz. In this section, the geometry of the network is varied, keeping constant the total number of nodes, M ) mxmymz. In the first set of simulations, no internal or external pores were removed and the value of REI changed, according to the network geometry (Figure 7a). The intrusion curves show that the differences between them are related to the change in REI, as observed in the previous section. The same observation is true for the extrusion curves and the residual mercury. In the second set of simulations, some external pores were removed randomly for each geometry in order to keep REI ) 0.11 (Figure 7b). No significant differences are observed now, leading to the conclusion that the porosimetry simulations are not affected by the network geometry, as long as accessibility is maintained constant. From the previous studies, it is clear that, if a sufficiently large network is used, the porosimetry simulations are mainly dependent only on two topology parameters, the mean coordination number, C, and the network accessibility parameter, REI. Although the intrusion and extrusion processes are impaired when REI decreases, this effect is less pronounced than the effect of C. For this reason, in the remainder of this work, the value of REI is kept constant and only C will be varied. Sensitivity to the Initial PSD In the previous examples, the same reference InitialPSD was used. In this section, the influence of the Initial- PSD is addressed by studying the effect of the spread (standard deviation), the number of modes (bimodal distribution), and the symmetry (non-Gaussian). Spread. Figure 8 shows the results of simulations done with Initial-PSDs obtained from Gaussian distribution functions with a mean value, µ ) 750 µm, and standard deviations σ ) 250, 150, and 50 µm, truncated from a minimum value of dmin ) 0 and a maximum value of dmax ) 1125 µm. Figure 8a shows that narrower Initial-PSDs, that is, lower values of σ, result in narrower intrusion curves, narrower extrusion curves, and lower values of the residual mercury. This result may be explained by the fact that, as the distributions become narrower, the range of pore diameters decreases, and interconnectivity effects are less important. Again, these results are in good qualitative agreement with those of Tsakiroglou and Payatakes25 and Portsmouth and Gladden25 and also with those of Tsetsekou et al.,6 who used a probabilistic model applied to a onedimensional model of nonintersecting corrugated pores. Figure 8b shows that the apparent PSD, f ap N , becomes closer to the Initial-PSD, f in N , as σ decreases, a result that was also reported by Conner and Lane.5 Number of Modes. In real porous media, it is common that the PSD presents more than one maximum value.6 To illustrate this case, the initial PSD was based on a bimodal distribution composed of two Gaussian distributions with mean values of µ1 ) 750 µm and µ2 ) 2250 µm and the same standard deviation value, σ ) 250 µm, and truncated between dmin ) 0 and dmax ) 3000 µm. Figure 9 shows the results of the porosimetry simulations for two values of the mean coordination number, C ) 6 and 3. The intrusion curves show the characteristic shapes for a bimodal intrusion curve with

Figure 8. Simulations with a Gaussian Initial-PSD. Effect of the standard deviation, σ, on (a) the intrusion and extrusion curves and (b) the number diameter probability density function, f ap N. Dotted lines show the corresponding p-p-m or Initial-PSD curves.

two plateaus in the volume diameter probability distribution function curve, FV. The extrusion curves show that the mercury starts to flow out of the network, only after the capillary pressure corresponding to di ) µ2 is reached, and as before, the residual mercury increases as C decreases. When Figures 5a and 9a are compared, bimodal distributions seem to give rise to higher values of Vres Hg than single-mode distributions and this result is in agreement with Tsetsekou et al.6 Analysis of the Apparent-PSD curves shows totally different behavior for each value of the mean coordination number. For C ) 6 the Apparent-PSD closely follows the initial PSD for small diameters (0 < di < 2000), but the fraction of larger diameters (2400 < di < 3000) is underestimated, resulting in an overestimation of diameters in the range 2000 < di < 2400. This means that the pores with small diameters are directly connected to the external mercury but that some large pores are only accessible to the external mercury when smaller pores have been invaded. For C ) 3 the differences between the Apparent and Initial-PSDs are even more accentuated. The fraction of diameters in the region of the distribution’s second mode (1500 < di < 3000) is underestimated, the fraction of diameters in the range 0 < di < 600 is relatively well represented, but the fraction of diameters in the range 600 < di < 1500 is overestimated. These results indicate that the pores with larger diameters are accessible to the external mercury only through the smaller pores, showing clearly a large interconnectivity effect. A possible explanation is that for C ) 6 there are more percolation pathways to reach any given branch and so the probability of a branch with a large pore being hidden by

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 3519

Figure 9. Simulations with a bimodal Initial-PSD. Effect of the mean coordination number, C, on (a) the intrusion and extrusion curves and (b) the number diameter probability density function, f ap N . Dotted lines show the corresponding p-p-m or Initial-PSD curves.

small pores is lower than that in the case of C ) 3. The result is that in this last case a relatively high number of the larger pores is only connected to the bulk mercury through the smaller pores and is only filled for large pressures. Using the p-p-m implies that the intrusion volume of these large pores these pores is associated to small pores and will be counted as a high number of small pores rather than a few large pores. Symmetry. In many porous media, the PSD covers a wide range of diameters, sometimes over several orders of magnitude, and in this case the distribution is no longer symmetric. The software used in most commercial porosimetry equipment often uses a logarithmic scale to represent the diameter, suggesting the use of a log-normal distribution function. This function has the advantage of being defined only for positive pore diameter values, in contrast with the Gaussian distribution that ranges from -∞ to +∞. This function has been used to analyze membrane PSDs.29 Here an upperlimit log-normal (ULLN) distribution function is used (see the Appendix for details) for the Initial-PSD. This function has the distinct advantage of having an upper value, and it has been successfully used to represent data in other fields of study, namely, drop size distributions in annular gas-liquid flow.30,31 Figure 10 shows the number and volume probability density functions used as the Initial-PSD for the porosimetry simulations shown in Figure 11. The behavior of these curves is similar to the simulations with a Gaussian initial distribution, when the porosimetry curves are plotted on a logarithmic scale. As the coordination number decreases, the intrusion curves get further apart from

Figure 10. ULLN-type Initial-PSD showing the number diameter probability density function, f in N , and the corresponding volume diameter probability density function, f in V : (a) linear diameter scale; (b) logarithmic diameter scale.

the p-p-m curve and the amount of residual mercury increases. Again, it is clear that the Apparent-PSD curves do not match the Initial-PSD, and in this case the differences are more accentuated in the range of smaller diameters. These results are in agreement with Conner and Lane,5 who analyzed the effect of distributions skewed toward both the small and large radii. Conclusions A porosimetry simulator associated with a 3D pore network model has been developed. Although the pore network model is versatile and general, there are a couple of restrictions that must be pointed out. First, the pores are assumed to be of equal length, resulting in the fact that the actual pore volume and the porosity of the porous material cannot be matched. Another problem related to this issue is that possible correlations between the pore length and pore diameter cannot be used. This is likely to be important in the case of porous media with pore structures ranging over several orders of magnitude. Another restriction imposed on the network formation is that at least two pores are connected to each node, which implies that there are no dead-end pores in the network. This is a useful assumption if one is interested in analyzing only the percolating part of the pore space. Nevertheless, in porosimetry simulations, the presence of dead-end pores will not significantly affect the intrusion curves, because the mercury can reach all pores. On the other hand, extrusion simulations may be affected by the presence of deadend pores, because new rules will have to be assumed

3520

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001

was found to have a major effect on the shape of both the intrusion and extrusion curves. The accessibility of the internal pores to the external mercury was also analyzed, and again both the intrusion and extrusion curves are affected but to a lesser degree than the mean coordination. The simulated curves were found to be independent of the network geometry as long as the other parameters, such as the mean coordination number and the accessibility, were kept invariant. The sensitivity analysis to the Initial-PSD was done by changing the spread of a Gaussian distribution, by using a bimodal distribution, and by using the ULLN nonsymmetric distribution. It was found that, as the distributions become narrower, the interconnectivity effects become less important. On the other hand, these interconnectivity effects become very important for bimodal distributions and nonsymmetric curves. This simulator is used in a joint paper to infer the actual PSD of a catalyst particle by means of a new iterative method. From the sensitivity analysis, it becomes clear that the simulations in the iterative method have to be implement on networks with a minimum size of 20 × 20 × 20. The mean coordination number, C, is the parameter that most affects the simulations, but although the accessibility parameter, REI, affects the simulated curves on a much lesser degree, this parameter should be kept constant in all simulations, even when C varies. Acknowledgment Financial support for this work was in part provided by the research program PRAXIS XXI (PRAXIS/2/2.1/ QUI/07/94 and PRAXIS/2/2.1/CEG/2564/95) for which the authors are thankful. V.G.M. acknowledges her Ph.D. scholarship by JNICT/PRAXIS. Appendix The Gaussian or normal distribution function is defined as

f(ξ) )

Figure 11. Simulations with a ULLN Initial-PSD. Effect of the mean coordination number, C, on (a) the intrusion and extrusion curves, (b) the volume diameter probability density function, f ap V , and (c) the number diameter probability density function,f ap N. Dotted lines show the corresponding p-p-m or Initial-PSD curves.

as to what happens to the mercury in these dead-end pores; the mercury may be allowed to flow out of these pores, or the mercury will be trapped if the dead pores are connected to the bulk mercury by a relatively large pore. A study of the sensitivity of the simulator to the model parameters has been performed. These parameters were divided into two main groups: the network parameters (including the number of pores, the mean coordination number, the accessibility, and the network geometry) and the Initial-PSD parameters (including the PSD spread, symmetry, and number of modes). The network size was chosen so that the invasion curves were representative of a large network and the required number of pores was found to be dependent on the mean coordination number. The mean coordination number

1 1ξ-µ2 exp 2 σ σx2π

[ (

)]

-∞ < ξ < +∞

(8)

where µ is the mean value and σ is the standard deviation. This function is often used as a model to random continuous variables that describe physical phenomena by adjusting the values of µ and σ. It is symmetric around µ, and σ measures the spread. The ULLN distribution function is obtained from the Gaussian distribution function by means of a transformed variable, η, and is given by

f(η) )

δ exp[-δ2η2] xπ

-∞ < η < +∞

(9)

where δ is a parameter and η is related to ξ by

(

η ) ln

)

aξ -ξ

d/max

0 < ξ < d/max

(10)

/ being two additional parameters. The with a and dmax ULLN distribution function, expressed in terms of ξ, can be obtained by means of a statistics theorem for a function of random variables that relates the one-toone correspondence of η and ξ, that is

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001 3521

dξ f(η) ) f(ξ) | | dη

(11)

Using eqs 9-11 results in

f(ξ) )

δd/max

xπξ(d/max

- ξ)

[

exp -δ2 ln2

(

)]

aξ -ξ

d/max

0 < ξ < d/max (12)

List of Symbols a ) parameter of the ULLN distribution function C ) mean coordination number di ) pore diameter dmax ) distribution maximum diameter dmin ) distribution minimum diameter d/max ) parameter of the ULLN distribution function fN ) number(-based) diameter probability density function fV ) volume(-based) diameter probability density function FN ) number(-based) diameter probability distribution function FV ) volume(-based) diameter probability distribution function i ) pore number max Ik ) estimate of the integral ∫ddmin ξ2fN(ξ) dξ j ) discretization interval J ) number of discretization intervals l ) pore length mx ) number of nodes in the x direction my ) number of nodes in the y direction mz ) number of nodes in the z direction M ) total number of nodes k ) estimated number for Ik NExt ) number of external pores NInt ) number of internal pores NT ) total number of pores P ) pore mercury pressure REI ) NExt/NInt S ) square root of the sum of the squares of the residuals of each curve relative to the mean value V ) pore volume Vres Hg ) dimensionless residual mercury volume VT ) total pore volume x ) Cartesian coordinate y ) Cartesian coordinate z ) Cartesian coordinate Greek Letters δ ) parameter of the ULLN distribution function η ) parameter of the ULLN distribution function θ ) contact angle ΘExt ) number of external pores removed ΘInt ) number of internal pores removed µ ) distribution mean value ξ ) dummy variable (diameter) σ ) distribution standard deviation σmv ) interfacial tension Superscripts ap ) apparent in ) initial

Literature Cited (1) Payatakes, A. C.; Dias, M. M. Immiscible microdisplacement and ganglion dynamics in porous media. Rev. Chem. Eng. 1984, 2, 85-174. (2) McGreavy, C.; Andrade, J. S., Jr.; Rajagopal, K. Size exclusion chromatography in pore networks. Chromatographia 1990, 30, 639-644.

(3) Mann, R. Developments in chemical reaction engineering: issues relating to particle pore structure and porous materials. Trans. Inst. Chem. Eng. 1993, 71, 551-562. (4) Sahimi, M.; Gavalas, G. R.; Tsotis, T. T. Statistical and continuum models of fluid-solid reactions in porous media. Chem. Eng. Sci. 1990, 45, 1443-1502. (5) Conner, W. C.; Lane, A. M. Measurement of the morphology of high surface area solids: effect of network structure on the simulation of porosimetry. J. Catal. 1984, 89, 217-225. (6) Tsetsekou, A.; Androutsopoulos, G. P.; Mann, R. Mercury porosimetry hysteresis and entrapment predictions based on a corrugated random pore model. Chem. Eng. Commun. 1991, 110, 1-29. (7) Fatt, I. The Network model of porous media. I. Capillary pressure characteristics. AIME Pet. Trans. 1956, 207, 144-159. (8) Fatt, I. The network model of porous media. II. Dynamic properties of a single size tube network. AIME Pet. Trans. 1956, 207, 160-163. (9) Fatt, I. The network model of porous media. III. Dynamic properties of networks with tube radius distribution. AIME Pet. Trans. 1956, 207, 164-181. (10) 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. (11) Mann, R.; Golshan, H. Application of a stochastic network pore model to a catalyst pellet. Chem. Eng. Commun. 1981, 12, 377-391. (12) Dias, M. M.; Payatakes, A. C. Network models for twophase flow in porous media. Part 1. Immiscible microdisplacement of nonwetting fluids. J. Fluid Mech. 1986, 164, 305-336. (13) Dias, M. M.; Payatakes, A. C. Network models for twophase flow in porous media. Part 2. Motion of oil ganglia. J. Fluid Mech. 1986, 164, 337-358. (14) Constantinides, G. N.; Payatakes, A. C. Network simulation of steady-state two-phase flow in consolidated porous media. AIChE J. 1996, 42, 369-382. (15) Hollewand, M. P.; Gladden, L. F. Modelling of diffusion and reaction in porous catalysts using a random three-dimensional network model. Chem. Eng. Sci. 1992, 47, 1761-1770. (16) Andrade, J. S., Jr.; Rajagopal, K.; McGreavy, C. Chromatography in pore networks IIsThe role of structure and adsorption in the band broadening. Chromatographia 1991, 32, 345-349. (17) Andrade, J. S., Jr.; Shibusa, Y.; Arai, Y.; McGreavy, C. A network model for diffusion and adsorption in compacted pellets of bidispersed grains. Chem. Eng. Sci. 1995, 50, 1943-1951. (18) 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. (19) Washburn, E. W. The dynamics of capillary flow. Phys. Rev. 1921, 17, 273-283. (20) 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. (21) Androutsopoulos, G. P.; Mann, R. Evaluation of mercury porosimeter experiments using a network pore structure model. Chem. Eng. Sci. 1979, 34, 1203-1212. (22) Cox, M. A. A. Mercury injection measurements used in the prediction of rock pore dimensions employing a crystalline lattice of capillaries. Chem. Eng. Sci. 1991, 46, 167-172. (23) Portsmouth, R. L.; Gladden, L. F. Determination of pore connectivity by mercury porosimetry. Chem. Eng. Sci. 1991, 46, 3022-2036. (24) Portsmouth, R. L.; Gladden, L. F. Mercury porosimetry as a probe of pore connectivity. Trans. Inst. Chem. Eng. 1992, 70, 63-70. (25) 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. (26) Tsakiroglou, C. D.; Payatakes, A. C. Effects of pore-size correlations on mercury porosimetry curves. J. Colloid Interface Sci. 1991, 146, 479-494. (27) Constantinides, G. N.; Payatakes, A. C. A three-dimensional network model for consolidated porous media. Basic studies. Chem. Eng. Commun. 1989, 81, 55-81. (28) Mata, V. G. Characterization of porous media: porosimetry, 3D simulation and serial tomography. Application to catalytic support. Ph.D. Thesis, University of Porto, Porto, Portugal, 1998.

3522

Ind. Eng. Chem. Res., Vol. 40, No. 16, 2001

(29) Zydney, A. L.; Aimar, P.; Meireles, M.; Pimbley, J. M.; Belfort, G. Use of log-normal probability density function to analyze membrane pore size distributions: functional forms and discrepancies. J. Membr. Sci. 1994, 91, 293-298. (30) Tatterson, D. F.; Dallman, J. C.; Hanratty, T. J. Drop sizes in annular gas-liquid flows. AIChE J. 1977, 23, 68-76.

(31) Lopes, J. C. B.; Dukler, A. E. Droplet Dynamics in GasLiquid Vertical Annular Flow. AIChE J. 1987, 33, 1013-1023.

Received for review February 1, 2001 Accepted May 22, 2001 IE010112E