Regular Article
Application of near critical behavior of equilibrium ratios to phase equilibrium calculations
^{1}
CNRS UMR 5150, Laboratoire des Fluides Complexes et leurs Réservoirs, Université de Pau et des Pays de l’Adour, BP 1155, 64013 Pau Cedex, France
^{2}
Centre Scientifique et Technique Jean Feger, Total, avenue Larribau, 64018 Pau Cedex, France
^{*} Corresponding author: dnichita@univpau.fr
Received:
14
June
2019
Accepted:
18
September
2019
We examine the asymptotic behavior of the equilibrium ratios (K_{i}) near the convergence locus in the pressuretemperature plane. When the Equation of State (EoS) is analytical, which is the case of most EoS of engineering purpose, K_{i} tends towards unity or, equivalently, its logarithm lnK_{i} tends to zero, according to a power ½ of the distance to this locus. As a consequence, if lnK_{i} is expressed as a linear combination of pure component parameters with coefficients only depending on mixture phase properties (i.e., reduction parameters), these coefficients obey a similar power law. Deviations from the ½ power law are thus fairly limited for lnK_{i} and for the reduction parameters (at least in the negative flash window between the convergence locus and the phase boundaries), which can be exploited to speed up flash calculations and for quickly determining approximate saturation points and convergence pressures and temperatures. The chosen examples are representative synthetic and natural hydrocarbon mixtures, as well as various injection gashydrocarbon systems.
Key words: critical point / asymptotic behavior / convergence pressure / negative flash / reduction method / equation of state / Binary Interaction Parameters (BIPs)
© D.V. Nichita et al., published by IFP Energies nouvelles, 2019
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_{i} : Component parameter in the EoS
B_{i} : Component parameter in the EoS
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
M: Number of reduction parameters
q_{αi} : Elements of the reduction matrix
x_{i} : Liquid mole fraction, component i
y_{i} : Vapor mole fraction, component i
z_{i} : Feed mole fraction, component i
Greek letters
φ_{i} : Fugacity coefficient of component i in a mixture
Subscripts
Superscripts
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 multicomponent fluids with a dozen or more components or pseudocomponents, 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}; highquality 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 socalled “negative flash”, as an alternative to phase stability testing.
Twophase vaporliquid 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 Kvalues 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 tielinebased 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 twoparameter 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 computerassisted 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 multicomponent 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 dewpoint 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 dewpoint pressure and P_{conv} corresponds to “negative flashes” (Whitson and Michelsen, 1990): the flash calculation gives a nontrivial 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 dewpoints). 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 twophase 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 RachfordRice 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, Kvalues 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 nearcritical behavior in the twophase liquidvapor 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 Taylorexpanded 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 meanfield) 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 Taylorexpansion (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 thirdorder 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 = T_{c}, in the twophase region, and
(2)on the critical isobar, P = P_{c}, in the twophase 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 lnK_{i} near a critical point (where K_{i} − 1 are small), but on smaller pressure intervals, since
At this point it is important to remind that the above meanfield or Landautype approach provides only an approximation of the nearcritical behavior of fluids, which is characterized by scaling exponents that differ from the meanfield (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 multicomponent mixture. A widely occurring circumstance is when the EoS is a twoparameter cubic EoS, such as the PengRobinson (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 C_{0}, C_{1} and C_{2} 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 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_{0}, C_{1} and C_{2} 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) 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 twoparameter cubic EoS with nonzero 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 (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 C_{k}, and equation (5) can be rewritten in the more compact form,
(6)where q_{ki} (, , and , 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 pseudocomponent 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 twophase (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 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 twophase 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 twoparameter 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 quasilinearity 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 twophase region near phase boundaries where flash calculations are usually more difficult (as compared to flashes well inside the twophase 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 twophase flash near a phase boundary. Third, it can be used to inter and extrapolate Kvalues from known values. Fourth, it allows a rapid approximate determination of saturation pressures.
Another advantage of this formulation is that precalculated tables of Kvalues 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
(7)the equation of a straight line between a certain reference pressure, P* (where the results of a negative or twophase flash calculation are available), and P_{conv} is
At P = P_{conv}, equation (8) reads
(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_{0} (note that its expression does not depend on the number of nonzero BIPs).
The convergence pressure (its approximation considering a linear variation) is obtained from equations (9) and (11) written for k = 0:
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
(12)then the coefficients C_{k} can be calculated as functions of pressure from
(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 twophase region or in the negative flash region) are available, one can rapidly estimate the phase boundary location. Using Kvalues 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 Kvalues tables (some compositional simulators have an option for using Kvalues tables for flash calculations). For domains where the coefficients C_{k} are nonlinear, a procedure similar to that presented by Chien and Lee (1983) for Kvalues can be set up. Note that C_{k} generally have a “smoother” variation with pressure than Kvalues, and less storage is required.
5 Results
The linear dependence of the Kvalues and of C_{k} with , which we call hereafter quasilinearity, is tested for a model gascondensate system (Y8) taken from the literature (Yarborough, 1972), and for mixtures of this gascondensate with nitrogen (Y8/N_{2}): our purpose is to study the influence of an injection gas on convergence pressures and Kvalues 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 normalalkanes, 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.
Fig. 1 Phase envelope and convergence locus of the Y8 mixture. 
Composition and component properties of the Y8 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_{i}vs. (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 quasilinearly when pressure approaches the critical pressure. Interestingly, lnK_{i} obeys this quasilinear 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 at P = P_{c} (Fig. 4b) over a large pressure or temperature interval (comparable to that of quasilinear behavior for lnK_{i}): 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.
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. 
Fig. 3 (a) Equilibrium ratios vs. (1 − T/T_{c})^{0.5} on the critical isobar (P = 210.66 bar) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. (1 − T/T_{c})^{0.5} on the critical isobar (P = 210.66 bar) for the Y8 mixture. 
Fig. 4 (a) Coefficients C_{k} (k = 0, 1, 2) vs. (1 − P/P_{c})^{0.5} on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Coefficients C_{k} (k = 0, 1, 2) vs. (1 − T/T_{c})^{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 ≠ T_{c}. On the isotherm of T = 335 K (in the region where retrograde condensation occurs), the Kvalues (Fig. 5) and coefficients C_{k} (Fig. 6) are plotted vs. pressure up to the convergence pressure. Note again the quasilinear 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} = 0.23876, and its derivative with respect to pressure is . 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 at any pressure between P_{sat} and P_{conv}, and for a certain pressure interval inside the twophase 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).
Fig. 5 Kvalues vs. pressure for the Y8 mixture at T = 335 K. 
Fig. 6 Coefficients C_{k}vs. (1 − P/P_{conv})^{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 Kvalues 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 Kvalues (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 Kvalues 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 widelyused Wilson (1969) ideal Kvalues relationship (which are poor initial estimates for these conditions). It should be noted that AFC is an iterative method (it solves a nonlinear 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 Kvalues 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 Kvalues; the relative errors in Kvalues for each component are given for comparison in Figure 7. In fact, excellent approximations of Kvalues can be obtained for a wide pressure interval.
Fig. 7 Relative errors on Kvalues for different (various) initializations. Y8 mixture, flash at T = 335 K and P = 215 bar. 
Exact flash results and different initial guess for Kvalues for Y8 mixture at T = 335 K and P = 215 bar (the reference pressure is P* = P_{sat} = 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 nonlinear 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 twophase region or in the negative flash region) and we want a rapid estimate of the phase boundary location. For example, using P* = 215 bar (twophase 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 Kvalues (Fig. 8) and coefficients C_{k} (Fig. 9) vs. 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.
Fig. 8 Kvalues vs. temperature for the Y8 mixture at P = 100 bar. 
Fig. 9 Coefficients C_{k}vs. (1 − T/T_{conv})^{0.5} for the Y8 mixture at P = 100 bar. 
5.2. Y8/Nitrogen mixtures
In order to investigate the case of nonzero BIPs, different amounts of nitrogen are combined with the Y8 mixture. The BIPs between N_{2} and hydrocarbon components are taken equal to 0.1 (m = 1). For a mixture made up of the Y8 fluid + 25% moles N_{2}, 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 nonzero 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_{2}, 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.
Fig. 10 Phase envelope and convergence locus of Y8 + 25% N_{2} mixture. 
Fig. 11 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the Y8 + 25% N_{2} mixture at T = 300 K. 
Fig. 12 Phase envelope and convergence locus of Y8 + 50% N_{2} mixture. 
Fig. 13 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the Y8 + 50% N_{2} mixture at T = 300 K. 
5.3. Reservoir fluid
The reservoir fluid composition is described by 29 components; eight components have nonzero BIPs with the remaining ones (m = 8). Mixture composition and component properties are given in Table 3 and the nonzero 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.
Fig. 14 Phase envelope and convergence locus of the reservoir fluid. 
Reservoir fluid composition and component properties.
Nonzero BIPS for the reservoir fluid.
There are 18 (2m + 2) reduction parameters, and 11 coefficients C_{k} in equation (6). Figure 15 depicts the coefficients C_{k}vs. 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, quasilinearity holds for all coefficients C_{k} over a large pressure interval.
Fig. 15 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the reservoir fluid at T = 500 K. 
5.4. CO_{2}/reservoir fluid mixture
In order to study the influence of an injection gas on the new regularity, CO_{2} is added in various amounts to the above reservoir fluid. The phase envelope and the convergence locus for the reservoir fluid + 50% moles CO_{2} mixture are presented in Figure 16. Figure 17 plots the coefficients C_{k}vs. for the reservoir fluid + 50% moles CO_{2} mixture at T = 500 K; some of the C’s (i.e., C_{0}, C_{1}, C_{2}) are somewhat deviating from linearity, while for the Cparameters 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 Cparameters.
Fig. 16 Phase envelope and convergence locus for the reservoir fluid + 50% CO_{2} mixture. 
Fig. 17 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the reservoir fluid + 50% CO_{2} 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 Kvalues 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 Kvalues 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 twophase pressure/temperature domain (except for systems containing a hydrocarbon mixture and a high concentration of a nonhydrocarbon 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 Kvalues 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 twoparameter 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 onefluid mixing rules are used for the energy, A and for the volume, B parameters of the EoS
(A3)where and . Ω_{a}, Ω_{b} and m(ω) are specific to each EoS.
The logarithm of equilibrium constants is
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, SPE12243MS, Reservoir Simulation Symposium, November 15–18, San Francisco, CA. [Google Scholar]
 Dalton B.J. (1970) Twophase 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 twocomponent system from an empirical equation of state, including liquidvapor equilibria data, Technical Report BMRI7076, 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 pressuretemperature 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 equationofstate 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 analyticalcomponent delumping procedure for equations of state with nonzero 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 twoconstant 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) Tielinebased Kvalue 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 PengRobinson programs, Research Report, Gas Processors Association, Tulsa, Okla, RR28. [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) Noniterative 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, SPE37979MS 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 RedlichKwong 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) Vaporliquid equilibrium data for multicomponent mixtures containing hydrocarbon and nonhydrocarbon 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 EoSbased and Kvaluesbased methods for threephase thermal simulation, Transp. Porous Media 116, 663–686. [Google Scholar]
All Tables
Exact flash results and different initial guess for Kvalues for Y8 mixture at T = 335 K and P = 215 bar (the reference pressure is P* = P_{sat} = 224.39 bar).
All Figures
Fig. 1 Phase envelope and convergence locus of the Y8 mixture. 

In the text 
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 
Fig. 3 (a) Equilibrium ratios vs. (1 − T/T_{c})^{0.5} on the critical isobar (P = 210.66 bar) for the Y8 mixture. (b) Logarithm of equilibrium ratios vs. (1 − T/T_{c})^{0.5} on the critical isobar (P = 210.66 bar) for the Y8 mixture. 

In the text 
Fig. 4 (a) Coefficients C_{k} (k = 0, 1, 2) vs. (1 − P/P_{c})^{0.5} on the critical isotherm (T = 293.78 K) for the Y8 mixture. (b) Coefficients C_{k} (k = 0, 1, 2) vs. (1 − T/T_{c})^{0.5} on the critical isobar (P = 210.66 bar) for the Y8 mixture. 

In the text 
Fig. 5 Kvalues vs. pressure for the Y8 mixture at T = 335 K. 

In the text 
Fig. 6 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the Y8 mixture at T = 335 K. 

In the text 
Fig. 7 Relative errors on Kvalues for different (various) initializations. Y8 mixture, flash at T = 335 K and P = 215 bar. 

In the text 
Fig. 8 Kvalues vs. temperature for the Y8 mixture at P = 100 bar. 

In the text 
Fig. 9 Coefficients C_{k}vs. (1 − T/T_{conv})^{0.5} for the Y8 mixture at P = 100 bar. 

In the text 
Fig. 10 Phase envelope and convergence locus of Y8 + 25% N_{2} mixture. 

In the text 
Fig. 11 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the Y8 + 25% N_{2} mixture at T = 300 K. 

In the text 
Fig. 12 Phase envelope and convergence locus of Y8 + 50% N_{2} mixture. 

In the text 
Fig. 13 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the Y8 + 50% N_{2} mixture at T = 300 K. 

In the text 
Fig. 14 Phase envelope and convergence locus of the reservoir fluid. 

In the text 
Fig. 15 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the reservoir fluid at T = 500 K. 

In the text 
Fig. 16 Phase envelope and convergence locus for the reservoir fluid + 50% CO_{2} mixture. 

In the text 
Fig. 17 Coefficients C_{k}vs. (1 − P/P_{conv})^{0.5} for the reservoir fluid + 50% CO_{2} mixture at T = 500 K. 

In the text 