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

List of symbols

A: EoS parameter

A_i : Component parameter in the EoS

a: EoS parameter

B: EoS parameter

B_i : Component parameter in the EoS

b: EoS parameter

C_k : Coefficients in lnK_i expression

f_i : Fugacity of component i in the mixture

K_i : Equilibrium ratio (= y_i/x_i)

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

L: liquid mole fraction

m: Number of nonzero BIPs

m_k : Critical amplitudes

M: Number of reduction parameters

n: Number of components

n_i : Number of moles

P: Pressure

q_αi : Elements of the reduction matrix

R: Universal gas constant

T: Temperature

V: Vapor mole fraction

x_i : Liquid mole fraction, component i

y_i : Vapor mole fraction, component i

z_i : 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/P_conv)^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 K_i = y_i/x_i (x_i and y_i are the equilibrium compositions of component i in the liquid and vapor phases), starting from “initial” values of the composition ratios K_i0. The process is more rapid if the initial guess K_i0 are closer to the equilibrium ratios K_i; high-quality initial estimates are important especially for “difficult” regions of the phase envelope (critical points and convergence points, where the equilibrium ratios K_i 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 K_i0 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” K_i 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 P_c or temperature T_c are one particular convergence pressure and temperature: at T = T_c the equilibrium ratios converge to unity when pressure P is increased and approaches P_c. For other temperatures T ≠ T_c, the equilibrium ratios converge to unity for a pressure that exceeds the bubble- or dew-point pressure: the convergence pressure P_conv falls in a region in the T–P plane where the fluid is in the single phase state. The interval between the bubble- or dew-point pressure and P_conv corresponds to “negative flashes” (Whitson and Michelsen, 1990): the flash calculation gives a non-trivial solution (the K_i’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 T_conv is defined as the temperature where the K_i’s converge to unity; T_conv lies outside the two-phase domain except when P = P_c (and then T_conv = T_c). 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 $1/\left[1-\max \left(K_i\right)\right]$ and $1/\left[1-\min \left(K_i\right)\right]$, 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 z_ci are different from feed compositions z_i, except indeed when the convergence point is a “true” critical point. However, those mole fractions z_ci are those of a mixture which has its “true” critical point equal to the convergence point. Moreover, the results of isothermal flashes for P < P_conv, or of isobaric flashes for T < T_conv, are very close (identical for binary mixtures) for the feed composition z_i and for the “critical composition” z_ci (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 P_conv or T_conv, 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 K_i’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 − T_c or ΔP = P − P_c, or of their dimensionless counterparts ΔT = (T − T_c)/T_c and ΔP = (P − P_c)/P_c, 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 T_c and P_c), 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

$$\ln K_i\approx A_{\mathrm{iP}}{\left(\frac{P_c-P}{P_c}\right)}^{1/2}$$(1)in the vicinity of the critical point along the critical isotherm, T = T_c, in the two-phase region, and

$$\ln K_i\approx A_{\mathrm{iT}}{\left(\frac{T_c-T}{T_c}\right)}^{1/2}$$(2)on the critical isobar, P = P_c, in the two-phase region, where $A_{\mathrm{iP}}=f\left(D_P,z_i,u_i\right)$ and $A_{\mathrm{iT}}=f\left(D_T,z_i,u_i\right)$ 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 lnK_i near a critical point (where K_i − 1 are small), but on smaller pressure intervals, since

$$\ln K_i=\ln \left[1+\left(K_i-1\right)\right]=K_i-1-1/2(K_i-1{)}^2+\dots$$

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):

$$\ln K_i=C_0+C_1\sqrt{A_i}+C_2{B}_i,$$(3)where C₀, C₁ and C₂ 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 $\sqrt{A_i}$ and B_i are related to pure component parameters (see Appendix) and therefore are non singular near the critical point. One may assume that the coefficients C₀, C₁ and C₂ are exhibiting a similar asymptotical behavior as equilibrium ratios and they are obeying the same scaling law, i.e.,

$$C_k=m_k(P_c-P{)}^{1/2};\hspace{1em}k=\mathrm{0,1},2,$$(4)where the prefactors (critical amplitudes) m_k 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 k_ij) provided lnK_i 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 ${\gamma }_{\mathrm{ki}}=\sqrt{A_i\left(k_{\mathrm{ki}}\right)}$ (k = 1, …, m; i = k + 1, …, n: the first m components have nonzero BIPs with the remaining ones) such that:

$$\ln K_i=C_0+C_1\sqrt{A_i}+C_2{B}_i+\sum _{k=3}^{m+2}{C}_k{\gamma }_{\mathrm{ki}},$$(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 C_k, and equation (5) can be rewritten in the more compact form,

$$\ln K_i=\sum _{k=0}^{m+2}{C}_k{q}_{\mathrm{ki}};\hspace{1em}i=1,...,n,$$(6)where q_ki ($q_{0i}=1,q_{1i}=\sqrt{A_i}$, $q_{2i}=B_i$, and $q_{\mathrm{ki}}=\sqrt{A_i}(1-k_{\mathrm{ki}});k=3,...,m+2$, for i = 1,…, n) are the elements of the reduction matrix (Hendriks, 1988), and the coefficients C_k depend only on the reduction parameters (directly and via the compressibility factors). Expressions for the coefficients C_k 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 C_k (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 z_i and for the critical composition z_ic (corresponding to the current temperature) are very close, therefore equation (1) holds when P_c is replaced by P_conv and z_i is replaced z_ic. 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 K_i tend asymptotically to unity (or, equivalently, lnK_i tend to 0) according to the “universal” scaling law $(T_{\mathrm{conv}}-T{)}^{1/2}$ or $(P_{\mathrm{conv}}-P{)}^{1/2}$, 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 lnK_i (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 m_k in equation (4) must be stored for given temperature and composition conditions).

Introducing the notation

$$\xi \left(p\right)={\left(\frac{P_{\mathrm{conv}}-P}{P_{\mathrm{conv}}}\right)}^{0.5},$$(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 P_conv is

$$C_k=C_k\left(P\mathrm{*}\right)+{{\left(\frac{\mathrm{\partial }{C}_k}{\mathrm{\partial }\xi }\right)}_T|}_{P\mathrm{*}}\left[\xi \left(P\right)-\xi \left(P\mathrm{*}\right)\right],$$(8)with

$$\frac{\mathrm{\partial }{C}_k}{\mathrm{\partial }\xi }=\frac{\mathrm{\partial }{C}_k}{\mathrm{\partial }P}\frac{\mathrm{d}P}{\mathrm{d}\xi }\hspace{1em}k=0 ,\dots ,\hspace{0.5em}m+2,$$(9)where $\mathrm{d}P/\mathrm{d}\xi =-2\xi P_{\mathrm{conv}}$, since $P=P_{\mathrm{conv}}(1-{\xi }^2)$.

At P = P_conv, equation (8) reads

$$0=C_k\left(P\mathrm{*}\right)-{{\left(\frac{\mathrm{\partial }{C}_k}{\mathrm{\partial }\xi }\right)}_T|}_{P\mathrm{*}}\xi \left(P\mathrm{*}\right).$$(10)The partial derivatives (∂C_k/∂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 C₀ (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:

$$P_{\mathrm{conv}}=P\mathrm{*}-\frac{1}{2}\frac{C_0\left(P\mathrm{*}\right)}{{{\left(\frac{\mathrm{\partial }{C}_0}{\mathrm{\partial }P}\right)}_T|}_{P\mathrm{*}}}.$$(11)

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

Once P_conv is calculated, the amplitudes m_k (slopes) are obtained from

$$m_k=C_k\left(P\mathrm{*}\right)/\xi \left(P\mathrm{*}\right),$$(12)then the coefficients C_k can be calculated as functions of pressure from

$$C_k\left(P\right)=m_k\xi \left(P\right);\hspace{1em}k=0,...,m+2;$$(13)and finally the dependences K_i(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 $\left(\sum _{i=1}^n{z}_i/K_i=1\right)$ is satisfied (on the bubble point side, the equation is $\sum _{i=1}^n{z}_i{K}_i=1$). 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 C_k are non-linear, a procedure similar to that presented by Chien and Lee (1983) for K-values can be set up. Note that C_k 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 C_k with $(P_c-P{)}^{1/2}$, 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/N₂): 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.

The calculated critical point is T_c = 293.78 K and P_c = 210.66 bar. Plots of equilibrium ratios K_i and of their logarithms lnK_ivs. (1 – P/P_c)^0.5 on the critical isotherm (indeed at T = T_c 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)), K_i tends to unity and lnK_i tends to zero quasi-linearly when pressure approaches the critical pressure. Interestingly, lnK_i obeys this quasi-linear behavior and it appears to be closer to linearity in a much larger pressure interval than does K_i (it can be observed in Fig. 2b for the two heaviest components that lnK_i behave almost linearly, while an important curvature can be clearly observed in Fig. 2a for these K_i on the same interval). The same trends are observed along the critical isobar for the equilibrium ratios K_i (Fig. 3a) and for the logarithms lnK_i (Fig. 3b) vs. (1 – T/T_c)^0.5. The coefficients C_k (k = 0, 1, 2) also vary linearly with (1 − P/P_c)^0.5 at T = T_c (Fig. 4a) and with $(1-T/T_c{)}^{0.5}$ at P = P_c (Fig. 4b) over a large pressure or temperature interval (comparable to that of quasi-linear behavior for lnK_i): the coefficients of determination for the three coefficients are in both cases R ² > 0.999 over the intervals represented in Figures 4a and 4b, about 10 bar and 10 K, respectively.

Let us now analyze the same behavior for a temperature T ≠ T_c. On the isotherm of T = 335 K (in the region where retrograde condensation occurs), the K-values (Fig. 5) and coefficients C_k (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 P_sat = 224.39 bar, and the convergence pressure is P_conv = 229.09 bar. At P_sat, the value of the first coefficient is C₀ = 0.23876, and its derivative with respect to pressure is ${\left(\mathrm{\partial }{C}_0/\mathrm{\partial }P\right)}_T^{\mathrm{sat}}=-0.025$. From equation (11), taking P* = P_sat, an excellent approximation of the convergence pressure is obtained: P_conv = 229.16 bar (0.07 bar absolute error; the relative error is about 0.03%). Note that similar results are obtained using ${\left(\mathrm{\partial }{C}_0/\mathrm{\partial }P\right)}_T$ at any pressure between P_sat and P_conv, 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).

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* = P_sat. 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.

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 P_sat = 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 C_k (Fig. 9) vs. $(1-T/T_{\mathrm{conv}}{)}^{0.5}$ at P = 100 bar shows a similar behavior (with T_conv = 554.32 K and T_sat = 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.

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 N₂ and hydrocarbon components are taken equal to 0.1 (m = 1). For a mixture made up of the Y8 fluid + 25% moles N₂, the phase envelope and the convergence locus are plotted in Figure 10. The coefficients C_k (Fig. 11) are plotted vs. pressure at T = 300 K. Note that in this case we have four coefficients C_k, 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 N₂, the phase envelope and the convergence locus are plotted in Figure 12 and the coefficients C_k are plotted vs. pressure at T = 300 K in Figure 13. Note that this mixture has no critical points.

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.

There are 18 (2m + 2) reduction parameters, and 11 coefficients C_k in equation (6). Figure 15 depicts the coefficients C_kvs. $(1-P/P_{\mathrm{conv}}{)}^{0.5}$ 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 C_k over a large pressure interval.

5.4. CO₂/reservoir fluid mixture

In order to study the influence of an injection gas on the new regularity, CO₂ is added in various amounts to the above reservoir fluid. The phase envelope and the convergence locus for the reservoir fluid + 50% moles CO₂ mixture are presented in Figure 16. Figure 17 plots the coefficients C_kvs. $(1-P/P_{\mathrm{conv}}{)}^{0.5}$ for the reservoir fluid + 50% moles CO₂ mixture at T = 500 K; some of the C’s (i.e., C₀, C₁, C₂) 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.

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 C_k 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 C_k 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 C_k

The general form of two-parameter cubic EoS is

$$p\left(T,v\right)=\frac{\mathrm{RT}}{v-b}-\frac{a\left(T\right)}{\left(v+{\delta }_{1}b\right)\left(v+{\delta }_{2}b\right)}.$$(A1)Equation (A1) includes the SRK EoS (Soave, 1972, for δ₁ = 1 and δ₂ = 0) and the PR EoS (Peng and Robinson, 1976, for ${\mathit{\delta }}_{\mathbf{1,2}}=\mathbf{1}\pm \sqrt{\mathbf{2}}$). With Z = pv/RT, A = ap/(RT)², 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

$$A=\sum _{i=1}^n\sum _{j=1}^n{x}_i{x}_j\sqrt{A_i}\sqrt{A_j}(1-k_{\mathrm{ij}}),$$(A2)

$$B=\sum _{i=1}^n{x}_i{B}_i,$$(A3)where $A_i={\mathrm{\Omega }}_{\mathrm{ai}}{p}_{\mathrm{ri}}/T_{\mathrm{ri}}^2{\left[1+m_i\left({\omega }_i\right)\left(1-T_{\mathrm{ri}}^{0.5}\right)\right]}^2$ and $B_i={\mathrm{\Omega }}_{\mathrm{bi}}{p}_{\mathrm{ri}}/T_{\mathrm{ri}}$. Ω_a, Ω_b and m(ω) are specific to each EoS.

The logarithm of equilibrium constants is

$$\ln K_i=C_0+C_1\sqrt{A_i}+C_2{B}_i+\sum _{k=3}^{m+2}{C}_k{\gamma }_{k-2,i};\hspace{1em}i=1,n,$$(A4)where

$$C_0=-\ln \left(Z_L-B_L\right)+\ln \left(Z_V-B_V\right),$$(A5a)

$$C_{1p}=\frac{2}{({\delta }_1-{\delta }_2)}\frac{{\lambda }_p^{\left(i\right)}}{B_p}\cdot \ln \left(\frac{Z_p+{\delta }_1{B}_p}{Z_p+{\delta }_2{B}_p}\right);p=L,V,$$(A5b)

$$C_1=C_{1V}-C_{1L},$$(A5c)

$$C_2=\frac{Z_L-1}{B_L}-\frac{Z_V-1}{B_V}-\frac{1}{2}\left(\frac{C_{1L}}{{\lambda }_L^{\left(i\right)}}\frac{A_L}{B_L}-\frac{C_{1V}}{{\lambda }_V^{\left(i\right)}}\frac{A_V}{B_V}\right)$$(A5d)and

$$C_k={\mathrm{\chi }}_k^{\left(i\right)}\left(\frac{C_{1L}}{{\lambda }_L^{\left(i\right)}}{\theta }_{\mathrm{kL}}-\frac{C_{1V}}{{\lambda }_V^{\left(i\right)}}{\theta }_{\mathrm{kV}}\right);k=3,m+2,$$(A5e)where

$$A_p={\alpha }_p^2+2\sum _{k=1}^m{\theta }_{\mathrm{kp}}{\gamma }_{\mathrm{kp}}+\sum _{k=1}^m{\theta }_{\mathrm{kp}}^2;p=L,V,$$(A6)

$$\alpha =\sum _{j=m+1}^{\mathrm{nc}}{y}_j\sqrt{A_j},$$(A7)

$${\theta }_k=y_k\sqrt{A_k};\hspace{1em}k=1,m,$$(A8)

$${\gamma }_k=\sum _{i=k+1}^n{y}_i{\gamma }_{\mathrm{ki}};\hspace{1em}k=1,m,$$(A9)

$${\gamma }_{\mathrm{ki}}=\sqrt{A_i}\left(1-k_{\mathrm{ki}}\right);\hspace{1em}k=1,m;\hspace{1em}i=k+1,n,$$(A10)

$${\lambda }^{\left(i\right)}=\{\begin{array}{l}{\theta }_i+{\gamma }_i\\ \alpha \end{array}\begin{array}{ll};& i=1,m\\ ;& i=m+1,n\end{array}$$(A11)and ${\mathrm{\chi }}_j^{\left(i\right)}=0;i\le j;{\mathrm{\chi }}_j^{\left(i\right)}=1;i>j$.

Finally, the elements of the reduction matrix in equation (6) are: ${\mathit{q}}_{\mathbf{0}\mathit{i}}=\mathbf{1},{\mathit{q}}_{\mathbf{1}\mathit{i}}=\sqrt{{\mathit{A}}_{\mathit{i}}}$, ${\mathit{q}}_{\mathbf{2}\mathit{i}}={\mathit{B}}_{\mathit{i}}$, and ${\mathit{q}}_{\mathit{ki}}={\mathit{\gamma }}_{\mathit{ki}};\mathit{k}=\mathbf{3},\mathit{m}+\mathbf{2}$, for i = 1, …, n.

References

All Tables

All Figures