Znd. Eng. Chem. Res. 1995,34, 3981-3994
3981
A Noncubic Equation of State for the Treatment of Hydrocarbon Fluids at Reservoir Conditions Giorgio S. Soave Eniricerche SPA, Via Maritano 26, 1-20097 S u n Donato Milanese, Italy
A modified version of the Benedict-Webb-Rubin
equation of state (EOS) has been developed. Attention was focused on near-critical conditions, where good behavior was obtained by forcing a tight fit to the experimental critical isotherms, including the critical density. The new model requires component critical constants T, and P, and acentric factor w . Generalized expressions have been developed for the calculation of the EOS parameters from the reduced temperature and the acentric factor. Pure-component vapor-pressure and density data were reproduced accurately in a wide range of conditions (reduced temperatures down to 0.4 and pressures up t o 1000 bar, both for liquids and gases, including the critical region) and w values up to 0.6. Mixtures are treated like pure compounds. Their T,, P,, and w are calculated by new mixing rules derived from those of cubic EOS’s. With these rules, experimental VLE data were reproduced accurately, with no need of binary interaction constants for alkane-alkane pairs and one constant for the other pairs. Procedures were developed for a fast and safe solution of the equation of state.
Introduction Conditions prevailing in oil and gas reservoirs imply high pressures and are very often close to the critical point of the system. The phase and volumetric behavior of fluids under such conditions can be hardly reproduced with accuracy by current calculation methods. They are generally based on cubic equations of state, whose inability to describe correctly the volumetric properties of fluids in general, and particularly in the critical field, is well-known and cannot be overcome for the intrinsic limitations imposed by their cubic form. This is seen clearly when trying to reproduce experimental critical isotherms of pure compounds over wide pressure ranges. Best, but still unsatisfactory, results are obtained with too high values of the critical compressibility factor and seriously inaccurate densities near the critical point. The aim of this work was to develop a new equation of state with marked ability to describe the phase and volumetric behavior of hydrocarbon fluids in a wide range of conditions and, in particular, in the nearcritical region. A second objective was t o give the new equation good predictive features through the estimation of its parameters from limited data (e.g., the component critical constants and the acentric factor) and the use of mixing rules with few empirical binary parameters (hopefully, no binary constant for hydrocarbon pairs and one only for pairs involving COZ,HzS, and Nz), as is common with cubic EOS’s. In addition, robust solution techniques were to be developed, to make the solution of the equation of state as easy and safe as possible and t o reduce the computation times.
The New Equation of State The starting point was the old, well-known fiveparameter Benedict-Webb-Rubin (BWR) equation of state (Benedict et al., 1940). Neglecting the temperature dependence of the parameters, it can be written as:
Z = PvIRT= 1 + B e
+ Ce2 + De5 +
Ee2(1+ F e z ) exp(-F& (1) where B, C,D , E , and F are the five, temperature- and composition-dependent, parameters.
Many modifications of the BWR equation have been proposed in the literature (e.g., Starling et al.,1972; Lee and Kesler, 1975), whose analytical form is the same or very similar to that of eq 1,differing in the parameter expressions or the mixing rules. The strong point of BWR-type equations is the presence of the exponential term, which is active at intermediate densities and improves drastically the calculated pure-compound critical isotherms with respect t o cubic EOS’s. Some residual inaccuracy is left anyway, which becomes evident when imposing the critical constraints (i.e., forcing the critical isotherm to pass through the critical point with zero slope and zero curvature). An initial consideration was that, if we replace the exponential term in eq 1by a series expansion:
+ + (C+ E)$ + De5 - EF2e6/2+
Z = 1 Be
EF3e8/3 - EFle”I8
+ ...
an “anomalous”distribution of powers of e (1,2,5, 6 , 8 , 10, ...) is evidenced. A more “regular” sequence (1,2,4, 6,8, ...) would be obtained by replacing e5with e4. The intuition that this may give better results will prove true further on. For the moment, let us write eq 1 in the more general form:
Z = PvIRT= 1+ B e + Ce2 +De” + Ee2(1+Fe2)exp(-Fe2) (2) with a general exponent n in place of the original 5. The proposed equation of state (like most) is analytic at the critical point and therefore cannot describe properly the nonclassical limiting behavior that real fluids display. This will cause slight inaccuracies near the critical point, which are of limited extent anyway (Michels and Meijer, 1984) and were considered as acceptable, t o avoid the need of a much more complicated model.
Critical Density In order t o ensure a correct behavior of eq 2 in the critical region, the first step was to force it to reproduce the experimental pure-compound critical isotherms by
0888-5885/95/2634-3981$09.Q~IQ 0 1995 American Chemical Society
3982 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995
imposing the critical constraints and also the value of the critical density. Unfortunately, critical density data are rather scarce in the literature, especially for heavy compounds, and the available predictive methods (cf. Reid et al., 1987) require the knowledge of the molecular structure or seem to cover too limited ranges of the acentric factor w (usually adopted for nonpolar substances). So, a different approach was followed, by replacing the critical compressibility factor 2, with the Rackett factor (Rackett, 1970), which can be easily determined from liquid-density data and is normally within 1%of 2,. Details about data collection and correlation are given in Appendix A. The values of calculated from the critical constants and the liquid densities by the Rackett equation (eq A-1) were correlated to w , and the following expression was obtained:
2, = 0.2908 - 0.0990
+ 0.040~
(3)
where the quadratic term was required t o ensure a good accuracy up to w = 0.8 (see Figure 1). The values of 2%estimated by eq 3 will be assumed for 2, from now on.
Critical Isotherm If the value of 2, and the critical constraints are
ZRa(exp.)-ZRa(eq.3)
0.006
Onn
0
U
O
x
b
+ + c + d + e(1 + f ) exp(-f) = Z, 1 + 2 b + 3c + (n + l)d + e ( 3 + 3f - 2 f 2 ) x
f
U
X
x
A
