Two-fluid mathematical model for compressible flow in fractured porous media

A. L. Khlaifat¹

Department of Chemical Engineering, Faculty of Engineering, Mutah University, Mutah 61710, Jordan
khlaifat@mutah.edu.jo
¹ Director of Prince Faisal Center for Dead Sea, Environmental and Energy Research

Abstract — A three-dimensional isothermal transient compressible two phase flow of liquid and gas in a low permeability fractured porous media model has been developed from the general mass and momentum balances using volume averaging techniques. The emphasis of the paper is on the available analytical and semi-empirical correlations and closure relations. The interfacial mass and momentum transfer terms in the averaged formulation are explained. The assumptions used in the model's formulation are that a condition of capillary equilibrium exists throughout the media, momentum transfer between mobile fluid phases is negligible in the porous media and exists in the fracture, and fluid phase change has been neglected. The pore size distribution of the studied porous media was represented by three mean pore diameters: liquid pore network, gas pore network, and fractures. An average pore diameter and length for each network were determined using the gamma distribution function. The model was validated and solved for a hypothetical porous media and fractured domain.

Keywords — Two-Fluid. Fractured. Porous Media. Capillary Pressure. Transient. Compressible Flow. Isothermal. Mass. Momentum Transfer.

I. INTRODUCTION

The increased sophistication of oil/gas recovery technologies has brought with it increased operating and material costs and therefore a greater demand for sound process designs. Mathematical models of fluids flow in petroleum and gas reservoirs have become key tools by which reservoir engineers develop and implement these designs. Using mathematical models together with various characterizations of the rock-fluid system being modeled, the engineer can test various operating strategies, compare different recovery technologies, and formulate hypotheses in diagnosing the performance of ingoing projects.

It is impossible to specify the microstructure in a realistic porous medium completely. Porous media can be characterized without specifying the porous geometry in all its details. The well-known approach of doing so is by constructing a simplified geometric model, like bundle of capillary tubes, grain models, network models, and percolation model, for each specific porous medium of interest.

Most of the models describing the flow through fractured porous media, such as (Thomas et al., 1983; Evans, 1982; Guzman and Khaled., 1992; Fung, 1993; Warren and Root, 1963; Zekai, 1988; Nitao, 1990; Arbogast, 1993; Celis et al., 1992; Lee and Tan, 1987; Alder and Thovert, 1999; Panfilov and Panfilov, 2000) were based on Zheltov's model (Zheltov et al., 1960), and used his concept of dual-porosity dual-permeability region. They differentiate between two flow regions, one representing the discrete matrix where the other represents the continuous fracture network. Based on that, they considered mainly two types of equations: the equations describing flow in fracture system and equations describing flow in the matrix.

The main objective of the manuscript is the development of a mathematical model which predicts the flow phenomena in low permeability fractured porous media. It is rather impossible to apply conservation laws directly to each pore in porous medium and to find a possible numerical or analytical solution. Instead a semi-empirical approach has been used where some constitutive equations are needed, and a set of equations are developed based on flow through the entire medium.

It is known in reservoir engineering that the productive strata (producing formation or beds) could be made up of rocks described by the following characteristics and properties; porosity, permeability, granulemetric composition, elasticity, resistance to rupture, compression, deformation and saturation. Using these characteristics and properties, the conditions of oil and gas fields' development can be determined.

Pore structure of reservoir rocks is very complex. In the case of fractured porous media, the porosity can be placed into two classes:

1. Primary porosity, which is highly interconnected and can be correlated with permeability, this could be the porosity of the homogeneous rocks. This region always has high resistance (low permeability) to flow.
2. Secondary porosity, formed by a fracture, usually this porosity is the product of geological movements, hydraulic fracturing and chemical processes. Although it does not contain a fraction of fluid reserve as large as that of the first class, it greatly affects the flow. Moreover, it has low resistance to flow compared to the primary porosity region.

For the case of low permeability fractured porous media a network model is developed, the simultaneous existence of two phases of water and gas is considered. The wetting fluid occupies the smaller pores while the non-wetting fluid occupies the larger ones, and when these two fluids exist in the same pore diameters, they move in plug flow fashion. Simultaneous two-phase flow exists in the fractures and similar flow patterns as in the case of two-phase flow in a pipe with momentum transfer across phase boundaries are assumed. Further, there is continuous movement of the rigid porous media due to changes in stress distribution, particularly during production.

The pore size distribution in the medium is represented by three mean pore diameters: gas pore network, liquid pore network, and fractures. This representation along with the movement of the rock matrix assumption result in a model with five different phases: 1) gas in the porous matrix (gas network), 2) liquid in the porous matrix (liquid network), 3) gas in the fracture, 4) liquid in the fracture, and 5) rock matrix. An average pore diameter and length for each network were calculated using the gamma distribution function.

II. GOVERNING EQUATIONS

In the formulation of the mathematical model the following assumptions are made:

1) The reservoir matrix is homogeneous and isotropic. Furthermore, there is inaccessible pore volume

2) The reservoir rock and the working fluids are slightly compressible

3) Because of the absence of both thermal action and chemical reaction, the process is considered to be isothermal

4) We assume a constant mass generation of phase π, where π can be either solid (β) or fluid (α) phase

5) Gas and liquid are immiscible; furthermore we assume that there is no slip condition at the interface between the phases

The continuity equation for the Cartesian system of coordinate may be written, in tensor notation, as:

(1)

where: is the rate of accumulation of mass per unit volume at specific point (say p); is the net flow rate of mass out of p per unit volume; is mass generation of phase π; and the momentum equation as:

(2)

where is the rate of π momentum increase at the fixed point (say p), is the net rate of π momentum carried into p by the fluid or solid flow , is the net π pressure force at p, is the net π stress force, where is the viscous stress tensor for phase π, is π body force at p, where if gravity is the only body force we consider and is momentum generation at point p.

III. AVERAGING PROCEDURE

The governing equations were averaged using the volume averaging technique. This technique has been used widely for flow in porous media. It was demonstrated by William and O'Neill (1976), that averaging technique can be applied to various transport processes in porous media. William and O'Neill (1976), William and Sampath (1982) and Slattery (1972), were among the first to use this technique to develop the governing equations for flow in porous media. Other researcher such as (Dullien, 1975; Le Gallo, 1990; Pakdel, 1994; Jacob and Yehuda, 1992, and Al-Khlaifat, 1998) has applied this technique to multiphase flow system.

Volume averaging technique results in achieving manuscript initial objective which is to associate with every point in a porous medium a local volume average of the differential equations of continuity (1) and momentum (2). When we say every point, we include the solid phase as well as the fluid phase and the solid-fluid phase interface.

The phase average of the continuity (1) and the momentum (2) equations can be written as:

(3)

and

(4)

Applying the averaging rules (William and O'Neill, 1976; William and Sampath, 1982), to each term of Eq. (3) and (4). When assuming the no slip condition at the interface and a constant mass generation for α phase, the final volume averaged continuity equation will take the following form:

(5)

After accounting for the no-slip condition at the interface, the final form of averaged momentum equation will be:

(6)

IV. CONSTITUTIVE EQUATIONS

In general a fluid flow problem is governed by several equations. First there are the basic equation of continuity, three momentum equations and an energy relation. Second, there are the constitutive equations, which are not basic, but they do apply to group of substances. Various transport coefficients are introduced in the constitutive relations. They are quasi-thermodynamic properties that depend on the composition of the fluid and its thermodynamic state. Third, the thermodynamics of fluid must be specified, this may be done through the fundamental equation s(ρ,e) for the substance, or, more commonly, through two equation of state P(ρ,T) and e(ρ,T). All of these equations are required to give a well-posed problem for a general flow situation.

For the considered isothermal case, in order to use Eq. (5) and Eq. (6) we must insert expressions for the density, pressure, and stress of each phase. Also mass and momentum transfer between the phases should be defined and inserted.

A. Pore Geometry andVolume Fractions

It has been believed that the smaller pores in the consolidated material of tight sands cause the segregation of mobile fluids such that the wetting fluid occupies the smaller pores. The non-wetting fluid, and possibly an immobile wetting fluid layer lining the walls, would therefore occupy the large pores of this water wet material. This is a result of the large difference between the gas and liquid pressure caused by the small radius of curvature of the gas-liquid interface.

In addition to the assumption that the mobile fluids are separated on a pore level, it has been assumed that two networks for the fractured porous media (Fig. 1) exist.

Based on this network assumption, five different phases exist in the model as specified bellow:

1. Solid phase:
   1. gas in low permeability porous media (gp)
   2. liquid in low permeability porous media (lp)
   3. rock matrix (r)
2. Fracture
   1. gas in the fracture (gf)
   2. liquid in the fracture (lf),

Governing equations for all these phases should be developed from Eq. (5) and Eq. (6).

The volume fraction for each phase can be defined as:

(7)

where V^α(t,x) is the volume of phase α, which is a function of the position of the volume element (V) in the reservoir and also a function of time. Using equation (7), the volume fraction for all five phases is:

ε^qp + ε^lp + ε^r + ε^lf = 1 (8)

as it can be seen from Eq. (8) gas and liquid phases exist simultaneously in both pore and fracture networks, then one can determine the gas and liquid volume fractions, respectively, as:

(9)

and

(10)

then, Eq. (8) can be rewritten as:

ε^g + ε^l + ε^r = 1 (11)

but from the definition of porosity , which is defined as:

Figure 1: Network Model of Fractured Porous Media

(12)

equation (11) becomes:

(13)

B. Pressures

The pressure under reservoir conditions is very high and depends on the location of production formation. According to Huinink and Michels (2002) and Holditch (1989), the pressure for some reservoirs can reach several thousand psi (10⁶ and 10⁸ Pa). This requires taking the compressibility of moving phases into account. The α phase pressure P^α in the volume element was defined by the intrinsic average phase pressure , in what follows and for simplicity the following notation is used:

(14)

In order to reduce the number of independent variables in the model, the equation of state to relate the density to the pressure for each phase is used.

Gas Equation of State

The equation of state for compressible gas is:

P^gV^g = n^gZ^gR T, (15)

substitute for V^g=m^g/ρ^g and n^g=m^g/M^g into the above equation, then, in both fracture and pore networks the gas density is:

(16)

where M^g is the gas molecular weight, p^g is the pressure of the gas, Z^g is the compressibility factor of the gas, R is the universal gas constant, and T is the temperature.

Liquid equation of State

The liquid phase in this study is water. Under reservoir conditions water usually moves with gas, and it contains organic material and salt in dissolved form. Pakdel (1994) had done an analytical analysis of traveling water into and out of the porous media, her chemical analysis showed significant difference in the concentration of some minerals. This means that the liquid phase contains some salt materials in dissolved form that should be taken into consideration in the model.

The solution of most of the technological problems in gas reservoir engineering requires determination of volume coefficient, coefficients of volume and heat expansions, viscosity and density of water and reservoir condition. The necessary data to find all these physical characteristics of water are pressure, temperature, and mineralization (saltiness).

The density of water under reservoir conditions was found by Mishinko et al. (1984), in the form:

(17)

where is water density under standard conditions (kg/m³), and is the volume coefficient of water under reservoir conditions. Both and are discussed in more details in Mishinko et al. (1984), Shirkovski (1987) and Al-Khlaifat (2005).

Rock Equation of State

In order to validate the assumption of deformable porous medium, the relationship between pressure and density of matrix rock must be taken into account in the model.

Petrographic thin-section analysis of pore geometry and grain size for low permeability porous media have been done by William and Sampath (1982) and Holditch (1989), their analyses showed that all sandstone samples in staged field experiments are fine- or very-fine-grained (0.06 to 0.25 mm). Most silt is coarse silt, between 0.031 and 0.063 mm. The samples are classified texturally as sandstone, silty mudstone and silty clay-stone. The porous media for which the model is developed belongs to texturally (sandstone) group, this kind of stone is classified as granular rock.

A porous granular rock was modeled by an aggregate of identical, randomly stacked spherical particles for which the effective density was found, according to Digby (1981), to be:

(18)

where ρ^s is the grain density; is the local porosity of the volume element which is assumed to be equal to the flowing porosity of the rock defined by Eq. (13); and δ* is half of the distance by which the centers of two adhering spheres approach one another when a purely normal force Y^s acting through the center of each sphere is applied, δ* can be determined as:

(19)

where a is an average radius (a>b); b is the radius of adhesion region typically between 0.0 and 0.05R^s; and R^s is the grain's radius.

The normal force acting through the center of two adhering spheres is:

(20)

where μ^s is the grain shear modulus (38x10⁹ Pa), and v^s is the grain Poisson's ratio (0.2).

The force Y^s for a set of spheres, where the pressure can reach the pressure of granular rock in reservoir P^r, is purely hydrostatic loading (overburden) and determined as:

(21)

Equation (21) agrees exactly with the similar equation discussed in Brandt (1955), for the special case of a "dry" packing of spheres, where in Brandt's paper the average number of contact points is n=8.84.

From Eqs. (18), (19) and (21), one obtains the following equation for the normalized contact radius, a/R^s, as:

(22)

From Eq. (22) we can find the value of (a).x/R^s satisfies the following cubic equation:

(23)

C. Capillary Equilibrium

When two immiscible fluids are in contact in the interstices of a porous medium, a discontinuity in pressure exists across the interface separating them. Its magnitude depends on the interface curvature at the point. Here "point" is the microscopic point inside the void space. The difference in pressure is called capillary pressure P^c, and determined as:

P^c=P^nw-P^w (24)

where P^nw, P^w are the pressures in the nonwetting and wetting phases, respectively. Because of the two regions assumption in the model, Eq. (24) can be applied to both pore networks and fractures in the form:

P^cp=P^gp-P^lp (25)

and

P^cf=P^gf-P^lf (26)

where P^cp, P^cf are the capillary pressures in the pore and fracture networks, respectively; P^gp, P^gf are the gas pressures in the pore and fracture networks, respectively; and P^lp, P^lf are the liquid pressures in the pore and fracture networks, respectively.

Capillary effect takes place usually when wetting phase becomes in contact with porous matrix. If the fluid pressure in the fracture is high enough to overcome capillary pressure in the matrix, the fluid penetration into the pore networks occurs. Pressure required to overcome capillary forces at the entrance of the pore is called pore entry pressure. Capillary pressure value depends on the diameter of the pores in the pore networks and on the fracture openings in the fracture network.

Because of the two networks' assumption in the matrix, Eq. (24) can be written simultaneously for gas and liquid pore networks as the following:

P^cgp=P^gp-P^lf (27)

and

P^clp=P^gf-P^lp (28)

where P^cgp is the capillary pressure between the gas in the pore network and the liquid in the fracture; and P^clp is the capillary pressure between the liquid in the pore network and the gas in the fracture.

As long as Reynolds number in porous media is small, all of the above capillary equilibrium equations are valid under both static and dynamic conditions (Semrau, 1986). This pressure equilibrium across various interfaces of pores results in the equality between the sums of Eqs. (25-26) and Eqs. (27-28).

In the current model, under constant hydrostatic (overburden) pressure, the decrease of fluid pressure because of its production causes an increase which leads to closing the smaller diameter (water network) and decrease the number of interconnection between water and gas networks, this causes significant change in the rock matrix porosity, which modifies the mean pore sizes. Because of the deformable medium assumption, the average pore size for both liquid and gas networks will be found by using statistical approach (Section IV.G). The above effect does not take place in the fracture, because the later is composed of large pores that lead to minimal volume variations of the fracture network. Sources of useful information for the dynamic response of porosity to changes in rock pressure are limited and the information itself tends to be more qualitative than quantitative. Thus in the current model a constant fracture volume is acceptable in the considered volume element, therefore ε^gf +ε^lf =const.

D. Gas and Liquid Mass Transfer

The characteristic feature of an unsteady-state motion of a fluid in fractured rocks is the fluid transfer between the rock matrix and the fractures. Therefore, in investigating the flow of fluids in fractured porous medium it is necessary to take into consideration the outflow of fluids from the matrix blocks into the fractures.

The process of fluid transfer from the pores takes place essentially under a sufficiently smooth change of pressure, and, therefore, it can be assumed that this pressure is quasi-stationary, i.e. it is, explicitly, independent of time. It is obvious in such a case that during homogeneous fluid flow in the fractures, the volume of the fluid V^f, which flows from the matrix blocks into the fractures per unit of time and unit of volume of the rock, depends on the following (Zheltov, 1960),: (a) viscosity of the fluid μ^f; (b) pressure difference between the pores and the fractures P^fp, P^ff; and (c) certain characteristic of the rock, which can only be geometrical one, i.e. they might be the dimensions of length, area, volume, etc. The fluid volume was found to be:

(29)

where η is some new dimensionless characteristics of the fractured rock; and P^fp, P^ff are the fluid pressures in the pore and fracture, respectively. Thus for the mass of the fluid which flows from the pores into the fracture per unit of time, per unit volume of the rock, the following equation is valid (Zheltov, 1960):

(30)

where ρ^f is the density of the fluid. Because of the two phase flow in the fracture, the mass transfer between either gas or liquid network and the fracture is determined by Eq. (30) as:

(31)

and

(32)

where c^g and c^l are mass transfer coefficients which characterize the gas and liquid transfer from the pore network to the fracture, these coefficients may not be directly related to the permeability of the porous medium. There is a similarity in the form between the above equations and Darcy's law (Petkovic et al., 2004).

E. Viscous Stress Tensors

The model of fractured porous media is characterized by three volume fractions: ε^r, ε^p and ε^f corresponding to the solid phase, fluid phase in the pores and fluid phase in the fractures, respectively. Newton's viscosity law implies that the fluid has the following properties (Bird et al., 2002): (1) Stress is a linear function of strain rate; (2) The coefficients in the expression for the stress are functions of the thermodynamic state; (3) When the fluid is stationary, the stress is the thermodynamic pressure; (4) The fluid is isotropic. The porous medium is isotropic; (5) The stress tensor is symmetric; (6) The mechanical and thermodynamic pressures are equal. The average normal viscous stress is zero. Stokes assumption applies λ= -2μ/3.

The averaged viscous stress expression derived by Al-Khlaifat (2005) is:

(33)

In the above equation μ^α has been assumed approximately constant within an averaging volume but may vary globally, especially for rock phase. The last two terms in this equation involve surface integrals of the fluid velocity along the α-β phase interface. Knowing that the porous medium deforms but not rapidly, the gradient of these surface integrals will be negligible because of the no-slip condition over the fluid solid interface. Taking what stated above into consideration and applying averaging property to the left side of Eq. (33), then Eq. (33) reduces to:

(34)

By further manipulation, Eq. (34) becomes:

(35)

For simplicity, terms containing ∂_j μ^α and ∂_j μ^λ will be assumed to be negligible in comparison with the second derivative terms, then Eq. (35) becomes:

(36)

F. Interfacial Momentum Transfer

Average momentum transfer across the interface A*^αβ between phases α and β in the volume element V is represented by the last term of the right hand side of Eq. (6). The interfacial momentum transfer term is expressed in terms of body forces, Q_i^βα. In the pore network, both fluids will transfer momentum only with the rock (or with connate water), in this case Darcy and non-Darcy equations will be used to express the momentum transfer terms for liquid and gas phases, respectively, in the following forms:

(37)

and

(38)

where Q_i^lr and Q_i^rg are the momentum exchanges from liquid and gas in both pores and fractures to the rock, respectively; k^l and k^g are the effective phase permeabilities of the liquid and gas phases, respectively, during the simultaneous filtration of multiphase systems. This type of permeability depends on the physico-chemical properties of both porous medium and each phase taken separately, the percentage of phases in the system, and the actual pressure gradients. These permeabilities can be defined by:

k^ζ=k^ABSζk^RELζ (39)

where k^ABSζ ≡ k^Aζ is the absolute permeability of a porous medium, which measures the ability of porous medium to transmit only one of the phases, either gas or water, and always determined experimentally. It should be mentioned that there is no physicochemical interaction between the porous medium and the one phase fluid; k^RELζ ≡ k^rζ is the relative permeability, which usually defined as the ratio between the effective and the absolute permeabilities. Relative permeability for granular rock has been determined by Gimattodinov (1974) as a function of phase saturation by the following formulas:

- Relative permeability of gas phase k^rg, where S^g is the coefficient of the gas saturation is:

k^rg (S^gp)=0 for 0≤S^gp ≤0.1 (40)

and

(41)

- Relative permeability of water phase k^rw is:

(42)

and

k^rw (S^gp)=0 for 0.8≤S^gp ≤1 (43)

where saturation of ζ phase S^ζ is the fraction of the pore volume occupied by phase ζ. Obviously, for the considered volume element the sum of S^ζ is equal to one, and for the two network assumption, one can write:

S^gp+ S^gf + S^lp + S^lf = S^g + S^l = 1 (44)

The saturation of ζ phase S^ζ, where ζ can be either gas or liquid is determined as:

(45)

where V^ζ(t,x) is the volume of phase ζ that is a function of the position of the volume element V in the reservoir, and also a function of time; V^υ is the total volume of pores (voids) in the volume element.

Next, Eq. (44) is expressed in terms of volume fraction and porosity, where the porosity is defined as the property of the rock to contain voids (pores, coverns, and fractures (fissures)). It thus determines the ability of a rock to hold gas and liquid and can be found by the following formula:

(46)

substituting the V^υ value from Eq. (46) into Eq. (45) and using the definition of volume fraction, Eq. (7), to obtain:

(47)

using Eq. (47), the relative permeability Eqs. (40-43) will be expressed as function of ε^gp/ instead of S^gp, and Eq. (44) becomes:

ε^gp +ε^gf +ε^lp +ε^lf = . (48)

Because of the deformable porous medium assumption gas and liquid volume fractions will be found using statistical approach (Section IV.G).

As the fracture dimensions (fracture opening) is much larger than the pore diameters, the momentum transfer across the interfaces of gas-liquid, gas-rock, and liquid-rock, should be considered differently from that one in the pore.

As long as water is the wetting fluid in this study, it tends to be in touch with porous medium and wet it, leaving the gas phase in the fracture in-between two layers of water moving in the same direction. Interfacial momentum transfer for annular flow with liquid on the wall for flow in a round tube was found by Lyczkowski et al. (1982) as:

(49)

Because the cross section of the fracture is not circular the hydraulic radius r^h is used instead of the tube diameter D, where the former is defined as r^h=s/p, in which S is the cross section of the stream (fracture) and p is the wetted perimeter. Replacing the diameter D of circular tube by 4r, Eq. (49) for two-phase flow in the fracture becomes:

(50)

In the case when gas phase flows on the surface of the porous medium, ε^gp is replaced by ε^lf.

G. Statistical Model of Porous Medium

Parameters to be Statistically Averaged

Previous statistical models (Haring and Greenkorn, 1970; Fatt, 1956; Dullien, 1975; Koplik and Lasseter, 1985; Brutsaert, 1966; Tsang and Tsang, 1987) included non-uniformity of porous medium, expressed by the change of pore radius. Real porous media even though seemingly homogeneous and isotropic are most often non-uniform, and this non-uniformity may affect flow phenomena.

A parametric statistical model is used usually to calculate macroscopic properties of porous medium such as volume fractions, capillary pressure, and permeability in order to relate these properties and others to the structure of the model. As it was stated in Sections IV.D and IV.F, the average pore size for both liquid and gas networks should be found using statistical approach. Moreover the pores can be described in terms of their, orientation, radius, length, and number, all these terms appear simultaneously in volume fraction expression Eq. (7), which can be expressed as:

(51)

in which ζ can be either gas or water, is the average ζ network pore diameter, and is an averaged parameter defined as:

(52)

where is the average number of pores occupied by ζ phase, and is the average ζ network pore length.

It is obvious from Eq. (51) that the current model is constructed of two parameter distributions functions, where the orientation is considered to be random with all directions allowed.

Classical method of estimating and was described by Dullien (1979) as follows:

(53)

where is an averaged quantity for ζ phase, r=0,2,3 corresponds to number average, surface average, and volume average, respectively, and n(D^ζp) is the density function of the number distribution.

The purpose of the current paper is to develop a mean of estimating of two parameters and of random network. The number average (r=0; Herro, 1973) is used. As we have two networks assumption in the pore, there is an intermediate diameter d*, at which we cannot say exactly, which phase occupies the pore, based on that, the liquid phase diameter in pore network will be d^min≤d^lp≤d^*, and the gas phase diameter varies d^*<d^gp<d^max. When the diameter is smaller than d^min, we consider no movement occurs for water phase, and when it is greater than the maximum diameter the flow will take place in the fracture where we have both phases flow simultaneously. Using Eq. (53) the average liquid and gas pore diameters are expressed as:

(54)

and

(55)

where f(d^lp) and f(d^gp) are liquid and gas density functions, respectively.

Let us now go back to the parameter defined by the multiplication of the number of pores by the length . We consider the volume element to be consisted of for both phases in pore networks, so:

(56)

Because of water wettability character, we consider . Moreover, in Eq. (56) will be arbitrarily assigned, which reduces our effort just to find either or . Also we assign some floating parameter to be defined as:

(57)

where this floating parameter is the upper bound for water phase in the considered porous medium. When , the porous medium is equally saturated by both phases, and when , the porous medium is fully saturated by water phase, no gas phase exists.

In order to determine the quantity using Eq. (53), we have to assign another parameter , as the minimum pore parameter filled with liquid phase, which does not contribute to the flow. Then using equation (53) and taking into consideration the above assumptions, the average of the number of pores times length parameter for liquid phase will be:

(58)

where g(L^lp) is the liquid pore density function for parameter L^lp.

Density Distribution Functions

Most of the distribution functions, were applied to porous medium studies, have given different results depending on the values of parameters involved in each function. Pakdel (1994) was the first researcher, who applied asymptotic distribution of extremes to such kind of study. The asymptotic theory works well in case of knowledge of initial distribution, sample size, and parameters involved in initial distribution. However, in general, the initial distribution is unknown, which makes the problem more complicated to be solved.

Our main purpose of choosing a distribution function is to be able to draw inferences about pore size of the considered porous medium sample. Usually an empirical sampling distribution will be more or less irregular, but its form may suggest that the pore size distribution is closely normal, Poisson or some other well-known mathematical type. As long as we do not have the empirical distribution function for the considered porous medium, the model will be made non-uniform by distributing the two parameters d^ζp and L^ζp, according to the gamma function. The gamma function is a parametric distribution, which gives a variety of skew, and symmetric shapes where its density function (Abramowitz and Stegun, 1972) is:

(59)

The above density function is valid for 0<x<∞, otherwise its value is zero. In Eq. (59) α* and β* are parameters that determine the specific shape of the curve. It is obvious that the gamma density reduces to the exponential form when α*=1. The constant Γ(α*)(β*)^α* is simply the correct one to make the integral of f(x) from zero to infinity equal to one. The parameter α* need not to be an integer, but α* and β* must both be positive. The gamma function, Γ(α*), is a kind of generalized factorial, which can be found by taking the improper convergent integral, for n*>0, as:

(60)

to be a gamma function of n*. Integrating by parts yields:

(61)

If we let n*=1 in Eq. (60), we obtain Γ(1)=1, then in general for any positive integral value of n*, we have:

(62)

then, Γ(α*) is determined as:

(63)

For the purpose of parametric study, one of the shape parameters can be kept constant while changing the other.

Averaged Parameters

Using Eq. (59), the liquid and gas density functions in terms of pore diameter d^ζp and L^ζp parameters, are determined, respectively, as:

(64)

and

(65)

substitute Eq. (64) into Eqs. (54) and (55), respectively, we have:

(66)

and

(67)

by substituting Eq. (65) into Eq. (58), we have

(68)

Integrands in the last three equations are similar. To find the value of these integrals, denote the numerator of Eq. (66) by Iⁿ and the denominator of the same equation by I^d, as follows:

(69)

and

(70)

These two integrals can be solved using gamma function. Use the following transformation; t^lp=d^lp/β*, then dd^lp=β*dt^lp, and, in Eq. (69), let α*=α'-1. After making the substitution the last two equations above become:

(69)

and

(72)

If we split the gamma function into three integrals with limits from 0 to d^min/β*, d^min/β* to d^*/β*, and from d^*/β* to ∞, the above integrals become:

(73)

and

(74)

Because of the extensions of the gamma function to integrals of variable limits, we obtained two new integrals that usually called incomplete gamma functions, they can be written for Eq. (73) as:

	(75)
	(76)

and for Eq. (74) as:

	(77)
	(78)

Substitute the last four incomplete gamma functions back into Eqs. (73) and (74), respectively, taking into consideration that the first term in the right side of these equations is gamma function, we obtain

(79)

and

(80)

The first term in the right side of the last two equations can be defined by Eq. (62), while the second and third terms are found according to Abramowitz and Stegun (1972), in general, as:

(81)

and

(82)

where

(83)

Applying Eqs. (81) through (83) to Eqs. (79) and (80), we obtain:

(84)

and

(85)

Using Eq. (62) to find Γ(α) and substitute for α'=α*+1 into Eq. (84), then substitute them back into Eq. (66) to get the final form of the averaged diameter as:

(86)

Similarly, we can find the gas network averaged diameter and the parameter to be, respectively:

(87)

and

(88)

V. FINAL FORM OF GOVERNING EQUATIONS

After substituting all constitutive equations into the continuity Eq. (5) and momentum Eq. (6) for all five phases a set of continuity and momentum equations is obtained as follows:

Continuity Equations:

Liquid phase in the pore network

(89)

Gas phase in the pore network

(90)

Liquid phase in the fracture network

(91)

Gas phase in the fracture network

(92)

Solid (rock) phase

(93)

Momentum Equations

Liquid phase in the pore network

(94)

Gas phase in the pore network

(95)

Liquid phase in the fracture network

(96)

Gas phase in the fracture network

(97)

Solid phase

(98)

VI. MODEL VALIDATION EXAMPLE

A. Physical Domain

The finite volume technique incorporated in the CFX-Flow3D code that was used by Al-Khlaifat and Arastoopour (2001) was used to obtain the numerical values for the flow parameters for the geometry shown in Figure 2. A 42 cm long low permeability porous physical domain with a height of 10 cm and width of 10 cm was considered.

A single vertical fracture was represented by the small blocks in the middle of larger matrix blocks. The fracture thickness was chosen as 0.2 cm with an absolute permeability of 6μm² and a porosity of 80% and negligible capillary pressure. Uniform matrix blocks with low permeability of 2x10^-5μm² and porosity of 10% were assumed. The matrix block size increased from the fracture towards the outer blocks. The pressure drop was assumed to be in x-direction only.

Figure 2: Grid Representation (19x5x5)

All the simulations carried out were two-phase flow (water and gas), three-dimensional, isothermal, laminar, transient, and compressible flow. The physical domain was divided horizontally into 19 blocks, with the fractured layer located at the tenth block from both the most-left and the most-right blocks. The matrix block size was chosen in a way that avoids large size contrast between neighboring blocks, while minimizing the number of blocks. Block size is as follows: 2 cm in z-direction; 2 cm in y-direction; and 6.3078, 3.6900, 2.6178, 2.0310, 1.6590, 1.4028, 1.2150, 1.0728, 0.8488 and 0.2000 cm in x-direction from either the most-left or the most-right block to the fracture.

B. Initial and Boundary Conditions

A uniform initial gas saturation of 0.6 was assumed in the physical domain. The inlet gas saturation was set to be 0.4 at x=0. The initial water pressure in the physical domain was assigned lower than gas pressure by the capillary pressure which was a function of initial water saturation distribution. At the walls (y=0, y=10cm; z=0, z=10cm) the no-slip conditions were assumed. Inlet and outlet pressures were assumed 4.1 MPa and 0.1 MPa, respectively.

C. Results and Discussion

Figure 3 shows the water saturation distribution throughout the physical domain. It is clear from this figure that water saturation increases from its initial value of 0.4 to become constant at steady state of 0.6 which is the water saturation value assumed at the inlet.

Water pressure profile is shown in Fig. 4. It behaves almost linearly, except for the change of slope in the fracture of high porosity and permeability.

Water and gas velocity profiles are shown in Fig. 5 and 6, respectively. They have almost similar trends but the gas phase velocity is almost twice larger than water phase velocity. This increase is caused by higher gas permeability compared to water. At steady state, one can see clearly that the fracture has higher gas and water velocities than anywhere else in the matrix blocks.

Figure 3: Water Saturation Distribution

Figure 4: Water Pressure Distribution

Figure 5: Water Velocity Distribution

Outlet water flow rates for different permeabilities are shown in Fig. 7.

Low water flow, for lower permeability, during the initial period was caused by low water saturation. The movement of water-gas front toward the exit (x=42cm) results in increasing water flow rate to reach steady state value. Also, Fig. 7 shows the effect of different fracture permeabilities on water flow rates. Higher fracture permeability results in higher exit water flow rate during the transient condition. Under steady state conditions, variation in fracture permeability, does not have pronounced influence on water flow rate, this means that water flow rate at steady state is mainly controlled by low permeability region. An increase in the fracture permeability results in a slight increase of gas flow rate for both transient and steady state regimes as shown in Fig. 8.

Figure 6: Gas Velocity Distribution

Figure 7: Outlet Water Flow Rate. (Fracture Permeability-I = 6μm²; Fracture Permeability-II = 7μm²; Fracture Permeability-III = 8μm²)

Figure 8: Outlet Gas Flow Rate(Fracture Permeability-I = 6μm²; Fracture Permeability-II = 7μm²; Fracture Permeability-III = 8μm²)

VII. CONCLUSIONS

A transient three-dimensional isothermal two-phase flow model has been developed for the flow of gas and liquid through low permeability fractured porous media using the general continuity and momentum equations. The fundamental assumptions made in the development of the model are that the momentum transfer between fluids inside the porous structure is negligible and exists in the fracture, and a condition of capillary equilibrium exists throughout the media.

The developed two-phase flow model is very important because of its detailed description of the occurring hydraulic phenomena and phase interactions. Furthermore, the model is most flexible because it uses separate conservation equations for each phase in both fracture and matrix networks.

Moreover, the formulation of the constitutive and closure relations for the two-fluid model is at highly acceptable level. The major problems in the modeling of various interfacial transfer terms have been addressed in detail. The interfacial transfer of mass and momentum is proportional to the interfacial area and driving force. The geometrical effects on the flow phenomena are taken into account by the pore size distribution. The pore size distribution of the studied porous media was considered to be represented by three mean pore diameters: liquid pore network, gas pore network, and fractures. An average pore diameter and length for each network were calculated using the gamma distribution function. Furthermore, the movement of rock matrix was considered, this results in a model with five different phases: 1) gas in the porous matrix (gas network), 2) liquid in the porous matrix (liquid network), 3) gas in the fracture, 4) liquid in the fracture, and 5) rock matrix.

The model was validated using a hypothetical fractured domain with dual porosity, dual permeability for the pressure and capillary pressure driven flow. Sensitivity analysis for different fracture permeabilities was conducted as well.

REFERENCES
1. Abramowitz, M.,I. and Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publication (1972).         [ Links ]
2. Alder, P.M. and J.F. Thovert, Fractures and Fracture Networks, Kluwer (1999).         [ Links ]
3. Al-Khlaifat, A., Two Phase Flow through Low Permeability Fractured Tight Sand Porous Media, Ph.D. Thesis, IIT (1998).         [ Links ]
4. Al-Khlaifat, A., and H. Arastoopour, "Simulation of two phase flow through anisotropic porous media", J. Porous Media, 4, 275-281 (2001).         [ Links ]
5. Al-Khlaifat, A., "Mathematical Modeling of Two Phase Flow in Low Permeability Porous Media", V Jordan International Conference, Amman-Jordan, 12-15 Sept. (2005).         [ Links ]
6. Arbogast, T., "Gravitational Forces in Dual-Porosity Systems: II. Computational Validation of the Homogenized Model", Transport in Porous Media, 13, 205-220 (1993).         [ Links ]
7. Bird, R.B., W.E. Stewart and E.N. Lightfoot, Transport Phenomena, John Wiley & Sons, 2^nd Edition (2002).         [ Links ]
8. Brandt, H., "A Study of the Speed of Sound in Porous Granular Media", ASME Journal of Applied Mechanics, 22, 479-486 (1955).         [ Links ]
9. Brutsaert, W., 'Probability Laws for Pore Size Distributions', Soil Science,101, 85-93 (1966).         [ Links ]
10. Celis, V., J. Guerra and G.D. Prat, "A new Model for Pressure Transient Analysis in Stress Sensitive Naturally Fractured Reservoirs", SPE Advance Technology Series, 2, 126-135 (1992).         [ Links ]
11. Digby, P.J., "The Effective Elastic Moduli of Porous Granular Rocks", Journal of Applied Mechanics, 47, 803-808 (1981).         [ Links ]
12. Dullien, F.A., "Single Phase Flow through Porous Media and Pore Structure", Chemical Engineering Journal, 10, 1-34 (1975).         [ Links ]
13. Dullien, F.A., Porous Media Fluid Transport and Pore Structure, Academic Press, New York (1979).         [ Links ]
14. Evans, R., "A Proposed Model for Multiphase Flow through Naturally Fractured Reservoirs", SPE Journal, 669-680 (1982).         [ Links ]
15. Fatt, A., "The Network Model of Porous Media", Petroleum Transaction AIME, 207, 144-159 (1956).         [ Links ]
16. Fung, L.S., "Numerical Simulation of Naturally Fractured Reservoirs", 8th Middle East Oil Show & Conference, 2, 203-213 (1993).         [ Links ]
17. Gimattodinov, S.K., Manual of Oil and Gas Production, Nedra, Moscow (1974).         [ Links ]
18. Guzman, R.E. and A. Khaled, "Fine Grid Simulation of Two-Phase Flow in Fractured Porous Media", SPE Journal, 605-615 (1992).         [ Links ]
19. Haring, E. and R. Greenkorn, "A Statistical Model of a Porous Medium with Non uniform Pores", AIChE Journal, 16, 477-483 (1970).         [ Links ]
20. Herro, J.J., Estimation of the Average Frequency of a Random Process, Ph.D. Thesis, IIT (1973).         [ Links ]
21. Holditch, S.A., Staged Field Experiment No. 2. Application of the Advanced Geological, Petrophysical and Engineering Technologies to Evaluate and Improve Gas Recovery from Low Permeability Sandstone Reservoirs, Report No. GRI-89/0140, Gas Research Institute, CER Corporation (1989).         [ Links ]
22. Huinink, H.P. and M.A.J. Michels, "The influence of buoyancy on drainage of a fractal porous medium", Physical Review,E 66, 064301 (2002).         [ Links ]
23. Jacob, B. and B. Yehuda, "Transport Phenomena in Porous Media -Basic Equations", Transport in Porous Media, 7, 3-61 (1992).         [ Links ]
24. Koplik, J., and T.J. Lasseter, "Two-Phase Flow in Random Network Models of Porous Media", SPE Journal, 89-100 (1985).         [ Links ]
25. Le Gallo, Y., Two Phase Flow of Gas-Liquid in Pipes and Low Permeability Porous Medium, Ph.D. Thesis, IIT (1990).         [ Links ]
26. Lee, B.Y.Q. and T.B.S. Tan, "Application of Multiple Porosity/Permeability Simulator in Fractured Reservoir Simulation", SPE Journal, 181-192 (1987).         [ Links ]
27. Lyczkowski, R.W., D. Gidaspow and C.W. Solbrig, "Multiphase Flow Models for Nuclear Fossil and Biomass Energy Production", Advances in Transport Processes, 2, 198-351 (1982).         [ Links ]
28. Mishinko, I.T., B.A. Sakharov, B.G. Gron and G.I. Bogomolny, Collection of Problems in the Technology of Oil Production, Nedra, Moscow (1984).         [ Links ]
29. Nitao, J.J., Theory of Matrix and Fracture Flow Regimes in Unsaturated Fractured Porous Media, Lawrence Livermore National Laboratory (1990).         [ Links ]
30. Pakdel, P., Analysis of Two Phase Flow through Low Permeability Tight Sand Porous Media, Ph.D. Thesis, IIT (1994).         [ Links ]
31. Panfilov, M. and M.B. Panfilov, Macroscale Models of Flow through Highly Heterogeneous Porous Media, Kluwer (2000).         [ Links ]
32. Petkovic, J., H.P. Huinink, L. Pel and K. Kopinga, "Diffusion in porous building materials with high internal gradients", Journal of Magnetic Resonance, 167, 97-106, 2004.         [ Links ]
33. Semrau, J.T., Analysis of Two Phase Flow through Low Permeability Media, Ph.D. thesis, IIT (1986).         [ Links ]
34. Shirkovski, A.I., Production and Exploitation of Gas and Gas-condensate fields, Nedra, Moscow (1987).         [ Links ]
35. Slattery, J.C., Momentum, Energy, and Mass Transfer in Continua, McGraw Hill (1972).         [ Links ]
36. Thomas, L.K., N.D. Thomas and R.G. Pierson, "Fractured Reservoir Simulation", SPE Journal, 42-54 (1983).         [ Links ]
37. Tsang, Y.W., and C.F. Tsang, "Channel Model of Flow through Fractured Media", Water Resources Research, 23, 467-479 (1987).         [ Links ]
38. Warren, J.E. and P.J. Root, The Behavior of Naturally Fractured Reservoirs, Gulf Research & Development Co. (1963).         [ Links ]
39. William, G.G. and K. O'Neill, "On the General Equations for Flow in Porous Media and their Production to Darcy's Law", Water Resources Research, 12, 184-194 (1976).         [ Links ]
40. William, C. and K. Sampath, "Evaluation of Pore Geometry of Some Low-Permeability Sandstones-Unita Basin", Journal of Petroleum Technology, 65-70 (1982).         [ Links ]
41. Zekai, S., "Fractured Media Type Curves by Double-Porosity Model of Naturally Fractured Rock Aquifers", Bull. Tech. Univ. Istanbul, 41, 103-133 (1988).         [ Links ]
42. Zheltov, I.P., G.I. Barenblatt and I.N. Kochina, "Basic Concepts in the Theory of Seepage of Homogeneous Liquids in Fissured Rocks (Strata)", PMM, 24, 1286-1303 (1960).         [ Links ]

Received: February 12, 2007.
Accepted: October 30, 2007.
Recommended by Subject Editor José Pinto.