- Email: [email protected]
Network modelling of capillary pressure curves, permeability, and diffusivity
Chemical Engineering Science 62 (2007) 5117 – 5122 www.elsevier.com/locate/ces
Network modelling of capillary pressure curves, permeability, and diffusivity a,∗ ˇ P. Capek , V. Hejtmánek b , L. Brabec c , A. Zikánová c , M. Koˇciˇrík c a Department of Organic Technology, Institute of Chemical Technology, Technická 5, 16628 Prague, Czech Republic b Academy of Sciences of the Czech Republic, Institute of Chemical Process Fundamentals, Rozvojová 135, 16502 Prague, Czech Republic c Academy of Sciences of the Czech Republic, J. Heyrovský Institute of Physical Chemistry, Dolejškova 3, 18223 Prague, Czech Republic
Received 15 June 2006; received in revised form 21 December 2006; accepted 10 January 2007 Available online 25 January 2007
Abstract A pore network model is presented to predict permeability and diffusivity in porous bodies with a relatively high porosity. The model application is exemplified on a macroporous sample of -alumina with the porosity of ≈ 0.4. The network model is constructed on the basis of proximity of the computed data on its total porosity, pore-size function, and simulated mercury intrusion curve to the respective experimental data. The experiment related pore-size function is extracted from a 3-D stochastic replica of the -alumina sample obtained by stochastic reconstruction. The reconstruction technique employs morphological information based on a set of 2-D cuts through the porous medium. The only free parameter in the network construction is connectivity. The impact of the connectivity adjustment on the transport properties is studied. If the calculated mercury intrusion curve is forced to fit the experimental one, changes in connectivity are counterbalanced by the presence of slightly wider or narrower throats in the network. This compensation effect decreases the span of calculated values of permeability and diffusivity, which agree well with experiment. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Diffusion; Laminar flow; Microstructure; Porous media; Pore network model; Stochastic reconstruction
1. Introduction The prediction and estimation of macroscopic properties of porous media are long-standing problems of great interest for oil recovery, hydrology, catalysis, membrane separations, and the intelligent design of adsorbents, membranes, and catalysts. The dependence of macroscopic properties of porous media, such as drainage and imbibition capillary curves, permeability, and diffusivity, on the microscopic pore structure makes the use of pore structure models essential (Dullien, 1992; Sahimi, 1995). A pore structure model must involve those geometrical and topological characteristics of porous media, which have the major effect on the equilibrium or transport properties under consideration, while retaining the relative mathematical simplicity. It is generally accepted that the most realistic models of porous media are those that approximate the pore space as a
∗ Corresponding author. Tel.: +420 220444213; fax: +420 220444340.
ˇ E-mail address: [email protected] (P. Capek). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.01.011
network of pores, in which relatively large pores (chambers or cavities) are connected through relatively narrow ones (necks or throats). For the reliable prediction of macroscopic properties of a real porous medium, its pore space has to be reconstructed. The stochastic reconstruction is one of advanced experimental techniques which may be used for the reconstruction of the pore space in the three dimensions. It requires only photomicrographs of many 2-D cuts through a porous medium. Given that the photomicrographs of the pore space are available, the key problem in the stochastic reconstruction is to identify microstructural descriptors that suffice to reproduce the pore space of the selected porous medium (Adler et al., 1990; Torquato, 2002). The macroscopic transport properties of porous solids, such as permeability, are strongly dependent on the throat size distribution (TSD). In this context Tsakiroglou and Payatakes (2000) indicated that the estimation of the TSD with the aid of the 3-D reconstruction had some shortcomings. They proposed the combination of the 3-D reconstruction with mercury
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5118
porosimetry for the purpose of determination of the “actual” TSD. The principal aim of this study is therefore to employ Tsakiroglou and Payatakes’ (2000) idea for the development of a pore network model, whose parameters would be partially extracted from the stochastic replica of a macroporous solid (an -alumina sample) and also adjusted on the basis of mercury intrusion. The other aim of this work is to formulate a simulator with the aid of which the permeation and diffusion flows will be predicted in the pore network. 2. Stochastic reconstruction The structure of a two-phase digitized porous medium is completely defined (Torquato, 2002) in terms of the indicator function for the void phase, I( x ), 1 if x belongs to the pore space, I( x) = (1) 0 otherwise, where x is the position vector. In the 2-D representation of a porous medium, the position vector x takes discrete values x = (ia, ja) determined by a square lattice with i = 0, . . . , lx − 1 and j = 0, . . . , ly − 1 and the size lx × ly . Similarly, in the 3-D digitized media the discrete values of the position vector x = (ia, ja, ka) are defined by a simple cubic lattice and the size lx × ly × lz . The lattice constant, a, corresponds to the pixel or voxel resolution. The focus of this study is on isotropic porous media characterized by the following microstructural descriptors: (i) the total porosity, , (ii) the autocorrelation function, (u), (iii) the lineal-path function, L(u), and (iv) the pore-size probability density function, P(). The definitions of these descriptors can be found, for instance, in the book by Torquato (2002). For the isotropic porous media, the total porosity is a constant and both the autocorrelation function as well as the lineal-path function depend only on the distance u. Moreover, the descriptors , (u), and L(u), can be obtained from any planar cut through the isotropic medium. On the other hand, the pore-size function, P(), is an intrinsically 3-D descriptor and, therefore, cannot be obtained from planar cuts.
2.1. Evaluation of reference microstructural descriptors Cylindrical pellets of -alumina (Carborundum Electrite, CZ, 5 mm height × 5 mm diameter) were impregnated with an epoxy resin under vacuum. The hardened epoxy resin blocks were cut, ground, and polished to achieve the smooth surface. Series of back-scatter scanning electron images of crosssections through the medium were recorded with an appropriate size (1280 × 960 pixels) and resolution (1 × 1 m2 ). Ten images selected for analysis were subjected to a sequence of operations, which involved filtering and image segmentation, i.e., all grey pixels were separated into black (the pore space) and white (the solid matrix). The segmented images were sampled in the two principal directions. Both the two-point probability function and the lineal-path function were nearly independent of the direction chosen. The reference functions (Fig. 1) for the purpose of the 3-D stochastic reconstruction were determined by averaging individual functions obtained from the selected images. 2.2. Simulated annealing Simulated annealing, an optimization technique originally developed for combinatorial optimization problems (Press et al., 1992), has been employed by many researchers (Torquato, 2002; Manwart et al., 2000; Talukdar et al., 2002) for reconstruction of 3-D random media. The reconstruction consists in the gradual transformation of a high-energy configuration of pore and solid voxels into a minimum energy configuration. The objective function or “energy” for the purpose of this work was defined as the sum of squared deviations between the ref˜ erence functions, ˜ (u) and L(u), and the functions i (u) and Li (u) calculated for 3-D replicas of the porous medium. These functions were calculated in the three principal directions applying periodic boundary conditions. The deviations between the reference and calculated functions beyond the distance of 130 m were considered to be negligible because the reference autocorrelation function first got its long-range value of 0 at the distance of 86 m and the lineal-path function did not deviate
1.0
(u)/φ
~
χ(u)
1.0
0.5
0.0
0.5
0.0 0
50
100 u, μm
0
50
100 u, μm
˜ (u), divided by the total porosity, . Fig. 1. Reference microstructural descriptors: the autocorrelation function, ˜ (u) and the lineal-path function, L
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
significantly from 0 beyond the distance of 130 m, see Fig. 1. A 3-D replica of -alumina was realized on the simple cubic lattice of size lx × ly × lz = 256 × 256 × 256. An approximation to P() was computed by placing a random point in the reconstructed pore space and measuring its distance to the nearest point on the phase interface, which was modelled by faces of cubic voxels. This procedure was repeated for a large number of random points in the pore space. The pore-size probability density function was then obtained by discretizing the distances and dividing by the total number of random points placed into the pore space. Integration of the resulting function gave an approximation ∞ of the cumulative pore-size distribution function, F(0 ) = 0 P() d. 3. Pore network model For the sake of simplicity, we model the pore structure as a cubic network of cubic chambers interconnected through throats of rectangular cross-section, see Fig. 2. The general procedure of the cubic network construction can be very briefly described as follows. After assignment of chamber sizes to network nodes and throat sizes to network bonds, the distance between the centres of gravity of all pairs of chambers was adjusted so that the total porosity of the network could match that of prototype; both the chambers and the throats contributed to the total pore volume. If, during the network construction, two adjacent chambers were found to intersect, one of them was randomly interchanged with another chamber. A similar interchange was made if a throat was found to be larger than one of the chambers at its ends. If the chamber size distribution (CSD) and TSD overlap, there is a certain degree of correlation between sizes of chambers and throats. The simple cubic network has regular topology with the coordination number of 6. However, if the network was to be of irregular topology,
Fig. 2. Schematic representation of a cubic chamber of volume, Vc , with three throats of rectangular cross-section. Throat dimensions, at and bt , are limited 1/3 by the chamber dimension, 2bc =Vc : at bc and bt < bc (a self-consistency limitation). The throat aspect ratio, = at /bt , satisfies the inequality: 1.
5119
the mean coordination number, , ¯ of the network had to be reduced by removing bonds selected at random. Then, the coordination numbers of individual chambers, , were equal to random integers between 1 and 6. In this work the network connectivity was the only free parameter characterized by themean coordination number. Since the mean pore size, ∞ = 0 F() d, was equal to 5.17 m, the network size of mx × my × mz was set up so that the outer particle size, dp , could match the size of the real particles without excessive demands on computer resources.
3.1. Simulation of the mercury intrusion The mercury porosimetry is a widely used method of measuring pore-size distributions. The common model of the mercury porosimetry, also accepted in the present work, assumes that the mercury intrusion is solely controlled by the throat diameters (Sahimi, 1995). The critical capillary pressure for the mercury intrusion into a pore throat of rectangular cross-section was calculated readily using the approximate method proposed by Mason and Morrow (1984) and the intrusion contact angle of 2/9. A mercury intrusion algorithm which was accepted for the purpose of this work is usually called invasion percolation (Sahimi, 1995). In all the mercury porosimetry simulations the simulated cubic particle was immersed in mercury so that mercury could invade the pore space through its six faces.
3.2. Estimation of the pore network parameters The influence of the network topology and the TSD on mass transport was investigated for three networks differing in the mean coordination number. Its typical values of 4, 4.5, and 5.0 were considered for the purpose of this work. The construction of the chamber-and-throat network was performed iteratively: (1) The domain of the CSD was derived from the domain of the cumulative pore-size distribution function (Fig. 4). (2) An initial approximation of the CSD, based on the previous ˇ experience (Capek and Hejtmánek, 2004), was accepted. (3) A first estimate of the TSD was obtained by differentiating the mercury intrusion curve. (4) The pore space of the network model was sampled to get the approximation of F(). This function was compared with the reference function characterizing the reconstructed pore space (Fig. 4). (5) Mercury intrusion was simulated and the resulting capillary pressure curve was compared with the experimental (reference) one. (6) If either the calculated F() or the simulated capillary pressure curve did not match the reference one, parameters of the CSD and TSD were altered and a new realization of the pore network was generated. (7) This sequence of the steps (4)–(6) was repeated until the calculated characteristics did fit well the reference ones.
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5120
ϖ = 4.0 ϖ = 4.5 ϖ = 5.0
1 1 -1
TSD
ϕ (bc)×10, μm
ϕ (bt)×10, μm
-1
CSD
0
0 0
10
20 bc, μm
30
0
10
20 bt, μm
30
Fig. 3. The number density function of CSD, (bc ), and the number density functions of TSDs, (bt ), for the final realizations of the networks differing in the mean coordination numbers.
Fig. 4. Comparison of the calculated (−) and experimental (·) curves. Points in the right-hand figure represent the cumulative pore-size distribution function, ¯ = 4.5. F(0 ), of the 3-D replica of -alumina. Both the capillary pressure curve and pore-size function were calculated by using the pore network with
Since the volume of individual pores is proportional to the third power of their characteristic sizes, a large fraction of the total pore volume has to be located in chambers of the chamberand-throat network. Therefore, the CSD was particularly expected to account for the course of the cumulative pore-size distribution function, F(0 ). The preliminary construction revealed that the maximum of the CSD was skewed strongly, with a long tail extending out towards more positive bc . This CSD depicted in Fig. 3 with an asymmetric tail was modelled by a combination of two beta distributions. A similar combination of two beta distributions was found to be useful for modelling the TSD, which had the positive skewness as well. The choice of the mean coordination number had a certain effect on the TSD. Since the simulated capillary pressure curves had to fit the experimental curve in the same way, three “optimum” TSDs were obtained, see Fig. 3. The worse accessibility of the internal pores due to lower connectivity was counterbalanced by the presence of slightly wider throats in the TSD.
The pore network of mx × my × mz = 92 × 92 × 92 nodes provided the cubic particle whose external size, dp , was very close to the real particle diameter and height (5 mm × 5 mm). This network size also ensured the excellent statistical stability of simulated quantities, which was confirmed by a few realizations of each network model. Fig. 4 illustrates how the calculated pore-size function and the mercury intrusion curve fit the experimental data. The long tail of the TSD was necessary to get the good fit to a low-pressure part (P < 48 kPa) of the intrusion curve. For the pore networks of the other mean coordination numbers, fits of similar quality were obtained. 4. Mass transport in pore networks The analysis was focused on an isothermal steady state regime of mass transport involving bulk diffusion, Knudsen diffusion, viscous flow, Knudsen flow, and also the transport in the transition region (Dullien, 1992; Keil, 1999). The
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
network description of mass transport was formulated via molar flows of gases and their binary mixtures through “single” pores. Each single pore consisted of three segments, a throat of rectangular cross-section and a half of each adjacent cubic chamber. The molar flow rates in the single pores were calculated on the basis of state variables located in chamber centres. These state variables (unknown node pressures or mole fractions) were determined by simultaneous solution of 92 × 92 × 92 = 778 688 non-linear equations, which resulted from application of the mass conservation law to each chamber, namely
n˙ ij = 0,
j = 1, . . . , 778 688,
(2)
i
where n˙ ij are the molar flow rates between the jth chamber and the adjacent chambers denoted by appropriate subscripts i. The simulated porous particle had the shape of a cube with the two opposite faces open to the gaseous phase and the remaining four faces closed in order to resemble our experimental setup. Boundary conditions of the Dirichlet type were used in the direction of macroscopic transport. The system of equations (2) was solved iteratively by a Newton-type method (Press et al., 1992). For the sake of simplicity, expressions for the molar flow rates, n˙ ij , were derived for a single pore between two arbitrary nodes 1 and 2. 4.1. Permeation flow in a single pore
w
(1) () ct
(1) n˙ 12 = ⎣gk + 16 + + 3 3 4 bc1 at bt3 bc2 ×
−1
⎤ P1 + P 2 ⎦ 2
P1 − P2 , Rg T
4.3. Comparison of experimental and simulated transport properties Experimental pressure dependences of the effective permeability, B(P¯ ), on the mean pressure in the particle, P¯ =[P (0)+ P (dp )]/2, were determined by the quasi-stationary permeation method applying small pressure differences P =P (0)−P (dp ) across the porous body (Fott and Petrini, 1982) for H2 , He, N2 , and Ar. The experimental dependences were fitted by the following correlation:
¯ 2 P B(P¯ ) = w + 3 8
(8)
8464 Rg T n˙ 0j , dp P
(9)
j =1
(3)
(4)
The dependence of the shape factor, (), on the throat aspect ratio, published by Perry et al. (1984) in the tabular form, was approximated by the power expansion in −1
() = 12 + 7.37043−1 + 5.47478−2 + 3.60699−3 .
(6)
where y1 and y2 are the mole fractions of A in the chambers 1 and 2, respectively, gk is given by (4), DAB √ is the binary diffusion coefficient of the pair A–B, and = [MA /MB ] − 1. For the single-pore conductance, gd , the following expression was derived 1 ct 1 −1 gd = 4 + + . (7) bc1 at bt bc2
B(P¯ ) =
where P1 and P2 stand for the pressures in the nodes 1 and 2, respectively, w is the mean arithmetic velocity of molecules (atoms) of the gas, is its dynamic viscosity, Rg is the gas constant, and T is the thermodynamic temperature. The quantity gk (single-pore conductance for Knudsen flow) was defined by √ −1 a t bt 3 ct gk = 4at bt 1 − + √ . bc1 bc2 16 at bt
expressed in this case as 4 1 y1 + y2 −1 + + n˙ 12 = g k wA gd DAB gd DAB 2 P (y1 − y2 ), × Rg T
and the transport parameters and were evaluated and averaged over the set of those four gases. The effective permeability of the pore networks was defined as follows:
In this case the molar flow rate n˙ 12 was given by ⎡
5121
(5)
where n˙ 0j refers to the molar flow rate between the bulk gaseous phase and the node j adjacent to the external particle surface (92 × 92 = 8464 nodes). Increasing the connectivity of the network increased the effective network conductance in the Knudsen region ( ↑) as well as in the region of continuum ( ↑). The predicted values agreed quite well with the experimental data. Specifically, the choice of ¯ = 4.50 provided the pore network revealing the transport properties, which almost matched the experimental ones (Table 1). Simulations of diffusion flow in the binary mixture of hydrogen (component A) and nitrogen (component B) were repeated for 20 discrete levels of the total pressure P in the wide interval 0.01, 500 kPa. By analogy with (9), the effective flux of hydrogen, N, was defined as N = 8464 ˙ 0j /dp2 . The resultj =1 n ing dependence of N on the total pressure was approximated by a functional relationship, which is commonly used (Dullien, 1992) for the reduction of experimental data k + y(0) 1 + DAB /DA P DAB ln , k + y(d ) Rg T dp 1 + DAB /DA p
4.2. Diffusion flow in a single pore
N (P ) =
Diffusion molar flow of a gas A in a binary mixture A–B which takes place under isobaric conditions (P = const.) was
where y(0) and y(dp ) are the mole fractions of hydrogen at k = 2 w ) the particle boundaries. The parameters and (DA A 3
(10)
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5122
Table 1 Calculated and experimental transport properties of -alumina estimated from single gas flow (S) and binary diffusion flow (B). Flow
Parameter
Experiment
Pore network
–
¯
–
4.00
4.50
5.00
S
(m)
(m 2 )
1.46 ± 0.14 21.4 ± 0.28
1.30 16.3
1.50 19.9
1.72 23.4
B
(m)
0.124 –
0.099 1.06
0.120 1.22
0.141 1.41
characterize the transport properties of the porous medium in the regimes of ordinary and Knudsen diffusion, respectively. Diffusion experiments were performed in the single-pellet Wicke–Kallenbach cell, which was operated at the atmospheric pressure. Since the Knudsen number was small under these conditions (prevailing bulk diffusion in macropores), the only parameter could be estimated reliably from a large set of experimental data including binary diffusion in various mixtures of H2 , He, N2 , and Ar. Again, the network calculations were in a good agreement with the experiments (Table 1). 5. Conclusions The stochastic reconstruction of a macroporous -alumina from limited morphological information provided the basis for the construction of a 3-D random pore network. The pore network model and the simulator of the mercury intrusion enabled the fine adjustment of the TSD with regard to the experimental mercury intrusion curve. Since connectivity was not extracted from the 3-D replica of the porous medium, the mean coordination number was the only free parameter in the network construction. It was shown that the mercury intrusion curve could be reproduced very well and independently of the mean coordination number. The change in connectivity was compensated by small changes in the TSD. Due to this correlation of the TSD and connectivity, the pore networks of the three different mean coordination numbers revealed only a small span of the calculated permeability and diffusivity. The predicted permeability and diffusivity agreed well with the experimental data. Further extension of this work will focus on the extraction of topological and geometrical information from 3-D replicas of porous solids. Acknowledgement The authors gratefully acknowledge the support by the Grant Agency of the Czech Republic (Research Grant # 203/05/0347).
References Adler, P.M., Jacquin, C.G., Quiblier, J.A., 1990. Flow in simulated porous media. International Journal of Multiphase Flow 16 (4), 691–712. ˇ Capek, P., Hejtmánek, V., 2004. A relationship between capillary pressure and permeability as revealed by a pore network model. In: 16th International Congress of Chemical and Process Engineering CHISA 2004. Prague, Czech Republic, p. P3.194. Dullien, F.A.L., 1992. Porous Media: Fluid Transport and Pore Structure. second ed. Academic Press, San Diego. Fott, P., Petrini, G., 1982. Determination of transport parameters of porous catalysts from permeation measurements. Applied Catalysis 2, 367–378. Keil, F., 1999. Diffusion und Chemische Reaktionen in der Gas/Feststoff–Katalyse. Springer, Berlin. Manwart, C., Toquato, S., Hilfer, R., 2000. Stochastic reconstruction of sandstones. Physical Review E 62 (1), 893–899. Mason, G., Morrow, N.R., 1984. Meniscus curvatures in capillaries of uniform cross-section. Journal of the Chemical Society—Faraday Transactions 1 80, 2375–2393. Perry, R.H., Green, D.W., Maloney, J.O., 1984. Perry’s Chemical Engineers’ Handbook. sixth ed. McGraw-Hill, New York. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in C: The Art of Scientific Computing. second ed. Cambridge University Press, Cambridge. Sahimi, M., 1995. Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches. VCH, Weinheim, Germany. Talukdar, M.S., Torsaeter, O., Ioannidis, M.A., Howard, J.J., 2002. Stochastic reconstruction of chalk from 2D images. Transport in Porous Media 48, 101–123. Torquato, S., 2002. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Tsakiroglou, C.D., Payatakes, A.C., 2000. Characterization of the pore structure reservoir rocks with the aid of serial sectioning analysis, mercury porosimetry and network simulation. Advances in Water Resources 23, 773–789.
Network modelling of capillary pressure curves, permeability, and diffusivity a,∗ ˇ P. Capek , V. Hejtmánek b , L. Brabec c , A. Zikánová c , M. Koˇciˇrík c a Department of Organic Technology, Institute of Chemical Technology, Technická 5, 16628 Prague, Czech Republic b Academy of Sciences of the Czech Republic, Institute of Chemical Process Fundamentals, Rozvojová 135, 16502 Prague, Czech Republic c Academy of Sciences of the Czech Republic, J. Heyrovský Institute of Physical Chemistry, Dolejškova 3, 18223 Prague, Czech Republic
Received 15 June 2006; received in revised form 21 December 2006; accepted 10 January 2007 Available online 25 January 2007
Abstract A pore network model is presented to predict permeability and diffusivity in porous bodies with a relatively high porosity. The model application is exemplified on a macroporous sample of -alumina with the porosity of ≈ 0.4. The network model is constructed on the basis of proximity of the computed data on its total porosity, pore-size function, and simulated mercury intrusion curve to the respective experimental data. The experiment related pore-size function is extracted from a 3-D stochastic replica of the -alumina sample obtained by stochastic reconstruction. The reconstruction technique employs morphological information based on a set of 2-D cuts through the porous medium. The only free parameter in the network construction is connectivity. The impact of the connectivity adjustment on the transport properties is studied. If the calculated mercury intrusion curve is forced to fit the experimental one, changes in connectivity are counterbalanced by the presence of slightly wider or narrower throats in the network. This compensation effect decreases the span of calculated values of permeability and diffusivity, which agree well with experiment. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Diffusion; Laminar flow; Microstructure; Porous media; Pore network model; Stochastic reconstruction
1. Introduction The prediction and estimation of macroscopic properties of porous media are long-standing problems of great interest for oil recovery, hydrology, catalysis, membrane separations, and the intelligent design of adsorbents, membranes, and catalysts. The dependence of macroscopic properties of porous media, such as drainage and imbibition capillary curves, permeability, and diffusivity, on the microscopic pore structure makes the use of pore structure models essential (Dullien, 1992; Sahimi, 1995). A pore structure model must involve those geometrical and topological characteristics of porous media, which have the major effect on the equilibrium or transport properties under consideration, while retaining the relative mathematical simplicity. It is generally accepted that the most realistic models of porous media are those that approximate the pore space as a
∗ Corresponding author. Tel.: +420 220444213; fax: +420 220444340.
ˇ E-mail address: [email protected] (P. Capek). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.01.011
network of pores, in which relatively large pores (chambers or cavities) are connected through relatively narrow ones (necks or throats). For the reliable prediction of macroscopic properties of a real porous medium, its pore space has to be reconstructed. The stochastic reconstruction is one of advanced experimental techniques which may be used for the reconstruction of the pore space in the three dimensions. It requires only photomicrographs of many 2-D cuts through a porous medium. Given that the photomicrographs of the pore space are available, the key problem in the stochastic reconstruction is to identify microstructural descriptors that suffice to reproduce the pore space of the selected porous medium (Adler et al., 1990; Torquato, 2002). The macroscopic transport properties of porous solids, such as permeability, are strongly dependent on the throat size distribution (TSD). In this context Tsakiroglou and Payatakes (2000) indicated that the estimation of the TSD with the aid of the 3-D reconstruction had some shortcomings. They proposed the combination of the 3-D reconstruction with mercury
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5118
porosimetry for the purpose of determination of the “actual” TSD. The principal aim of this study is therefore to employ Tsakiroglou and Payatakes’ (2000) idea for the development of a pore network model, whose parameters would be partially extracted from the stochastic replica of a macroporous solid (an -alumina sample) and also adjusted on the basis of mercury intrusion. The other aim of this work is to formulate a simulator with the aid of which the permeation and diffusion flows will be predicted in the pore network. 2. Stochastic reconstruction The structure of a two-phase digitized porous medium is completely defined (Torquato, 2002) in terms of the indicator function for the void phase, I( x ), 1 if x belongs to the pore space, I( x) = (1) 0 otherwise, where x is the position vector. In the 2-D representation of a porous medium, the position vector x takes discrete values x = (ia, ja) determined by a square lattice with i = 0, . . . , lx − 1 and j = 0, . . . , ly − 1 and the size lx × ly . Similarly, in the 3-D digitized media the discrete values of the position vector x = (ia, ja, ka) are defined by a simple cubic lattice and the size lx × ly × lz . The lattice constant, a, corresponds to the pixel or voxel resolution. The focus of this study is on isotropic porous media characterized by the following microstructural descriptors: (i) the total porosity, , (ii) the autocorrelation function, (u), (iii) the lineal-path function, L(u), and (iv) the pore-size probability density function, P(). The definitions of these descriptors can be found, for instance, in the book by Torquato (2002). For the isotropic porous media, the total porosity is a constant and both the autocorrelation function as well as the lineal-path function depend only on the distance u. Moreover, the descriptors , (u), and L(u), can be obtained from any planar cut through the isotropic medium. On the other hand, the pore-size function, P(), is an intrinsically 3-D descriptor and, therefore, cannot be obtained from planar cuts.
2.1. Evaluation of reference microstructural descriptors Cylindrical pellets of -alumina (Carborundum Electrite, CZ, 5 mm height × 5 mm diameter) were impregnated with an epoxy resin under vacuum. The hardened epoxy resin blocks were cut, ground, and polished to achieve the smooth surface. Series of back-scatter scanning electron images of crosssections through the medium were recorded with an appropriate size (1280 × 960 pixels) and resolution (1 × 1 m2 ). Ten images selected for analysis were subjected to a sequence of operations, which involved filtering and image segmentation, i.e., all grey pixels were separated into black (the pore space) and white (the solid matrix). The segmented images were sampled in the two principal directions. Both the two-point probability function and the lineal-path function were nearly independent of the direction chosen. The reference functions (Fig. 1) for the purpose of the 3-D stochastic reconstruction were determined by averaging individual functions obtained from the selected images. 2.2. Simulated annealing Simulated annealing, an optimization technique originally developed for combinatorial optimization problems (Press et al., 1992), has been employed by many researchers (Torquato, 2002; Manwart et al., 2000; Talukdar et al., 2002) for reconstruction of 3-D random media. The reconstruction consists in the gradual transformation of a high-energy configuration of pore and solid voxels into a minimum energy configuration. The objective function or “energy” for the purpose of this work was defined as the sum of squared deviations between the ref˜ erence functions, ˜ (u) and L(u), and the functions i (u) and Li (u) calculated for 3-D replicas of the porous medium. These functions were calculated in the three principal directions applying periodic boundary conditions. The deviations between the reference and calculated functions beyond the distance of 130 m were considered to be negligible because the reference autocorrelation function first got its long-range value of 0 at the distance of 86 m and the lineal-path function did not deviate
1.0
(u)/φ
~
χ(u)
1.0
0.5
0.0
0.5
0.0 0
50
100 u, μm
0
50
100 u, μm
˜ (u), divided by the total porosity, . Fig. 1. Reference microstructural descriptors: the autocorrelation function, ˜ (u) and the lineal-path function, L
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
significantly from 0 beyond the distance of 130 m, see Fig. 1. A 3-D replica of -alumina was realized on the simple cubic lattice of size lx × ly × lz = 256 × 256 × 256. An approximation to P() was computed by placing a random point in the reconstructed pore space and measuring its distance to the nearest point on the phase interface, which was modelled by faces of cubic voxels. This procedure was repeated for a large number of random points in the pore space. The pore-size probability density function was then obtained by discretizing the distances and dividing by the total number of random points placed into the pore space. Integration of the resulting function gave an approximation ∞ of the cumulative pore-size distribution function, F(0 ) = 0 P() d. 3. Pore network model For the sake of simplicity, we model the pore structure as a cubic network of cubic chambers interconnected through throats of rectangular cross-section, see Fig. 2. The general procedure of the cubic network construction can be very briefly described as follows. After assignment of chamber sizes to network nodes and throat sizes to network bonds, the distance between the centres of gravity of all pairs of chambers was adjusted so that the total porosity of the network could match that of prototype; both the chambers and the throats contributed to the total pore volume. If, during the network construction, two adjacent chambers were found to intersect, one of them was randomly interchanged with another chamber. A similar interchange was made if a throat was found to be larger than one of the chambers at its ends. If the chamber size distribution (CSD) and TSD overlap, there is a certain degree of correlation between sizes of chambers and throats. The simple cubic network has regular topology with the coordination number of 6. However, if the network was to be of irregular topology,
Fig. 2. Schematic representation of a cubic chamber of volume, Vc , with three throats of rectangular cross-section. Throat dimensions, at and bt , are limited 1/3 by the chamber dimension, 2bc =Vc : at bc and bt < bc (a self-consistency limitation). The throat aspect ratio, = at /bt , satisfies the inequality: 1.
5119
the mean coordination number, , ¯ of the network had to be reduced by removing bonds selected at random. Then, the coordination numbers of individual chambers, , were equal to random integers between 1 and 6. In this work the network connectivity was the only free parameter characterized by themean coordination number. Since the mean pore size, ∞ = 0 F() d, was equal to 5.17 m, the network size of mx × my × mz was set up so that the outer particle size, dp , could match the size of the real particles without excessive demands on computer resources.
3.1. Simulation of the mercury intrusion The mercury porosimetry is a widely used method of measuring pore-size distributions. The common model of the mercury porosimetry, also accepted in the present work, assumes that the mercury intrusion is solely controlled by the throat diameters (Sahimi, 1995). The critical capillary pressure for the mercury intrusion into a pore throat of rectangular cross-section was calculated readily using the approximate method proposed by Mason and Morrow (1984) and the intrusion contact angle of 2/9. A mercury intrusion algorithm which was accepted for the purpose of this work is usually called invasion percolation (Sahimi, 1995). In all the mercury porosimetry simulations the simulated cubic particle was immersed in mercury so that mercury could invade the pore space through its six faces.
3.2. Estimation of the pore network parameters The influence of the network topology and the TSD on mass transport was investigated for three networks differing in the mean coordination number. Its typical values of 4, 4.5, and 5.0 were considered for the purpose of this work. The construction of the chamber-and-throat network was performed iteratively: (1) The domain of the CSD was derived from the domain of the cumulative pore-size distribution function (Fig. 4). (2) An initial approximation of the CSD, based on the previous ˇ experience (Capek and Hejtmánek, 2004), was accepted. (3) A first estimate of the TSD was obtained by differentiating the mercury intrusion curve. (4) The pore space of the network model was sampled to get the approximation of F(). This function was compared with the reference function characterizing the reconstructed pore space (Fig. 4). (5) Mercury intrusion was simulated and the resulting capillary pressure curve was compared with the experimental (reference) one. (6) If either the calculated F() or the simulated capillary pressure curve did not match the reference one, parameters of the CSD and TSD were altered and a new realization of the pore network was generated. (7) This sequence of the steps (4)–(6) was repeated until the calculated characteristics did fit well the reference ones.
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5120
ϖ = 4.0 ϖ = 4.5 ϖ = 5.0
1 1 -1
TSD
ϕ (bc)×10, μm
ϕ (bt)×10, μm
-1
CSD
0
0 0
10
20 bc, μm
30
0
10
20 bt, μm
30
Fig. 3. The number density function of CSD, (bc ), and the number density functions of TSDs, (bt ), for the final realizations of the networks differing in the mean coordination numbers.
Fig. 4. Comparison of the calculated (−) and experimental (·) curves. Points in the right-hand figure represent the cumulative pore-size distribution function, ¯ = 4.5. F(0 ), of the 3-D replica of -alumina. Both the capillary pressure curve and pore-size function were calculated by using the pore network with
Since the volume of individual pores is proportional to the third power of their characteristic sizes, a large fraction of the total pore volume has to be located in chambers of the chamberand-throat network. Therefore, the CSD was particularly expected to account for the course of the cumulative pore-size distribution function, F(0 ). The preliminary construction revealed that the maximum of the CSD was skewed strongly, with a long tail extending out towards more positive bc . This CSD depicted in Fig. 3 with an asymmetric tail was modelled by a combination of two beta distributions. A similar combination of two beta distributions was found to be useful for modelling the TSD, which had the positive skewness as well. The choice of the mean coordination number had a certain effect on the TSD. Since the simulated capillary pressure curves had to fit the experimental curve in the same way, three “optimum” TSDs were obtained, see Fig. 3. The worse accessibility of the internal pores due to lower connectivity was counterbalanced by the presence of slightly wider throats in the TSD.
The pore network of mx × my × mz = 92 × 92 × 92 nodes provided the cubic particle whose external size, dp , was very close to the real particle diameter and height (5 mm × 5 mm). This network size also ensured the excellent statistical stability of simulated quantities, which was confirmed by a few realizations of each network model. Fig. 4 illustrates how the calculated pore-size function and the mercury intrusion curve fit the experimental data. The long tail of the TSD was necessary to get the good fit to a low-pressure part (P < 48 kPa) of the intrusion curve. For the pore networks of the other mean coordination numbers, fits of similar quality were obtained. 4. Mass transport in pore networks The analysis was focused on an isothermal steady state regime of mass transport involving bulk diffusion, Knudsen diffusion, viscous flow, Knudsen flow, and also the transport in the transition region (Dullien, 1992; Keil, 1999). The
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
network description of mass transport was formulated via molar flows of gases and their binary mixtures through “single” pores. Each single pore consisted of three segments, a throat of rectangular cross-section and a half of each adjacent cubic chamber. The molar flow rates in the single pores were calculated on the basis of state variables located in chamber centres. These state variables (unknown node pressures or mole fractions) were determined by simultaneous solution of 92 × 92 × 92 = 778 688 non-linear equations, which resulted from application of the mass conservation law to each chamber, namely
n˙ ij = 0,
j = 1, . . . , 778 688,
(2)
i
where n˙ ij are the molar flow rates between the jth chamber and the adjacent chambers denoted by appropriate subscripts i. The simulated porous particle had the shape of a cube with the two opposite faces open to the gaseous phase and the remaining four faces closed in order to resemble our experimental setup. Boundary conditions of the Dirichlet type were used in the direction of macroscopic transport. The system of equations (2) was solved iteratively by a Newton-type method (Press et al., 1992). For the sake of simplicity, expressions for the molar flow rates, n˙ ij , were derived for a single pore between two arbitrary nodes 1 and 2. 4.1. Permeation flow in a single pore
w
(1) () ct
(1) n˙ 12 = ⎣gk + 16 + + 3 3 4 bc1 at bt3 bc2 ×
−1
⎤ P1 + P 2 ⎦ 2
P1 − P2 , Rg T
4.3. Comparison of experimental and simulated transport properties Experimental pressure dependences of the effective permeability, B(P¯ ), on the mean pressure in the particle, P¯ =[P (0)+ P (dp )]/2, were determined by the quasi-stationary permeation method applying small pressure differences P =P (0)−P (dp ) across the porous body (Fott and Petrini, 1982) for H2 , He, N2 , and Ar. The experimental dependences were fitted by the following correlation:
¯ 2 P B(P¯ ) = w + 3 8
(8)
8464 Rg T n˙ 0j , dp P
(9)
j =1
(3)
(4)
The dependence of the shape factor, (), on the throat aspect ratio, published by Perry et al. (1984) in the tabular form, was approximated by the power expansion in −1
() = 12 + 7.37043−1 + 5.47478−2 + 3.60699−3 .
(6)
where y1 and y2 are the mole fractions of A in the chambers 1 and 2, respectively, gk is given by (4), DAB √ is the binary diffusion coefficient of the pair A–B, and = [MA /MB ] − 1. For the single-pore conductance, gd , the following expression was derived 1 ct 1 −1 gd = 4 + + . (7) bc1 at bt bc2
B(P¯ ) =
where P1 and P2 stand for the pressures in the nodes 1 and 2, respectively, w is the mean arithmetic velocity of molecules (atoms) of the gas, is its dynamic viscosity, Rg is the gas constant, and T is the thermodynamic temperature. The quantity gk (single-pore conductance for Knudsen flow) was defined by √ −1 a t bt 3 ct gk = 4at bt 1 − + √ . bc1 bc2 16 at bt
expressed in this case as 4 1 y1 + y2 −1 + + n˙ 12 = g k wA gd DAB gd DAB 2 P (y1 − y2 ), × Rg T
and the transport parameters and were evaluated and averaged over the set of those four gases. The effective permeability of the pore networks was defined as follows:
In this case the molar flow rate n˙ 12 was given by ⎡
5121
(5)
where n˙ 0j refers to the molar flow rate between the bulk gaseous phase and the node j adjacent to the external particle surface (92 × 92 = 8464 nodes). Increasing the connectivity of the network increased the effective network conductance in the Knudsen region ( ↑) as well as in the region of continuum ( ↑). The predicted values agreed quite well with the experimental data. Specifically, the choice of ¯ = 4.50 provided the pore network revealing the transport properties, which almost matched the experimental ones (Table 1). Simulations of diffusion flow in the binary mixture of hydrogen (component A) and nitrogen (component B) were repeated for 20 discrete levels of the total pressure P in the wide interval 0.01, 500 kPa. By analogy with (9), the effective flux of hydrogen, N, was defined as N = 8464 ˙ 0j /dp2 . The resultj =1 n ing dependence of N on the total pressure was approximated by a functional relationship, which is commonly used (Dullien, 1992) for the reduction of experimental data k + y(0) 1 + DAB /DA P DAB ln , k + y(d ) Rg T dp 1 + DAB /DA p
4.2. Diffusion flow in a single pore
N (P ) =
Diffusion molar flow of a gas A in a binary mixture A–B which takes place under isobaric conditions (P = const.) was
where y(0) and y(dp ) are the mole fractions of hydrogen at k = 2 w ) the particle boundaries. The parameters and (DA A 3
(10)
ˇ P. Capek et al. / Chemical Engineering Science 62 (2007) 5117 – 5122
5122
Table 1 Calculated and experimental transport properties of -alumina estimated from single gas flow (S) and binary diffusion flow (B). Flow
Parameter
Experiment
Pore network
–
¯
–
4.00
4.50
5.00
S
(m)
(m 2 )
1.46 ± 0.14 21.4 ± 0.28
1.30 16.3
1.50 19.9
1.72 23.4
B
(m)
0.124 –
0.099 1.06
0.120 1.22
0.141 1.41
characterize the transport properties of the porous medium in the regimes of ordinary and Knudsen diffusion, respectively. Diffusion experiments were performed in the single-pellet Wicke–Kallenbach cell, which was operated at the atmospheric pressure. Since the Knudsen number was small under these conditions (prevailing bulk diffusion in macropores), the only parameter could be estimated reliably from a large set of experimental data including binary diffusion in various mixtures of H2 , He, N2 , and Ar. Again, the network calculations were in a good agreement with the experiments (Table 1). 5. Conclusions The stochastic reconstruction of a macroporous -alumina from limited morphological information provided the basis for the construction of a 3-D random pore network. The pore network model and the simulator of the mercury intrusion enabled the fine adjustment of the TSD with regard to the experimental mercury intrusion curve. Since connectivity was not extracted from the 3-D replica of the porous medium, the mean coordination number was the only free parameter in the network construction. It was shown that the mercury intrusion curve could be reproduced very well and independently of the mean coordination number. The change in connectivity was compensated by small changes in the TSD. Due to this correlation of the TSD and connectivity, the pore networks of the three different mean coordination numbers revealed only a small span of the calculated permeability and diffusivity. The predicted permeability and diffusivity agreed well with the experimental data. Further extension of this work will focus on the extraction of topological and geometrical information from 3-D replicas of porous solids. Acknowledgement The authors gratefully acknowledge the support by the Grant Agency of the Czech Republic (Research Grant # 203/05/0347).
References Adler, P.M., Jacquin, C.G., Quiblier, J.A., 1990. Flow in simulated porous media. International Journal of Multiphase Flow 16 (4), 691–712. ˇ Capek, P., Hejtmánek, V., 2004. A relationship between capillary pressure and permeability as revealed by a pore network model. In: 16th International Congress of Chemical and Process Engineering CHISA 2004. Prague, Czech Republic, p. P3.194. Dullien, F.A.L., 1992. Porous Media: Fluid Transport and Pore Structure. second ed. Academic Press, San Diego. Fott, P., Petrini, G., 1982. Determination of transport parameters of porous catalysts from permeation measurements. Applied Catalysis 2, 367–378. Keil, F., 1999. Diffusion und Chemische Reaktionen in der Gas/Feststoff–Katalyse. Springer, Berlin. Manwart, C., Toquato, S., Hilfer, R., 2000. Stochastic reconstruction of sandstones. Physical Review E 62 (1), 893–899. Mason, G., Morrow, N.R., 1984. Meniscus curvatures in capillaries of uniform cross-section. Journal of the Chemical Society—Faraday Transactions 1 80, 2375–2393. Perry, R.H., Green, D.W., Maloney, J.O., 1984. Perry’s Chemical Engineers’ Handbook. sixth ed. McGraw-Hill, New York. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in C: The Art of Scientific Computing. second ed. Cambridge University Press, Cambridge. Sahimi, M., 1995. Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches. VCH, Weinheim, Germany. Talukdar, M.S., Torsaeter, O., Ioannidis, M.A., Howard, J.J., 2002. Stochastic reconstruction of chalk from 2D images. Transport in Porous Media 48, 101–123. Torquato, S., 2002. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Tsakiroglou, C.D., Payatakes, A.C., 2000. Characterization of the pore structure reservoir rocks with the aid of serial sectioning analysis, mercury porosimetry and network simulation. Advances in Water Resources 23, 773–789.