The Similarity Theory Applied to the Analysis of ... - ACS Publications

1592

Energy & Fuels 2007, 21, 1592-1600

The Similarity Theory Applied to the Analysis of Multiphase Flow in Gas-Condensate Reservoirs Luis F. Ayala* and Jean P. Kouassi The PennsylVania State UniVersity, UniVersity Park, PennsylVania, and British Petroleum GoM Deepwater DeVelopment, Houston, Texas ReceiVed October 9, 2006. ReVised Manuscript ReceiVed January 22, 2007

An important part of the world’s hydrocarbon reserves is constituted of gas-condensate deposits, and as such, they have been the subject of intensive research over the years. Gas-condensate reservoirs are complex systems, and their performance analysis involves the handling of highly nonlinear partial differential equations and sophisticated thermodynamic models. For the case of noncondensing natural gas reservoirs, analytical solutions based on the exponential integral function are readily available for the analysis of infinite-acting behavior in production and injection wells. No such tool is available for the analysis of infinite-acting behavior in gas-condensate reservoirs, and the associated analysis of well deliverability is a topic of continuous research. This paper investigates the implementation of the similarity variable theory in the context of the analysis of depletion behavior in natural gas reservoirs. The main objective of this work is to present a novel application of the similarity theory for the case of the analysis of complex systems such as gas-condensate reservoirs, without introducing any simplifications to the inherent nonlinearities present in the associated governing equations.

Introduction The use of similarity variables is common-place in the solution of the differential equations typically encountered in engineering problems. Classic engineering textbooks utilize the concept of the similarity variable in order to arrive at analytical solutions of various problems of interest. Carslaw and Jaeger1 discussed the use of the Boltzmann’s transformation η ) r/xt for the solution of heat conduction problems and employed several other variable transformations that made possible the analytical solution of other problems. Bird et al.2 present several problems of fluid flow, heat, and mass transfer where the use of the method of similarity solutions (also known as the method of combination of variables) is instrumental for the transformation of partial differential equations (PDEs) into one or more ordinary differential equations (ODEs) and the attainment of an analytical solution. The importance of the similarity variable theory in reservoir engineering problems cannot be overstated. The field of well-test analysis in petroleum and natural gas engineering has been built around the concept of the exponential integral solution, which is the result of applying the similarity theory to the partial differential equations for single-phase flow in porous media in a radial geometry.3-5 Nevertheless, in the area of reservoir engineering and over the years, the application * To whom correspondence should be addressed. 122 Hosler Building, University Park, PA 16802. Phone: (814) 235-1386. Fax: (814) 865-3248. E-mail: [email protected]. (1) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solids, 2nd ed.; Oxford University Press: Oxford, U.K., 1959. (2) Bird, R. B.; Stewart, W.; Lightfoot, E. Transport Phenomena, 2nd ed.; John Wiley & Sons, Inc.: New York, 2002. (3) Theory and Practice of the Testing of Gas Wells, 3rd ed.; Energy Resources and Conservation Board: Calgary, Canada, 1975. (4) Horne, R. Modern Well Test Analysis, 2nd ed.; Petroway Inc.: Palo Alto, CA, 2002; 5th reprint. (5) Lee, J.; Rollins, J. B.; Spivey, J. P. Pressure Transient Testing; Henry Doherty Memorial Fund of AIME: Richardson, Texas, 2003; SPE Textbook Series, Vol. 9.

of the similarity solution theory has been limited to the analysis of relatively simple and linear (or linearized) problems. In 1981, O’Sullivan6 recognized that the similarity variable method can also be used in a rigorous way for the solution of highly nonlinear problems without neglecting any of the nonlinearities. O’Sullivan proposed the use of the similarity method for the analysis of the highly nonlinear problem of fluid flow in geothermal reservoirs. Later on, Doughty and Pruess7 utilized O’Sullivan’s approach to model the two-phase fluid flow and heat-transfer phenomena taking place in geologic repositories of nuclear waste such as those in Yucca Mountain, Nevada. Shortly after this publication, Doughty and Pruess8 improved their formulation and enhanced the mathematical representation of the problem by incorporating a new component in their system (air). Moridis and Reagan9 have recently utilized the similarity method to obtain analytical solutions for the highly complex and nonlinear problem of fluid flow in gashydrate reservoirs. Following O’Sullivan6 and Doughty and Pruess,7 the present work explores the applicability of the similarity variable concept in the formulation of ODEs capable of capturing highly nonlinear flow phenomena encountered during the analysis of the depletion behavior in retrograde gascondensate reservoirs. A gas condensate is a fluid mixture (mainly formed of hydrocarbon components) encountered at near-critical conditions and at reservoir temperatures higher than the fluid’s critical point. (6) O’Sullivan, M. A Similarity Method for Geothermal Well Test Analysis. Water Resour. Res. 1981, 17 (2), 390-398. (7) Doughty, C.; Pruess, K. A Similarity Solution for Two-Phase Fluid and Heat Flow Near High-Level Nuclear Waste Packages Emplaced in Porous Media. Int. J. Heat Mass Transfer 1990, 33 (6), 1205-1222. (8) Doughty, C.; Pruess, K. A Similarity Solution for Two-Phase Water, Air, and Heat Flow Near a Linear Heat Source in a Porous Medium. J. Geophys. Res. 1992, 97 (B2), 1821-1838. (9) Moridis, G.; Reagan, M. A Similarity Solution for Gas Production from Dissociating Hydrates in Geologic Media. Transp. Porous Media 2006, in press.

10.1021/ef060505w CCC: $37.00 © 2007 American Chemical Society Published on Web 03/13/2007

Multiphase Flow in Gas-Condensate ReserVoirs

Energy & Fuels, Vol. 21, No. 3, 2007 1593

Gas condensates are also called retrograde gas condensatess or simply, retrograde gasessbecause of the distinct (retrograde) behavior that they exhibit upon pressure depletion. In a gascondensate reservoir, and as a result of isothermal depletion and vapor expansion, a hydrocarbon liquid phase drops out of the wholly gaseous mixture once the pressure of the system falls below the fluid’s dewpoint pressure. This behavior is characterized as “retrograde” because the appearance of a liquid phase upon vapor expansion would never be anticipated for the case of simple (pure) fluid systems. Condensate dropout at reservoir conditions during production operations in a gas-condensate field is not desirable. The reason is two-fold: it causes an appreciable loss of highly marketable surface condensate and adversely affects gas well deliverability. The analyses of the effect of condensate dropout on depletion performance, well deliverability, and pressure transient analysis of gas-condensate reservoirs have been under the spotlight over the years.10-15 Bøe et al.11 and Vo12 were the first to recognize that the similarity variable could be applied to the numerical analysis of gas-condensate reservoirs of uniform permeability; but no attempt was made at generalizing and solving the corresponding ordinary differential equation system. In this study, we explore this avenue and develop a method, outlined in the following sections, that has the capability of rigorously solving the highly nonlinear two-phase flow equations in radial geometry without introducing simplifications in terms of nonlinearities found among the fluid phase behavior, fluid properties, rock properties, and rock/fluid interaction for this type of system. The method applies to the analysis gas-condensate reservoirs as long as the reservoir is infinite-acting, and it is treated as having a uniform permeability.

) 0. By adding eqs 1 and 2, one obtains the combined continuity equation for the analysis of fluid flow in gas-condensate reservoirs:

(1)

-

1 ∂ ∂ (rF V ) + mog ) (φSoFo) r ∂r o o ∂t

(2)

These are the statements of continuity (mass conservation) for the gas and condensate phases. The terms mgo and mog account for the mass transfer that occurs below dewpoint conditions between the phases, and mgo) -mog. When the condensate is immobile (Vo ) 0), the divergence of the mass flux on the left-hand side of eq 2 vanishes. Under such conditions, the presence of the mass-transfer term (mog) allows for the condensate saturation (So, right-hand side) to build upsdriven by condensation of the gas phaseseven if Vo (10) Fussel, D. D. Single-Well Performance Predictions for GasCondensate Reservoirs. JPT, J. Pet. Technol. 1973, July, 860-870. (11) Bøe, A.; Skjaeveland, S. M.; Whitson, C. H. Two-Phase Pressure Test Analysis. SPE Form. EVal. 1989, December, 601-610; SPE paper 10224. (12) Vo, D. T. Well Test Analysis for Gas Condensate Systems. Ph.D. Dissertation, The University of Tulsa, OK, 1989. (13) Fevang, O.; Whitson, C. H. Modeling Gas-Condensate Well Deliverability. SPE ReserVoir Eng. 1996, NoVember, 221-230. (14) Xu, S.; Lee, J. Two-Phase Well Test Analysis of Gas-Condensate Reservoirs, SPE paper 56483 presented at the SPE Annual Technical Conference, Houston, Texas, Oct. 3-6, 1999. (15) Gringarten, A. C. Well Test Analysis in Gas-Condensate Reservoirs, SPE paper 62920 presented at the SPE Annual Technical Conference, Dallas, Texas, Oct. 1-4, 2000.

F ) mass flux term ) - r(FgVg + FoVo)

(3b)

A ) mass accumulation ) φSgFg + φSoFo

(3c)

and

As it is customary in reservoir engineering analysis, the depletion process of the reservoir is assumed to be isothermal, and thus no energy equation is coupled with the above equations. The mass flux term (F) represents the flow of mass per unit reservoir thickness, and it is considered to be positive when the flow goes toward the center of the radial coordinate system. Darcy’s law governs the flow of fluid in porous media and relates the velocity of each phase to the gradient of the reservoir pressure field, as shown below: Vg ) -

kkrg ∂p µg ∂r

(4a)

where velocities, flow rates, and mass fluxes are related through kkro ∂p µo ∂r

(4b)

Vg ) qg/2πrh

(5a)

Vo ) qo/2πrh

(5b)

F ) -(Fgqg + Foqo)/2πh

(5c)

Vo ) the following equations:

The mass conservation equations that govern the simultaneous flow of gas and condensate in a gas-condensate reservoir in 1-D radial-cylindrical coordinates are written as 1 ∂ ∂ (rF V ) + mgo ) (φSgFg) r ∂r g g ∂t

(3a)

where

Governing Equations

-

1 ∂F ∂A ) r ∂r ∂t

For the sake of brevity and simplicity, the contribution of capillary pressuresswhich could be otherwise included without any fundamental change to the methodology hereby discussedshas been neglected in the presentation of Darcy’s law. By incorporating Darcy’s law into eq 3b, the mass flux term can be expressed in terms of the pressure gradient as follows: F ) krλt

∂p ∂p ) rτ ∂r ∂r

(6a)

where λt ) total mobility ) kro

Fo Fg + krg µo µg

τ ) transmissibility ) kλt

(6b) (6c)

For the purposes of this paper, the classical pressure transient analysis boundary conditions are examined, that is, a single well, located at the center of an infinite reservoir, producing at a constant specified rate Fsp. The traditional line-source wellbore representation is implemented (i.e., wellbore radius approaches a zero value and the mass flux is specified at r ) 0) while the pressure remains unchanged at the outer boundary. Symbolic, initial, and boundary conditions can be stated as F f Fsp as r f 0

(7a)

1594 Energy & Fuels, Vol. 21, No. 3, 2007

Ayala and Kouassi

For the analysis of gas-condensate reservoirs, eqs 3a-6c need to p f pi as r f ∞

(7b)

F ) 0 and p ) pi ∀r at t ) 0

(7c)

be solved within the domain 0 < r < ∞ and 0 < t < ∞, along with the stated boundary and initial conditions.

Implementing the Similarity Variable O’Sullivan6 demonstrated that the nonlinear PDEs that govern fluid flow in geothermal reservoirs can be collapsed to a set of ODEs by introducing the similarity variable η ) r/xtsBoltzmann’s similarity variable. An equivalent and common choice in classical well-testing theory is the use of the transformationη ) r2/t, which leads more neatly to the classical exponential integral solution in the case of single-phase equations.3-5 For the case of gas-condensate reservoirs, the introduction of the similarity variable η ) r/xt in eq 3a yields 1 ∂F ∂η ∂A ∂η ) r ∂η ∂r ∂η ∂t

(8)

1 dF 1 dA )- η η dη 2 dη

(9)

or

Following the same procedure, eq 6a can be rewritten as dp F ) ητ dη

the liquid and gas phases,17 in order to ensure that both values approach each other at near-critical conditions. Thermodynamic equilibrium is imposed by invoking the criterion of equal fugacities for all components both in the gas and in the condensate phases. Equations 11a and 11b, along with their boundary and initial conditions, need to be integrated within the domain 0 < η < ∞. Since the right-hand sides of the ODEs presented in eqs 11a and 11b are functions of η, p, and F, the system of equations takes the form dp ) f1(η,p,F) dη

(12a)

dF ) f2(η,p,F) dη

(12b)

F f Fsp as η f 0

(12c)

p f pi as η f ∞

(12d)

where eqs 12c and 12d represent the corresponding boundary and initial conditions expressed in terms of η. As a matter of convenience, no mobile water was considered in the system, but the flow of a third phase can be accommodated by modifying the mass flux and accumulation terms accordingly. The presence of immobile water does not modify the procedure hereby described, as such a system could be treated as an analogous reservoir with a reduced porosity when rock and compressibility effects are neglected.

Noncondensing Gases: Analytical Treatment (10)

Equations 9 and 10 represent the system of ODEs that governs the flow of fluids in gas-condensate reservoirs, expressed in terms of the similarity variable. If one chooses pressure (p) and mass flux (F) as the dependent variables, this system of ODEs can be expressed as dp F ) dη ητ

(11a)

dF η2 dA )dη 2 dη

(11b)

At this point, it is clear that the similarity method does collapse the original system of PDEs into a set of first-order ODEs in terms of a single independent variablesthe similarity variable, η. Note that the right-hand sides of the ODEs presented in eqs 11a and 11b are sole functions of the independent variable of the system (η) and the selected dependent variables themselves (p and F). When thermodynamic equilibrium is invoked at all points and times, as routinely done in compositional reservoir simulation, the values of mass accumulation (A in eq 3c) can be calculated as a sole function of pressure. A thermodynamic flash calculation at a given pressure and prescribed reservoir temperature would reveal hydrocarbon saturations (So and Sg) at any given reservoir fluid composition. In addition, the values of fluid properties (Fo, Fg, µo, and µg) would be uniquely defined at those conditions. Relative permeabilities (kro and krg), needed in the calculation of transmissibility τ, are calculated as functions of fluid saturation by means of the appropriate relative permeability curves. For the calculation of Fo, Fg, µo, µg, So, and Sg, a hydrocarbon phase behavior and property prediction model is required. In the present study, a thermodynamic model based on the Peng-Robinson equation of state16 is implemented. Viscosities are calculated using the same model for both (16) Peng, D.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59.

When no condensate is present in the reservoir (p > pdew, i.e., single-phase gas conditions), the system of ODEs presented above does accept the well-known exponential integral analytical solution. In the case of this study, the analytical solution for noncondensing (dry or wet) gases can be used as a yardstick to assess the accuracy of the proposed numerical scheme for undersaturated conditions in a gas condensate (p > pdew). For the case of single-phase flow of noncondensing gases, it can be demonstrated that the proposed set of ODEs and initial/boundary conditions can be collapsed to the exponential integral-type solution. For such conditions, the ODE system would be written as µg dp F ) dη ηkFg

(13a)

dF η2 dFg )- φ dη 2 dη

(13b)

F f Fsp as η f 0

(13c)

p f pi as η f ∞

(13d)

where the porous medium has been considered to be incompressible and krg ) 1 (Sg ) 1) and kro ) 0 (So ) 0). In order to obtain an analytical solution, the governing equations above are linearized by introducing the concept of pseudo-pressure,18 defined below: ψ)

∫

p

o

p dp µgZ

(14)

When the real gas law and the concepts of pseudo-pressure and gas isothermal compressibility (cg) are introduced, the above set (17) Lorentz, J.; Bray, B.; Clark, C. Calculating Viscosities of Reservoir Fluids from Their Compositions. JPT, J. Pet. Technol. 1964, October, 1171-1176; SPE Paper 915. (18) Al-Hussainy, R.; Ramey, H. J.; Crawford, P. B. The Flow of Real Gases Through Porous Media. JPT, J. Pet. Technol. 1966, May, 624-636.



of ODEs is recast in the following form: RT F dψ ) dη kMW η

(15a)

dF η2 φMW dψ )µc dη 2 RT g g dη

(15b)

F f Fsp as η f 0

(15c)

ψ f ψi as η f ∞

(15d)

An analytical expression for the mass flux (F) can thus be obtained, by neglecting the remaining (usually weak) nonlinearity µgcgsits value being evaluated at initial conditionssand assuming that the porous medium is homogeneous:

∫

F

Fsp

φµgicgi dF )F 2k

(

F ) Fsp exp -

∫

η

0

ηdη

φµgicgi 2 η 4k

)

(16)

An alternative approach for the handling of the nonlinearity introduced by the term “µgcg” is to define an additional pseudovariable (“pseudo-time”),4-5 which can be necessary at low-pressure operating conditions when compressibility variations with pressure can be significant. When eq 15a and the previous result are utilized, the classical exponential integral solution that governs the flow of a real gas through a porous medium is readily obtained: ψ ) ψi +

(

RTFsp φµgicgi 2 E η 2kMW i 4k

)

(17)

(

)

( )

6

Ui(X + ∆X) ) Ui(X) +

where Ei is the exponential integral function. Equation 17 describes the flow of a real gas through a porous medium, under the stated boundary and initial conditions and using the concept of pseudopressure. It is customarily used in the analysis of pressure transient data obtained from single-phase natural gas reservoirs. This equation can also be expressed in a more conventional form as a function of gas volumetric flow rate specifications at wellbore or surface conditions, using the corresponding equivalences shown below: pwbMW pscMW qg ) 2πhFsp ) qsc ZwbRT RTsc

This set of equations, boundaries, and initial conditions is to be integrated within the domain X ) -∞ (long times, close to the wellbore) to X ) +∞ (early times, far away from the wellbore). The set of equations 19a-19d represents the proposed set of governing equations for fluid flow in gas-condensate reservoirs in terms of the similarity variable. Please note that the implementation of the similarity variable is carried out without introducing any linearization or simplification to the highly nonlinear terms in the equations and that this implementation renders a system of equations that can be solved with traditional ODE solvers, as discussed below. In contrast, the typical linearizations traditionally used for dry (i.e., noncondensing) gases were examined in the previous section. For the case of the forward solution of the proposed well-testing problem, the pressure response is unique for each combination of reservoir and fluid properties and production specification. The backward solution of this problem, however, might accept multiple numerical solutions; that is, different combinations of reservoir and fluid properties might generate the same pressure response for a given production specification. This latter case is not examined in this manuscript. In terms of numerical integration variable X, and for most gascondensate problems, single-phase conditions would exist at X ) +∞ (r ) ∞ and t ) 0) and two-phase conditions would be very likely to exist at X ) -∞ (r ) 0 and t ) ∞). Even though this set of ODEs is not readily amenable to an analytical solution, it can be straight forwardly integrated via the Runge-Kutta method. This is a much less expensive alternative to traditional finite-difference numerical simulators used to solve the original system of PDEs. Press et al.19 presented a detailed explanation of an efficient RungeKutta method, by means of which a fifth-order accurate solution of the system can be generated as follows:

∑ aG j

ij

i ) 1, 2

(20a)

j)1

where U is the two-component vector containing the unknowns or dependent variablessU1 ) p and U2 ) Fs, aj and bj are constants, and the quantities Gij are computed as j-1

Gij ) ∆Xfi′[Ui(X) +

∑d F] jk k

1eje6

(20b)

k)1

(18)

Equations 16 and 17 can be utilized as benchmark equations to evaluate the soundness of the numerical solution of ODE equations proposed previously when single-phase flow conditions prevail, such as it is in the case of natural gases found above or far away from dewpoint conditions.

Condensing Gases: Numerical Treatment A final manipulation of the equations is needed in order to avoid the singularity at η ) 0 in the ODE system of eqs 13a-13d. In order to avoid this singularity, the proposed system of ODEs can be recast in terms of “ln η” in lieu of “η”.6,7 By introducing the proposed change of variables, X ) ln η and d/dX ) ηd/dη, the proposed set of ordinary differential equations with its initial/ boundary conditions becomes dp F ) f1′(X,F,p) ) dX T

(19a)

dF e2X dA ) f2′(X,F,p) ) dX 2 dX

(19b)

F f Fsp as X f -∞

(19c)

p f pi as X f +∞

(19d)

where aj’s and djk’s are constants; f ′ is the vector containing the right-hand-side functions f1′ and f2′ defined by eqs 19a and 19b. For numerical purposes, finite values for these lower and upper limits of integration must be selected (X ) -∞ and X ) +∞). In practice, these infinite limits must be replaced by a set of appropriate finite boundary valuesssay, Xmin and Xmax. For the conditions investigated in this study, the values Xmin ) -15 and Xmax ) +3 proved to be reasonable and convenient choices. In general, we seek the smallest absolute values of Xmin and Xmax that could still properly characterize the boundary conditions represented in eqs 19c and 19d. Clearly, the final solution must be insensitive to the final choice of Xmin and Xmax values. Because of the nature of the boundary conditions stated in eqs 19c and 19d, the Runge-Kutta method cannot “freely” march from beginning (X ) -∞) to end (X ) +∞) of the domain. The reason for that is that the system is being required to simultaneously satisfy very specific boundary conditions at both ends of the domain. This is a two-point or mixed boundary value problem19 for which an “initial guess” must be made for pressure at the boundary X ) Xmin where F has been specified, and then the Runge-Kutta solution should be marched forward up to the other end (X ) Xmax). If the (19) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: New York, 1994.


Figure 1. ψ-p curve for pure methane at T ) 355 K.

Figure 2. Pressure predictions for the dry gas reservoir.

Ayala and Kouassi

Figure 3. Mass flux predictions for the dry gas reservoir.

Figure 4. Pressure gradient profilesdry gas reservoir.

Table 1. Reservoir and Fluid Properties, Single-Phase Case initial pressure, pi reservoir temperature, T permeability, k porosity, φ initial pseudo-pressure, ψi initial gas viscosity, µgi initial gas isothermal compressibility, cgi

24.48 MPa (3550 psia) 355 K (180 F) 1 × 10-13 m2 (100 md) 0.15 1.97 × 1019 Pa/s 2.02 × 10-5 Pa-s 3.6 × 10-8 Pa-1 (1/Pa)

value of pressure obtained at X ) Xmax agrees with the boundary condition in eq 19d (within a tolerance), the actual solution has been attained. If not, the guess for pressure at X ) Xmin must be refined. A single-variable Newton Raphson search is typically implemented to speed up the refinement process.

Discussion of Results: Benchmarking and Single-Phase Flow Conditions In order to assess the suitability of the proposed method to model single-phase gas flow, the proposed method was tested for the case of a noncondensing (dry) reservoir natural gas for which an analytical solution is available. Figure 1 displays the values of pseudo-pressure (ψ) as a function of pressure (p) for a pure methane reservoir fluid at T ) 355 K. Reservoir and initial fluid properties and listed in Table 1. Mass production flow rates were specified at Fsp ) 0.05, 0.15, and 0.25 kg/s m (qsc ) 0.431, 1.295, and 2.157 MMSCFD/ft, respectively, for the case of methane). Figures 2 and 3 reveal the pressure and mass flux numerical predictions (represented by dots) for the case of the dry gas reservoir under consideration, with pressure values converted to pseudo-pressure via Figure 1. These figures

Figure 5. Mass flux gradient profilesdry gas reservoir.

display the behavior of the system as calculated by the proposed ODE method and under the prescribed conditions. The analytical solution of the problem, based on the analytical exponential integral solution, has been included in both figures for comparison purposes. In these figures, it is clear that both methods are in excellent agreement. Figures 2 and 3 also demonstrate the expected behavior prescribed by the boundary conditions. As X f +∞, the reservoir pressure approaches initial pressure (Figure 2), and the mass flux approaches its prescribed value near at the sandface as X f -∞ (Figure 3). In addition, Figures 4 and 5 illustrate the behavior of the two major players in the proposed numerical solution: pressure and mass flux gradients along the


Figure 6. Gas-condensate fluid phase envelope.

integration domain X. Pressure gradients take a zero value at conditions of X > 2, at which reservoir pressure is constant and equal to initial pressure, which is the expected behavior at early times and/or conditions far away from the wellbore. The largest pressure gradients are found as X f -∞, in other words, for large times or conditions within the wellbore vicinity. Figure 5 indicates that the mass flux does not change appreciably throughout the domain; except for the region -2 < X < 2, which experiences the largest local mass flux variation. These figures also serve as verification of the validity of the numerical solution. Once the method has generated the results, the lefthand side and right-hand side of eqs 19a and 19b must be identical at all points of the domain. That is exactly the case in Figures 4 and 5, where the values of the left- and right-hand sides of these equations are plotted using a continuous and discontinuous line, respectively. The agreement is excellent. It is important to note the generality of the solution presented in Figures 2 and 3. These figures implicitly provide the behavior of pressure and mass flux at any particular time and at any location of this reservoir. Embedded within the ratio r/xt, a myriad of pressure and flux profiles versus the radius (at constant t) and pressure and flux profiles versus time (at a given reservoir location) can be readily constructed using these characteristic curves. This would become especially useful in the analysis of two-phase flow in gas-condensate reservoirs, where no analytical solution is available to undertake such analysis. Discussion of Results: Multiphase Flow in Gas-Condensate Reservoirs In this study, the gas-condensate reservoir fluid described by Kenyon and Behie20 is studied. Figure 6 exhibits the fluid’s phase envelope and highlights the initial reservoir conditions. Initial reservoir conditions (pi ) 24.48 MPa and T ) 355 K) lie outside the phase envelope but close to saturation conditions at the prevailing reservoir temperature. Unlike the dry gas reservoir studied in the previous section, this reservoir fluid will experience retrograde liquid dropout and multiphase flow upon fluid withdrawal once the reservoir pressure reaches dewpoint conditions (p < psat ) 23.974 MPa). In this work, and for illustration purposes, the van Genuchten function21 and a (20) Kenyon, D.; Behie, G. Third SPE Comparative Project: Gas Cycling of Retrograde Condensate Reservoirs. JPT, J. Pet. Technol. 1987, August, 981-997; SPE Paper 12278. (21) van Genuchten, M. T. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 1980, 44, 892-898.


Figure 7. Pressure predictionssgas-condensate reservoir using proposed ODE model. Table 2. Reservoir and Fluid Properties, Two-Phase Gas-Condensate Case initial pressure, pi reservoir temperature, T permeability, k porosity, φ van Genuchten parameter, λ critical So, Soc

24.48 MPa (3550 psia) 355 K (180 F) 1 × 10-13 m2 (100 md) 0.15 1 0.20

modified Corey22 relative permeability function were used for the prediction of the relative permeability values (kro and krg, respectively) needed for the study of two-phase flow. These correlations are presented below:

kro ) xSo*{1 - [1 - (So*)1/λ]λ}2

(21a)

krg ) (1 - So)2(1 - So2)

(21b)

In these correlations, So* ) (So - Soc)/(1 - Soc), Soc is the condensate critical saturation, and λ is the parameter in the van Genuchten relative permeability curves that controls the steepness of the kro-relative permeability curve. In general, the larger the value of λ, the steeper the rise of kro with increasing So for So > Soc. Any other choice of (continuous) relative permeability functions can be made with the proposed method. Table 2 shows the values of the various properties selected for this study. Figures 7-9 display the pressure, mass flux, and condensate saturation predictions, respectively, for the gas-condensate reservoir under consideration via the proposed ODE method. For comparison purposes, the numerical and analytical solution for an analogous single-phase condition is presented. In this latter case, a fictitious single-phase condition was artificially created by forcing the reservoir fluid to remain as a stable single phase throughout the simulation. Figure 7 demonstrates that much higher pressure losses are to be expected during multiphase flow conditions, especially as flow-rate specification increases and as X f -∞. Figure 8 shows that the mass flux term behavior follows the same qualitative trends studied for single-phase flow, while Figure 9 illustrates the evolution of condensate buildup throughout the domain as the production specification (Fsp) changes. It is important to stress the universal character of the similarity variable characteristic curves presented here. Once the characteristic curve is available in terms (22) Corey, A. T. The Interrelation Between Gas and Oil Relative Permeabilities. Prod. Mon. 1954, 19, 38.


Ayala and Kouassi

Figure 8. Mass flux predictionssgas-condensate reservoir using proposed ODE model.

Figure 10. CMG pressure vs radius prediction, Fsp ) 0.23 kg/s m.

Figure 9. Condensate saturation predictionssgas-condensate reservoir using proposed ODE model.


of the similarity variable, no further calculations are needed to generate the myriad of possible pressure and saturation profiles for various times and reservoir locations for the flow rates under consideration. Therefore, a single set of results is sufficient to describe the behavior of the system at any point and at any time, as long as the assumptions inherent in the model remain valid. Results: Comparison with Standard Finite-Difference Simulators The simplicity of the system of ODEs presented above enables the quick calculation of gas-condensate reservoir performance and circumvents the need of standard finitedifference simulators for the problem under consideration. In standard commercial simulators, the strongly nonlinear porous media flow equations are solved simultaneously by implementing finite-difference discretizations. In this section, a standard finite-difference commercial simulator will be used to further verify the applicability of the similarity solution during twophase flow in gas-condensate reservoirs. Figures 10-12 display the pressure-radius profile predicted by the CMG-GEM commercial simulator23 for three different specified flow rates (Fsp ) 0.23, 0.46, and 0.65 kg/m s, respectively), for which multiphase flow conditions (gas + condensate) prevail at locations where p < psat. The CMG (23) GEM: AdVanced Compositional ReserVoir Simulator - User’s Manual; Computer Modeling Group Ltd: Calgary, Canada, 2005.


simulator uses the same fluid and rock characterization described in the preceding section and implements a wellbore of finite radius placed at the center of a very large cylindrical reservoir. Since the problem admits a similarity solution, pressure distributions numerically calculated by the simulator should be invariant when plotted against the ratio r/xt (or r2/t, equivalently). Figure 13 demonstrates that the use of the similarity variable as the independent variable in the x axis indeed collapses numerical profiles onto a single characteristic line for each of the flow rates under study, and regardless of the significant



Even though the Boltzmann transformation remains valid for the cases under consideration (multiphase flow, infinite acting reservoir, homogeneous and isotropic systems), the single-phase exponential integral solution itself is no longer applicable to multiphase systems, and thus the proposed ODE model can be instrumental in the analysis of reservoir performance. Conclusions

Figure 13. CMG pressure prediction in terms of the similarity variable.

nonlinearities associated with multiphase flow. As pointed out by Doughty and Pruess,7 numerical solutions calculated using a finite grid spacing have only approximate invariance with respect to the similarity variable, and thus some small spread between the profiles should be expected. Figure 13 also displays the solution as predicted by the proposed system of ODEs, in the form of a continuous line. In this figure, both models basically predict the same behavior, while the proposed ODE model is much simpler in nature and in terms of conceptualization and implementation. It is to be expected that these two predictions would agree as long as the background assumptions of both models are compatible. For example, large discrepancies between the models can typically occur in the region closest to wellbore conditions (i.e., as X f -∞). This discrepancy can be explained because of the use of line-source (zero-radius) wellbore approximation in the proposed ODE model versus the use of an actual finite-radius well in standard commercial simulators. The implementation of skin values other than zero around the wellbore would also make these two predictions differ at X f -∞. It is important to point out that the imposed boundary conditions (constant rate line-source well and infinite reservoir, the standard assumptions of well-test analysis) are important prerequisites for the validity of the similarity theory in this problem because they are amenable to being expressed in terms of the Boltzmann variable. Other sources of discrepancies between the two models can occur in the region far away from the wellbore X f -∞ for the cases when the infinite acting assumption becomes no longer valid. In such a case, a CMG simulation would predict depletion at the reservoir external boundary while the proposed ODE model will maintain the external boundary at a constant pressure equal to initial pressure. Bøe et al.11 and Vo12 also showed numerical experiments that demonstrated that the similarity variable can be successfully used to express simulator-generated pressures and saturations in terms of r/xt for gas-condensate reservoirs. They did not attempt, however, to postulate the associated ODE model and its solution. Vo12 additionally demonstrated some conditions where the Boltzmann transformation would no be longer valid for the analysis, namely, those conditions for which the reservoir is considered heterogeneous (e.g., when wellbore skin effects are important) and/or when the influence of the boundaries does influence pressure transient behavior (e.g., small reservoirs). In this paper, it has been shown that the method of searching similarity solutions by looking for a scaling-invariance (namely, r/xt or Boltzmann transformation) can indeed be valuable in the analysis multiphase numerical data based on ODE equations.

This study has presented an ODE model that has the capability of rigorously solving the highly nonlinear governing equations for gas-condensate multiphase flow in porous media without the need of simplifications in terms of fluid phase behavior, fluid properties, rock properties, and rock/fluid interaction modeling. The simplicity of the system of ODEs presented here enables the quick calculation of reservoir performance and circumvents the need of standard finite-difference simulators for the problem under consideration. The method applies to the analysis of infinite-acting gas-condensate reservoirs as long as the reservoir is treated as having a uniform permeability. The model provides a useful tool for the quick testing of the influence of various parameters, such as relative permeability and other rock and fluid properties, on reservoir performance. Work is underway to capture the influence of non-Darcian behavior, compositional effects and proximity to critical point conditions, and the effect of relative permeability shapes on the overall performance of a gas-condensate system using the proposed model. Acknowledgment. The authors are grateful to Dr. George Moridis of the Lawrence Berkeley National Laboratory for his helpful introduction and discussion about the applicability of the similarity variable to highly nonlinear reservoir problems in radialcylindrical coordinates.

Nomenclature A ) accumulation term (kg m-3) a ) Runge-Kutta constant coefficient b ) Runge-Kutta constant coefficient c ) isothermal compressibility (Pa-1) d ) Runge-Kutta constant coefficient e ) Runge-Kutta constant coefficient F ) mass flux term (kg s-1 m-1) f, f′ ) right-hand side functions of the ODE system G ) Runge-Kutta variable coefficient k ) absolute permeability (m2) kr ) relative permeability MW ) molecular weight (kg/kg mol) p ) pressure (Pa ) q ) volumetric flow rate (m3 s-1) R ) universal gas constant (8.314 × 103 m3 Pa kg mol-1 K-1) r ) radial distance (m) S ) saturation T ) temperature (K) t ) time (s) V ) fluid velocity (m s-1) X ) numerical integration variable, X ) ln η Z ) gas compressibility factor φ ) porosity λ ) mobility (kg m-3 Pa-1 s-1) η) similarity variable, η ) r/xt µ ) viscosity (Pa s) π ) ratio circle circumference to diameter, 3.12159265(...) F ) density (kg m-3) τ ) transmissibility (s) ψ ) gas pseudo-pressure, (Pa s-1)

1600 Energy & Fuels, Vol. 21, No. 3, 2007 Subscripts g ) gas i ) initial conditions or running subscript (Runge-Kutta) j ) running subscript (Runge-Kutta) k ) running subscript (Runge-Kutta) max ) maximum

Ayala and Kouassi o ) condensate wb ) wellbore sc ) standard conditions sp ) specified EF060505W

The Similarity Theory Applied to the Analysis of ... - ACS Publications

Recommend Documents