1962
Ind. Eng. Chem. Res. 1994,33, 1962-1970
An Equation-of-State-BasedReservoir Model Incorporating Continuous-Mixture Thermodynamics Dwarkaprasad S. Chakravarty and Michael A. Matthews' Chemical Engineering Department, University of Wyoming, P.O. Box 3295, University Station, Laramie, Wyoming 82071
The application of continuous-mixture thermodynamics for representing the hydrocarbon phases in a compositional reservoir simulator is considered as opposed to the conventional approach which uses several discrete pseudocomponents to characterize the phases. A method of formulating the continuous material balance and phase equilibrium relationships within an implicit pressure-explicit framework is demonstrated. A one-dimensional reservoir simulator was written, saturation (IMPES) and continuous-mixture thermodynamics was successfully employed. The results obtained from the continuous and the discrete versions of the simulator are in good agreement for the test cases considered. Introduction Compositional reservoir simulation is required to characterize mass transfer during miscible gas displacement. However, most reservoir fluids are complex ensembles of hydrocarbons (continuous mixtures) containing so many components that specifying each individual component is impossible. Hence, current methods use pseudocomponents to represent the hydrocarbon phases. Each pseudocomponent is in reality a lumped version of many hydrocarbons which boil over a narrow range of temperatures. Physical properties are then assigned to these lumped pseudocomponents so that each may be treated as a discrete component. The pseudocomponent method has been well developed and documented (e.g., Fussell and Fussell, 1979; Coats, 1980; Nghiem et al., 1981; Acs etal., 1985;Collinsetal., 1986;Changetal., 1990). Though accuracy increases with the number of pseudocomponents, typically only 5-10 pseudocomponents are employed in reservoir simulation due to CPU time and storage constraints. Also, many analytical procedures force the engineer into using lumped species because a practical implementation of continuous-mixture thermodynamics for modeling stagewise equilibrium processes has been lacking. Potentially valuable thermodynamic and chemical information is sacrificed in such lumping procedures. Continuous-mixture thermodynamics describes a complex mixture by means of one or more continuous distribution functions. The concept has been developed by several authors (e.g., Gaultieri et al., 1982; Kehlen and Ratzch, 1983; Prausnitz et al., 1985; Prausnitz and Cotterman, 1985a; Behrens and Sandler, 1988; Wang and Whiting, 1988; Haynes and Matthews, 1991). The application of continuous-mixture thermodynamics for simulating phase equilibria for enhanced oil recovery has been addressed by Prausnitz and Cotterman (1985b) and more recently by Galeana and Ramirez (1990). A convenient computational procedure is inplemented by using numerical quadrature procedures to define a set of fixed components at the quadrature roots. These "quadrature" components are treated in the same manner as pseudocomponents selected by more arbitrary means. The simplest implementation of the continuous ap-
* To whom correspondence should be addressed. Present address: ChemicalEngineeringDepartment, Universityof South Carolina, Columbia, SC 29208. E-mail:
[email protected] 0888-588519412633-1962$04.50/0
proach would be to select quadrature components based on the initial distribution of the continuous fluid, and to retain these quadrature components unchanged throughout the course of a multistage equilibrium simulation. Such an approach would constitute a mathematically rigorous way of selecting pseudocomponenta, and requires no changes to existing computer simulation strategy. However, such a strategy may not utilize the full potential of the continuous-mixtureconcept. In a miscible reservoir flood,for example, it is likely that the composition distributions would change shape drastically over the course of time as components are successively stripped and removed from the reservoir. Maintaining a fixed set of pseudocomponents (or quadrature components) does not allow one to track the complete distributions and changes that occur. Matthews et al. (1991) proposed to recompute the complete composition distribution of the various equilibrium phases in a multiple contacting process. To incorporate this proposal in a reservoir simulator, however, one cannot easily replace the discrete thermodynamics formulation with the proposed continuous approach. The purpose of this work was to incorporate a continuous thermodynamics approach to address a demanding multistage problem-the compositional petroleum reservoir simulator. An equation of state simulator is needed to account for the high pressures and asymmetric nature of the fluids. The objective of this paper is to explore the use of continuous-mixturethermodynamics, including the calculation of complete composition distributions, in a compositional reservoir simulator. To achieve this, continuous-mixture mass balance equations are required to generate the fluid distributions in each reservoir cell before phase equilibrium calculations can be carried out. Rather than attempting to modify an existing commercial simulator, we developed a one-dimensional, threephase (oil, water, gas) simulator based on the implicit pressure-explicit saturation (IMPES) solution procedure of Nghiem et al. (1981). Simplifying assumptions (constant rock properties, no capillary pressure, and no gravity effects) were made so that the computed results would emphasize the effects of the thermodynamic routines. The simulator developed would be suitable for modeling onedimensional core flooding experiments. For comparison, a conventional pseudocomponent simulator was also written so that comparisons between the continuous and discrete thermodynamics strategies can be made. 0 1994 American Chemical Society
An Approach to Multistage Contacting of Continuous Phases While there have been many papers addressing phase equilibrium calculations for single-stageprocesses (bubble point, dew point, various flash and cloud point calculations), there are few detailed demonstrationsof continuousmixture thermodynamics for multistage processes. The mathematical problems of maintaining mass balances and the limitations of assuming closed form distribution functions are well-known (Luks et al., 1990; Sandler and Libby, 1991). These mathematical difficulties may be overcome, but at the expense of a number of physical/ chemical assumptions. For example, Kehlen and Ratzch (1987) proposed an approach to steady-state distillation of a continuous mixture. However, their development was valid at low pressure, with vapor-phase fugacitycoefficients independent of composition and liquid activity coefficients of unity. In addition, the method of moments was invoked for flow distribution functions, which were of a fixed (Gaussian) form. The shortcomings of the method of moments have been discussed (Luks et al., 1990; Sandler and Libby, 1991). Most previous work on continuous-mixture thermodynamics assumes that the composition distribution assumes a closed functional form (Gamma distribution, Gaussian distribution, or truncated exponential, for instance). It has been reported that no single functional form adequately represents the composition of an arbitrary petroleum (Peng et al., 1987). In addition, closed functional forms are not preserved, even for a single-stageflash calculation (Sandler and Libby, 1991). Sandler and Libby (1991) proposed to fix the form of the distribution for one of the phases and calculate the distribution of the second phase based on mass balances and the results of the phase split calculation. An alternative is to use an experimentally-determined composition distribution of arbitrary shape. In the continuous version of the reservoir simulator, the Haynes and Matthews (1991) procedure is employed. This uses cumulative weight or mole fraction as the dependent variable and true boilingpoint as the independent variable. This continuous boiling point curve is fit to a cubic spline, and no fixed functional form for the distribution is assumed. Representative components are computed using the quadrature method as outlined in any standard book on numerical analysis (Haynes and Matthews, 1991).For petroleum, a boiling point distribution is standard experimental information, and it is sufficient to represent the distribution with a cubic spline, the coefficients of which are determined from the analytical data. The boiling point curve alone is insufficient information; either specificgravity or the Watson characterization factor (Watson et al., 1933) is needed also. It is desirable to have specific gravity data on each of the boiling-point cuts established by the distillation analysis, but lacking this a single specific gravity of the whole oil can be used to calculate a Watson characterization factor for the whole oil, and the specific gravity of the individual cuts can then be calculated by assuming the Watson factor is a constant for each cut. With boiling point and specific gravity defined over the entire distribution, one uses well established petroleum correlations like Riazi-Daubert (Riazi and Daubert, 1980)andLeeKesler (Lee and Kesler, 1976) for computing critical properties and acentric factors. These correlations were developed using data for pure components and real petroleum fractions, and give critical properties and acentric factors as power-law functions of boiling point and specific gravity. Once the quadrature
1.00 7
0.00
1
Ind. Eng. Chem. Res., Vol. 33, No. 8, 1994 1963 ,
, , , , , I , , ,
200
::
, I I I I I I , ~ , , , , , t, I , , , , , , Tm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
300
400
500
600
700
800
900
True boiling point ( K )
Figure 1. Initial distribution. pi = pseudocomponent i; qj = quadrature component i.
component mole fractions and their physical properties are computed, phase equilibrium calculations are carried out in the same manner as for a set of pseudocomponent mole fractions. The procedure is applicable to any equation of state whose parameters are derived from critical properties and acentric factor. After a single vapor-liquid equilibrium calculation, the mole fractions of the most recent quadrature components in each phase are available,as well as the fraction vaporized. As shown by Matthews et al. (1991), this is sufficient information to reconstruct the true boiling point curves of both vapor and liquid phases, without prescribing a fixed functional form for either phase. In this latter respect, the procedure differs from that suggested, for example, by Sandler and Libby (1991) or by Luks et al. (19931, who specify a declining exponential curve with one adjustable parameter to represent the liquid phase. Our method is very similar to that of Hendriks (19871, who demonstrated how to reconstruct the equilibrium vapor and liquid phases when a two-parameter equation of state was adequate to describe the mixture. The procedure we employ is outlined in the Appendix.
Pseudocomponent vs Continuous-Mixture Description of Fluid In the conventional approach, the same pseudocomponents are retained throughout the simulation. Some of the pseudocomponents may become wasteful as lighter ends from the oil are stripped by the gas and heavier ends are retained in the oil phase. The present continuousmixture approach involves choosing discrete species not arbitrarily, but by doing quadrature on the continuous distribution curve such that the quadrature components are representative of the current state of the fluid (see Figures 1and 2). These figures are an indication of how the distributions at a particular location (cell) may change during the simulation. In Figure 1, four pseudocomponents are used to describe the initial distribution, and each pseudocomponent pi represents a significant weight fraction in the boiling point distribution. Figure 2 shows a different distribution, derived from the feed, that has shifted to a lighter boiling range as a result of a separation process. Figure 2 shows that now pseudocomponents p3
1964 Ind. Eng. Chem. Res., Vol. 33, No. 8, 1994
following equilibrium constraints are employed i n each cell j . 0.80
1
/
fl-fi=O
0.60
/ 2 0.40
E X i = 1.0
(4)
CYi= 1.0
(5)
The phase saturations are then computed as follows:
s;? = 1.0-soj1+1 -23,
Figure2. New distribution. pi = peeudocomponenti;qi = quadrature component i.
and p4 represent negligible fractions of the oil. The importance of these components depends on the type of phase equilibrium calculation being considered. PseudoComponentsp3 and p4 would strongly influence a dew point calculation, but would have little influence on a bubble point calculation,for example. With the present approach, however, the four quadrature components are reassigned on the basis of the new distribution and hence always represent a significant mass fraction of the overall distribution.
Method Development In this section, the pseudocomponent IMPES method for compositional reservoir simulation involving three phases, i.e., oil, gas, and water, is outlined first, followed by a development of the continuous solution method. Pseudcomponent Solution Outline. Briefly, the IMPES equations and constraints for an immobile water phase are as follows. For details, the reader is referred to Nghiem et al. (1981). The overall flow equation is
(3)
(7)
Then, the process is repeated until convergence of pressures results, upon which we move on to time step n + 1. Continuous Solution Method. We may have one or more discrete components to represent the oil along with a heavier, continuous fraction. In miscible gas displacement, the flooding species are generally lighter hydrocarbons or carbon dioxide and hence also discrete components. First of all, one needs the distillation curve and the specificgravity of the inplace oil, data which are readily obtainable for petroleum. As in the discrete case, the pressure equation is solved first. For a cell j , at any pressure iteration step, say E$++',we can find the mole fractions of the discrete components using eq 2. The mole fractions of the quadrature components representing the continuous distribution are calculated. The mole fractions (discrete and continuous) are then normalized. Now, an isothermal flash calculation is performed on this normalized set of mole fractions. The critical properties and acentric factors of the quadrature components are calculated using the Riazi-Daubert and Lee-Kesler correlations referred to earlier. Thus, we can obtain the oil and gas fractions, mole fractions of each component in the oil and gas phases, phase densities, and saturations as in the pseudocomponent case. The mole fraction of continuous oil is given by XC':
= 1.0 - X X d f ; '
The mole fraction of continuous gas is given by At a particular time step n,the flow equation is discretized using fiiite differenceswith upstream weighting and solved for pressure in each grid block. Once the grid block pressures are obtained, the mole fractions of the various pseudocomponents are then computed for each grid block from the discretized form of the component mass balance equations (eq 2).
Knowing the component mole fractions, an isothermal flash calculation is carried out in each cell to compute the oil and gas fractions, mole fractions of each component in the oil and gas phases, and the phase densities. The
The overall continuous mole fraction is given by
where the summation is over the discrete component mole fractions. From the results of the phase equilibrium calculations, we can reconstruct the true boiling point curves of the vapor and liquid streams in cell j to obtain Fuj, Fl! as per the method outlined in the Appendix (Angelos et a!.,1992; Matthews et al.,1991). Figure 3 shows the feed, liquid (oil), and vapor (gas) distributions representing a typical continuous isothermal flash. Again it is emphasized that no a priori assumption about the functional form of the distribution curves is required. The difficulty really lies in obtaining the true boiling point curve of the recombined reservoir fluid in cell j for the next iterative step ( 1 + 2) since fluids (vapor and liquid
Ind. Eng. Chem. Res., Vol. 33, No. 8,1994 1966 1.00 1
F(T) =
M;. = V&(Sg, + Sgpg)j"ZCj
0.80 1 C
.-c
