Calculation of the Permeability of Porous Media from the Navier-Stokes Equation M. I. S. Azzam and F. A. L. Dullien' Department of Chemical Engineering, University of Waterloo, Waterloo, Ontario, Canada
A set of cubic networks consisting of tubes with periodic step changes in diameter is used to model the flow in porous media. The structural information used in the model is the bivariate pore volume distribution function. The model is used to calculate the permeability of tightly consolidated sandstones. The converging-diverging character of the flow channels is taken into consideration. Pressure drop in excess of that predicted by Hagen-Poiseuille type equations is estimated by solving the complete Navier-Stokes equation numerically. An unsteadystate gas permeametry technique is used to measure the permeability of sandstones. The theoretically calculated and experimentally measured permeabilities agree to within about f20 % . The present calculations yield lower predicted permeabilities by about 159'0, on the average, than the values obtained from the Hagen-Poiseuille equation; the greatest single difference was 70 % .
Introduction The study of flow through porous media has a wide spectrum of applications in industry and a strong relevance to many natural processes. The prediction of various properties of porous materials has been, and still is, a basic task of this research area. For a long time most of the models of porous media have been variations on the simple model consisting of a bundle of uniform capillary tubes. Most of the shortcomings of these models have been discussed critically by Scheidegger (1960) and by Dullien (1975a),among others. I t has been realized, e.g. Scheidegger (1953), that using tubes with variations in the cross-sectional area along their length (so-called serial models) will account for some effects that are ignored in the models using tubes of uniform cross section along their entire length (so-called parallel models). The idea of using periodically constricted tubes, in particular, was apparently introduced by Petersen (1958) and by Houpeurt (1959). Payatakes et a1 (1973) presented a model of unconsolidated porous media consisting of periodically constricted tubes and a solution of the complete Navier-Stokes equation for the problem of flow through such tubes. Their structural analysis, and structural model, seems to be applicable only to granular media of nearly monosized grains. This limitation is also inherent in most of the cell models, for example, Marmur and Rubin (1972);for an extensive review see Happel and Brenner (1965). Very recently several contributions in the field of pore structure analysis have been made by Dullien and his coworkers. An important result of that research is the introduction of bivariate pore size distribution function, determined experimentally in a few sandstones by Dullien and Dhawan (1975), and of the network model advanced by Dullien (1975a,b) which is based on the bivariate pore size distribution concept. This model consists of a set of independent networks each of which is built of tubes with variations in the cross-sectional area along their length. When calculating the pressure drop across these tubes Dullien (1975a) assumed Poiseuille flow in each tube segment. In the present work this model is examined and the extra pressure loss arising from the convergent-divergent character of flow is estimated by solving the complete Navier-Stokes equation, retaining the inertia terms, for the model. The Bivariate Pore Size Distribution Dullien and Dhawan (1975) have recently presented a unique approach to pore structure characterization by means
of a bivariate pore volume distribution function f(D,, D ) dD,dD, giving the fraction of pore volume which is accessible to the non-wetting phase through controlling constrictions (entry pores) in the diameter range De - D e dD, and that D dD. Because the exhave diameters in the range D perimental technique is very tedious and time consuming, it was used, by Dullien and Dhawan (1975) only for a few porous samples. Dullien (1975a) used a computational scheme to estimate the bivariate pore size distributions of a large number of sandstones, based on the knowledge of the mercury porosimetry and photomicrographic pore size distributions; see, for example, Dullien and Mehta (1971). An illustration of such pore size distribution curves is given in Figure 1. Ilullien and Dhawan (1975) and Dullien (1975a) have presented the bivariate distributions in the form of histograms as shown in Figure 2 . In the matrix shown, the index i denotes an entry pore and the index J denotes any pore; thus, the diameter D, indicates a pore entry diameter and D, is the diameter of any pore. The percent pore volume, in the sample, associated with a diameter D., controlled by a neck with diameter D , is denoted by V,,. The entries in the matrix corresponding to i = J give percentage pore volumes occupied by the pores of entry.
-
+
+
Dullien's Network Permeability Model In order to appreciate the process of building upon the analysis presented by Dullien (1975a), a very brief account of Dullien's network model is given here. The model assumes a porous medium to consist of a set of three-dimensional networks. Each class (i) of flow capillaries, characterized by a pore entry with diameter D,, forms an independent network in the model. A network (i) is assumed to he built up of identical capillary elements. Each capillary element, in turn, consists of several capillary segments with different diameters D, controlled by a neck with diameter Di, Figure 3. I t is assumed, moreover, that each element has the same length, the same size distribution and, thus, the same fluid conductivity as the rest of the elements in a particular network. The elements in different networks, on the other hand, have different size distributions, lengths, and hydraulic conductivities. Cubic networks having an arbitrary orientation with respect to the macroscopic flow direction have been assumed by I>ullien (1975a), who showed that a cubic network has isotropic flow properties; in other words, it has the same permeability in any arbitrary direction as the permeability measured in any one of the principal directions. Ind. Eng. Chem., Fundam., Vol. 15, No. 4 , 1976
281
SAMPLE
W
-
u
z
I20
C L E A R CREEK 118 I F W P I
MERCURY POROSIMETRY CURVE PEAK 8 microns PHOTOMICROGRAPHY ( S P H E R E MODEL)
CURVE I
b--
’
CAPILLARY ELEMENT OF NETWORK
, ’
Figure 3. Schematic presentation of capillary element of network (Dullien, 1975).
a 0
D or D e , prn
Figure 1. Comparison of pore volume distributions for clear creek sandstone.
+
1,,,.
,,,,,,, ,,,,,,,, /j
Figure 4. Equivalent capillary element E.
changed also in the equivalent element; however, all the different diameters Dl>, in the original network element were replaced with the equivalent diameter D,+l, defined by the following moment equation Figure 2. Bivariate pore volume distribution.
Through a detailed derivation, accounting for the geometry of the model and assuming the validity of Hagen-Poiseuille equation in each capillary segment, Dullien (1975a) arrived a t the following expression for the permeability of the model
where V,+1 is the volume of the equivalent segment of diameter D,+1, Le., the volume of the bulge in Figure 4. We have also the condition on V i + ,
VI+, =
x v,,
(3)
I>,
because the volume of the network element obviously has to remain unchanged. From eq 1the permeability of network i of Dullien’s model ( K , ) l )is evidently All the quantities appearing on the right-hand side of eq 1are obtainable from the bivariate pore volume distributions histograms, and e is the porosity of the medium.
t
(Ki)n = 96
Further Work on the Network Permeability Model As has been mentioned earlier, Dullien (1975a), when calculating the pressure drop across various capillary segmens of the model, based his formulation on the validity of the Hagen-Poiseuille equation. In other words, extra pressure loss due to the convergent-divergent nature of flow through network elements has been ignored. The significance of extra momentum losses due to a converging-diverging flow path, however, has been demonstrated by Batra et al. (1969),Payatakes et al. (19731, Dullien and Azzam (1973), and Azzam and Dullien (1975). It was felt necessary, therefore, to estimate the magnitude of this effect in the case of the network model by solving the complete Navier-Stokes equation for a similar geometry as that of the capillary elements of the model. A numerical technique suitable for solving the NavierStokes equation for such geometries has been developed by Azzam (1975) (also, Azzam and Dullien, 1975). Examination of the geometry of a network element (Figure 3), however, reveals that the great number of grid points required to cover this field of flow would result in a prohibitively long computer time. Therefore, it was decided to replace the network element shown in Figure 3 with the equivalent element in Figure 4 which has the same hydraulic conductivity or permeability as the original element. The crucial controlling segment or neck diameter D, in the original element was retained un282
Ind. Eng. Chem., Fundam., Vol. 15,No. 4, 1976
‘‘ (
1
[ V~i/DI’+ [ [ V,,/DL6+
VLJ/Dj2)]2 (4)
V,,/D,6]]
whence, by substituting in the numerator the left-hand side of eq 2, the expression for the permeability of the equivalent network i is obtained
where V I V I ,was introduced for the sake of a consistent notation for the equivalent network. In the numerators of eq 4 and 5 the quantities in the braces are exactly equal. In the denominators of eq 4 and 5, (V,,/D,6 , is the dominant term so that the approximation used in the term VL+1/DL+l6 is good enough. Therefore, eq 5 represents a very good approximation of eq 4 as it has been also verified by numerical calculations. Having established the volumes and diameters of the two capillary segments, there remains to calculate their respective lengths. The total length of all “bulges” of diameter D,+l associated with the volume V I +1 is given by
where is the (average) length of a bulge and ni+l is the total number of the bulges in the equivalent i network element. Similarly, the total length of all necks, associated with
Table I. Model Geometry, Calculated Excess Momentum Loss Factors, and Permeabilities of Two Sandstones
4
Dl
Sandstone
3 2 1
25 21 21
0
1,+1
0
tl
Pm
i
Bartlesville 3 31 e = 0.27 2 24 KeXp= 1307.5md 1 15 Berea BE1 e = 0.22 K e x p= 400md
Dl+1
67 53 38
24 54 62 18 9 42.5 10 13.5 30
X
62 1.14 42.5 1.89 30 1.30
/
/
I
206 166 24
(7)
where n,= n,+1 because in the equivalent i network element the number of “bulges” is equal t o the number of “necks”. T o obtain the average length of a “bulge”, 1 , + 1 and that of a “neck”, I,, the major assumption is made t h a t the length of a “bulge” is equal t o its diameter, based on the picture that voids in packed beds of more-or-less isometric particles also tend t o be isometric, on the average. Examination of the photomicrographs of many sandstone samples, Batra (1973) is also in support of this assumption, particularly in the case of large pores. Thus = Dl+i
=0.5
I
‘
67 ‘1.47 848 53 1.43 252 38 1.25 64
volume V, is
L+i
=043
f I / X
K , , md
0.1
n,l, = ~ V , / T D , “
D,/X
(8)
Then, by applying eq 6 and 7
or, using eq 8 (10)
Equations 8 and 10 are the relations used for estimating the length of the “bulge” and the “neck”, respectively.
The Excess Momentum Loss Factor Excess momentum losses due to convergent-divergent flows through tubes described by the present model are ignored where applying the Hagen-Poiseuille equation. T o account for this effect an excess momentum loss factor ( has been introduced by Azzam (1975) and Azzam and Dullien (1975), defined as the ratio of the pressure drop obtained by solving the complete Navier-Stokes equation numerically, 2, to that calculated using the Hagen-Poiseuille equation, APH p
AP (11) IPH P For low Reynolds numbers (Re 5 1) a linear relationship exists between the pressure gradient and the flow rate. This can be written as
E=-
(12) Using eq 12 both for the case when the pressure drop is estimated by the Hagen-Poiseuille equation and when it is calculated by solving the Navier-Stokes equation, along with the expression of the permeability in the H.P. approximation, given by eq 5, the following expression is obtained for the permeability which includes‘the additional pressure drop due to convergent-divergent flow.
1.0
10 Re
100
1000
Figure 5 . Deviation from linearity of the pressure drop vs. Reynolds number relationship (A = I I t 12 = length of equivalent capillary segment).
Here (, is the excess momentum loss factor of the network i , as given by the definition in eq 11 (see Table I). The geometry of each network element is completely defined by the dimensionless ratios Dl+l/DLand l,/D,, as, by assumption, in every case l,+l/Df+lwas equal t o one. T h e Navier-Stokes equation was solved numerically as described in detail by Azzam and Dullien (1975), for each different network element. From the solution the dimensionless pressure drop A p ( p = Pplu2)was obtained and this was used to calculate the excess momentum loss factor, (, by the equation
The numerical solution was performed a t low values of Reynolds number (Re I1;Re = [u(l, l , + ~ ) p ] / pwhere , the relation between l p and Re is linear on a log-log scale and can be represented by the equation
+
A p = C/Re
(Re 5 1.0)
(15)
In such a case, the value of the excess momentum loss factor is constant, independent of the Reynolds number. Accordingly, excess momentum loss is uniquely determined by the geometry of tubes. At higher values of Reynolds number gradual deviation from eq 15 starts to occur as exemplified in Figure 5 for a typical case. This deviation is attributed to inertial effects which become increasingly important a t higher values of Reynolds number. Under these conditions additional excdss momentum losses take place. For length-to-diameter ratios of the neck, l,/D! > 5 , the excess momentum loss factor was practically unity in every case. The smaller ll/Dl and the bigger Df+l/DL,the greater was E . The highest value of ( encountered in the calculation was 2.33, obtained for the case where ll/Dl = 0.45 and D,+1/DI = 3.4. The geometrical parameters, calculated excess momentum loss factors and permeabilities of the various network elements of the two sandstones for which the bivariate pore size distributions were determined experimentally by Dullien and Dhawan (19751, are listed in Table I. It can be seen that the l,/D, ratios are realistic values for each network element, indicating the basic soundness of the approach used in this work. There is additional, and completely independent, evidence for this from a comparison of pore size distributions and displacement experiments with miscible liquids (Dullien and Azzam, 1976).
Flow Measurements In order to test the predictions of the model presented in this work, a n experimental study has been undertaken to Ind. Eng. Chem., Fundam., Vol. 15,
No. 4, 1976
283
COMPRESSED AIR
R
BEREA 108
caso4
-
TO VAC
PUMP
PRESS RELEASE
HASSLER CELL
B
RESERVOIR
Figure 6. Schematic presentation of the experimental setup.
determine the permeability of a wide range of sandstones. For reasons of reliability and convenience, an unsteady-state technique has been adopted, and gas was used as the flowing fluid. This unsteady gas permeametry technique parallels in several aspects the unsteady-state methods used by Barrer and Grove (1950) and by Jenkins and Roberts (1962). The apparatus designed and constructed for the purpose of the present experiments is versatile and suitable for a wide range of porous samples. The Apparatus. A schematic representation of the apparatus used in this work is shown in Figure 6. The apparatus consists of the following major components: (1) A pressure vessel packed with calcium sulfate to ensure moisture-free ingoing stream of air. (2) A pressure gauge for measuring the inlet pressure, P1.Two pressure gauges with different pressure ranges were used depending on the flow conditions. (3) A Hassler-type cell for holding the sample. The Hassler sleeve core holder consisted of a brass cylinder lined with a rubber tube and fitted with flanges a t both ends. A hole in the middle of the brass cylinder provided a connection to a nitrogen pressure-cylinder. Under operating conditions, nitrogen gas pressure of about 400 psig was applied to the outside of the rubber sleeve. This ensured a tight fit of the rubber tube to the surface of the core thus preventing peripheral flow of air. (4) A reservoir vessel which provided extra volume that could be brought into the circuit in case of high permeability samples. This extra volume was found necessary sometimes to allow reasonable rate of change of downstream pressure, P?. Also, knowing the volume of this vessel aided in determining the volume of the associated pipework on the downstream side of the sample. ( 5 ) A pressure transducer the signal from which was displayed on a D.V.M. and could be recorded when required. (6) A vacuum pump. Experimental Procedure. Various constant pressures PI were applied upstream of the sample with the aid of a compressed air line. At the same time, the space downstream of the sample was evacuated. When steady conditions were reached, as detected by the constancy of the D.V.M.’s signal, the vacuum pump was isolated whereupon the pressure P2 started to rise. The rate of pressure rise was recorded. It has been found necessary to keep Pz > 1 wm) that would not normally allow slip to occur under the flow conditions of the present experiments, it was felt necessary to test this effect.
lieved that it is equally applicable to unconsolidated materials such as packed beds of particles.
I B E L T SERIES 120 2 WHETSTONE 141 3 BANDERA 136
>
k
m i
SOUIRREL50
r
NOXIE 1 2 8
Nomenclature
9 TORPEDO 43
/
a 100 W I
4
I1 B E R E A I08 12 B A R T L E S W L L E
W Dc
n 0 W l-
a
O K
3
Y
a
Kg
1 l .0o / ,
0
0.1 0
1
10 10 100 1000 MEASURED PERMEABILITY
Figure 9. Comparison of
measured a n d calculated permeabilities.
Thus, the experimental data were fitted using both eq 17 and 20; statistical significance tests showed that slippage may be ignored in all the cases studied. Only in few of the cases investigated did the relation between the pressure drop and flow rate display a nonlinear behavior. Accordingly, the experimental observations were fitted using the Forchheimer equation (Figure 8 ) , which in case of gases is written in the form
P,AP p -- - -k ppu2 P ~ u ~ KL It is noted that the Ergun equation (e.g., Bird et al., 1960) is a special form of the Forchheimer equation. The analysis of this work, expressed ultimately by eq 13, is tested by comparing the predicted permeability K with the measured permeability Kexpin case of 12 widely different sandstones (Figure 9). The permeability value K Dcalculated by Dullien (1975a), based on the Hagen-Poiseuille approximation, is also compared with value K in the figure. In several cases considerable improvement over the predictions of Dullien (1975a) has been achieved. For example, in case of the “Belt Series 120” sample improvement of about 70% was obtained, over 40% improvement was obtained in case of “Berea BEl”, and over 30% in case of “Whetstone i41”. In seven cases the present calculations showed improvement over the original predictions. In case of “Bandera 136” the calculations resulted in a permeability value K that is equal to KD.In four cases the uncorrected predictions agreed better with the experimental values. However, the difference between K and K D in these cases did not exceed 14%. On a collective basis, the average deviation of the K Dvalues from the 45’ line is over +lo%,while the average deviation of the K values is less than -5%. The average scatter of the Kll values around the mean is about f31%, while the average scatter of the K values around the mean is about f 2 1 % . It is apparent from the foregoing discussion of the results that taking the effect of the converging-diverging nature of flow into consideration resulted, on the average, in a better prediction of the permeabilities of sandstones. In some cases, appreciable improvement in the predictions was obtained. In some other cases, however, the Hagen-Poiseuille approximation appeared quite reasonable. Although the analysis of this work was applied to consolidated porous media, it is be-
A = cross-sectional area of porous samples C = constant D = diameter of capillary segment, or a pore K = permeability I, = length of porous sample 1 = length of capillary segment n = number of capillary segments P = pressure p = dimensionless pressure Q = volumetric flow rate Re = Reynoldsnumber S = slip term t = time u = velocity in the narrow segment (neck) of the equivalent capillary element V = percentage pore volume in the bivariate distribution V , = volume on the vacuum side of porous sample u = superficial velocity Greek Symbols = inertial resistence coefficient of porous media t = porosity p = viscosity = excess momentum loss factor p = density Subscripts 1 = inlet of porous sample 2 = outlet of porous sample D = pertains to Dullien’s calculations based on the HagenPoiseuille equation E = refers to the element E e = entry pore i = specifies a particular network or neck size ij = specifies both the size of neck and bulge j = refers to a bulge m = mea.n value L i t e r a t u r e Cited Azzam, M. I. S., Ph.D. Tbesis, Uriiversity of Waterloo, Waterloo, Ont., Canada, 1975. Azzam, M. I. S., Dullien, F. A. L., paper submitted for publication to Chem. Eng. Sci., (197!5). Barrer. 2 . M., Grove, D. M., Trans. Faraday Soc., 47, 826 (1951). Batra, V. I., Fulford, G. D., Dullien, F. A. L., Can. J. Chem. Eng., 48, 622 (1970). Bird, R. B., Stewart, W. E., Lightfoot, E. N., “Transport Phenomena”, Wiley, New York, N.Y., 1960. Carman, P. C., “Flow of Gases through Porous Media”, Butterworths, London, 1950. Dullien, F. A. L., AlChE J., 21, 299 (1975a). Dullien, F. A. L., AlChEJ., 21, 820 (1975b). Dullien, F. A. L., Azzam, M. I. S., AIChE J., 19, 222 (1973). Dullien, F. A. L., Azzam, M. I. S., lnd. Eng. Chem.. Fundarn., 15, 147 (1976). Dullien, F. A. L., Dhawan. G. K.. J. Colloid lnterface Soc., 52, 129 (1975). Dullien, F. A. L., Mehta. P. N., Powder Techno/., 5, 179 (1971/72). Happel. J., Brenner, H.. “Low Reynolds Number Hydrodynamics”, Prentice-Hall. Englewoocl Cliffs, N.J., 1965. Houpeurt. A., Rev. lnst. Fr. Pet., 14, 1468 (1959). Jenkins. T. R., Roberts, Fifth Carbon Conference, 335 (1962). Marmur, A., Rubin, E., lnd. Eng. Chem., Fundam., 11, 497 (1972). Payatakes, A . C., Tien, C.. Turian, R. M., AlChE J., 19 (Parts I & II), 58 (1973). Petersen. E. E., AlChE J., 4, 343 (1958). Scheidegger, A. E.,Producers Monthly, 17, 17 (1953). Scheidegger, A. E., “The Physics of Flow Through Porous Media,” University of Toronto Press, Toronto, 1960.
Receiued for reuieu>October 13, 1975 Accepted M a y 17,1976
Ind. Eng. Chern., Fundarn., Vol. 15, No. 4, 1976
285
