5010
Ind. Eng. Chem. Res. 1997, 36, 5010-5014
Reconstruction of Nonhomogeneous Porous Media Alessandra Adrover*,†,§ and Massimiliano Giona‡,§ Dipartimento di Ingegneria Chimica, Universita` di Roma “La Sapienza”, Via Eudossiana 18, 00184 Roma, Italy, Dipartimento di Ingegneria Chimica, Universita` di Cagliari, Piazza d’Armi, 09123 Cagliari, Italy, and Centro Interuniversitario sui Sistemi Disordinati e sui Frattali nell’Ingegneria Chimica, Universita´ di Roma, Via Eudossiana 18, 00184 Roma, Italy⊥
An analytical method for the generation of model porous structures (lattice structures) exhibiting a given nonuniform void fraction distribution and a given two-point (pore-pore) correlation function is presented. The method proposed for the generation of correlated and nonhomogeneous porous models is based on the application of a nonlinear filter to a stochastic correlated process. This process originates from the linear superposition of a suitable set of stochastic processes possessing Gaussian decay in the second-order correlation function. 1. Introduction The generation of model porous structures possessing the same porosity as well as the same two-point (porepore) correlation function as real porous materials belongs to the class of inverse problems and is referred to as the reconstruction of porous media. This topic has been extensively studied by Joshi (1974), Adler et al. (1990), Quiblier (1984), and Giona and Adrover (1996) for isotropic and homogeneous porous materials. In particular, Giona and Adrover (1996) (henceforth referred to as the GA approach) propose a closed-form solution for the reconstruction of isotropic and homogeneous porous media by reducing the problem to a linear functional equation which can be solved by means of a Laplace inversion. The reduction of the reconstruction problem to the form of a linear functional equation is made possible by the choice of a suitable set of stochastic basis processes possessing Gaussian decay in the second-order correlation function. The purpose of this paper is to extend the GA approach to solve the problem posed by the reconstruction of isotropic and nonhomogeneous porous materials, i.e., porous structures exhibiting a nonuniform void fraction distribution. This problem has many practical applications, especially in the field of fluid-solid noncatalytic reactions in which solid particles with a nonuniform distribution of solid reactant are generally encountered (Dudukovic, 1984; Sohn and Xia, 1986, 1987). Indeed, the reconstruction of porous media is a typical inverse problem aimed at obtaining a way to generate lattice structures with the same statistical geometrical features as the original sample. Many of the ideas developed in order to solve the reconstruction problem can, however, also be applied to percolation models with correlations. As noted by Sahimi (1994), correlated percolation models are very useful in application to porous media and transport because they mimic the structure of real materials better than the usual (uncorrelated) percolation schemes. Many authors have undertaken a separate analysis of the effects of spatial correlation properties and pore* To whom correspondence should be addressed. † Universita ` di Roma “La Sapienza”. ‡ Universita ` di Cagliari. § Centro Interruniversitario sui Sistemi Disordinati e sui Frattali nell’Ingegneria Chimica. ⊥ Phone: +39-6-44585892. Fax: +39-6-44585339. E-mail:
[email protected]. S0888-5885(97)00165-6 CCC: $14.00
space nonuniformities in convection/diffusion transport (Adler et al., 1990; Adrover and Giona, 1996; Sahimi, 1994) and in reaction/diffusion problems (Sohn and Xia, 1986; Dudukovic, 1984), but no studies focus on the interplay between these two features. The possibility of generating model structures with assigned nonuniform void-fraction distribution (or equivalently solid-reactant distribution) and spatial correlation properties enables us to analyze these coupled textural properties (Adrover and Giona, 1997), which are encountered in real porous media. 2. Problem Let us consider an isotropic, radially symmetric, nonhomogeneous porous particle. Under the assumption of isotropy, a generic axial cross section of the particle is representative of the entire structure. Let us analyze a two-dimensional experimental section of the particle by means of its digital image (blackand-white picture). This leads to identification of the pore space P in terms of its characteristic function χP(r):
χP (r) )
{
1 r∈P 0 elsewhere
}
(1)
Knowledge of the characteristic function χP (r) makes it possible to evaluate the void fraction (porosity) (r)
(r) )
∫02πχP (r,θ)r dθ
1 2π
(2)
which depends exclusively on the radial coordinate r and the normalized pore-pore correlation function C(d) 2χ (x)
C(d) 2χ (x) )
〈(χP (r) - (r))(χP (r + x) - (|r + x|))〉 〈(r)〉 - 〈2(r)〉
(3)
where the symbol 〈‚〉 stands for spatial average
〈f(r,θ)〉 )
∫02πdθ∫0Rdr f(r,θ)p(d)(r)
1 2π
p(d)(r) ) d
rd-1 Rd (4)
R being the particle radius and p(d)(r) the d-dimensional radial density function, (d ) 2 for two-dimensional structures, d ) 3 for the three-dimensional ones). © 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997 5011
The aim of the reconstruction is to generate two- and/ or three-dimensional lattice models possessing the same (r) and the same C(d) 2χ (x) as the original porous structure. The solution is presented in the form of a linear filter acting on Gaussian processes by means of a superposition of elementary correlated processes with prescribed correlation properties, as developed by Giona and Adrover (1996). Transformation from the resulting Gaussian correlated process y(r) to the binary (0/1) process representing the reconstructed medium can be achieved by means of a pointwise nonlinear filter depending on the porosity function (r). Before developing the method, let us examine in some detail the differences between the isotropic, homogeneous, and isotropically nonhomogeneous porous media and the relevant implications. By definition, a porous medium is homogeneous if it is statistically invariant under translation and isotropic analogously if it is invariant under rotation Tθ, where Tθ indicates the rotation operator of angles θ (Adler, 1992). In two dimensions, θ reduces solely to the angle θ of a cylindrical coordinate system and in three dimensions to the angle θ ) (θ,φ) of a spherical coordinate system. For homogeneous structures, translational invariance ensures that spatial (volume) and ensemble averages coincide. The equivalence between volume and ensemble averages does not hold for nonhomogeneous structures. For this reason, the definition of the porepore correlation function given by eq 3 implies a spatial average. Let us further analyze the properties of nonhomogeneous isotropic structures, which are the main objects of investigation in this article. Examples of such a structure occur for radially symmetric pellets, for which (r) depends on the radial coordinate. Since we are mainly interested in correlation properties up to the second order, let us consider the statistical (ensemble) description of the isotropic porous structures up to this order by introducing (a) the one-point probability function p1(x) such that χP (x) ) 1 at x and (b) the two-point probability function p2(x1,x2) such that the random field χP (x) attains the value of 1 at both x1 and x2. Under the condition of isotropy (radial symmetry), it follows that
p1(Tθ(x)) ) p1(x) for all θ and x
(5)
p2(Tθ(x1),Tθ(x2)) ) p2(x1,x2) for all θ and x1, x2 (6) Equation 5 implies that p1(r) depends exclusively on the radial coordinate (r); i.e., p1(r) ) (r). C(d) 2χ (x), as given by eq 3, can be expressed as the volume average of p2(x1,x2). To simplify the notation, let us consider the two-point term 〈χP (r)χP (r+x)〉. By definition,
〈χP (r)χP (r+x)〉 )
1 V(d)
∫V
(d)
p2(r,r+x) dV(d)(r) ) G(x) (7)
where V(d) is a d-dimensional reference volume. By
enforcing isotropy, i.e., eq 6, it follows that
1
∫
p (r,r+x) (d) V(d) 2
V
dV(d)(r) )
Tθ(x)) dV(d)(r) )
1 V
1 V(d)
∫V
(d)
p2(Tθ(r),Tθ(r) +
∫
p (r′,r′+Tθ(x)) (d) V(d) 2
dV(d)(r′) ) G(Tθ(x)) (8)
In deriving eq 8, we have made use of a coordinate change r′ ) Tθ(r) of the integration (dummy) variable r. From eqs 7 and 8, it follows that G(Tθ(x)) ) G(x) for all θ and x, and the function G in eq 7 therefore depends exclusively on x ) |x|. A similar calculation can be extended to all the other terms, linear in χP, entering into the definition of the normalized two-point correlation function, eq 3. This proves that C(d) 2χ (x) depends exclusively on x ) |x| for isotropic (radially symmetric) structures. To conclude this analysis, it should be observed that while a generic two-dimensional thin section is representative of the entire structure for isotropic and homogeneous porous media, in the case of isotropic and nonhomogeneous structures, the only representative two-dimensional sections are those passing through the origin, i.e., through the center of symmetry of the structure. 3. Method The method develops in the following steps: Step I involves definition of a set of Gaussian correlated basis processes {y(r,λ)} by convoluting a system of Gaussian uncorrelated processes ξλ(r) (with zero mean and unit variance) with a Gaussian kernel,
y(r,λ) )
∫E a(u,λ)ξλ(u+r) du ) d/4 (4λπ) ∫E e-2λu ξλ(u+r) du d
2
d
(9)
where Ed is the Euclidean d-dimensional space, Ed ) {r|-∞ < ri < ∞ (i ) 1, ..., d)}. Each basis process y(r,λ) remains Gaussian with zero mean and unit variance (Doob, 1953) but exhibits a Gaussian decay of the correlation function 2
C2y(x,λ) ) e-λx
(10)
(The basis processes {y(r,λ)} for different λ are uncorrelated with each other; i.e., 〈y(r,λ1)y(r′,λ2)〉 ) 0 for λ1 * λ2.) Step II involves expansion of a generic correlated Gaussian process y(r) as a linear superposition of the basis processes {y(r,λ)}:
y(r) )
∫0∞p(λ)y(r,λ) dλ
(11)
The definition, eq 11, of the y(r) process implies
C2y(x) )
∫0∞π(λ)C2y(x,λ) dλ ) ∫0∞π(λ)e-λx
2
dλ
(12)
where π(λ) ) p2(λ) is the weight function. Step III involves transformation from the y process to the binary process representing the characteristic function of the reconstructed porous structure χRP (r) by means of a nonlinear filter G, depending on the Gauss-
5012 Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997
ian distribution function Fy of y(r) and on the (r) of the original porous medium:
χRP(r) ) G(y(r),(r)) )
{
1 Fy(y(r)) < (r) 0 Fy(y(r)) > (r)
(13)
The definition, eq 13, of G ensures statistically that the reconstructed porous structure admits the correlation function of the reconstructed lattice (r). The only further condition to be imposed is that the porosity (d) function C2χR (x) coincides with C(d) 2χ (x). (d) C2χR(x) is related to C2y(x) through the integral relation (d) C2χR (x) )
∫0Rdr1∫0Rdr2∫-∞∞dy1∫-∞∞[H(y1,y2,r1,r2)] × g(y1,y2,x)p(d)(r2/r1;x)p(d)(r1) dy2 (14)
with
H(y1,y2,r1,r2) )
(G(y1,(r1)) - (r1))(G(y2,(r2)) - (r2)) 〈〉 - 〈2〉
(15)
where g(y1,y2,x) is a bivariate Gaussian density function
g(y1,y2,x) )
1 × 2π(1 - C2y2(x))1/2
[
exp -
]
y12 + y22 - 2C2y(x)y1y2 2(1 - C2y2(x))
(16)
and F(d)(r2/r1;x) ) ∫or2p(d)(r/r1;x) dr is the probability distribution function of finding a point possessing the radial coordinate r2 at a relative distance x from a point at a radial coordinate r1 in a d-dimensional radial symmetric structure. By defining
γ(r,r1,x) )
r2 - (r12 + x2) 2xr1
(17)
the probability distribution function F(d)(r2/r1;x) attains the form
{
F(d)(r2/r1;x) ) 1 1 - arccos (γ(r2,r1,x)) π if (x + r1) e R and |γ(r2,r1,x)| < 1 π - arccos(γ(r2,r1,x))
(18) π - arccos(γ(R,r1,x)) if (x + r1) > R and |γ(r2,r1,x)| < 1 1 if γ(r2,r1,x) g 1 0 if γ(r2,r1,x) e -1
(It should be observed that F(2)(r2/r1;x) ) F(3)(r2/r1;x).) The integral equation (14), expressing the functional (d) relationship between C2χR (x) and C2y(x), can be integrated numerically for the given porosity function. In this way, for each x, starting from a set of values {C2yi} ∈ [-1,1], it is possible to evaluate the corresponding set (d) of values {C2χRi }, which can be used as a calibration curve. (For homogeneous porous media, the integral (d) equation relating C2χR (x) and C2y(x) does not depend
(2) Figure 1. Calibration curves C2χR (x) vs C2y(x) > 0 for a linear porosity function, eq 15 (0 ) 0.2, R ) 0.8, R ) 150 lattice units), and for three different values of x/R.
explicitly on x, and consequently, there is only one calibration curve for each value of the average (and uniform) porosity 〈〉 ) .) (d) Figure 1 shows three calibration curves C2χR (x) vs C2y(x) > 0 for a linear porosity function
(r) ) 0 + (R - 0)r/R
(19)
and for three different values of x/R. The integration of eqs 14-18 was performed by means of a Monte Carlo algorithm. Step IV involves the evaluation of the weight function (d) π(λ). By imposing C2χR (x) ) C(d) 2χ (x) and by making use of the relations eqs 14-18 (or equivalently of the calibration curves), it is possible to obtain the correlation function C2y(x). From the knowledge of C2y(x), the weight function π(λ) can then be evaluated from eq 12 either numerically or, where possible, analytically, following the same procedure as developed in the GA approach in connection with homogeneous porous media. The key point in the GA approach is that, based on the condition that C2y(x) g 0, putting z ) x2 and C ˜ 2y(z) ) C2y(x), it follows that
C ˜ 2y(z) )
∫0∞π(λ)e-λz dλ
(20)
and therefore, π(λ) is the inverse Laplace transform of the analytic continuation on the complex plane C ˜ (p) 2y (z) of the correlation function of the y(x) process C ˜ 2y(z). The analytic continuation C ˜ (p) ˜ 2y(z), valid for all 2y (z) of C complex z (whose restriction to real z coincides with C ˜ 2y(z)), can be achieved by considering rational approximations (e.g., Pade´ approximants) or by means of other methods such as orthogonal polynomial expansion. In all the cases in which closed-form solutions for the inverse Laplace transform of C ˜ (p) 2y (z) cannot be obtained, the weight function π(λ) can be evaluated numerically. Step V involves generation of the reconstructed structure. Once the weight function π(λ) has been evaluated, the generation method consists in the successive application of steps I-III. 4. A Numerical Example The validity of the method proposed can be checked directly and simply by analyzing the correlation properties and the porosity function of a two-dimensional lattice model generated from a single basis process y(r) ) y(r,λ) with a given nonlinear filter G, and ascertaining
Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997 5013
5. Conclusions
Figure 2. Two-dimensional lattice models generated directly from two different basis processes y(r,λ) (for two different values of λ) by applying a nonlinear filter G with a linear porosity function (r), eq 19 (0 ) 0.2, R ) 0.8, R)150 lattice units). (a, left) λa2 ) 0.1; (b, right) λa2 ) 0.01. a is the unit lattice site size.
Figure 3. Comparison of the porosity function (r) of the twodimensional lattice structure of Figure 2a (points) and the theoretical linear behavior, eq 19 (continuous line).
The GA approach to the solution in closed form of the problem of reconstructing porous media, originally developed for homogeneous and isotropic structures, has been extended to nonhomogeneous porous structures. In order to take into account a nonuniform porosity (r), we have introduced a modified definition of the pore-pore correlation function, eq 3, a pointwise nonlinear filter depending on (r), and an integral relation (d) between C2χR and C2y depending on (r), the radial density function p(d)(r1), and the probability distribution function F(d)(r2/r1;x). This makes it possible to reduce the problem of reconstruction to a linear functional equation, eq 12 or 20, which can be solved by means of Laplace inversion. The method developed is of great interest in the lattice simulation of fluid-solid noncatalytic reactions, diffusion processes, and transport in porous particles, in order to study the reaction/diffusion properties of porous solids displaying a nonuniform porosity as well as correlation properties in the pore-network structure. It has been shown by Adrover and Giona (1997) that the effects of correlations in nonuniform reacting pellets modify the overall reaction evolution (conversion-time plots), and correlation properties should be taken into account in all models which attempt to analyze the reaction evolution of solid particles starting directly from structural data on the porous pellets. From an experimental point of view, the information required to apply the method proposed is essentially optical image analysis of thin sections of the structures, which has been already successfully applied in the study of transport and permeability of porous media up to the length scales of about 0.1 µm. Acknowledgment This work was supported in part by a MURST grant for bilateral research exchange. Nomenclature
Figure 4. Comparison of the pore-pore correlation function (2) C2χR (x) of the lattice structure of Figure 2a (points) and the theoretical behavior (continuous line) obtained from the numerical integration of eqs 14-18, with the linear porosity function, eq 19 and with C2y(x) ) C2y(x,λ) ) exp(-λx2) (λa2 ) 0.1). x is expressed in lattice units.
their level of agreement with the theoretical predictions of eq 13 and eqs 14-18. Figure 2 shows two two-dimensional lattice models generated from two different basis processes y(r,λ) (for two different values of λ) by applying a G with a linear porosity function eq 19 (0 ) 0.2, R ) 0.8, R ) 150 lattice units). Figure 3 shows the good level of agreement attained between the porosity distribution function of the twodimensional lattice structure of Figure 2a (points) and the theoretical linear behavior, eq 19 (continuous line). Figure 4 shows the good level of agreement attained between the pore-pore correlation function of the lattice structure of Figure 2a (points) and the theoretical behavior (continuous line) obtained from the numerical integration of eqs 14-18, with the linear porosity function, eq 19, and with C2y(x) ) C2y(x,λ) ) exp(-λx2).
a ) unit lattice site size a(r,λ) ) kernel defining the basis processes; see eq 9 y(r,λ) ) stochastic basis process, see eq 9 (d) ) normalized pore-pore correlation function of the C2χ d-dimensional porous structure C2y ) normalized second-order correlation function of y(r) (d) C2χR ) normalized pore-pore correlation function of the d-dimensional reconstructed porous structure d ) space dimension Fy ) distribution function of y(r) G ) nonlinear filter, see eq 13 r ) radial coordinate r ) position vector x ) |x| [lattice units] y(r) ) continuous correlated stochastic process, see eq 11 y(r,λ) ) stochastic basis processes, see eq 9 Greek Letters (r) ) radial porosity function λ ) Gaussian decay parameter, see eq 10, [1/a2] ξλ(r) ) normalized Gaussian random process π(λ) ) weight function, see eq 12 χP (r) ) characteristic function of the pore space P
5014 Ind. Eng. Chem. Res., Vol. 36, No. 11, 1997 χR P (r) ) characteristic function of the reconstructed porous structure
Literature Cited Adler, P. M. Porous MediasGeometry and Transports; Butterworth-Heinemann: Boston, 1992. Adler, P. M.; Jacquin, C. G; Quiblier, J. A. Flow in Simulated Porous Media. Int. J. Multiphase Flow 1990, 16, 691. Adrover, A.; Giona, M. A Predictive Model for Permeability of Correlated Porous Media. Chem. Eng. J. 1996, 64, 7. Adrover, A.; Giona, M. A versatile lattice simulator for fluid-solid noncatalytic reactions in complex media. Ind. Eng. Chem. Res 1997, 36, 4993. Doob, J. L. Stochastic Processes; John Wiley: New York, 1953. Dudukovic, M. P. Reactions of Particles with Nonuniform Distribution of Solid Reactant. The Shrinking Core Model. Ind. Eng. Chem., Process Des. Dev. 1984, 23, 330. Giona, M.; Adrover, A. Closed-form solution for the reconstruction problem in porous media. AIChE 1996, 42, 1407. Joshi, M. Y. A class of stochastic models for porous media. Ph.D. Dissertation, University of Kansas, Lawrence, 1974.
Quiblier, J. A. A new three-dimensional modeling technique for studying porous media. J. Coll. Interface Sci. 1984, 98, 84. Sahimi, M. Long-Range Correlated Percolation and Flow and Transport in Heterogeneous Porous Media. J. Phys. I Fr. 1994, 4, 1263. Sohn, H. Y.; Xia, Y. N. Effects of Nonuniform Distribution of Solid Reactant on Fluid-Solid Reactions. 1. Initially Nonporous Solids. Ind. Eng. Chem., Process Des. Dev. 1986, 25, 386. Sohn, H. Y.; Xia, Y. N. Effects of Nonuniform Distribution of Solid Reactant on Fluid-Solid Reactions. 2. Porous Solids. Ind. Eng. Chem., Process. Des. Dev. 1987, 26, 246.
Received for review February 24, 1997 Revised manuscript received July 15, 1997 Accepted July 22, 1997X IE970165P
X Abstract published in Advance ACS Abstracts, September 15, 1997.
