Application of near critical behavior of equilibrium ratios to phase equilibrium calculations

. We examine the asymptotic behavior of the equilibrium ratios ( K i ) near the convergence locus in the pressure-temperature 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 ln K i tends to zero, according to a power ½ of the distance to this locus. As a consequence, if ln K i is expressed as a linear combination of pure component parameters with coef ﬁ cients only depending on mixture phase properties ( i


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 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 lawcharacterized 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(Zaydullin et al., , 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).

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 6 ¼ 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 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 Rachford-Rice function 1= 1 , is wider (and in many cases much wider) than the physical interval V 2 [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 ?±1, 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 DT = T À T c or DP = P À P c , or of their dimensionless counterparts DT = (T À T c )/T c and DP = (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 |DT| or |DP|, 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 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 in the vicinity of the critical point along the critical isotherm, T = T c , in the two-phase region, and on the critical isobar, P = P c , in the two-phase region, where Þ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 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 b % 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): 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 ffiffiffiffiffi A i p 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., 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 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 . ., n: the first m components have nonzero BIPs with the remaining ones) such that: 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, . ., 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., liquidvapor) region: the equilibrium ratios K i tend asymptotically to unity (or, equivalently, lnK i tend to 0) according to the "universal" scaling law ðT conv À T Þ 1=2 or ðP 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).

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 precalculated 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 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 where dP=dn ¼ À2nP conv , since The partial derivatives (oC k /oP) 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 non-zero 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 then the coefficients C k can be calculated as functions of pressure from 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 P n i¼1 z i =K i ¼ 1 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 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.

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

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.78K and P c = 210.66bar.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 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 2 > 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 6 ¼ 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.39bar, and the convergence pressure is P conv = 229.09bar.At P sat , the value of the first coefficient is C 0 = 0.23876, and its derivative with respect to pressure is oC 0 =oP ð Þ sat T ¼ À0:025.From equation ( 11), taking P* = P sat , an excellent approximation of the convergence pressure is obtained: P conv = 229.16bar (0.07 bar absolute error; the relative error is about 0.03%).Note that similar results are obtained using oC 0 =oP ð Þ 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.36bar, 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 conv Þ 0:5 at P = 100 bar shows a similar behavior (with T conv = 554.32K 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.

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

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  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 k vs. ð1 À P=P 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.

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. ð1 À P=P conv Þ 0:5 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 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.

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.The general form of two-parameter cubic EoS is Equation (A1) includes the SRK EoS (Soave, 1972, for d 1 = 1 and d 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 where Ã 2 and B i ¼ X bi p ri =T ri .X a , X b and m(x) are specific to each EoS.The logarithm of equilibrium constants is where and where and v ðiÞ j ¼ 0; i j; v ðiÞ j ¼ 1; i > j.Finally, the elements of the reduction matrix in equation ( 6) are: q 0i ¼ 1; q 1i ¼ ffiffiffiffiffi A i p , q 2i ¼ B i , and q ki ¼ c ki ; k ¼ 3; m þ 2, for i = 1, . .., n.

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* = P sat = 224.39bar).
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

Table 3 .
Reservoir fluid composition and component properties.

Table 4 .
Non-zero BIPS for the reservoir fluid.Phase envelope and convergence locus of the reservoir fluid.D.V.Nichita et al.: Oil & Gas Science and Technology -Rev.IFP Energies nouvelles 74, 77 (2019)