Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 74, 2019
Numéro d'article 77
Nombre de pages 13
DOI https://doi.org/10.2516/ogst/2019049
Publié en ligne 25 octobre 2019

© D.V. Nichita et al., published by IFP Energies nouvelles, 2019

Licence Creative Commons
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

List of symbols

A: EoS parameter

Ai : Component parameter in the EoS

a: EoS parameter

B: EoS parameter

Bi : Component parameter in the EoS

b: EoS parameter

Ck : Coefficients in lnKi expression

fi : Fugacity of component i in the mixture

Ki : Equilibrium ratio (= yi/xi)

kij : Binary Interaction Parameter (BIP) between components i and j

L: liquid mole fraction

m: Number of nonzero BIPs

mk : Critical amplitudes

M: Number of reduction parameters

n: Number of components

ni : Number of moles

P: Pressure

qαi : Elements of the reduction matrix

R: Universal gas constant

T: Temperature

V: Vapor mole fraction

xi : Liquid mole fraction, component i

yi : Vapor mole fraction, component i

zi : Feed mole fraction, component i

Z: Compressibility factor

Greek letters

α: Reduction parameter

β: Critical exponent

γ: Reduction parameter

φi : Fugacity coefficient of component i in a mixture

ξ = (1 − P/Pconv)0.5 :

θ: Reduction parameter

Subscripts

c: Critical

conv: Convergence

i,j: Component index

L: Liquid phase

V: Vapor phase

sat: Saturation

Superscripts

T: Transposed

*: Reference conditions

1 Introduction

Flash calculations represent a large fraction of the computer time in the simulation of many industrial processes. This is for instance the case of compositional reservoir simulations, which consist in the simulation of phase and flow behavior of oil and gas through heterogeneous porous media for reservoir engineering purposes (Ben Gharbia and Flauraud, 2019, He et al., 2019, Luo et al., 2019, Montel, 1998). In these simulations, a huge number (up to billions) of flash calculations (one for each grid block representing a homogeneous region of the reservoir) are carried out for multi-component fluids with a dozen or more components or pseudo-components, and these flash calculations must be repeated for every time step. A flash calculation is an iterative computational process for finding the equilibrium ratios Ki = yi/xi (xi and yi are the equilibrium compositions of component i in the liquid and vapor phases), starting from “initial” values of the composition ratios Ki0. The process is more rapid if the initial guess Ki0 are closer to the equilibrium ratios Ki; high-quality initial estimates are important especially for “difficult” regions of the phase envelope (critical points and convergence points, where the equilibrium ratios Ki tend to unity, see next section, as well as the immediate vicinity of phase boundaries). Some compositional reservoir simulators (Wang et al., 1997) take advantage of the so-called “negative flash”, as an alternative to phase stability testing.

Two-phase vapor-liquid phase equilibrium calculations can be initialized using the results of phase stability testing (Michelsen, 1982a) or from ideal equilibrium constants; The latter initial estimates are very poor, except at low pressures, while the former require the resolution of a nonlinear system of equations (twice if the trial phase is not known a priori to be vapor or liquid).

One route to obtain “good” values of Ki0 is to extrapolate the equilibrium ratios obtained at the previous time step (Mehra et al., 1982; Nghiem and Li, 1990; Nichita et al., 2007a; Wang and Stenby, 1994) at close pressure, temperature and composition. Another possible route, which is pursued in this paper, consists in exploiting a regularity related to the analyticity of the Equation of State (EoS) and thermodynamic potentials: upon approaching a convergence point, these ratios tend to unity according to a simple scaling law – characterized by an exponent ½ – of the temperature or pressure. The EoS of concern here are analytical, which is the case of most EoS used for engineering applications.

It appears that there is an increasing interest to implement simplified rapid phase behavior calculations procedures in compositional reservoir modeling, to achieve significant performance improvement without losing prediction accuracy. Gaganis and Varotsis (2014) proposed an integrated approach using automatically generated classification and regression models (fully replacing conventional routines in the simulator) to provide direct answers to both the phase stability and phase split problems; the dimensionality of the model can be lowered by using the reduced variables framework.

Gaganis (2018) presented a new phase stability method applicable when repeated phase behavior calculations are required, based on simple off line generated discriminating functions. The K-values based methods (interpolating equilibrium constants as functions of pressure and composition) are widely used in thermal compositional simulation (Zaydullin et al., 2014, 2016). Rannou et al. (2013) used a tie-line-based equilibrium constants method that captures the compositional dependence of the phase behavior.

The outline of the paper is as follows. After reminding of the properties of convergence points and the concept of negative flash, we present the general expressions for the asymptotic behavior of the natural logarithms of equilibrium constants when a convergence point is approached. Several examples are then considered, including synthetic and natural hydrocarbon fluids and their mixtures with injection gases, which are examined by using a general form of two-parameter cubic EoS (given in an Appendix together with the associated reduction parameters).

2 A reminder on convergence points

The concept of convergence pressure dates back to the early 1950’s when it was used in conjunction with charts of equilibrium ratios for the determination of equilibrium phase compositions (Kaliappan and Rowe, 1971; Kazemi et al., 1978; Rowe, 1967) until the development in the 1970’s of equations of state of engineering purpose and computer-assisted flash calculations. In these calculations, the “initial” Ki values required for starting the iterative computational process are an essential input, which is in many cases were obtained by using approximate methods based on the concept of convergence.

For a given multi-component mixture, the convergence pressure (or temperature) is the pressure (or temperature) where the equilibrium ratios of all mixture components converge toward unity when pressure (or temperature) is increased, the other parameter (i.e., temperature or pressure) being fixed. The critical pressure Pc or temperature Tc are one particular convergence pressure and temperature: at T = Tc the equilibrium ratios converge to unity when pressure P is increased and approaches Pc. For other temperatures T ≠ Tc, the equilibrium ratios converge to unity for a pressure that exceeds the bubble- or dew-point pressure: the convergence pressure Pconv falls in a region in the TP plane where the fluid is in the single phase state. The interval between the bubble- or dew-point pressure and Pconv corresponds to “negative flashes” (Whitson and Michelsen, 1990): the flash calculation gives a non-trivial solution (the Ki’s are different from unity) but the liquid and gas mole fractions lie outside the interval of physical solutions [0,1] (the limits of this interval correspond to the bubble- and dew-points). Likewise, for fixed P the convergence temperature Tconv is defined as the temperature where the Ki’s converge to unity; Tconv lies outside the two-phase domain except when P = Pc (and then Tconv = Tc). The locus in the T–P plane of convergence pressures and temperatures is termed the Convergence Locus (CL).

The “negative flash” domain is comprised between the CL and the phase boundaries. The mathematical domain of the vapor mole fraction V, comprised between two adjacent asymptotes of the Rachford-Rice function and , is wider (and in many cases much wider) than the physical interval V ∈ [0, 1]. In a negative flash calculation, the vapor (or liquid) mole fraction is allowed to go out of this interval. When a convergence point that is not a critical point is approached, V → ±∞, where the sign depends on which of the two asymptotes is approached (whereas V → 1/2 when a “true” critical point is approached). The main advantage of the negative flash is that all phase properties are continuously derivable at phase boundaries; this property is useful to safely treat phase appearance and disappearance at the crossing of phase boundaries.

On the CL, K-values are all equal to unity but component mole fractions zci are different from feed compositions zi, except indeed when the convergence point is a “true” critical point. However, those mole fractions zci are those of a mixture which has its “true” critical point equal to the convergence point. Moreover, the results of isothermal flashes for P < Pconv, or of isobaric flashes for T < Tconv, are very close (identical for binary mixtures) for the feed composition zi and for the “critical composition” zci (Rowe, 1967) except for the vapor and liquid mole fractions V and L = 1 − V. In the latter case, these fractions remain in the physical domain [0,1] and tend to ½ on approach to Pconv or Tconv, which are “true” critical pressure or temperature. This feature is exploited in the next section, which goes along with a very peculiar asymptotic behavior of Ki’s on approach to a convergence point and several examples will be given in Section 4.

3 Asymptotic behavior of equilibrium ratios and their logarithm near convergence points

The near-critical behavior in the two-phase liquid-vapor region was studied by Dalton and Barieau (1968) and Dalton (1970) for the equilibrium constants of binary mixtures and by Fleming and Vinatieri (1979) in ternary mixtures. It relies on an EoS or a thermodynamic potential that is analytical near the critical point, which allows it to be Taylor-expanded in powers of ΔT = T − Tc or ΔP = P − Pc, or of their dimensionless counterparts ΔT = (T − Tc)/Tc and ΔP = (P − Pc)/Pc, depending on whether pressure or temperature is varied in the process of interest. The departure from unity of the equilibrium constants turns out to vary according to a “classical” (or mean-field) scaling law as a function of the distance to the critical point, namely |ΔT| or |ΔP|, with an exponent equal to ½. This feature is directly related to the Taylor-expansion (or Landau expansion, Landau and Lifshitz, 1959) of the appropriate thermodynamic potential (in this case the Gibbs free energy) around the critical point (of coordinates Tc and Pc), and on vanishing at the critical point of the second- and third-order derivatives of the thermodynamic potential with respect to composition. A similar behavior was proved by Michelsen (1984) for multicomponent mixtures.

The natural logarithms of the equilibrium constants can be approximated as

(1)in the vicinity of the critical point along the critical isotherm, T = Tc, in the two-phase region, and

(2)on the critical isobar, P = Pc, in the two-phase region, where and depend on mixture composition, on the eigenvector corresponding to the minimum eigenvalue of the Hessian at the critical point and on third and fourth order partial derivatives of the thermodynamic potential with respect to composition and pressure or temperature. While the calculation of these derivatives is a difficult task and it requires the prior calculation of the critical point it is interesting to note that the critical amplitudes need not be directly calculated; they can be approximated using the information available from a previously performed flash calculation by assuming the ½ power law from equations (1) and (2). Note that the equilibrium constants follow the same asymptotic behavior as lnKi near a critical point (where Ki − 1 are small), but on smaller pressure intervals, since

At this point it is important to remind that the above mean-field or Landau-type approach provides only an approximation of the near-critical behavior of fluids, which is characterized by scaling exponents that differ from the mean-field (or Landau) values: for instance, the exponent ½ in equations (1) and (2) should be replaced by β ≈ 0.33 (for a review see Levelt Sengers et al., 1983). The latter behavior is not considered further in this paper, which is focused on EoS and thermodynamic potentials of engineering purpose, all of which are analytical functions of thermodynamic and composition variables.

Until now, no hypothesis has been made as to the specific analytical EoS or thermodynamic potential used for describing the phase properties of the multi-component mixture. A widely occurring circumstance is when the EoS is a two-parameter cubic EoS, such as the Peng-Robinson (PR) EoS (Peng and Robinson, 1976; Robinson and Peng, 1978), used in next section. In the simplest case where all Binary Interaction Parameters (BIPs) are equal to zero, the equilibrium ratios can be written as (Michelsen, 1986):

(3)where C0, C1 and C2 depend only on bulk quantities (EoS parameters and compressibility factors of both phases at given T, P and composition conditions, cf. Appendix) and the coefficients and Bi are related to pure component parameters (see Appendix) and therefore are non singular near the critical point. One may assume that the coefficients C0, C1 and C2 are exhibiting a similar asymptotical behavior as equilibrium ratios and they are obeying the same scaling law, i.e.,

(4)where the prefactors (critical amplitudes) mk depend on the fluid system and the EoS used. As will be seen later, the above assumption was confirmed by all examples examined in this paper (as well as in many others not reported here). This feature can be extended to two-parameter cubic EoS with non-zero BIPs (hereafter denoted kij) provided lnKi is properly decomposed (i.e., the mixing rules must be linear forms or decomposable into linear forms). One example (Nichita and Minescu, 2004) is given by the reduction parameters (k = 1, …, m; i = k + 1, …, n: the first m components have nonzero BIPs with the remaining ones) such that:

(5)where the term under summation is for the first m components having nonzero BIPs with the remaining ones. In this case there are m additional coefficients Ck, and equation (5) can be rewritten in the more compact form,

(6)where qki (, , and , for i = 1,…, n) are the elements of the reduction matrix (Hendriks, 1988), and the coefficients Ck depend only on the reduction parameters (directly and via the compressibility factors). Expressions for the coefficients Ck are given in the Appendix. Equation (6) remains valid for any reduction procedure (Hendriks and van Bergen, 1992; Nichita, 2006; Nichita and Minescu, 2004), provided appropriate elements of the reduction matrix are used. It is worth noting that equation (6) is a key equation for flash calculations using reduction methods (Nichita and Graciaa, 2011; Petitfrere and Nichita, 2015) and for pseudo-component delumping (Nichita and Leibovici, 2006). Expressions for the coefficients Ck (k > 2) are given in the Appendix.

As emphasized in the previous section, the results of an isothermal flash (i.e., the equilibrium constants) performed for the fluid composition zi and for the critical composition zic (corresponding to the current temperature) are very close, therefore equation (1) holds when Pc is replaced by Pconv and zi is replaced zic. More generally, any fluid system that is described by an analytical EoS or thermodynamic potential behaves similarly when approaching a critical (or convergence) point from the two-phase (i.e., liquid – vapor) region: the equilibrium ratios Ki tend asymptotically to unity (or, equivalently, lnKi tend to 0) according to the “universal” scaling law or , depending on whether the critical (or convergence) point is approached by varying the temperature (at constant pressure) or pressure (at constant temperature).

It must be noted that the conditions near the convergence locus are extremely difficult ones for flash calculations. From a computational point of view, a negative flash near the convergence locus is even more difficult than a regular flash near critical points, since the negative flash solution corresponds to a saddle point of the Gibbs free energy hypersurface, while at two-phase conditions the solution is at the global minimum of the Gibbs free energy G. The Successive Substitution Iterations (SSIs) scheme only guarantees the convergence to a local minimum of G (Michelsen, 1982b), and at certain conditions in the negative flash domain the SSI method converges to the trivial minimum of the G surface (a more detailed explanation can be found in Whitson and Michelsen, 1990). Thousands of successive substitution iterations may in fact be needed to ensure convergence, and more efficient methods are needed to estimate the convergence pressure without repeatedly performing flash calculations under very difficult conditions.

In the following examples, the negative flash routine and the algorithm for convergence pressure calculation presented by Nichita et al. (2007b) are used in the “exact” flash calculations. A general form of two-parameter cubic EoS is used (see Appendix), and numerical results are obtained by using the Peng–Robinson EoS. Using the reduction parameters limits the scope of the paper to cubic EoS (the most widely used in petroleum engineering). However, a similar methodology can be used for any EoS by exploiting the quasi-linearity of lnKi (on narrower intervals than in the case of reduction parameters, as will be shown in the next section).

4 Numerical procedure

Surprisingly, as will be seen below from numerical examples, little deviations from the asymptotic behavior, equations (1) and (4), are generally observed in the entire negative flash region (between the CL and the saturation curve), and also inside the two-phase region near phase boundaries where flash calculations are usually more difficult (as compared to flashes well inside the two-phase region). These equations can be exploited in various manners that are illustrated in the next section with particular fluid examples. First, it can be exploited to set up a fast algorithm (as compared to the existing ones, Jensen and Michelsen, 1990; Nichita et al., 2007b) for approximate calculations of convergence pressure. Second, this property can be used to extrapolate the results of previously performed flash calculations and provide a good initialization for a negative flash or a two-phase flash near a phase boundary. Third, it can be used to inter- and extra-polate K-values from known values. Fourth, it allows a rapid approximate determination of saturation pressures.

Another advantage of this formulation is that pre-calculated tables of K-values are not required, because by using equations (4) and (5) only very limited storage is required (essentially only the slopes mk in equation (4) must be stored for given temperature and composition conditions).

Introducing the notation

(7)the equation of a straight line between a certain reference pressure, P* (where the results of a negative or two-phase flash calculation are available), and Pconv is

(8)with

(9)where , since .

At P = Pconv, equation (8) reads

(10)The partial derivatives (∂Ck/∂P)T have rather simple expressions, and are easily calculated analytically (as described by Nichita, 2008); in fact, in practice it suffices to calculate only one of these partial derivatives, say that of C0 (note that its expression does not depend on the number of non-zero BIPs).

The convergence pressure (its approximation considering a linear variation) is obtained from equations (9) and (11) written for k = 0:

(11)

A similar equation written in terms of temperature can be used for an approximate convergence temperature calculation at given pressure.

Once Pconv is calculated, the amplitudes mk (slopes) are obtained from

(12)then the coefficients Ck can be calculated as functions of pressure from

(13)and finally the dependences Ki(P) are obtained from equation (6) for i = 1 ,…, n.

If the results of a flash calculation at a certain pressure (which can be in the two-phase region or in the negative flash region) are available, one can rapidly estimate the phase boundary location. Using K-values estimated from equation (6), saturation pressure can be rapidly approximated. Few iterations are required until the dewpoint locus equation is satisfied (on the bubble point side, the equation is ). A numerical application will be presented in the next section.

The new regularity can also be used to efficiently generate K-values tables (some compositional simulators have an option for using K-values tables for flash calculations). For domains where the coefficients Ck are non-linear, a procedure similar to that presented by Chien and Lee (1983) for K-values can be set up. Note that Ck generally have a “smoother” variation with pressure than K-values, and less storage is required.

5 Results

The linear dependence of the K-values and of Ck with , which we call hereafter quasi-linearity, is tested for a model gas-condensate system (Y8) taken from the literature (Yarborough, 1972), and for mixtures of this gas-condensate with nitrogen (Y8/N2): our purpose is to study the influence of an injection gas on convergence pressures and K-values in the negative flash region. Then, a reservoir fluid and its mixtures with different amounts of carbon dioxide are studied. The PR EoS is used in all calculations.

5.1 Y8 synthetic mixture

Test calculations have been performed first on a synthetic mixture of six normal-alkanes, known in the literature as the Y8 mixture (Yarborough, 1972). Mixture composition and component properties are given in Table 1. All BIPs are set to zero in the PR EoS. Figure 1 depicts the phase envelope and the convergence locus (calculated with the procedure described in Nichita et al., 2007b) of this mixture.

thumbnail Fig. 1

Phase envelope and convergence locus of the Y8 mixture.

Table 1

Composition and component properties of the Y8 mixture.

The calculated critical point is Tc = 293.78 K and Pc = 210.66 bar. Plots of equilibrium ratios Ki and of their logarithms lnKivs. (1 – P/Pc)0.5 on the critical isotherm (indeed at T = Tc the critical pressure is a convergence pressure) are presented in Figures 2a and 2b, respectively. As expected from a general argument (see Sect. 3, Eq. (1)), Ki tends to unity and lnKi tends to zero quasi-linearly when pressure approaches the critical pressure. Interestingly, lnKi obeys this quasi-linear behavior and it appears to be closer to linearity in a much larger pressure interval than does Ki (it can be observed in Fig. 2b for the two heaviest components that lnKi behave almost linearly, while an important curvature can be clearly observed in Fig. 2a for these Ki on the same interval). The same trends are observed along the critical isobar for the equilibrium ratios Ki (Fig. 3a) and for the logarithms lnKi (Fig. 3b) vs. (1 – T/Tc)0.5. The coefficients Ck (k = 0, 1, 2) also vary linearly with (1 − P/Pc)0.5 at T = Tc (Fig. 4a) and with at P = Pc (Fig. 4b) over a large pressure or temperature interval (comparable to that of quasi-linear behavior for lnKi): the coefficients of determination for the three coefficients are in both cases R 2 > 0.999 over the intervals represented in Figures 4a and 4b, about 10 bar and 10 K, respectively.

thumbnail Fig. 2

(a) Equilibrium ratios vs. on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. on the critical isotherm (T = 293.78 K) for the Y8 mixture.

thumbnail Fig. 3

(a) Equilibrium ratios vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture.

thumbnail Fig. 4

(a) Coefficients Ck (k = 0, 1, 2) vs. (1 − P/Pc)0.5 on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Coefficients Ck (k = 0, 1, 2) vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture.

Let us now analyze the same behavior for a temperature T ≠ Tc. On the isotherm of T = 335 K (in the region where retrograde condensation occurs), the K-values (Fig. 5) and coefficients Ck (Fig. 6) are plotted vs. pressure up to the convergence pressure. Note again the quasi-linear behavior for a pressure range of about 10 bar. At T = 335 K, the saturation (dew point) pressure is Psat = 224.39 bar, and the convergence pressure is Pconv = 229.09 bar. At Psat, the value of the first coefficient is C0 = 0.23876, and its derivative with respect to pressure is . From equation (11), taking P* = Psat, an excellent approximation of the convergence pressure is obtained: Pconv = 229.16 bar (0.07 bar absolute error; the relative error is about 0.03%). Note that similar results are obtained using at any pressure between Psat and Pconv, and for a certain pressure interval inside the two-phase domain. Convergence pressures have been calculated using this approximation scheme for the entire temperature range up to the cricondentherm (where the difference between convergence and saturation pressure is of the order of 100 bar); the results turn out to be very close to the exact ones (at the scale of Fig. 1, the approximated CL turns out to be undistinguishable from the exact CL).

thumbnail Fig. 5

K-values vs. pressure for the Y8 mixture at T = 335 K.

thumbnail Fig. 6

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 mixture at T = 335 K.

Suppose a flash calculation is needed at T = 335 K and P = 215 bar. One can use as initial guess the K-values from a phase stability calculation, or those obtained by an extrapolation technique based on previous flash calculations at some reference conditions. Table 2 lists the exact K-values (flash results), those obtained in this work, and those obtained using two extrapolation techniques: the Approximate Flash Calculations (AFC) proposed by Nghiem and Li (1990) and Direct Flash Calculations (DFC) of Wang and Stenby (1994) (earlier suggested by Mehra et al., 1982; it appears to be the most widely used), for a reference pressure P* = Psat. The K-values obtained by phase stability testing are also listed, as well as those at the dew point (which give a reasonable initialization for these conditions) and those calculated from the widely-used Wilson (1969) ideal K-values relationship (which are poor initial estimates for these conditions). It should be noted that AFC is an iterative method (it solves a non-linear system with the Jacobian matrix available from the reference conditions; for this specific calculation it requires six iterations) and DFC requires the resolution of a linear system of equations (also using the Jacobian matrix at the reference conditions); phase stability testing requires six Newton iterations using reduced variables (Nichita and Petitfrere, 2013). The third column in Table 1 gives the K-values calculated from equations (11), (12) and (6), using the approximate value of the convergence pressure calculated above. Obviously, the method proposed here provides the best initial guess for K-values; the relative errors in K-values for each component are given for comparison in Figure 7. In fact, excellent approximations of K-values can be obtained for a wide pressure interval.

thumbnail Fig. 7

Relative errors on K-values for different (various) initializations. Y8 mixture, flash at T = 335 K and P = 215 bar.

Table 2

Exact flash results and different initial guess for K-values for Y8 mixture at T = 335 K and P = 215 bar (the reference pressure is P* = Psat = 224.39 bar).

In practice the difference between current and reference pressure are smaller (in compositional reservoir simulators, the pressure step is bounded via the time step restriction imposed by the stability requirements of the iterative procedure for solving the non-linear system of flow equations); the above example was chosen to enlighten the capabilities of the proposed method.

Suppose now we have the results of a flash calculation at a certain pressure (which can be in the two-phase region or in the negative flash region) and we want a rapid estimate of the phase boundary location. For example, using P* = 215 bar (two-phase flash), we obtain Psat = 224.36 bar, which is very close to the exact one.

Let us now look at what’s happening at isobaric conditions. The plot of K-values (Fig. 8) and coefficients Ck (Fig. 9) vs. at P = 100 bar shows a similar behavior (with Tconv = 554.32 K and Tsat = 435 K). Note that in this case the linearity holds for hundreds of K. A similar method can be readily set up to approximate the convergence temperature at given pressure.

thumbnail Fig. 8

K-values vs. temperature for the Y8 mixture at P = 100 bar.

thumbnail Fig. 9

Coefficients Ckvs. (1 − T/Tconv)0.5 for the Y8 mixture at P = 100 bar.

5.2. Y8/Nitrogen mixtures

In order to investigate the case of non-zero BIPs, different amounts of nitrogen are combined with the Y8 mixture. The BIPs between N2 and hydrocarbon components are taken equal to 0.1 (m = 1). For a mixture made up of the Y8 fluid + 25% moles N2, the phase envelope and the convergence locus are plotted in Figure 10. The coefficients Ck (Fig. 11) are plotted vs. pressure at T = 300 K. Note that in this case we have four coefficients Ck, and their linearity holds over the entire negative flash domain. This observation is important because in the description of reservoir fluids by cubic EoS non-zero BIPs are usually required. The results of CL calculations with the exact (Nichita et al., 2007b) and approximate (this work) methods are very close. For a mixture made up of the Y8 fluid + 50% moles N2, the phase envelope and the convergence locus are plotted in Figure 12 and the coefficients Ck are plotted vs. pressure at T = 300 K in Figure 13. Note that this mixture has no critical points.

thumbnail Fig. 10

Phase envelope and convergence locus of Y8 + 25% N2 mixture.

thumbnail Fig. 11

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 + 25% N2 mixture at T = 300 K.

thumbnail Fig. 12

Phase envelope and convergence locus of Y8 + 50% N2 mixture.

thumbnail Fig. 13

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 + 50% N2 mixture at T = 300 K.

5.3. Reservoir fluid

The reservoir fluid composition is described by 29 components; eight components have non-zero BIPs with the remaining ones (m = 8). Mixture composition and component properties are given in Table 3 and the non-zero BIPs are listed in Table 4 (only elements in the lower triangular part of the symmetric BIPs matrix are given). The PR EoS is used. The phase envelope, the spinodal, and the convergence locus of the reservoir fluid are plotted in Figure 14.

thumbnail Fig. 14

Phase envelope and convergence locus of the reservoir fluid.

Table 3

Reservoir fluid composition and component properties.

Table 4

Non-zero BIPS for the reservoir fluid.

There are 18 (2m + 2) reduction parameters, and 11 coefficients Ck in equation (6). Figure 15 depicts the coefficients Ckvs. at T = 500 K; at this temperature, there is a wide negative flash region (the dew point pressure is 352 bar, and the convergence pressure is 379 bar). As expected, quasi-linearity holds for all coefficients Ck over a large pressure interval.

thumbnail Fig. 15

Coefficients Ckvs. (1 − P/Pconv)0.5 for the reservoir fluid at T = 500 K.

5.4. CO2/reservoir fluid mixture

In order to study the influence of an injection gas on the new regularity, CO2 is added in various amounts to the above reservoir fluid. The phase envelope and the convergence locus for the reservoir fluid + 50% moles CO2 mixture are presented in Figure 16. Figure 17 plots the coefficients Ckvs. for the reservoir fluid + 50% moles CO2 mixture at T = 500 K; some of the C’s (i.e., C0, C1, C2) are somewhat deviating from linearity, while for the C-parameters corresponding to small BIPs deviate from linearity negligibly over a wide pressure range (here to the order of hundreds bar, as the dew point pressure is 398 bar, and the convergence pressure is 628 bar). This example indicates that it is preferable to assume linearity for the latter C-parameters.

thumbnail Fig. 16

Phase envelope and convergence locus for the reservoir fluid + 50% CO2 mixture.

thumbnail Fig. 17

Coefficients Ckvs. (1 − P/Pconv)0.5 for the reservoir fluid + 50% CO2 mixture at T = 500 K.

6 Conclusion

The asymptotic behavior of the equilibrium ratios and their logarithm is examined near critical points and convergence points. For any analytical EoS, the K-values tend towards unity or, equivalently, their logarithms tend to 0, according to a power ½ of the distance to the critical point or convergence point.

Since the elements of the reduction matrix are never singular at critical points, the coefficients Ck for the multilinear expression of the logarithms of K-values also behave according to the same power law, which turn out to holds in a larger interval, as observed in a series of fluid examples using the conventional PR EoS. The coefficients Ck are frequently observed to behave according to this power law over the entire negative flash region, and a region inside the two-phase pressure/temperature domain (except for systems containing a hydrocarbon mixture and a high concentration of a non-hydrocarbon component, such as carbon dioxide or nitrogen).

This new regularity can be exploited in several ways; it is shown how to rapidly calculate approximate convergence pressures (avoiding costly flash calculations in the vicinity of the convergence locus), extrapolate the results of previously performed flashes to obtain a high quality initial guess for flash calculations, interpolate in rapidly generated K-values tables, rapidly approximate phase boundary location. Several numerical examples are given using a cubic EoS and fluid systems (synthetic and naturally occurring) taken from the literature, showing the extrapolating capabilities of the proposed method.

Acknowledgments

We thank Total SA for partial financial support and for the permission to publish this work.

Appendix

Cubic Equation of State and expressions of the coefficients Ck

The general form of two-parameter cubic EoS is

(A1)Equation (A1) includes the SRK EoS (Soave, 1972, for δ1 = 1 and δ2 = 0) and the PR EoS (Peng and Robinson, 1976, for ). With Z = pv/RT, A = ap/(RT)2, and B = bp/RT, the implicit form (in the compressibility factor Z) of the EoS is obtained.

The van der Waals one-fluid mixing rules are used for the energy, A and for the volume, B parameters of the EoS

(A2)

(A3)where and . Ωa, Ωb and m(ω) are specific to each EoS.

The logarithm of equilibrium constants is

(A4)where

(A5a)

(A5b)

(A5c)

(A5d)and

(A5e)where

(A6)

(A7)

(A8)

(A9)

(A10)

(A11)and .

Finally, the elements of the reduction matrix in equation (6) are: , , and , for i = 1, …, n.

References

  • Ben Gharbia I., Flauraud E. (2019) Study of compositional multiphase flow formulation using complementarity conditions, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 43. [CrossRef] [Google Scholar]
  • Chien M.C., Lee S.T. (1983) A new equilibrium coefficient correlation method for compositional simulators, SPE-12243-MS, Reservoir Simulation Symposium, November 15–18, San Francisco, CA. [Google Scholar]
  • Dalton B.J. (1970) Two-phase equilibria of analytical binary solutions near the critical point, U.S. Bureau of Mines, Information Circular 8486. [Google Scholar]
  • Dalton B.J., Barieau R.E. (1968) Equations for calculating various thermodynamic functions of a two-component system from an empirical equation of state, including liquid-vapor equilibria data, Technical Report BM-RI-7076, U.S. Bureau of Mines, Amarillo, TX. [Google Scholar]
  • Fleming P.D. III, Vinatieri J.E. (1979) Quantitative interpretation of phase volume behavior of multicomponent systems near critical points, AIChE J. 25, 493–502. [Google Scholar]
  • Gaganis V. (2018) Rapid phase stability calculations in fluid flow simulation using simple discriminating functions, Comput. Chem. Eng. 108, 112–127. [Google Scholar]
  • Gaganis V., Varotsis N. (2014) An integrated approach for rapid phase behavior calculations in compositional modeling, J. Petrol. Sci. Eng. 118, 74–87. [CrossRef] [Google Scholar]
  • He C., Mu L., Xu A., Zhao L., He J., Zhang A., Shan F., Luo E. (2019) Phase behavior and miscible mechanism in the displacement of crude oil with associated sour gas, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 54. [CrossRef] [Google Scholar]
  • Hendriks E.M. (1988) Reduction theorem for phase equilibrium problems, Ind. Eng. Chem. Res. 27, 1728–1732. [Google Scholar]
  • Hendriks E.M., van Bergen A.R.D. (1992) Application of a reduction method to phase equilibria calculations, Fluid Phase Equilib. 74, 17–34. [Google Scholar]
  • Jensen F., Michelsen M.L. (1990) Calculation of first contact and multiple contact minimum miscibility pressure, In Situ 14, 1–17. [Google Scholar]
  • Kaliappan C.S., Rowe A.M. (1971) Calculation of pressure-temperature phase envelopes of multicomponent systems, Soc. Petrol. Eng. J. 11, 243–251. [CrossRef] [Google Scholar]
  • Kazemi H., Vestal C.R., Shank D.G. (1978) An efficient multicomponent numerical simulator, Soc. Petr. Eng. J. 18, 355–368. [CrossRef] [Google Scholar]
  • Landau L., Lifshitz E. (1959) Fluid mechanics, Pergamon, Section 64. [Google Scholar]
  • Levelt Sengers J.M.H., Morrison G., Chang R.F. (1983) Critical behavior in fluids and fluid mixtures, Fluid Phase Equilib. 14, 19–44. [Google Scholar]
  • Luo E., Fan Z., Hu Y., Zhao L., Wang J. (2019) An evaluation on mechanisms of miscibility development in acid gas injection for volatile oil reservoirs, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 59. [CrossRef] [Google Scholar]
  • Mehra R.K., Heidemann R.A., Aziz K. (1982) Computation of multiphase equilibrium for compositional simulation, Soc. Petrol. Eng. J 22, 61–68. [CrossRef] [Google Scholar]
  • Michelsen M.L. (1982a) The isothermal flash problem. Part I. Stability, Fluid Phase Equilib. 9, 1–19. [Google Scholar]
  • Michelsen M.L. (1982b) The isothermal flash problem. Part II. Phase split calculation, Fluid Phase Equilib. 9, 21–40. [Google Scholar]
  • Michelsen M.L. (1984) Calculation of critical points and phase boundaries in the critical region, Fluid Phase Equilib. 16, 57–76. [Google Scholar]
  • Michelsen M.L. (1986) Simplified flash calculations for cubic equations of state, Ind. Eng. Chem. Proc. Des. Dev. 25, 84–188. [CrossRef] [Google Scholar]
  • Montel F. (1998) New tools for oil and gas reservoir fluid management, Revue de l’Institut Français du Pétrole 53, 9–11. [Google Scholar]
  • Nghiem L.X., Li Y.K. (1990) Approximate flash calculations for equation-of-state compositional models, Soc. Petrol. Eng. Res. Eng. 5, 107–114. [Google Scholar]
  • Nichita D.V. (2006) a reduction method for phase equilibrium calculation with cubic equations of state, Braz. J. Chem. Eng 23, 427–434. [CrossRef] [Google Scholar]
  • Nichita D.V. (2008) Phase envelope construction for mixtures with many components, Energy Fuels 22, 488–495. [Google Scholar]
  • Nichita D.V., Graciaa A. (2011) A new reduction method for phase equilibrium calculations, Fluid Phase Equilib. 302, 226–233. [Google Scholar]
  • Nichita D.V., Leibovici C.F. (2006) An analytical-component delumping procedure for equations of state with non-zero binary interaction consistent pseudoparameters, Fluid Phase Equilib. 245, 71–82. [Google Scholar]
  • Nichita D.V., Minescu F. (2004) Efficient phase equilibrium calculations in a reduced flash context, Can. J. Chem. Eng. 82, 1225–1238. [Google Scholar]
  • Nichita D.V., Petitfrere M. (2013) Phase stability analysis using a reduction method, Fluid Phase Equilib. 358, 27–39. [Google Scholar]
  • Nichita D.V., Broseta D., de Hemptinne J.-C., Lachet V. (2007a) Efficient phase equilibrium calculation for compositional simulation: the direct reduced flash, Petrol. Sci. Technol. 25, 315–342. [CrossRef] [Google Scholar]
  • Nichita D.V., Broseta D., Montel F. (2007b) Calculation of convergence pressure/temperature and stability test limit loci of mixtures with cubic equations of state, Fluid Phase Equilib. 261, 176–184. [Google Scholar]
  • Peng D.Y., Robinson D.B. (1976) A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15, 59–64. [CrossRef] [Google Scholar]
  • Petitfrere M., Nichita D.V. (2015) Multiphase equilibrium calculations using a reduction method, Fluid Phase Equilib. 401, 110–126. [Google Scholar]
  • Rannou G., Voskov D., Tchelepi H. (2013) Tie-line-based K-value method for compositional simulation, SPE J. 18, 1112–1122. [CrossRef] [Google Scholar]
  • Robinson D.B., Peng D.Y. (1978) The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs, Research Report, Gas Processors Association, Tulsa, Okla, RR-28. [Google Scholar]
  • Rowe A.M. (1967) The critical composition method – a new convergence pressure method, Soc. Petrol. Eng. J. 7, 54–60. [CrossRef] [Google Scholar]
  • Soave G. (1972) Equilibrium constants from a modified Redlich–Kwong equation of state, Chem. Eng. Sci. 27, 1197–1203. [Google Scholar]
  • Wang P., Stenby E.H. (1994) Non-iterative flash calculation algorithm in compositional reservoir simulation, Fluid Phase Equilib. 95, 93–108. [Google Scholar]
  • Wang P., Yotov I., Wheeler M., Arbogast T., Dawson C., Parashar M., Sepehrnoori K. (1997) A new generation EOS compositional simulator: Part I – Formulation and Discretization, SPE-37979-MS SPE Reservoir Simulation Symposium, 8–11 June, Dallas, Texas. [Google Scholar]
  • Whitson C.H., Michelsen M.L. (1990) The negative flash, Fluid Phase Equilib. 53, 51–72. [Google Scholar]
  • Wilson G. (1969) A modified Redlich-Kwong equation of state, application to general physical data calculations, Paper no. 15C presented at the AIChE 65th National Meeting, May 4–7, Cleveland, Ohio. [Google Scholar]
  • Yarborough L. (1972) Vapor-liquid equilibrium data for multicomponent mixtures containing hydrocarbon and non-hydrocarbon components, J. Chem. Eng. Data 17, 129–133. [Google Scholar]
  • Zaydullin R., Voskov D., James S.C., Lucia A. (2014) Fully compositional and thermal reservoir simulation, Comput. Chem. Eng. 63, 51–65. [Google Scholar]
  • Zaydullin R., Voskov D., Tchelepi H. (2016) Comparison of EoS-based and K-values-based methods for three-phase thermal simulation, Transp. Porous Media 116, 663–686. [Google Scholar]

All Tables

Table 1

Composition and component properties of the Y8 mixture.

Table 2

Exact flash results and different initial guess for K-values for Y8 mixture at T = 335 K and P = 215 bar (the reference pressure is P* = Psat = 224.39 bar).

Table 3

Reservoir fluid composition and component properties.

Table 4

Non-zero BIPS for the reservoir fluid.

All Figures

thumbnail Fig. 1

Phase envelope and convergence locus of the Y8 mixture.

In the text
thumbnail Fig. 2

(a) Equilibrium ratios vs. on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. on the critical isotherm (T = 293.78 K) for the Y8 mixture.

In the text
thumbnail Fig. 3

(a) Equilibrium ratios vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture.

In the text
thumbnail Fig. 4

(a) Coefficients Ck (k = 0, 1, 2) vs. (1 − P/Pc)0.5 on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Coefficients Ck (k = 0, 1, 2) vs. (1 − T/Tc)0.5 on the critical isobar (P = 210.66 bar) for the Y8 mixture.

In the text
thumbnail Fig. 5

K-values vs. pressure for the Y8 mixture at T = 335 K.

In the text
thumbnail Fig. 6

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 mixture at T = 335 K.

In the text
thumbnail Fig. 7

Relative errors on K-values for different (various) initializations. Y8 mixture, flash at T = 335 K and P = 215 bar.

In the text
thumbnail Fig. 8

K-values vs. temperature for the Y8 mixture at P = 100 bar.

In the text
thumbnail Fig. 9

Coefficients Ckvs. (1 − T/Tconv)0.5 for the Y8 mixture at P = 100 bar.

In the text
thumbnail Fig. 10

Phase envelope and convergence locus of Y8 + 25% N2 mixture.

In the text
thumbnail Fig. 11

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 + 25% N2 mixture at T = 300 K.

In the text
thumbnail Fig. 12

Phase envelope and convergence locus of Y8 + 50% N2 mixture.

In the text
thumbnail Fig. 13

Coefficients Ckvs. (1 − P/Pconv)0.5 for the Y8 + 50% N2 mixture at T = 300 K.

In the text
thumbnail Fig. 14

Phase envelope and convergence locus of the reservoir fluid.

In the text
thumbnail Fig. 15

Coefficients Ckvs. (1 − P/Pconv)0.5 for the reservoir fluid at T = 500 K.

In the text
thumbnail Fig. 16

Phase envelope and convergence locus for the reservoir fluid + 50% CO2 mixture.

In the text
thumbnail Fig. 17

Coefficients Ckvs. (1 − P/Pconv)0.5 for the reservoir fluid + 50% CO2 mixture at T = 500 K.

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.